Predictability of the solar cycle over one cycle

The prediction of the strength of future solar cycles is of interest because of its practical significance for space weather and as a test of our theoretical understanding of the solar cycle. The Babcock-Leighton mechanism allows predictions by assimilating the observed magnetic field on the surface. But the emergence of sunspot groups has the random properties, which make it impossible to accurately predict the solar cycle and also strongly limit the scope of cycle predictions. Hence we develop the scheme to investigate the predictability of the solar cycle over one cycle. When a cycle has been ongoing for more than 3 years, the sunspot group emergence can be predicted along with its uncertainty during the rest time of the cycle. The method for doing this is to start by generating a set of random realizations which obey the statistical relations of the sunspot emergence. We then use a surface flux transport model to calculate the possible axial dipole moment evolutions. The correlation between the axial dipole moment at cycle minimum and the subsequent cycle strength and other empirical properties of solar cycles are used to predict the possible profiles of the subsequent cycle. We apply this scheme to predict the large-scale field evolution from 2018 to the end of cycle 25, whose maximum strength is expected to lie in the range from 93 to 155 with a probability of 95\%.


Introduction
The various proxies of solar activity show a roughly 11-year cycle period with widely varying amplitudes. The cycle amplitude gives a rough idea of the frequency of space weather storms of all types, from radio blackouts to geomagnetic storms to radiation storms. For this reason there is a practical need for solar cycle prediction in our technological society. The longer the prediction extends with higher reliability, the more desirable it is for decision makers. Moreover, the predictions, especially the dynamo-based prediction, can also provide us the opportunity to examine our insight into the physical mechanisms underlying the solar cycle. Existing attempts to predict the level of solar activity can be broadly divided into two groups: (1) those predicting an ongoing cycle, and (2) those predicting the next or future cycle(s). The first group can be further divided into medium-term predictions (months in advance) and long-term predictions (years in advance). Three methods for short term prediction of ongoing cycles are the McNish-Lincoln method (McNish & Lincoln 1949), the standard method, and the combined method at Sunspot Index and Long-term Solar Observations (SILOS) 1 . For the long-term ongoing cycle predictions, people usually use the curve-fitting function (Waldmeier 1935;Hathaway et al. 1994;Li 1999;Du 2011;Li et al. 2017) or use the similarity of solar cycles, e.g., Li et al. (2002).
For the long-term predictions of the future cycle(s), existing attempts can be broadly divided into empirical methods and methods based on dynamo models. The empirical methods can be further divided into two subgroups. One group is based on extrapolating the past records using a purely mathematical or statistical analysis, which has wide applications in econometrics 2 . The other subgroup of the empirical methods are the precursor methods, which are based on correlations between certain measured quantities in the declining phase of a cycle and the strength of the next cycle, e.g., Thompson (1993); Svalgaard et al. (2005); Cameron & Schüssler (2007). Around the end of cycle 23, two group of people made the dynamo-based solar cycle prediction available for the first time Jiang et al. 2007;. The different observed data assimilated into prediction models and different understandings of the flux transport mechanisms cause markedly distinct predictions for cycle 24 (Karak et al. 2014). For a review about the solar cycle prediction, see Petrovay (2010) and references therein.
Although there have been numerous attempts to predict the solar cycle, the reliability of the different methods is still controversial. The predictions of cycle 24 using the different prediction methods, with a comparison to the true strength of solar cycle 24, are clearly presented in Pesnell (2012) who analysed 75 predictions for cycle 24. The proved most successful method is the precursor method based on the polar field (Schatten et al. 1978; The polar field evolution can be understood in terms of the Surface Flux Transport (SFT) model (Wang et al. 1989;Mackay & Yeates 2012;Jiang et al. 2014b). The tilt angle of the sunspot group with respect to the E-W direction, leads to each active region contributing to the global axial dipole moment of the Sun, and is an essential factor in the polar field generation. The tilt angle has a systematic component, which corresponds to the Joy's law, and a larger random component, i.e., the tilt angle scatter. Jiang et al. (2014a) found the observed scatter of the tilt angles causes a variation of 30%-40% in the resulting polar field around activity minima. By including detailed information about the individual tilt angles and magnetic polarities of the bipolar magnetic regions (BMRs) that emerged during cycle 23, the weak polar field around the end of cycle 23 is well reproduced by Jiang et al. (2015). Recently series of BL dynamo models successfully reproduce the variability of the solar cycle by introducing the random BMR tilts (Cameron & Schüssler 2017;Olemskoy & Kitchatinov 2013;Jiang et al. 2013;Nagy et al. 2017;Karak & Miesch 2017). The randomness of the sunspot emergence -part of the toroidal-to-poloidal part of the dynamo loop -makes the dynamo a stochastic process. Even if the poloidal-totoroidal part of the dynamo loop was fully deterministic, the intrinsic random features of flux emergence limits the scope of the solar cycle prediction. Furthermore, the randomness of the sunspot emergence due to the complex flux emergence dynamics (Weber et al. 2011(Weber et al. , 2013 also causes the uncertainty even for the prediction of an ongoing cycle. The randomness raises the question of whether the long-term predictions with reasonable levels of uncertainty are possible. A second question is the extent to which we can predict the time-dependence of the activity of a solar cycle beyond only predicting the cycle amplitude. Jiang et al. (2011) show the dependence of the statistical properties of sunspot emergence on the cycle phase and strength based on an empirical analysis of the historical sunspot number data (see also Muñoz-Jaramillo et al. 2015). Random realizations of sunspot group emergences during a cycle can be reconstructed for a given maximum sunspot number based on the statistical properties. The SFT model can be used to produce the large-scale field evolution over the Sun's surface. Cameron et al. (2016) and Jiang & Cao (2017) have applied the above scheme in a Monte-Carlo approach to predict the range of the polar field around the end of cycle 24 about 3-4 years before the minimum, including error bars. Similar approaches are also presented in Hathaway & Upton (2016) and Iijima et al. (2017). In this paper, we generalize and update the previous method to investigate the predictability of the solar cycle over one cycle at different phases of a cycle. The predictions include the time evolution of the monthly sunspot number of an ongoing cycle and the time evolution of the smoothed sunspot number over one cycle, all with estimated uncertainties. The investigation of the predictability beginning from different starting points distinguishes our scheme from the existing predictions of the solar cycle.
The paper is organized as follows. In Section 2, we study the properties of the solar cycle profiles, which will be adopted by Section 3 for the prediction of sunspot emergence for an ongoing cycle. The SFT model which will be used to describe and to predict the large-scale field evolution over the solar surface is presented in Section 4. The correlation between the dipole moment at cycle minimum and the subsequent cycle strength based on a homogeneous axial dipole moment dataset is given in Section 5. In Section 6, we present the results about the predictability of the subsequent cycle and its application to cycle 25. Our summary and discussion are given in Section 7.

Properties of solar cycle profiles
At present there are different time series of the sunspot number, which differ substantially before the cycle 12 (Clette et al. 2014). Hence in this study we only investigate the properties of solar cycle profiles during cycles 12-24, and have used the Sunspot Number Version 2.0 3 . The data are of the monthly mean total sunspot number (R mn ) and the corresponding 13-month smoothed sunspot numbers (R sm ) . The timing of each cycle minimum used throughout the paper is taken from the NGDC 4 . We first analyse the properties of the smoothed sunspot numbers. The monthly sunspot number will be analyzed at the end of this section. Figure 1 shows the observed smoothed sunspot number as a function of months from start of a cycle for cycles 12-24. We define the maximum phase of a cycle as the time period when the sunspot number surpasses 70% of maximum sunspot number of the cycle. The ascending phases, the maximum phases, and the declining phases are denoted in the dotted, the dashed, and the solid curves in Figure 1, respectively. There are two distinctive features about the shape of the solar cycle. The first is about the rising phase, which obeys the Waldmeier effect (Waldmeier 1955). It is that stronger cycles tend to show a faster rise of activity levels during their ascending phase than weaker cycles (Lantos 2000;Cameron & Schüssler 2008). The second is that once the solar cycle begins to decline, all cycles decline in a similar way (Ivanov & Miletsky 2014;Cameron & Schüssler 2016). A large scatter of the individual cycle about the means over all the cycles is also shown. During the declining phase, the profile of solar cycle shows noisier short-term variation than the early time period. These properties can be understood under the framework of the BL dynamo. A weaker polar field at the beginning of a cycle generates less toroidal field to form the sunspot groups than it would if the polar field was strong. This provides an explanation to the Waldmeier effect and also is demonstrated by Karak & Choudhuri (2011) using a BL dynamo model. Cameron & Schüssler (2016) interpret the similar decline phases in terms of oppositely directed toroidal flux bands in each hemisphere that diffuse and cancel across the equator when the distance of the center of the activity belts from the equator becomes about equal to their width. Stronger cycles show wider activity belts and thus start to decline earlier than weaker cycles. Nagy et al. (2017) and Kitchatinov et al. (2018) showed that peculiar BMRs with large tilt angles emerging during the rising phase of a cycle have large effects on the amplitude of the descending phase. Hence, the stochastic mechanism due to random emergence of the peculiar BMRs during the rising phase causes the later phases of the cycle to be noisy, as seen in Figure 5.
Several of these properties have been exploited in the development of a parameterization of the sunspot number during a cycle in terms of two simple paramters, the starting time (t 0 ) and amplitude (a) (as given in Hathaway et al. 1994, hereafter HWR94). The determination of the two parameters can be made a few years after the start of a cycle. The function captures the essence of the Waldmeier effect. The rising phase is presented in the form of t 3 , where t is the time in months from the start of a cycle. We use the function, i.e. Equation (1) of HWR94 to fit the cycles 12-24 individually. It is in the form of where b(a) = 27.12 + 25.15/(a × 10 3 ) 1/4 and c = 0.71. The solar cycle overlap, i.e., new cycle spots appearing at high latitudes while old cycles still seen at low latitudes, is about 3 years (Wilson et al. 1988;Harvey 1992). We set weights to the first 1.5 year of each cycle by a error function in the form of 1 2 [1 + erf(2 t−tovp tovp )], where t ovp =0.75yr to decrease the effect of the data during the first 1.5 year on the cycle fit. Figure 2 shows examples of the cycle fit overplotted with the observed sunspot number evolution from 1878 to the present, based on the fit either 4 years (upper panel), or 8 years (lower panel), after the start of the cycle. Since the function captures the essence of the Waldmeier effect, the fit and smoothed sunspot number are consistent with each other during the ascending phase. We measure the relative error s f2o between the fit and the observation during the maximum phases for different cycles when the fit is done at different timings t of cycles. Here s f2o is defined as where N is the number of months of the maximum phase for different cycles. The results are shown in Figure 3. The fits to cycles 12-24 are done from 1.5 years to 6 years into cycles with 0.5 year interval. The strength of a cycle defined by the maximum sunspot number is indicated by a given color shown by the color bar. The stronger cycles shown in colors between the green to the red tend to have smaller relative errors. They are usually less than 20% even the fits start from 2 years onwards. The weaker cycles shown in colors between the black to the cyan tend to match with the fits more poorly. From 3 years into a cycle, s f2o is always within 50%. This leads us to only make ongoing cycle prediction from at least 3 years into a cycle throughout the rest of this paper.
During the later phases, the fitting curves in Figure 2 shows larger deviations from the observations. BMRs during the descending phase tend to locate at low latitudes. Jiang et al. (2014a) demonstrated that the BMRs with the same initial axial dipole moment at lower latitudes have larger contributions to the polar field. Hence the BMRs during the descending phase play important roles in the polar field evolution. This motivates us to suggest an improvement in order to improve the predictions. We use two methods to measure the differences △f (t) between the observation R sm (t) and the fitted function f (t). The first is where where △f (t) will be added to function (1) to correct the systematic deviation from the cycle fits especially in the descending phase. Meantime, we estimate the fluctuation of the shape of cycles by measuring the standard deviation (σ △f ′ (t) ) of fits including the correction using △f (t) from observations. The △f (t) and σ △f ′ (t) are calculated by where the symbol i denotes different cycles 12-24 and N is equal to 13. The random component △f ′ (t) will be added to predict the sunspot number evolution as a set of random realizations with zero average and the standard deviation of σ △f ′ (t) .
The other method to correct the fit is where and where i-values are the same as that in Eq.(6), which are the 13 cycles from cycle 12 to cycle 24. The random component △r ′ (t) will be added to predict the sunspot number evolution as a set of random realizations with zero average and the standard deviation of σ △r ′ (t) . Figure 4 shows the time evolution of △f (t), σ △f ′ (t) , △r(t), and σ △r ′ (t) when the fits are done from 3 years or more into a cycle. The curves in each panel correspond to the different timings (from the 3rd years to the 9th years into cycles with one year cadence) to do the cycle fits. We see that △r(t) and σ △r ′ (t) increase for the later phase of cycles. When the fits are done at the 3rd year, the 4th year and the 5th year, △r(t) and σ △r ′ (t) have the almost same profiles. If the fits are done at the 6th year to the 8th year, there is a significant increase of the match between the observed values and the fitting function. And the △r(t)-values decrease accordingly. Close to the minimum time period, i.e., from the 9th year into a cycle, there are large deviations among different fits. The small sunspot number during the decline phase contributes to the large variations of the △r(t)-value. On the contrast, △f (t) and σ △f ′ (t) during this time period weakly depend on the fits which are done at different phases of a cycle.
The fitting functions for the time evolution of △r(t) and σ △r ′ (t) are listed in the following.

45.86
and z 2 = t−149.0 27.31 . That Eq.(11) or Eq. (12) is used depends on when the prediction is done. When it is less than 72 months into a cycle, we take Eq.(11). When it is more than 72 months and less than 108 months into a cycle, we take Eq.(12). The corresponding σ △r ′ (t) has a weak dependence on the fitting time and satisfies the following form: where z 3 = t−318.73 77.35 . During the later phases, the △f (t) obeys the following profiles, where where z 5 = t−108.16 1.455 . The two parameters △f (t) and σ △f ′ (t) are shown in the red curves in the lower left and the lower right panels, respectively. Now we analyze the monthly sunspot number R mn . The deviation of the R mn from the smoothed sunspot number R sm is regarded as a random component. We measure the relative deviations s mn2sm of R mn from R sm . It is defined as where i is the 13 cycles from cycle 12 to 24. The result is shown in Figure 5. We see during the maximum phase the relative deviations s mn2sm keep at a small value of about 0.2. Then the value increases with time during the decline phase, which means noisier later phases of solar cycle than earlier phases. The reason has been given in the second paragraph of this section. During the first 2 years, the large s mn2sm is due to the cycle overlap. The fitting functions are s mn2sm = 0.936 exp(−p 2 1 /2) + 0.196, t < 36 months, 0.375 exp(−p 2 2 /2) + 0.196, t ≥ 36 months, where p 1 = t+12.21 14.346 and p 2 = t−128.4 25.357 .

Prediction of the sunspot number evolution
With the features of the solar cycle shape listed in Section 2, we may predict the time evolution of the monthly sunspot number F mn (t) for the remainder of a cycle when it is some months, denoted as n, entering into the cycle. We separate the sunspot number into the systematic component and the random component. We first use the method described in Section 2 to fit the observed smoothed sunspot number to derive the function f (t). Then we get the predicted systematic part of smoothed sunspot number evolution F sm (t)=f (t) + △f (t) or f (t) + △r(t)f (t) for the different values of n according to Eqs.(11), (12), and (14). The random deviations from the systematic component measured by Eqs. (13) and (15) are determined by using Monte-Carlo simulations. The random realization of the sunspot emergence is denoted as F sm (t), which is equal to where △f ′ (t) and △r ′ (t) have zero averages, and their standard deviations satisfy Eqs. (13) and (15). The monthly sunspot number where the standard deviation of r sm satisfies Eq.(17) and the average of r sm is zero.
We take two timings, 4 years and 8 years into cycles, as examples to compare the differences between the predictions of ongoing cycles F sm (t) and the observations R sm (t). The difference is measured by the goodness-of-fit, which is given by where σ(t i ) corresponds to Eq.(6) and N corresponds to the numbers of the months from the timing of the prediction to the end of the cycle. When χ is equal to 1.0, it indicates that the fitting function passes within one standard deviation of the observed data points. Table 1 shows χ values of the predicted results of the sunspot number evolution during cycles 12-24 at 4 years and 8 years into cycles. The χ values are also calculated for the prediction just based on HWR94 method, i.e., fits using Eq.(1). We see that χ values are always much smaller than that of HRW94, especially when the predictions are at the earlier phase, i.e., 4 years into a cycle. This demonstrates the improvement of the predictive skill by the current strategy comparing with HWR94.
Another advantage of the method is to quantify the uncertainty of the predictions, which corresponds to the predictability of the sunspot number evolution. Tables 2 and  3 give the evaluation of the predictive skill of the smoothed sunspot number F sm (t) and the monthly sunspot number F mn (t), respectively. The values show the percentages of the observed sunspot numbers, which are out of the 1σ, 2σ and 3σ of the predicted variations.
For some cycles, e.g., cycle 12 and cycle 20, the percentages are larger than 32% and 4.6% for 1σ and 2σ fluctuations respectively. But all the observed sunspot numbers are within 3σ range of the prediction. The predictive skill for monthly sunspot number is worse than that of the smoothed sunspot number. Figure 6 shows the comparisons between the predicted and the observed sunspot number and the latitudinal location of sunspot groups (discussed in the next section). We take the following cases as examples: (1) cycle 20 (the first row) and cycle 23 (the second row) predicted at 4 years into cycles and showing smaller goodness-of-fits than HRW94, and (2) cycle 12 (the third row) and cycle 14 (the forth row) predicted at 8 years into the cycles and showing larger goodness-of-fits than HRW94. The first column is the smoothed monthly total sunspot number. Although the mean values of our predictions still always show the deviations from the observations, the observations are usually within our predicted 2σ uncertainty ranges denoted by the shaded regions. This also demonstrates the importance of giving the uncertainty range. The second column is the time evolution of the monthly sunspot number. The random realizations look very similar, given the uncertainty limits, to the observed one.

Prediction of the BMR emergence
Once we have the predicted sunspot number, the next step is to predict the emergence properties of the sunspot groups. Jiang et al. (2011) give the dependence of the statistical properties of sunspot group emergence, including the latitude, longitude, area and tilt angle of sunspot groups in the form of the BMRs on the cycle phase and strength. Those properties are used here. The following list summarizes the statistical properties aiming to convert the sunspot number time series into the time series of sunspot group emergence.
• The number of BMRs emerging per month was taken to be equal to F BM R (t) = 0.24F mn (t). The number of the daily emergence of the BMR is randomly realized satisfying the total number F BM R (t).
• The width of the latitude distribution obeys a Gaussian profile with a half width of σ, which is equal to (0.14 + 1.05x − 0.78x 2 )λ n (x). We exclude points deviating from the mean by more than 2.2 σ for the equatorward side. • The BMR emergence has the symmetric distribution in north and south hemispheres.
• The BMR emergence has the random longitudinal distribution.
• The number density function of sunspot group areas obey the following distributions, • The average areas of sunspot groups are cycle phase (x) dependent, which obeys • The umbra area of each BMR is A U = A S /5. The total area of BMR A = A S + A P , where A P is the plage area and is calculated by Chapman et al. 1997).
• The mean tilt angle, α n (λ), of the emerging BMRs for the cycle n is assumed to follow Joy's law in the form α n = T n |λ|, where T n = 1.72 − 0.0022S n .  • For the scatter of the tilt angle, ∆α, its averaged value is zero and the standard deviation is σ α . We use the empirical relation between σ α and A U , which is σ α = −11 log A U + 35 to derive σ α (Jiang et al. 2014a).
• A factor 0.7 is multiplied to the resulting tilt angle α n = α n + σ α to account for the effect of latitudinal inflow towards BMRs (Jiang et al. 2010;Martin-Belda & Cameron 2017).
The third column of Figure 6 shows examples of the time evolution of the sunspot latitudinal distribution, i.e., butterfly diagrams. The predicted one is consistent with the observed one.
The BMR emergence provides the source of the large-scale magnetic field over the solar surface. The evolution of the large-scale magnetic field can be described by the SFT model, which will be presented in the following section.  The relevant equation to describe the evolution of the large-scale magnetic flux distribution at the solar surface B(λ, φ, t) as a combined result of the emergence of BMRs, a random walk due to supergranular flows, and the transport by large-scale surface flows is as follows.
where λ and φ are heliographic latitude and longitude, respectively. The latitudinal differential rotation of Sun, Ω(λ), is taken from Snodgrass (1983) and the poleward meridional flow, υ(λ), is given by van Ballegooijen et al. (1998). It is consistent with the measurement by Hathaway & Rightmire (2011). The magnetic diffusivity η is 250 km 2 s −1 , which is within the range from the observational studies summarized by Schrijver et al. (1996). All the three transport parameters are time independent. S(λ, φ, t) is the time dependent flux source term, which will be obtained based on the predicted sunspot emergence using the empirical method in Section 3. The initial magnetic field configuration due to the BMR emergence based on its area, location and tilt is given by Baumann et al. (2004). The corresponding magnetic flux is determined by a single parameter, Bmax (=592G), which is calibrated by the total unsigned surface flux obtained from SOHO/MDI polar field corrected synoptic maps (Sun et al. 2011) after rebinning to the spatial resolution of the simulation (1 • in both latitude and longitude).
For the initial condition of each SFT simulation, we use the SOHO/MDI polar field corrected synoptic maps (available from 1996 June until 2010 May) and the radial synoptic maps with 3600 points in Carrington longitude and 1440 points equally spaced in sine latitude from HMI (2010 May until the present). The abnormal points which correspond to unobserved values or the absolute values over 50G above ±60 • latitudes are filled by the averaged value of neighbouring normal 3 points. Both the MDI and the HMI data are smoothed with the width of 7. Both data were reduced to a resolution of 1 • in latitude and longitude by the IDL CONGRID function and then were converted to the equal latitudes. Rightmost of Figure 8 shows that the cross-calibration of HMI and MDI based on the axial dipole moment during the overlap time period (Carrington Rotations, CRs2097-2104) is 1.3. This is consistent with Liu et al. (2012) who compared Line-of-Sight magnetograms taken by MDI and HMI. We multiply a factor 1.3 when the HMI synoptic map is used as the initial condition. The low resolution data of synoptic maps were used by Cameron & Schüssler (2016) and Jiang & Cao (2017).
We used the code originally developed by Baumann et al. (2004) to do the numerical calculations. The magnetic field is expressed in terms of spherical harmonics up to l = 63. A fourth-order Runge-Kutta method is used for time stepping with the interval of one day.

