# Predicting Time-to-Failure of Plasma Etching Equipment using Machine
Learning

Anahid Jalali,
Clemens Heistracher,
Alexander Schindler
and Bernhard Haslhofer  Austrian Institute of Technology,
Vienna, Austria
Email<EMAIL_ADDRESS>Tanja Nemeth,
Robert Glawar
and Wilfried Sihn  Fraunhofer Austria,
Vienna, Austria
Email<EMAIL_ADDRESS>Peter De Boer  Infineon Technologies
Austria AG,
Villach, Austria
Email<EMAIL_ADDRESS>
###### Abstract

Predicting unscheduled breakdowns of plasma etching equipment can reduce
maintenance costs and production losses in the semiconductor industry.
However, plasma etching is a complex procedure and it is hard to capture all
relevant equipment properties and behaviors in a single physical model.
Machine learning offers an alternative for predicting upcoming machine
failures based on relevant data points. In this paper, we describe three
different machine learning tasks that can be used for that purpose: (i)
predicting Time-To-Failure (TTF), (ii) predicting health state, and (iii)
predicting TTF intervals of an equipment. Our results show that trained
machine learning models can outperform benchmarks resembling human judgments
in all three tasks. This suggest that machine learning offers a viable
alternative to currently deployed plasma etching equipment maintenance
strategies and decision making processes.

##  I Introduction

Plasma etching is a key procedure in semiconductor wafer fabrication and can,
in case of equipment failures or breakdowns, lead to significant production
losses. Therefore, maintenance of plasma etching equipment, which aims at
minimizing unplanned breakdowns, has become a crucial task. Recently, there
has been a clear shift from re-active maintenance planning strategies such as
_Run-to-Failure_ or _Scheduled Maintenance_ to a more pro-active strategy,
which is called _Predictive Maintenance (PdM)_. The goal of that strategy is
to monitor the health state of an equipment and to predict upcoming failures
by estimating the _Time-to-Failure (TTF)_ before the next breakdown [1].

Known predictive maintenance approaches can roughly be categorized into
_model-based_ and _data-driven_ methods [2]. Model-based methods rely on
domain expertise and knowledge about the physical model of a system in order
to predict its degrading behavior. Data-driven approaches, on the other hand,
are used when it is not possible to draw a complete picture of a system’s
physical properties and behaviors. They usually employ machine Learning
techniques to model and detect changes in machine behavior. Their
effectiveness heavily depends on so called _Health Indicators_ [3], which are
quantitative features that were extracted from available sensor, product
quality, and production process data. A selection of relevant features is then
used to train a model that describes a machine’s health degradation to
eventually estimate its remaining lifetime.

A number of studies have already focused on prediction tasks in the plasma
etching context: Cheng et al. [4] developed a fault detection and isolation
system for plasma etching process chambers. Lou et al. [5] used information
about a chamber’s contamination and employed neural networks for predicting
the degradation in a semiconductor manufacturing processes. Puggini and
McLoone [6] applied _Extreme Learning Machines_ to predict the etch rate of
each wafer. Munirathinam et al. [7] predicted a machine’s state for
maintenance scheduling by focusing on product quality parameters. In summary,
existing work focuses on a specific fault type or specific physical properties
such as the thickness of the walls caused by particle contamination.

However, none of those approaches focus on predicting TTF of an entire plasma-
etching chamber without directly measuring its health status. Therefore, we
studied three different TTF prediction approaches in the context of plasma
etching equipment and can summarize our contributions as follows:

  1. 1.

_Task 1_ : We modeled TTF prediction as a _regression task_ and found that a
simple Linear Regression model can be trained to predict the TTF trend and
outperform a comparable benchmark resembling human judgment.

  2. 2.

_Task 2_ : We demonstrate that Task 1 can be transformed into a more effective
health state prediction task (_regression_) by converting the TTF target
variable.

  3. 3.

_Task 3_ : We modeled TTF prediction as a _classification task_ in order to
predict whether breakdowns will occur within defined intervals (e.g., upcoming
0-8h, 8-16h, etc.). We show that standard machine learning algorithms can
outperform comparable benchmarks.

All three tasks showed that alarm data, which is defined by the equipment
manufacturer, and sensor data limit violations, which are defined by experts
operating those equipment, provide the most informative features. We also
found that prediction effectiveness is higher in 50-200h intervals than
shortly before breakdowns (0-50h).

In the following, we briefly introduce related background information on
plasma etching and data-driven predictive maintenance approaches (Section II).
Then, in Section III, we present key findings of our exploratory analysis
before describing our experimental setup in Section IV and our results in
Section V.

##  II Background

###  II-A Plasma Etching

Semiconductor device production aims to create structures with specific
material properties on the surface of a silicon wafer. This can be achieved by
selectively adding or removing material and locally changing the chemical
structure of the wafer material (e.g., doting and oxidation).

Today, plasma etching is widely chosen for material removal as it provides the
high precision level required for the efficiency of modern devices. Prior to
the etching stage, areas of the wafer that should not be etched are masked.
Next, the wafer is placed in a low-pressure chamber where a plasma is ignited.
Ions in this plasma are accelerated towards the wafer and as they hit the
surface, they either remove material by mechanical impact, or they form a
chemical reaction with the material. The ions in the plasma must be chosen to
mostly react with the substrate and not the mask. Additionally, reaction
products should be volatile. Otherwise, they might cause deposits on the
wafer. To ensure a reliable etching process, control pressures, gas
temperatures, gas compositions, voltages, wafer cooling and the plasma
composition must be controlled.

Modern plasma etching equipment operate fully automated and process batches of
wafers with predefined recipes, where each recipe is a sequence of predefined
process parameters. In general, each equipment has a loading mechanism with
vacuum locks and several process chambers. Figure 1 illustrates the central
functional components of a modern multi-chamber plasma etching equipment.

Figure 1: Schematic illustration of chambers and loading mechanism of a modern
plasma etching equipment

###  II-B Data-Driven Predictive Maintenance Approaches

Most data-driven predictive maintenance approaches in the literature are based
on statistical, probabilistic or machine learning methods and can be
categorized into two groups [8]: the first group are prognostic models that
directly observe production state processes. They either apply regression-
based models [9], or Markovian-based models [10]. Furthermore, Wang et. al.
[11] proposed a generic model for probabilistic health condition estimation
and tested it on two scenarios: (i) electric cooling fans and (ii) an engine
dynamic simulation.

Machine learning algorithms such as Support Vector Machines (SVM) [12],
Support Vector Regression (SVR) [13] and binary logistic regression [14] form
the second group and have also proved to be an effective solution for
estimating TTF and other machine health descriptor variables. Deep Neural
Networks (DNNs) have recently shown strong performance on variety of complex
applications such as speech recognition [15], image classification [16] and
acoustic sound classification and detection [17]. However, as explainability
of trained models is a key requirement in the semiconductor industry, other
algorithms are typically preferred.

###  II-C Prediction Maintenance in the Semiconductor Industry

Outside the specific context of plasma etching, a number of data-driven
predictive maintenance approaches have been proposed for monitoring machine
health degradation in the semiconductor industry.

Munirathinam et al. [7] constructed a decision model for a semiconductor
fabrication plant by applying a variety of standard machine learning
algorithms (e.g., KNNs, SVM). In order to reduce data dimensionality, they
applied reduction approaches such as Principal Component Analysis (PCA),
Variable Importance Analysis (VIA) and Chi Square statistics. However,
although their models can predict whether a product passes a final quality
inspection, they were not designed to predict possible equipment breakdowns.

Lou et al. [5], proposed a two-step maintenance framework for degradation
prediction on a real case study in semiconductor manufacturing industries and
achieved 74.1% accuracy in predicting a machine’s health degradation. Their
work is divided into three stages: first, they adopted a back-propagation
Neural network (BPNN) to forecast the machine’s health. Second, as a backup
for the failed cases of the first stage, they have employed Restricted
Boltzmann Machine (RBM) and Deep Belief Networks (DBN), which have strong
inductive learning abilities. Finally, using multiple regression forecasting,
they checked the prediction accuracy in both stages. However, the model is
specialized to predict only the contamination of the process chamber. Failures
unrelated to the contamination are not considered in the model.