Prediction of the large-scale magnetic field evolution over surface
With 50 sets of random realizations of the magnetic flux emergence and the synoptic magnetograms as the initial condition, we may derive the possible large-scale magnetic field evolution over the surface using the SFT model. Since here our interest mainly concentrates on long-term predictions, we mainly show the results about the time evolution of the axial dipole moment D(t), which is defined as The uncertainties of the results originate from two ingredients. One is the uncertainty due to the scatter in the properties of the BMRs emergence, which are assumed to be reasonably represented in Sections 2 and 3. The other is the uncertainty due to measurement error in the magnetograms used as initial conditions. The imperfect measurements mainly result from the instrumental characteristics, e.g., spatial resolution, scattered light, and filter characteristics and from the inherent complexity of the solar magnetic field, e.g., the saturated factors for different spectral lines (Wang & Sheeley 1995). The properties of the magnetic field at lower latitudes have larger effects on the results (Jiang et al. 2014a). Since some imperfections of the measurements are inevitable, they are potential uncertainties to the predictive skill of the model. Here we only consider the averaged radial surface field over the surface, which deviates from zero due to measurement error. The same method described in Cameron et al. (2016) to estimate the error is used. The total error due to the magnetogram measurement and the random flux source emergence is determined by adding them quadratically.
We take the timing of 4 years and 8 years into cycles 23 and 24 as examples. Figure 7 shows the predicted time evolutions of the axial dipole moments for cycle 23 and cycle 24 at 4 years and 8 years into the cycles. Synoptic magnetograms CR1961 and CR2013 from MDI and CR2130 and CR2184 from HMI are used as the initial magnetic field, respectively. Solid green lines show the average of 50 SFT simulations with random sources. Dark and light red shading indicate the total σ and 2σ uncertainties, which correspond to the probabilities of 68% and 95.4% for the observed axial dipole moment to be within the error ranges, respectively. The dashed green lines give the 2σ range for the intrinsic solar contribution (source scatter). The errors from MDI measurements are significant and are larger than that from HMI measurements. This also demonstrates the big effects of observed magnetograms on the prediction. Figure 7 shows the following results. Firstly, the observed values are always within the 2σ uncertainty, although statistically it is possible that the observation is out of the 2σ range with the probability of 4.6%. This demonstrates the effectiveness of prediction method. Secondly, the error ranges increase with time. The error ranges of the predicted axial dipole moment at the end of the cycles are smaller when the prediction is made later. This is due to the accumulative effect of the random source emergence. Thirdly, the absolute values of the predicted mean axial dipole moments decrease when the prediction was made at a later time. Cycle 23 had a deep minimum with very low axial dipole moment. Jiang et al. (2015) found that the weakness of the cycle 23 minimum are mainly caused by a number of bigger bipolar regions emerging at low latitudes with a wrong (i.e., opposite to the majority for this cycle) orientation of their magnetic polarities in the north-south direction, which impaired the growth of the polar field. The abnormal emergences are supposed to be emerged during cycle 24 as well. Hence the mean value at cycle minimum decreases with time in both cycle 23 and cycle 24. The corresponding predictions of the following cycle strength decrease with time as well. But it is possible that the predicted mean values of the axial dipole moment increase with time if the peculiar BMR emergences (large areas and tilts at low latitudes) with positive contributions to the axial dipole moment surpass the peculiar BMR emergences with negative contributions to the axial dipole moment. In such case, a stronger subsequent cycle will be proceed.