Susto and McLoone [18] adopted an SVM classifier to predict mechanical faults
in semiconductor manufacturing. This type of failure is caused by usage and
stress of the equipment parts (e.g., filament breaks). In their work, they
classify a machine’s run as faulty when the SVM’s decision boundary falls
below a threshold. They specified this threshold by considering two factors:
unexpected breakdowns and unexploited lifetime. Their approach showed
robustness to cost changes associated with unexpected breaks and inefficient
lifetime of an equipment part. However, they do not discuss breakdowns of
plasma etching equipment. Additionally, they focus on specific parts (e.g.,
breakdown of a tungsten filament) and not on a whole plasma etching process
chamber.

##  III Exploratory Analysis

Before having focused on prediction model building, we first gathered data
from several nearly identically constructed plasma etching equipment and
process chambers. Then we computed TTF for each chamber, which represents the
target variable for all our prediction tasks and the ground-truth for
validating trained models. We also conducted an exploratory analysis of
available data points in order to identify possible correlations and to reduce
the dimensionality of our data. In order to minimize competitive intelligence
risks, all figures in this section are schematic and illustrate our findings
without exposing details of the underlying manufacturing process.

###  III-A Dataset Characteristics

Our dataset encompasses data recorded over a 6-month period and has been drawn
from the following sources:

  1. 1.

_Automatic Process Control (APC)_ : contains 492 distinct sensor data
recordings for a single wafer from the underlying plasma etching process
control system. This includes statistics of measurements, for instance
relevant gas flows and voltages, as well as process-related information such
as the time needed for a single etching step and the applied recipes.

  2. 2.

_APC Limit Violations_ : limits are upper- and lower-limit thresholds, which
were defined by domain experts for selected APC process control parameters.
Limit violations are categorized by severity (_error_ vs. _information_) and
can trigger actions ranging from automated equipment shutdowns to sending
informational emails to domain experts. In total, our dataset contains
recordings of 58 different limit violations.

  3. 3.

_APC Alarms_ : alarms are defined by the machine manufacturer and are
categorized into five classes; _warning_ , _information_ , _critical_ ,
_errors_ and _other_. Alarms can, analogous to limit violations, also trigger
a number of possible actions, including equipment shutdowns. In total, our
dataset contained recordings of 603 distinct alarms per process chamber.

  4. 4.

_Real Time Clock (RTC)_ : is a system that records state changes (e.g.,
standby, productive, breakdown) of plasma etching equipment and their
corresponding parts. We retrieved those state changes for all equipment over
the entire observation period.

  5. 5.

_Voltage Dips_ : describe the voltage reduction in the power supplies. The
occurrence of this feature is very sparse. However, domain experts assume that
the voltage dips are problematic and can cause equipment’s breakdown.

###  III-B TTF Computation

In order to build a ground truth for subsequent prediction tasks, we reverse
engineered the TTF for each process chamber by joining a machine’s state from
RTC with process durations extracted from APC. This metric yields a descending
counter of productive hours before the next recorded breakdown. Figure 2
illustrates the zigzag shape of an arbitrary TTF with four breakdowns.

Figure 2: Time-to-Failure (y-axis) of a single chamber over the observation
period (x-axis). Red marks indicate recorded breakdowns.

###  III-C Analysis of APC Limit Violations

Limit violations can, as discussed before, trigger further actions. An example
of a severe violation leading to an immediate shutdown are helium flows above
a certain threshold, which are typically caused by dust particles in a
chamber. Other violations are only recorded as warnings and do not cause
immediate actions.

To understand the impact of limit violations on future process chamber
breakdowns, we analyzed historical limit violation recordings and computed the
median TTF per defined violation. Figure 3 shows a selection of limit
violation occurrences (scatter plot) and computed median TTFs (bar plot) in
relation to TTF. The single box plots are arranged by their median TTF from
lowest to highest, which shows that some violations tend to occur closer to a
machine’s failure (TTF equals 0) than others. Limit violations with zero TTF
are bound to immediate shutdown actions. Intuitively, limit violations that
have a median TTF slightly above zero are the most interesting parameters as
they indicate upcoming breakdowns but do not immediately cause them.

Figure 3: Limit violation occurrences (scatter plot) in relation to TTF. Box
plots show medians, colors indicate severity (orange: error; green:
information).

###  III-D Analysis of APC Alarms

A plasma etching equipment can raise a number of alarms while processing a
wafer. Figure 4 illustrates APC alarm occurrences in relation to a chamber’s
TTF with each alarm colored by category. Analogous to the previous
illustration of limit violations, it arranges alarms in ascending order by
their median TTF. Thus, a lower median TTF indicates that an alarm was often
raised near an upcoming breakdown. Again, we see that some alarms are more and
some are less relevant for our prediction tasks and can assume that the most
informative alarms show a median TTF slightly above zero. We also see that
alarm severity does not necessarily correspond with the alarm categories
defined by the equipment manufacturer.

Figure 4: APC alarm occurrences (scatter plot) in relation to TTF. Box plots
show medians, colors indicate alarm categories (red: critical; yellow: error;
brown: warning; green: information; blue: other).

###  III-E Analysis of APC Sensor Data Correlations

We computed correlations between APC sensor data points in order to identify
linear relationships between parameters. This allowed us to reduce the
dimensionality of our dataset by 87.2% by discarding parameters that correlate
strongly with others and therefore add little information.

While conducting a principal component analysis, we also found that many
process control data recordings strongly depend on the recipe used in an
etching process run. This can also be observed when correlating two APC
process variables with each other, as shown in Figure 5. It shows that
correlations are clustered by some external factor, which in this case, is the
configured and applied recipe.

Figure 5: Illustration of APC process variable correlations. Arrows show
positive (green) or negative (red) correlations.

##  IV Experimental Setup

In our experimental setup, we followed the typical data science work-flow:
first, we cleansed and normalized our dataset, which had been aggregated from
a number of heterogeneous sources. Then, in collaboration with domain experts
and based on the observations made in the previous exploratory analysis phase,
we engineered a number of features that were potentially useful to the
prediction model. Next, in order to evaluate the effectiveness of our models,
we built a benchmark resembling human decision making. Finally, we built
predictive models for three different TTF prediction and interpretation
methods and evaluated their effectiveness against the respective benchmarks.
In the upcoming subsections we describe each of those steps in more detail.

###  IV-A Data Cleansing and Normalization

In this step, we removed non-relevant data points such as system-internal
identifiers, user names or null-values that do no contribute valuable
information to our models.

Given the recipe dependencies of sensor data correlations described in Section
III-E, we standardized all sensor features by the sample mean and standard
deviation of the corresponding recipe. This makes the values of each feature
in the data have zero-mean and unit-variance and allows comparison across
recipes. We also eliminated recipes used by non-productive runs (e.g.,
cleaning, experiments).

Filtering by relevant recipes reduced the number of missing sensor parameter
values to 53%. We set remaining missing values to zero when we could safely
assume that a certain process parameter was not used. Having a NaN in a helium
flow, for instance, means that no helium was used in a process step. This also
holds for voltage, intensity, and power values.

Finally, we computed the TTF target variable for each chamber, assigned unique
identifiers to each segment between breakdowns and normalized that variable to
zero mean and unit variance.

###  IV-B Feature Engineering

Given the cleansed and normalized datasets, we then engineered a number of
feature sets, which fall into three main groups: features derived from APC
sensor data (_APC_), alarm features (_AL_), and features derived from limit
violations (_LV_).

For alarms we constructed a cumulative feature that sums the number of alarm
occurrences (_counter alarms_) in each productive interval. Additionally, we
weighted each occurrence with a penalty value that captures the severity of an
alarm. The penalty for each alarm is computed by inversing its median TTF
($AP_{x}=1/median(TTF)_{x}$) and then added to an alarm occurrence. In analogy
to alarms, we created a similar feature (_counter violations_) for limit
violations. Figure 6 illustrates those weighted alarm and limit violation
occurrence features. In sub-figure (a) we can clearly see an increase of
violations during the run-time of a chamber and spikes shortly before
breakdowns. Sub-figure (b) shows similar behavior, which can be modeled by
computing the gradient of those features.

Figure 6: Limit violation and alarm counter features compared to Time-to-
Failure of one chamber. Sub-figure (a) depicts the cumulative weighted sum of
the limit violations over TTFs of multiple productive intervals. Sub-figure
(b) shows the cumulative weighted sum of the alarms within the same interval.

Table I summarizes the feature sets used in our experiments. We selected a
number of feature combinations based on feedback from domain experts who
defined limit violations and thereby indirectly pointed our importance of a
features based on past observations. In the following, we denote, for
readability purposes, our feature set combinations as follows:
$FS_{1}=APC_{V}+APC_{R}$
$FS_{2}=APC_{V}+LV_{P}$
$FS_{3}=APC_{V}+APC_{R}+LV_{P}+AL_{P}$
$FS_{4}=APC_{V}+APC_{R}+LV_{P}+AL_{P}+VoltageDips$
$FS_{5}=APC_{V}+LV_{P}+AL_{P}$
$FS_{6}=LV_{P}+AL_{P}$
$FS_{7}=AL_{P}$

Table I: Feature sets used in the model Feature Set |  Description
---|---
$APC_{V}$ |  Subset of APC process data containing only features with defined limit violations.
$APC_{R}$ |  Contains information on the recipe mix of x runs (e.g., gradient sum, gradient max).
$LV_{P}$ |  Engineered cumulative sum of limit violations with penalty (gradient sum, gradient max).
$AL_{P}$ |  Engineered Cumulative sum of alarms with penalty (gradient sum, gradient max).

###  IV-C Model Building

We have chosen a number of machine learning algorithms that are known (cf.
Section II) for their robustness among applications of classification and
regression for health monitoring and prognosis. For regression tasks we used
Linear Regression (LR), Support Vector Machines (SVM), Decision Trees (DT),
Random Forest (RF), and Multilayer Perceptrons (MLP). For the classification
task, we used the same algorithms except LR. Additionally, we also chose
Gradient Boosting Classifier (GBC), SVM with Stochastic Gradient Descent (SGD)
and K-Nearest Neighbours (KNN). For all our models we used implementations of
Scikit-learn [19].

The dataset is split in a way that ensures that recordings from the same
productive interval (the time between two breakdowns) are always assigned to
the same fold. This is necessary because data in a productive interval is
autocorrelated and would reveal information on the next breakdown, if it was
used for both training and testing.

To ensure the reliability of our results, we applied a 4-fold cross validation
and iteratively trained our models on three folds of data and tested it on one
other fold not included in the training data set.

###  IV-D Benchmark Definition

We evaluated the effectiveness of our models in comparison to three benchmarks
resembling human decision making: for the _naïve benchmark_ (B1), we assume
that the TTF is always constant and simply take the value of the mean TTF over
all runs. The _visionary benchmark (B2)_ serves as an idealized reference and
considers an adjusted mean TTF for each productive segment, which is
unrealistic because it requires knowledge on future breakdowns. The _realistic
benchmark (B3)_ is closer to reality and assumes that a breakdown occurs after
$x$ productive hours where $x$ denotes mean productive time between breakdowns
based on historic data. Figure 7 visualizes all three benchmarks.

Figure 7: Illustration of three benchmarks used for model evaluation. The red
line represents the first benchmark (mean TTFs over entire observation
period); the green line the second benchmark (adjusted mean TTF for each
productive interval) and the blue line represents the third benchmark (based
on historic means).

###  IV-E Model Evaluation

We evaluate our models using standard metrics used in machine learning. For
regression tasks we compute the root mean squared error (RMSE) of the
predicted and the actual value. Then we compute the relative difference to our
realistic benchmark B3 as follows: $x_{rel}=(x-B3)/B3$ with $x$ denoting the
RMSE of the prediction and $B3$ denoting the RMSE of B3. Thus, a negative $x$
characterizes an improvement to the realistic benchmark, whereas positive
values denote a prediction that is less effective than the benchmark.

Although the RSME is suitable metric for finding the prediction that is on
average closest to the actual TTF curve, it does not necessarily evaluate a
model’s practical applicability. A constant prediction of the mean value, for
instance, might result in a sound RMSE, but is less useful for maintenance
purposes as no degradations resembling the decay of TTF is shown. Therefore,
we define the best useful model as the model having the lowest RMSE and
showing degradation upon visual inspection of predictions.

Furthermore, we evaluate the effectiveness of the classification task with the
standard precision (P), recall (R) and F1-score (a trade off between precision
and recall). For a more detailed explanation of machine learning evaluation
metrics, we refer the reader to related literature such as Powers [20].

##  V Results

In the following, we present the results for three types of models we built
for predicting the TTF for plasma etching equipment.

###  V-A TTF Prediction (Regression)

The main goal of this task is to predict the remaining time to a machine’s
failure. In Section III-B we explained the calculation of the zigzag shaped
TTF curve, which is the target variable of this task. To find the most
informative set (or sets of features), we used different feature set
combinations as input of our regression models with the objective of
minimizing RMSE. Afterwards, we compared the results of each experiment with
our previously defined benchmarks.

Table II: Results TTF prediction with $B3_{RMSE}=223.96.$

Features | B1 | B2 | B3 | LR | SVM | MLP | Tree | RF
---|---|---|---|---|---|---|---|---
FS1 | -0.29 | -0.49 | 0 | -0.36 | -0.40 | -0.34 | -0.42 | -0.33
FS2 | -0.31 | -0.40 | -0.3 | -0.18 | -0.29
FS3 | -0.26 | -0.40 | -0.27 | -0.08 | -0.31
FS4 | -0.24 | -0.40 | -0.35 | -0.07 | -0.31
FS5 | -0.26 | -0.40 | -0.29 | -0.05 | -0.31
FS6 | -0.28 | -0.42 | -0.31 | -0.08 | -0.28
FS7 | -0.28 | -0.42 | -0.28 | -0.05 | -0.28

Table II presents the result of our experiments for one plasma etching chamber
after cleaning erroneous breakdown recordings. It shows that Support Vector
Machines (SVM) outperformed other models. However, when inspecting predictions
visually, we observed that Linear Regression (LR) had the lowest RMSE while
showing a degradation in its prediction. Trained with FS1 it outperformed B3
by 36%, showed 7% improvement over B1, and was 13% worse than the second,
visionary benchmark B2.

In summary, our first experiments on a single process chamber showed that a
trained regression model can predict the degrading trend of a TTF curve,
however, with relatively high RMSE. When analyzing the errors, we found that
some sequences of relatively short consecutive breakdowns were, according to
domain experts’ opinions, erroneous recordings caused by equipment starts and
almost immediate shutdowns within maintenance operations.

###  V-B Health State Prediction (Regression)

The goal of this task is to predict a machine’s health status. For this
purpose, we transformed the target value (TTFs) into a range between 0 and 1
(see Figure 8), where a health state of 1 is considered as healthy and 0 as
failure (breakdown). All other pre-processing steps were the same as in the
first task.

Figure 8: Illustration of time to failures in hours and the corresponding
machine’s health state.

Table III shows that both Multilayer Perceptron (MLP) and Linear Regression
(LR) models trained on all feature set combinations can outperform all
benchmarks and visual inspection of predictions indicates degradation in both
model families. Figure 9 illustrates health state prediction results of a
trained Linear Regression model from one cross-validation fold.

Table III: Results of Machine Health prediction with $B3_{RMSE}=0.36.$

Features | B1 | B2 | B3 | LR | SVM | MLP | Tree | RF
---|---|---|---|---|---|---|---|---
FS1 | -0.19 | -0.19 | 0 | -0.22 | -0.19 | -0.22 | -0.03 | -0.19
FS2 | -0.17 | -0.19 | -0.28 | -0.03 | -0.19
FS3 | -0.11 | -0.19 | -0.25 | 0.11 | -0.17
FS4 | -0.11 | -0.19 | -0.25 | 0.11 | -0.19
FS5 | -0.11 | -0.19 | -0.25 | 0.11 | -0.17
FS6 | -0.14 | -0.22 | -0.25 | 0.08 | -0.19
FS7 | -0.14 | -0.22 | -0.25 | 0.06 | -0.19

Figure 9: Machine’s health state prediction result of a Linear Regression
model from one cross-validation fold.

###  V-C TTF Prediction (Classification)

In this task, we consider TTF prediction as a binary classification task that
predicts whether a machine will face a breakdown within a predefined time
interval of $0$ and $x$, where
$x\in\\{8h,16h,24h,48h,72h,96h,120h,144h,168h,336h\\}$. For a given interval,
a run is labeled true when its TTF lies within that interval, otherwise the
run is labeled as false.

Figure 10 shows the most effective (highest F1-score) model and feature set
combination for each interval. We can clearly observe that trained models can
outperform the realistic benchmark (B3) in all intervals. Furthermore, we see
that prediction models trained for shorter intervals (e.g., 0-8h) have a lower
F1-score than those for longer intervals. This reflects our findings from the
previous regression tasks, which also yielded higher RMSE shortly before
breakdowns.