Correlation between the dipole moment at cycle minimum and the subsequent cycle strength
The prediction of the subsequent cycle depends on the correlation between the dipole moment at cycle minimum and the following cycle strength. The physical base has been clarified in Section 1. This correlation is a key ingredient to affect the prediction accuracy of the subsequent cycle. A long-term homogeneous dipole moment dataset, which is directly calculated based on the observed synoptic magnetograms, is expected. The earliest available synoptic magnetograms are from Mount Wilson Observatory (MWO) from July 1974 until December 2012 5 . Another long-term dataset which is widely used in the solar cycle prediction is the magnetic measurement from Wilcox Solar Observatory (WSO, continuously available from April 1976 until present 6 ). We use the axial dipole moment calculated based on MDI synoptic maps to cross-calibrate the other 3 axial dipole moments calculated based on MWO, WSO and HMI synoptic maps.  Figure 9 is the time evolution of the calibrated axial dipole moment after 13-month running averages in solid curves overplotted with sunspot number in dashed curve from 1974 onwards. We derive the time series of the averaged axial dipole moment D(t) in the following way. If there are more than one observations, we do the average over the different values. During the early time, we just take the values from MWO observations since only MWO observation is available. The 13-month running average of D(t) is regarded as the homogeneous axial dipole moment dataset to derive the correlation between the dipole moment at cycle minimum D n min and the following cycle strength S n+1 . To further reduce the random noise, we use the averaged value of 7 CRs around each minimum to get D n min . The cycle strength S n+1 is the average of 7 months around each maximum of the smoothed sunspot data. The correlation between the homogenous dipole moment at cycle minima D n min and the following cycle strength S n+1 is shown in Figure 10, which indicates S n+1 = 58.7 * D n min .
Although there are only 4 points, the correlation coefficient is r = 0.99 with the corresponding confidence level p = 0.045. The sudden jump of the WSO observations during 2016-2017 in Figure 9 is due to the reduced WSO polarization sensitivity 7 . This time period has no effect on the correlation.
6. Predictability of the subsequent cycle