Figure 10: F1-score values in multiple experiments. Each selected feature set
is defined with different color and each model is determined with a shape.

A summary of our experimental classification results is presented in Table IV
which compares the effectiveness of each interval-specific trained model to
the corresponding interval-specific benchmark. As in Task 1 and 2 we can
observe the importance of alarms and limit violations for shorter predictions
as they are included in feature sets FS4, FS5 and FS7 but not in FS1.

Table IV: Summary of best classification results per interval

Interval | Realistic Benchmark | Best Results
---|---|---
P | R | F1 | Features | Model | P | R | F1
0-8 h | 0.09 | 0.35 | 0.15 | FS4 | Tree | 0.25 | 0.36 | 0.29
0-16 h | 0.17 | 0.35 | 0.23 | FS5 | Tree | 0.36 | 0.53 | 0.42
0-24 h | 0.24 | 0.37 | 0.29 | FS7 | SVM | 0.39 | 0.75 | 0.51
2 days | 0.41 | 0.39 | 0.4 | FS7 | SVM | 0.59 | 0.76 | 0.66
3 days | 0.54 | 0.41 | 0.46 | FS1 | GBC | 0.68 | 0.91 | 0.78
4 days | 0.64 | 0.43 | 0.51 | FS1 | GBC | 0.74 | 0.96 | 0.83
5 days | 0.72 | 0.44 | 0.55 | FS1 | GBC | 0.79 | 0.98 | 0.88
6 days | 0.8 | 0.46 | 0.58 | FS1 | GBC | 0.84 | 0.99 | 0.91
7 days | 0.85 | 0.47 | 0.6 | FS1 | GBC | 0.88 | 0.99 | 0.93
14 days | 1 | 0.47 | 0.64 | FS6 | MLP | 1 | 1 | 1

##  VI Discussion

With the overall goal of predicting Time-to-Failure of plasma etching
equipment in the semiconductor industry, we experimented with three different
machine learning approaches. In all three tasks we were able to outperform
comparable benchmarks resembling human breakdown predictions based on historic
mean TTF observations. We highlighted the importance of precise breakdown
recordings and subsequent TTF computations for building accurate and effective
prediction models. We also emphasized the importance of including
manufacturer-defined alarms and expert-defined limit violations, which both
capture a high degree of domain knowledge.

One limitation of our approach lies in the relatively high manual effort
required for data cleansing, normalization and feature engineering. In order
to support genericity and applicability of our overall method on other
equipment types with similar data sources, we strongly focused on building
generic features that are derived from standard data sources (sensor data,
alarms, limit violations) and therefore applicable across equipment. However,
we believe that further boosts in effectiveness are possible by modeling
equipment-specific behaviors. A possible strategy is to investigate Deep
Learning methods, which support automated feature learning but still have the
drawback of being non-transparent and hardly explainable. Activation-based or
gradient-based methods for feature isolation are possible solutions for those
problems.

Finally, we can identify two main orthogonal challenges that need to be
tackled when implementing data-driven predictive maintenance in a production
setting. First, data quality, especially recordings of breakdowns, is a key
prerequisite for building effective prediction models. The most elaborate
machine learning method won’t provide more precision if the target variable
(TTF) is flawed. Second, data-driven breakdown predictions must be embedded
into existing maintenance management workflows and operations. This typically
requires formation of dedicated groups of professionals having both knowledge
of machine and process design as well as an understanding and intrinsic
interest in novel data-driven maintenance technologies.

##  VII Conclusion

In this paper, we described three different machine learning tasks that can be
used for predicting Time-to-Failure (TTF) or the health state of plasma
etching equipment in the semiconductor industry. Our results show that trained
prediction models provide acceptable effectiveness for periods exceeding
roughly 24 hours, which allows maintenance planners to react on those
predictions. We also highlighted the importance of alarms and limit
violations, which carry high degrees of domain knowledge. Since model
transparency was a key requirement, we restricted ourselves to manual feature
engineering and well-known machine learning techniques. In future, however, we
will proceed and investigate explainable Deep Learning techniques for those
prediction tasks.

##  VIII Acknowledgement

This work was partly funded by the Austrian Research Promotion Agency (FFG)
through the project DeepRUL (Project ID: 871357).

## References

  * [1] F. Ansari, R. Glawar. Knowledge Based Maintenance, Instandhaltungslogistik - Qualität und Produktivität steigern (Maintenance Logistics - Enhancing Quality and Productivity), K. Matyas (Ed.), 7th Edition, Carl Hanser, pp. 318-342, 2018.
  * [2] G. Susto,A. Schirru, S. Pampuri, S. and McLoone, A. Beghi. ”Machine learning for predictive maintenance: A multiple classifier approach”, IEEE Transactions on Industrial Informatics, vol. 3, 812–820, 11, 2015.
  * [3] L. Guo, N. Li, F. Jia, Y. Lei, J. Lin. ”A recurrent neural network based health indicator for remaining useful life prediction of bearings”, Neurocomputing, Bearing health indicator,Recurrent neural network,Related-similarity feature, pp. 98-109, 240, 2017.
  * [4] M.H. Cheng H. S. Lin, C. H. Liu, W. Y. Lee, C. H. Tsai. ”Fault detection and isolation for plasma etching using model-based approach”, IEEE: Advanced Semiconductor Manufacturing Conference and Workshop, 2003.
  * [5] M. Luo, H. C. Yan, B. Hu, J. H. Zhou,C. K. Pang. ”A data-driven two-stage maintenance framework for degradation prediction in semiconductor manufacturing industries”, Elsevier: Computers & Industrial Engineering, vol. 85, pp. 414-422, 2015.
  * [6] L. Puggini, S. McLoone. ”Extreme learning machines for virtual metrology and etch rate prediction”, IEEE: 2015 26th Irish Signals and Systems Conference (ISSC), 2015.
  * [7] Munirathinam, Sathyan and Ramadoss, Balakrishnan, ”Predictive models for equipment fault detection in the semiconductor manufacturing process”, IACSIT International Journal of Engineering and Technology, vol. 8, pp. 273-285, 2016.
  * [8] X. S. Si, W. Wang, C. H. Hu, M. Y. Chen, D. H. Zhou. ”A Wiener-process-based degradation model with a recursive filter algorithm for remaining useful life estimation”, 2012.
  * [9] D. Barraza-Barraza, V. G. Tercero-Gomez,M. G. Beruvides, J. and Limon-Robles. ”An adaptive ARX model to estimate the RUL of aluminum plates based on its crack growth”, Mechanical Systems and Signal Processing, vol. 82, pp. 519-536, 2017.
  * [10] Y. Xiang, C. R. and Cassady, E. A. Pohl. ”Optimal maintenance policies for systems subject to a Markovian operating environment”, Pergamon: Computers & Industrial Engineering, vol. 62, pp. 190-197, 2012.
  * [11] P. Wang, B. D. Youn, C. Hu. ”A generic probabilistic framework for structural health prognostics and uncertainty management”, Elsevier: Mechanical Systems and Signal Processing, vol. 28, pp. 622-637, 2012.
  * [12] M. A. Patil, P. Tagade, K. S. Hariharan, S. M. Kolake, T. Song, T. Yeo, S. Doo. ”A novel multistage Support Vector Machine based approach for Li ion battery remaining useful life estimation”, Elsevier: Applied Energy, vol. 159, pp. 285-297, 2015.
  * [13] H. Dong, X. Jin, Y. Lou, C. Wang. ”Lithium-ion battery state of health monitoring and remaining useful life prediction based on support vector regression-particle filter”, Elsevier: Journal of Power Sources, vol. 271, pp. 114-123, 2014.
  * [14] J. Phillips,E. Cripps, J. W. Lau, M. R. Hodkiewicz. ”Classifying machinery condition using oil samples and binary logistic regression”, Mechanical Systems and Signal Processing, pp. 316–325, 2015.
  * [15] D. Amodei, S. Ananthanarayanan, R. Anubhai, J. Bai, E. Battenberg, C. Case. ”Deep Speech 2 : End-to-End Speech Recognition in English and Mandarin”, International Conference on Machine Learning, vol. 48, pp. 173-182, 2016.
  * [16] E. H. Geoggrey, A. Krizhevsky, I. Sustskever. ”ImageNet Classification with Deep Convolutional Neural Networks”, Advances in Neural Information Processing Systems 25, 2012.
  * [17] A. Dang, T. H. Vu, J. C. Wang. ”Acoustic scene classification using convolutional neural networks and multi-scale multi-feature extraction”, 2018 IEEE International Conference on Consumer Electronics, ICCE 2018, pp. 1-4, 2018.
  * [18] G. A. Susto and S. McLoone and A. Schirru and S. Pampuri, A. Beghi. ”Prediction of Integral Type Failures in Semiconductor Manufacturing through Classification Methods”, IEEE. 18th Conference on Emerging Technologies & Factory Automation (ETFA), 2013.
  * [19] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg et al. ”Scikit-learn: Machine learning in Python”, Journal of machine learning researc, vol. 12, pp. 2825-2830, 2011.
  * [20] D. M. Powers. ”Evaluation: from precision, recall and F-measure to ROC, informedness, markedness and correlation”, Journal of Machine Learning Technologies, vol. 2, pp. 37-63, 2011.