Prediction of the subsequent cycle at different phases of a cycle
We can predict the possible sunspot emergence of an ongoing cycle by random realizations of the features of sunspot emergence given in Section 2 a few years after the start of a cycle. With the SFT simulations and the observed synoptic magnetograms as the initial state of the magnetic field, the possible large-scale field evolution over surface can be obtained. Based on the correlation between the dipole moment at cycle minimum and the following cycle strength constrained in Section 5 combined with the features of sunspot emergence given in Section 2, the predictability of the shape of the subsequent cycle can be given. Solar cycles also vary in length. Strong cycles tend to be shorter and vice versa, but with a large scatter (Solanki et al. 2002;Du 2006;Vaquero & Trigo 2008). Some of this is due to the overlapping of sunspots from adjacent cycles. In this paper we ignore the cycle length variation and assume that all the cycles have an 11-year cycle period. If we know the maximum of a cycle S n , we may derive the profile of the cycle based on the simplified formula of HRW94, which is just a function of a relevant to cycle amplitude S n . The parameter t 0 is not required. The correction function △r(t) is still required for the improvement of the later phase. We only predict the systematic part of the smoothed sunspot number F sm S (t) for the prediction of the subsequent cycle, since the standard deviation of the prediction is usually larger than the intrinsic short term variability. The uncertainty results from the uncertainty of the axial dipole moment at the end of the cycle.
The formula to describe the subsequent cycle is where and △r(t) corresponds to the formula of Eq.(11). The fits of cycles 12-24 using Eq.(9) indicate that the relationship between S n and a is S n = 9072.8a 0.706 .
The parameters b and c are the same as that in Eq.(1). The deviation of the maximum sunspot number given by Eq.(24) from S n is within 5. Figure 11 shows some examples illustrating the predictability of the subsequent cycle, starting from different phases of a cycle. The timings for the predictions are the same as the ones in Figure 7. From Figure 7, we can derive D n + 2σ D n , D n + σ D n , D n , D n − σ D n , and D n −2σ D n at the cycle n minimum. We derive the corresponding amplitude of the subsequent cycle strength, S n+1 +2σ S n+1 , S n+1 +σ S n+1 , S n+1 , S n+1 −σ S n+1 , and S n+1 −2σ S n+1 using Eq.(22). Still the σ S n+1 and 2σ S n+1 ranges correspond to the probabilities of 68% and 95.4% for the actual sunspot number to be within the error ranges. Based on Eqs.(23), (24), and (25), the profiles of the smoothed sunspot number in the subsequent cycle n + 1 can be obtained. We see that the observed sunspot number in cycle 24 is within σ to 2σ range of the predicted result whatever the prediction is at the early time, i.e., the 4th year or the later time, i.e., the 8th year. The error range decreases with time. The deviation between the observation and the prediction during the decline phase is expected since the prediction given here is only the systematic part. As shown in Table 1, the systematic part sometimes shows a deviation from the observed sunspot number evolution (but the observations still are within the given error range). The updated prediction of the cycle based on a later magnetogram and the random flux source is presented in the following subsection.