11institutetext: Helmholtz-Zentrum Dresden - Rossendorf, Bautzner Landstr.
400, 01328 Dresden, Germany
11email<EMAIL_ADDRESS><EMAIL_ADDRESS>22institutetext: Eawag,
Überlandstrasse 133, 8600 Dübendorf, Switzerland
22email<EMAIL_ADDRESS>
# No evidence for absence of solar dynamo synchronization

F. Stefani 11 J. Beer 22 T. Weier 11

(Received September 15, 1996; accepted March 16, 1997)

###### Abstract

Context. The old question of whether the solar dynamo is synchronized by the
tidal forces of the orbiting planets has recently received renewed interest,
both from the viewpoint of historical data analysis and in terms of
theoretical and numerical modelling.

Aims. We aim to contribute to the solution of this longstanding puzzle by
analyzing cosmogenic radionuclide data from the last millennium.

Methods. We reconsider a recent time-series of 14C-inferred sunspot data and
compare the resulting cycle minima and maxima with the corresponding
conventional series down to 1610 A.D., enhanced by Schove’s data before that
time.

Results. We find that, despite recent claims to the contrary, the 14C-inferred
sunspot data are well compatible with a synchronized solar dynamo, exhibiting
a relatively phase-stable period of 11.07 years, which points to a
synchronizing role of the spring tides of the Venus-Earth-Jupiter system.

###### Key Words.:

Sun: activity — Sun: dynamo

##  1 Introduction

The question of whether the solar dynamo might be “clocked” by the motion of
the planets traces back to early speculations by Wolf (1859), and has popped
up sporadically ever since (de la Rue et al. (1872), Bollinger (1952), Jose
(1965), Takahashi (1968), Wood (1972), De Jager & Versteegh (2005)). Recently,
new impetus was given to the issue by the exemplification of Hung (2007),
Scafetta (2012), Wilson (2013), and Okhlopkov (2016) that the 11.07-yr spring-
tide period of the tidally dominant planets Venus, Earth and Jupiter appears
to be in a phase-stable relation with the solar cycle. This finding turned out
to be in amazing agreement with the older results of Schove’s ambitious
“spectrum of time” project (Schove 1983) to determine the solar cycle maxima
and minima for the last two and a half millennia mainly from historical aurora
borealis sightings and naked-eye sunspot observations. Furthermore, the
identified 11.07-yr periodicity is also, within the error margin, well
compatible with the phase-stable 11.04-yr cycle as inferred by Vos et al.
(2004) utilizing two different algae data-sets from the early Holocene.

The key problem with those observed correlations is how they could be
substantiated by any kind of a viable causation. While the tidal forces of the
planets can easily be ridiculed by the minuscule tidal height of the order of
1 mm, several physical mechanisms have been invoked that could possibly lead
to noticeable effects, among them the extreme sensitivity of the storage
capacity for magnetic fields in the sub-adiabatic tachocline (Abreu et al.
2012, Charbonneau 2022), or the susceptibility of intrinsic helicity
oscillations of waves or instabilities (in particular, the Tayler instability)
to tidal forces (Weber et al. 2013, Weber et al. 2015, Stefani et al. 2016,
Stefani et al. 2019, Stefani et al. 2020a, Stefani et al. 2021).

Going beyond such mainly qualitative arguments, the recent work of Horstmann
et al. (2023) has shown that even weak tidal forces such as of Jupiter might
excite (magneto)-Rossby waves with typical velocity amplitudes of up to m/s. A
concurrent 2-dimensional simulation by Klevs et al. (2023) affirmed that
tidally triggered oscillations of the tachoclinic $\alpha$ effect of the order
of dm/s would be sufficient to synchronize an otherwise conventional
$\alpha-\Omega$ dynamo model. Together with the older argument of Öpik (1972)
that the “ridiculous” 1 mm tidal height corresponds energetically to a
velocity scale of 1 m/s, those recent results exemplify that tidal forces may
entail a serious potential for solar dynamo synchronization.

Still, a hard-to-solve problem of that kind would not even appear if the solar
dynamo was not phase-stable in the first place. Two recent papers (Nataf 2022,
Weisshaar et al. 2023) have seriously put into question the empirical evidence
for phase-stability. Ignoring the strong argument in favour of phase stability
coming from the algae-date in the early Holocene (Vos 2004), both papers focus
exclusively on the series of solar cycle extrema (i.e., minima or maxima)
during the last millennium. Nataf (2022) dismissed all the meticulous efforts
of Schove (1983) by claiming his cycle reconstruction to be construed by
simple rules. Strictly presupposing that the solar cycle is not clocked, he
argued that Schove’s “nine-per-century” rule would lead to a constrained, and
therefore wrong, series of extrema. What was not considered by Nataf (2022),
though, was the possibility that the solar cycle might indeed be clocked by a
11-07-yr trigger, in which case Schove’s auxiliary “nine-per-century” rule
would do absolutely no harm to an otherwise correctly inferred series of
extrema. In this respect it is interesting to note that Schove’s data actually
point to an average Schwabe cycle of 11.07 years rather than 11.11 years which
would result from a naive application of the “nine-per century” rule (see Fig.
1 of Stefani et al. (2020a)). For further critical remarks on Nataf (2022),
see the recent comment by Scafetta (2023).