Updated predictability of cycle 25
Here we take the synoptic magnetogram, CR2198 (December 2nd -December 31st, 2017), as the initial condition to investigate the large-scale field evolution over the solar surface during the rest of cycle 24 and the possible profiles of cycle 25. A possible (one from our Monte-Carlo ensemble) monthly sunspot emergence pattern during the rest of the ongoing cycle 24 is shown in Figure 12. The systematic part of the smoothed sunspot number is shown in the solid green curve. It is higher than the current observation. But the deviations are within 2σ range of the fluctuation. The thin green curve is one realization of the smoothed sunspot number. The corresponding time evolution of the latitudinal location of the sunspot groups (observed in black and predicted in red) is presented in the upper left panel of Figure 13. The detailed description of Figure 13 which is the predicted large-scale field evolution during the rest of cycle 24 is as follows.
The upper right panel of Figure 13 shows the time evolution of the axial dipole moment from MDI and HMI observations in blue and in black before the end of 2017 and the predicted one after 2017. The mean value will increase from 1.66 G at 2018.0 to 2.12G with 2σ range of 0.55G at 2020. We also calculated the averaged polar field over ±60 • to ±75 • latitudes to understand the different evolutions of the two hemispheres. The northern polar field is in dashed curve and the southern polar field is in solid curve. The green curves show the expected values and the dark and light red shades show the 1σ and 2σ range of the predicted polar field. During the first one year, the error range is very small. The polar field is determined by the initial conditions. The northern and the southern polar fields are almost in balance with values of 3.17G and -2.86G. The southern polar field will keep almost flat with the values varying to -2.94G, and the northern polar field will have a large increase to 4.56G at the end of 2018. The averaged northern polar fields keep increasing and the averaged southern polar fields keep stable until the end of the cycle. The mean values by then are -2.89G and 5.40G with 2σ range of 1.23G, respectively. The reason for the increasing northern polar field and the stable southern polar field can be explained by the lower right panel of Figure 13, which is the longitudinal averaged surface field evolution combined with the HMI observations (before the vertical line) and the simulation (after the vertical line) corresponding to the random realization in upper left panel. There are strong poleward positive plumes starting from second half of 2017. The positive plume in the northern hemisphere further increased the positive polar field. The positive plume in the southern hemisphere prevented the increase of southern negative polar field. Hence it keeps almost stable. The two notable plumes mainly result from two big ARs, i.e., AR12674 and AR12673 (Yang et al. 2017;Yan et al. 2018). They occurred on the solar disc with large tilts around September 5th, 2017 locating at northern and southern hemisphere respectively. According to Jiang et al. (2014aJiang et al. ( , 2015, the bigger bipolar regions emerging at low latitudes have significant effects on the large-scale field evolution. Their effects on the solar cycle will be studied in detail in a separated study.
The possible evolution of the smoothed sunspot number based on the axial dipole moment at the end of cycle 24 is given in Figure 12. The expected maximum amplitude of cycle 25 is 125, which is about 10% higher than current cycle 24. The 2σ range is 32, which means that the possibility of the amplitude of cycle 25 above 93 is 95.4%. Cycle 25 most probably is a normal cycle, rather than the Maunder minimum period. The possible profiles of cycle 25 denoted by the shaded region obey the empirical relations given in Section 2.
In order to avoid the effects due to the unexpected problems from synoptic magnetograms, we also repeated the predictions based on synoptic magnetograms, CR2196 and CR2197. The results are similar.

Summary and discussion
In this paper, we have developed a scheme to investigate the predictability of the solar cycle over one cycle. The scheme includes three steps. Firstly, empirical properties of the solar cycle are used to predict the possible sunspot emergence for an ongoing cycle. Then the SFT model is adopted to predict the possible large-scale field evolution over the surface, including the polar field at the end of the cycle. Finally, the correlation between the polar field and the subsequent cycle strength and empirical properties of the sunspot emergence are applied to get the possible profiles of the subsequent cycle. The scheme is verified by past cycles and is applied to predict the possible profiles of cycle 25. The results show the cycle 25 strength of 125 ± 32 (2σ uncertainty range), which is about 10% stronger than cycle 24 based on the mean value.
Comparing with the existing methods of the solar cycle prediction, the main progress of the current scheme is as follows. Firstly, the prediction scope is extended over one cycle. Secondly, the profiles of solar cycles during the whole prediction time period, rather than only the cycle amplitude, can be given. Thirdly, the uncertainty due to the randomness of the sunspot emergence and the data assimilated to the model during the whole prediction time period is given. Fourthly, not only the monthly and smoothed sunspot number but also the emergence properties, including the size, location, and tilt of sunspot groups, can be given. Fifthly, although the BL dynamo model is not directly used in the model, the empirical properties we used have the solid dynamo origin. These five properties distinguish the scheme from all the current existing attempts of the solar cycle predictions.
For one source of prediction uncertainties, the imperfect measurement of the initial magnetogram, we only estimate the error from the net flux density. The error due to other possible problems, like the center-to-limb correction of the magnetic saturation is not dealt with in this paper. Potential problems of the initial magnetograms can degrade the predictive ability. We have explicitly assumed that the scatter in the tilt angle is random in nature, consistent with the idea that it is due to turbulent buffeting by the convective motions. If there is a significant deterministic chaotic component in the dynamo process, which we currently do not see evidences for, the model can be improved. Conversely if there is some chaotic or random process operating on long timescales and which exceeds the randomness introduced by the scatter in the tilt angles, then our prediction can fail. Furthermore, we assumed all cycles have 11-year cycle period, and ignored the cycle overlap. So there remains room for improving the scheme, although these improvements are likely to be small compared to the uncertainties due to the tilt-angle scatter which we explicitly dealt with in this paper.
We are grateful to Robert Cameron, Sami Solanki, and  (1). The green solid curve is the predicted value using the scheme set out in this paper. The dark/light red shading gives the ±σ/±2σ variations of the predicted smoothed sunspot number. The three red curves below and above the mean values, respectively, show the boundaries of the ±σ, ±2σ and ±3σ variations of the prediction. The second column is similar except it shows unsmoothed data and shows both the expectation value (thick green) as well as an example of one realization (thin green curve). The third column shows the observed the sunspot time-latitudinal distribution, i.e., butterfly diagrams in black, and one realization from the Monte-Carlo ensemble in green.