Figure 1: Series of annual (pseudo) sunspot numbers (violet lines) for the
time interval between 975 and 1895, and inferred solar cycle maxima (red
circles) and minima (green circles) according to Usoskin et al. (2021),
together with a merger of the maxima (red squares) and minima (green squares)
of Schove (1983) and NOAA
(https://www.ngdc.noaa.gov/stp/solar/solardataservices.html). Each individual
panel shows a period of 250 years, with overlapping intervals of 25 years at
either margin. Red and green thin lines indicate putative correspondences
between the respective extrema (they are usually restricted to the central
200-yr interval of the corresponding panel, except in the first panel). The
size of the red and green circles mirrors the quality flag $q$ according to
Usoskin et al. (2021) (smallest circles denote very poor quality, $q=0$,
largest circles denote highest quality, $q=5$). Red and green question marks
indicate intervals with unclear correspondences, typically at times
characterized by low quality flags. The two dashed violet lines around 1425
and 1461 indicate corrections of obvious typing errors for the maxima in Table
1 of Usoskin et al. (2021). After these two trivial corrections, uninterrupted
one-to-one matches between the extrema of both series become visible between
1140 and 1620.

A second paper that claims to have finally debunked the clocking scenario was
recently published by Weisshaar et al. (2023). Based on 14C data of Brehm et
al. (2021) for the last millennium, it uses the series of cycle extrema
inferred by Usoskin et al. (2021) to show that this series points - with a
high statistical significance - to a random walk process instead of a clocked
process.

In the present paper we will reanalyze this series of cycle minima and maxima
and compare it with another series of extrema, for which we use a combination
of the standard Schwabe cycles for the later time interval starting at 1610
with Schove’s data (Schove 1983) for earlier times. We will show that - for
the most part of the interval - the extrema of both series can uniquely be
matched one-to-another with three exceptions. The latest of those shows up
around 1840 where Usoskin’s data exhibit two shallow minima at a place where
the telescopic data show only one. Given the shallowness of these two minima,
and the relatively low quality flag of the first one, we find it legitimate to
replace this pair by only one minimum, and to cancel the corresponding maximum
between them. A second ambiguity is found amidst the Maunder minimum around
1650 where all quality flags of Usoskin’s data are relatively low. Here again
we cancel one shallow minimum. The most problematic part appears in the
interval between 1040 and 1140 where Usoskin’s quality flags are generally
quite low. Specifically, we consider three different ways of correcting the
data which we all consider at least as plausible as the original selection of
Usoskin et al. (2021).

Then, we will analyze the considered times series with view on their phase
stability. At first we show the respective Observed-minus-Calculated (O-C)
diagrams of the residuals of the instants of the minima from a theoretical
linear trend with an alleged 11.07-yr period. While the Schove/NOAA data are
concentrated around a horizontal line with only slight ($\pm$ 4 years) upward
or downward deviations, Usoskin’s data show larger deviations exactly within
the three problematic intervals discussed above. We show how these deviations
are consecutively reduced by our corrections. Thereby we arrive at a one-to-
one matching of the 14C extrema with those of the Schove/NOAA series for the
entire 9 century long interval.

In the last step we compute - for the different time series - Dicke’s ratio
between the standard deviation of the residuals and the standard deviation of
the differences between neighbouring residuals. A closely related measure,
defined by Gough (1981), was used by Weisshaar et al. (2023) and Biswas et al.
(2023) to argue in favour of a random walk process. We show here that already
the two highly plausible corrections around 1840 and 1650 lead to a dramatic
move of Dicke’s ratio towards the corresponding theoretical curve for a
clocked process. Finally we show that a good deal of the remaining deviations
from a strictly clocked curve turns out to be due to the presence of a well-
expressed Suess-de Vries cycle in the data.

The paper closes with some conclusions.

##  2 Data sets and possible corrections

In the following, we will discuss two data sets. The first one is the annular
series of (pseudo) sunspot numbers as recently inferred by Usoskin et al.
(2021) from the 14C production rate data of Brehm et al. (2021). To start
with, we simply adopt the cycle minima and maxima as derived by Usoskin et al.
(2021).

The second sequence of minima and maxima consists of the standard data from
NOAA with the starting year 1610
(https://www.ngdc.noaa.gov/stp/solar/solardataservices.html), merged with the
earlier data as published by Schove (1983) and partly corrected in Schove
(1984).

In Fig. 1, these data are presented in five 200-yr intervals, with 25-yr
overlaps at each beginning and end. The size of the symbols of Usoskin’s
minima (green open circles) and maxima (red open circles) is scaled by the so-
called “quality flag” between $q=0$ (“cycle cannot be reliable identified”)
till $q=5$ (“clear cycle in both shape and amplitude”). The violet line
depicts the (pseudo) sunspot number according to Usoskin et al. (2021), where
in most cases the attribution to the minima and maxima is rather clear. Two
evident exceptions occur for the maxima at 1435 and 1468, which are obviously
due to typing errors in Table 1 of Usoskin et al. (2021), and which we
correct, according to the 14C-curve (violet dashed lines in the middle panel
of Fig. 1) to the years 1425 and 1461, respectively.

The red and green open squares at the upper and lower abscissa depict the
Schove/NOAA maxima and minima, respectively. After the two trivial corrections
mentioned above, we obtain a sequence of one-to-one matches between the 44
minima and maxima of Usoskin and Schove/NOAA that stretches uninterruptedly
over the time interval between 1140 and 1620.

Evidently, there are three distinctive segments where this one-to-one
correspondence fails or at least becomes problematic. The first one concerns
the long interval between 1040 and 1140, the second one is situated around
1650, the third one around 1840.

Figure 2: Plausible corrections of minima/maxima in two late time intervals.
Symbols and lines as in Fig. 1., except that two binomially filtered curves
bsm 9 (including 9 coefficients, dashed violet) and bsm 15 (including 15
coefficients, dotted violet) are added to the original 14C data (violet full
line). Lower panel: Interval from 1790 till 1890. The filtered curves
insinuate that the two minima at 1839 (with a low quality flag of $q=2$) and
1846 might indeed be only one minimum at 1843 which would also better fit to
the observational data. The merging of the minima is indicated by two
horizontal light-blue arrows, the cancellation of the maxima between them by a
corresponding light blue cross. Upper panel: Interval between 1600 and 1700.
Here, amidst the Maunder minimum, the quality flags of the minima and maxima
are typically low, making their unambiguous identification quite hard. A most
plausible cancellation of a flat minimum is indicated by a yellow cross.
Figure 3: As Fig. 2, but for the early interval between 1030 and 1160. The
vertical violet lines indicate the insertion of two minima/maxima pairs in the
later segment of that interval. The black horizontal arrows and crosses
indicate contractions and cancellations of various minima or maxima in the
early segment of the interval, of which we show three different permutations
indicated by full, dashed and dotted types of black crosses.

In Fig. 2 we consider two particular segments in more detail. Let us start
with the latest part, around 1840, which is shown in the lower panel of Fig.
2. Evidently, the NOAA data comprise only one minimum at 1843.5, whereas
Usoskin’s data show two minima here, at 1839 and 1847, the former of which
having a relatively low quality flag of $q=2$. In order to shed more light on
this issue we add to the original violet curve of sunspot data two further
curves representing two different binomial smoothings (bsm 9 and bsm 15) which
tend to smear out the two shallow minima into a single one centered at 1843,
which indeed corresponds to the NOAA value. Given the high validity of the
observationally constrained minimum in the middle of the 19th century, we
consider such contraction of two minima into a single one (indicated by the
two light blue arrows) and the corresponding cancellation (light blue cross)
of one maximum as highly plausible.

Next we turn to the situation around 1645 (upper panel in Fig. 2). Here,
amidst the Maunder minimum, we should have less trust into the NOAA data (as
for the problem of interpretation of naked-eye sunspot observations during
this time, see Carrasco et al. (2020)). Hence, it is not excluded that the
additional minimum/maximum pair of Usoskin’s data is indeed a real one. This
possibility will be discussed below. Still, it is also plausible that
Usoskin’s data show one minimum/maximum pair too much. From the most evident
variants to contract either the two (flat) maxima at 1638 and 1646 or those at
1646 and 1655, we show only the latter one, indicated with yellow arrows,
together with the cancelled minimum at 1650 (yellow cross).

Figure 4: O-C plot of various data sets of cycle minima. For the sake of
better visibility, the residuals (which all refer to a linear trend with
11.07-yr period), are differently shifted on the ordinate axis. The light-
blue, yellow, and black curves correspond to the different corrections in the
various panels of Figs. 2 and 3. The single violet full circle corresponds to
a possible shift of the 1109 year minimum as proposed in Fig. 3.

Finally, we treat, in Fig. 3, the long period between 1040 and 1140. This
interval, which strongly overlaps with the Oort minimum, is characterized by a
large number of low quality flags. While the total number of minimum/maximum
pairs in this segment is the same for the data of Usoskin and Schove, we
immediately notice the presence of two extremely long neighboring cycles
between the maxima at 1099 and 1119 (20 years) and 1119 and 1137 (18 years).
In either of those intervals we observe the existence of one (or two) local
minimum/maximum pair(s). Without overemphasizing the validity of Schove’s data
(at least the corresponding maxima were not labelled as uncertain by him, in
contrast to many others in Schove (1983)), we see at least that an insertion
of the two additional maximum/minimum pairs (indicated by the violet arrows)
leads to a one-to-one match of the extreme of both data sets. Yet, the
insertion of the first maximum leads to an even greater uncertainty for the
minimum just before it, which might well shift from its place at 1109
(according to Usoskin) to some position before. This variant is indicated by
the bent dotted violet curve ending in the alternative minimum at 1104 (violet
full circle).

Even more uncertain than in this late segment of the 1040-1140 interval is the
situation in the early segment which contains quite a number of low-q extrema,
leading to a significant number of possible permutations of minimum/maximum
pairs, in addition to those chosen by Usoskin et al. (2021). In order to “make
good” for the two insertions in the later segment, we opted here for a
compensating cancellation of two minima/maxima pairs in the early segment, a
choice which is, admittedly, strongly debatable. The different panels in Fig.
3 illustrate three most plausible permutations of minimum/maximum
cancellations in this early segment. The corresponding contractions or
cancellations are indicated by black arrows and three types of black crosses
(full, dashed, dotted). While in all three cases we obtain a one-to-one match
of the resulting minimum/maximum pairs with those of Schove, the first and
third variant (with the full and dotted crosses) lead to the most plausible
correspondences.

##  3 O-C plots, Dicke’s ratio, and the influence of the Suess-de Vries cycle

Figure 5: Dicke’s ratio for the residuals shown in Fig. 4, with corresponding
colors.

A first qualitative hint for phase-stability, or its absence, can be gained
from the so-called “Observed Minus Calculated” (O-C) diagram, which shows the
differences (or residuals) of a given data set from a linear trend for which
we use here an alleged 11.07-yr period of the Schwabe cycle. The lowermost
green curve (with open squares) in Fig. 4 shows the residuals for the combined
Schove/NOAA minima. Evidently, they are wiggling around a rather horizontal
line by typically not more than $\pm 4$ years, which - if confirmed - would
strongly speak in favour of a noise-perturbed, but nevertheless clocked
process. The corresponding residuals for Usoskin’s original minima are shown
as the uppermost green line with open circles (which is - for better
visibility - vertically shifted). Between 1140 and 1640 it exhibits already a
long horizontal segment, pointing again to phase stability in this interval.
Not surprisingly, however, it shows two steep downward-directed phase-jumps
around 1650 and 1840, where the two additional minima of Usoskin et al. (2021)
are intervening, as discussed earlier. Another remarkable feature is the
downward-pointing “nose” centered around 1090 which results from the two
additional minima in the early segment, combined with two missing minima in
the later segment of the 1040-1140 interval.

The light blue curve represents the residuals (again vertically shifted) after
having contracted the two minima at 1839 and 1846 into one at 1843. As a
consequence, the downward-directed phase jump disappears here. The next
correction, i.e. the cancellation of the minimum at 1650, leads to the yellow
curve which is already dominated by a long horizontal segment between 1140 and
1890.

We now turn to the most problematic “nose” between 1040 and 1140. All three
black curves in Fig. 4 rely on the insertion of two additional minima at 1115
and 1133 in the late segment, as was specified by the dashed violet arrows in
Fig. 3. They differ, however, by the specific combination of the two canceled
minima in the early segment, as shown in the three panels of Fig. 3. Full,
dashed, and dotted black lines in Fig. 4 correspond to the different types of
crosses in Fig. 3. In either case the previous downward-directed “nose” morphs
into a (mainly) upward-directed one which is significantly less pronounced,
though. Still, there remains one rather dominant peak at 1109. If we were to
shift the corresponding minimum to the not less plausible year 1104 (as shown
by the bent dotted violet line in Fig. 3), we would end up here at the violet
full circle in Fig. 4. Finally, with those corrections in the three intervals
we arrive at the rather horizontal black lines in Fig. 4 which are wiggling
around a horizontal by not more than $\pm 5$ years.

In Fig. 5 we show now the curves for Dicke’s ratio corresponding to all the
lines depicted in Fig. 4. This quantity is defined as the ratio
$\sum_{i}^{N}r_{i}^{2}/\sum_{i}(r_{i}-r_{i-1})^{2}$ between the mean square of
the residuals $r_{i}$ to the mean square of the differences $r_{i}-r_{i-1}$
between two consecutive residuals (Dicke 1978). For a random walk process,
Dicke’s ratio - for $N$ residuals taken into account - behaves as
$(N+1)(N^{2}-1)/(3(5N^{2}+6N-3))$ (dark blue line in Fig. 5) with its
asymptotic limit $N/15$, while the corresponding dependence for a clocked
process reads $(N^{2}-1)/(2(N^{2}+2N+3))$ (red line), with its asymptotic
limit $1/2$. Note that in Fig. 5 - in contrast to the residuals in Fig. 4 -
the period is not fixed to 11.07 years, but is separately computed for each
value $N$ of data points taken into account. The lowermost green curve shows -
again for the Schove/NOAA data - a close proximity to the curve for a clocked
process with its asymptotic limit 1/2. By contrast, Dicke’s ratio for
Usoskin’s original data (upper green curve) wiggles around the dark blue curve
for a random walk process, in good accordance with the observation by
Weisshaar et al. (2023) (who used, though, the slightly different ratio of
variances as derived by Gough (1981)). The light blue curve appears after the
first correction around 1840, the yellow curve after the second correction
around 1650. Hereby, we come already pretty close to the theoretical curve for
a clocked process (red). The additional corrections between 1040 and 1140
(only shown for the full black curve in Fig. 4) result then in the black curve
which is very close to the one for the Schove/NOAA data.

It is remarkable that the drastic difference between an apparently random-walk
like curve such as the light blue one and an apparently clocked-process like
curve such as the yellow one stems only from one single additional intervening
minimum at 1650. If we assume - for the sake of argument - that the additional
minimum at 1650 is indeed real, it would just correspond to one single phase
jump embedded into an otherwise nicely clocked process, just as discussed
recently (referring to other data sets) for two different phase-jump
candidates at 1565 and 1795 (Stefani et al. 2020b). Showing a strong
similarity to a random walk process, the shape of Dicke’s ratio would then be
completely misleading in this respect. In Fig. A1 of the Appendix we also
evidence a strong dependence of the shape of Dicke’s ratio on the very
position of such an intervening phase jump. It remains to be seen whether an
appropriately constructed statistical measure could be found to reliably
distinguish between random walk and clocked processes also in case of
intervening phase jumps. For the time being we advice to have always a
complementary glance on the O-C diagram (such as Fig. 4) which in some sense
is an even more telling device than Dicke’s ratio.

Apart from that problem, we also observe that even for the Schove/NOAA and for
the “optimally corrected” data of Usoskin Dicke’s ratio does not perfectly
approach the asymptotic limit 1/2. At this point we reiterate a pertinent
argument discussed already in Fig. 2 of Stefani et al. (2019) that a
significant share of the variance of the residuals is contained in a long-term
cycle of the Suess-de Vries type. This is again illustrated in Fig. 6a which
shows the residuals of the Schove/NOAA data - this time for the enlarged
interval from 980 to 2009, together with an optimal sinusoidal fit whose
period turns out to be 203 years. It is evident that this Suess-de Vries type
cycle entails quite a lot of the variance of the data. As a side remark: this
outcome speaks much in favour of the quality of Schove’s data who had - in the
construction of his series - never “put in” any long-term periodicity of this
kind (though he had well recognized it in hindsight, see p. 25 of Schove
(1983)). When subtracting the fitted 203-yr cycle from the original data, the
remaining residuals produce a Dicke ratio as shown as the lower green curve
(with triangles) in Fig. 6b, which now much better approaches the asymptotic
value of 1/2.

Figure 6: Illustration of the influence of the Suess-de Vries cycle on Dicke’s
ratio. (a) Residuals for the minima of the combined Schove/NOAA data for the
extended interval until 2009. The dashed thick line represents an optimal fit
of the residuals with a period of 203 years. (b) Dicke’s ratio for the
residuals of the data in (a), and for the corresponding data with the Suess-de
Vries trend being subtracted beforehand. Evidently, after subtraction of the
Suess-de Vries trend, the approachment of the curve towards the asymptotic
limit 0.5 (for clocking) becomes significantly closer. Note in particular that
the strong “overshooting” of the original Schove/NOAA curve for low $N$ is
widely suppressed by this subtraction.

##  4 Conclusions

In this paper, we have reanalyzed the series of annual (pseudo) sunspot
numbers from Usoskin et al. (2021) with particular view on a possible phase
stability of the minima and maxima of the solar cycle.

The corresponding sequence of extrema was compared with another sequence
comprising Schove’s data until 1609 and the standard Schwabe cycles after that
year. We have basically confirmed the outcome of Weisshaar et al. (2023) by
showing that the curve of Dicke’s ratio for the original sequence of Usoskin’s
minima looks formally similar to that of a classical random walk process. Yet,
we have also shown that this series comprises a very phase stable segment
interval between 1140 and 1640, with a one-to-one match of the corresponding
extrema with those of the Schove/NOAA series. Given that there is only one
possibility to get such a one-to-one match, in contrast to quite a number of
possibilities to infer less or more extrema from the original 14C data, this
is already a remarkable result that also reassures the plausibility of
Schove’s cycle reconstruction (which is by some “regarded as archaic” Usoskin
et al. (2017)).

A first correction of Usoskin’s series in form of a contraction of two minima
around 1840 into one seems highly plausible given the low quality flag of the
former of the two minima and also in view of the high observational validity
of only one minimum in this time span.

We have tried a second correction in form of a cancellation of one minimum at
1646. Here, amidst the Maunder minimum, the quality flags of all of Usoskin’s
minima are typically quite low. Admittedly, during this time the observational
validity of the standard Schwabe cycle is also not very high, so that the
justification of this cancellation remains doubtful. If we accept it for the
moment (also considering that two successive short cycles are not very likely
in this particular time of a very weak solar dynamo) we end up with a phase-
stable time interval between 1140 and 1890, whose Dicke ratio approaches
closely the curve for a clocked process. But even if the additional minimum at
1646 turned out to be real, it would just correspond to one single phase jump
embedded into a long period that is otherwise phase stable. This property is
best observed in the O-C diagram, while the curve of Dicke’s ratio makes the
data look like a random walk process.

The most ambiguous time interval is that between 1040 and 1140 (basically the
Oort minimum), where nearly all quality flags of Usoskin’s data are low. The
O-C diagram exhibits here a pronounced downward-pointing “nose”, stemming from
two very long cycles in the later segment and some correspondingly short
cycles in the earlier segment of this time interval. With two plausible
insertions of minima in the later segment, and two compensating
contractions/cancellations of minima in the earlier segment, we reach again a
reasonable one-to-one match with the minima of Schove. The resulting O-C
diagram is significantly smoothed and shows now a pretty horizontal line
between 970 and 1890, with not much more wiggling (approximately $\pm 5$
years) than in case of the Schove/NOAA data. In our view, this finding
strongly reinforces the validity of Schove’s data, and impugns Nataf’s
criticism of them as being simply construed by the “nine-per-century” rule.

Having been focused on a minimal number of corrections pointing (somewhat
biasedly) towards a clocked process, we admit that the high ambiguity of
Usoskins’s extrema data (with 29 of them having a quality flag $q=0$) also
entails the possibility that even more phase jumps might exist. Obviously,
with an increasing number of such events the entire notion of phase stability
would become more and more problematic.

At any rate, we conclude that, before entering into a statistical analysis of
the clocked - or non-clocked - character of the solar dynamo, the underlying
data should be carefully scrutinized. The data set of cycle extrema as
produced by Usoskin et al. (2021) entails quite a couple of intervals with low
quality flags where the specific extreme should be taken with a grain of salt.
While this was clearly expressed in Usoskin et al. (2021), a too uncritical
adoption of the data as by Weisshaar et al. (2023) and Biswas et al. (2023)
might lead to wrong conclusions. In this sense, their argument for a non-
clocked process appears premature since its allegedly high significance
depends crucially on the selection of the specific set of extrema according to
Usoskin et al. (2021). While we still refrain from claiming perfect evidence
for solar cycle synchronization, we argue that the work of Weisshaar et al.
(2023) does neither represent any conclusive evidence for its absence.

###### Acknowledgements.

This work has received funding from the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation programme (grant
agreement No 787544). We thank Carlo Albert and Antonio Ferriz Mas for
fruitful discussions about the solar dynamo and its possible synchronization.

## References

  * Abreu et al. (2012) Abreu, J.A., Beer, J., Ferriz-Mas, A., McCracken, K. G., Steinhilber, F. 2012, A&A, 548, A88
  * Biswas et al. (2023) Biswas, A., Karak, B. B. , Usoskin, I., Weisshaar, E. 2023, arXiv:2302.1484
  * Bollinger (1952) Bollinger, C. J. 1952, Proc. Okla. Acad. Sci., 33, 307
  * Brehm et al. (2021) Brehm, N., Bayliss, A., Christl, M., et al. 2021, Nature Geoscience, 14, 10
  * Carrasco et al. (2020) Carrasco, V. M. S., Gallego, M. C., Arlt, R., Vaquero, J. M. 2020, ApJ, 904, 60
  * Charbonneau (2022) Charbonneau, P. 2022, Front. Astron. Space Sci., 9, 853676
  * De Jager & Versteegh (2005) De Jager, C., Versteegh, G. 2005, Solar Phys., 229, 175
  * de la Rue et al. (1872) de la Rue, W., Stewart, B., Loewy, B. 1872, Proc. R. Soc. Lond. Ser., 21, 399
  * Dicke (1978) Dicke, R. H. 1978, Nature, 276, 676
  * Gough (1981) Gough, D. 1981, in NASA Conference Publication, Vol. 2191, 185-206
  * Horstmann et al. (2023) Horstmann, G. M., Mamatsashvili, G., Giesecke, A. Zaqarashvili, T. V, Stefani, F. 2023, ApJ, 944, 48
  * Hung (2007) Hung, C.-C. 2007, Apparent relations between solar activity and solar tides caused by the planets, NASA/TM-2007-214817.
  * Jose (1965) Jose, P. D 1965, Astron. J., 70, 193
  * Klevs et al. (2023) Klevs, M., Stefani, F., Jouve, L. 2023, arXiv:2301.05452
  * Nataf (2022) Nataf, H.-C. 2023, Solar Phys., 297, 107
  * Okhlopkov (2016) Okhlopkov, V. P. 2016, Mosc. Univ. Phys. B. 71, 440
  * Öpik (1972) Öpik, E. 1972, I. Astron. J., 10, 298
  * Scafetta (2012) Scafetta, N. 2012, J. Atmos. Sol.-Terr. Phys., 81-82, 27
  * Scafetta (2023) Scafetta, N. 2023, Solar Phys., 298, 24
  * Schove (1983) Schove, D. J. 1983, Sunspot Cycles (Hutchinson Ross Publishing Company, Stroudsburg)
  * Schove (1984) Schove, D. J. 1984, Chronologies of eclipes and comets AD 1-1000 (The Boydel Press, Woodbridge and Dover)
  * Stefani et al. (2016) Stefani, F., Giesecke, A., Weber, N., Weier, T. 2016, Solar Phys., 291, 2197
  * Stefani et al. (2018) Stefani, F., Giesecke, A., Weber, N., Weier, T. 2018, Solar Phys. 293, 12
  * Stefani et al. (2019) Stefani, F., Giesecke, A., Weier, T. 2019, Solar Phys., 294, 60
  * Stefani et al. (2020a) Stefani, F., Giesecke, A., Seilmayer, M., Stepanov, R., Weier, T. 2020a, Magnetohydrodynamics, 56, 269
  * Stefani et al. (2020b) Stefani, F., Beer, J., Giesecke, A., Gloaguen, T., Seilmayer, M., Stepanov, R., Weier 2020b, Astron. Nachr. 341, 600
  * Stefani et al. (2021) Stefani, F., Stepanov, W., Weier, T. 2021, Solar Phys., 296, 88
  * Takahashi (1968) Takahashi, K. 1968, Solar Phys., 3, 598
  * Usoskin et al. (2017) Usoskin, I. G., Living Rev. Sol. Phys., 14, 3
  * Usoskin et al. (2021) Usoskin, I. G., Solanki, S. K., Krivova, N. A., et al. 2021, A&A, 649, A141
  * Vos et al. (2004) Vos, H. et al. 2004, in Towards a Synthesis of Holocene Proxy Data and Climate Models, ed. H. Fischer, T. Kumke, G. Lohmann, G. Flöser, H. Miller, H., H. von Storch, J. F. Negendank, The Climate in Historical Times (GKSS School of Environmental Research, Springer) 293
  * Weber et al. (2013) Weber, N., Galindo, V., Stefani, F., Weier, T., Wondrak, T. 2015, New J. Phys., 15, 043034
  * Weber et al. (2015) Weber, N., Galindo, V., Stefani, F., Weier, T. 2015, New J. Phys., 17, 113013
  * Weisshaar et al. (2023) Weisshaar, E., Cameron, R. H., Schüssler, M. 2023, arXiv:2301.07469
  * Wilson (2013) Wilson, I. R. G. 2013, Pattern Recogn. Phys., 1, 147
  * Wolf (1859) Wolf, R. 1859, MNRAS, 19, 85
  * Wood (1972) Wood, K. 1972, Nature 240, 91

##  Appendix A Role of the position of a phase jump

In this appendix we illustrate the influence of an intervening additional
minimum - embedded into an otherwise clocked process - on the shape of the
curves of Dicke’s ratio. For that purpose, we utilize the Schove/NOAA minima
data, this time extended until 2009.

Fig. A.1 shows Dicke’s ratio for this data set, and for 6 further data where
phase jumps are artificially inserted at appropriate positions close to the
years 1000, 1200, 1400, 1600, 1800, and 2000. Obviously, the similarity of
Dicke’s ratio to that of a clocked process becomes most pronounced for phase
jumps inserted in the center of the time interval.

Figure 7: Dicke’s ratio in presence of the position of phase jumps, based on
the Schove/NOAA minimum data between 970 and 2009. Phase jumps are
artificially inserted at appropriate positions close to the years 1000, 1200,
1400, 1600, 1800, and 2000.