# Survival analysis for AdVerse events with VarYing follow-up times (SAVVY):
summary of findings and a roadmap for the future of safety analyses in
clinical trials

Kaspar Rufibach
Jan Beyersmann
Tim Friede
Claudia Schmoor
Regina Stegherr Methods, Collaboration, and Outreach Group, Product
Development Data Sciences (MCO), Hoffmann-La Roche Ltd, Basel,
SwitzerlandInstitute of Statistics, Ulm University, Ulm, GermanyDepartment of
Medical Statistics, University Medical Center GöttingenClinical Trials Unit,
Faculty of Medicine and Medical Center, University of FreiburgCharité -
Universitätsmedizin Berlin, Corporate Member of Freie Universität Berlin and
Humboldt-Universität zu Berlin, Institute of Biometry and Clinical
Epidemiology

###### Abstract

The SAVVY project aims to improve the analyses of adverse events (AEs) in
clinical trials through the use of survival techniques appropriately dealing
with varying follow-up times and competing events (CEs). This paper summarizes
key features and conclusions from the various SAVVY papers. Through
theoretical investigations using simulations and in an empirical study
including randomized clinical trials from several sponsor organisations,
biases from ignoring varying follow-up times or CEs are investigated. The bias
of commonly used estimators of the absolute and relative AE risk is
quantified. Furthermore, we provide a cursory assessment of how pertinent
guidelines for the analysis of safety data deal with the features of varying
follow-up time and CEs. SAVVY finds that for both, avoiding bias and
categorization of evidence with respect to treatment effect on AE risk into
categories, the choice of the estimator is key and more important than
features of the underlying data such as percentage of censoring, CEs, amount
of follow-up, or value of the gold-standard. The choice of the estimator of
the cumulative AE probability and the definition of CEs are crucial. SAVVY
recommends using the Aalen-Johansen estimator (AJE) with an appropriate
definition of CEs whenever the risk for AEs is to be quantified. There is an
urgent need to improve the guidelines of reporting AEs so that incidence
proportions or one minus Kaplan-Meier estimators are finally replaced by the
AJE with appropriate definition of CEs.

Keywords: Aalen-Johansen estimator; Cumulative incidence function; adverse
events; competing events; drug safety; incidence proportion; incidence
density; Kaplan-Meier estimator.

Running head: SAVVY overview.

##  1 Background

In randomized clinical trials (RCT), an essential part of the benefit-risk
assessment of treatments is the quantification of the risk of experiencing
adverse events (AEs), and comparing these risk between treatment arms. Methods
commonly employed to quantify absolute adverse event (AE) risk either do not
account for varying follow-up times, censoring, or for competing events (CEs),
although appreciation of these are important in risk quantification, see e.g.
O’Neill[1] and Procter and Schumacher[2].

Analyses of AE data in clinical trials can be improved through use of survival
techniques that account for varying follow-up times, censoring and CEs. The
AJE [3, 4] can be considered the non-parametric gold-standard method when
quantifying absolute AE risk. The reason is that the AJE is the standard (non-
parametric) estimator that accounts for CEs, censoring, and varying follow-up
times simultaneously, and, being non-parametric, does not rely on restrictive
parametric assumptions, such as constant hazards. Any other estimator of AE
probability, such as incidence proportion, probability transform incidence
density, or one minus Kaplan-Meier, delivers biased estimates in general.

To quantify that bias for all these methods in an ideal scenario Stegherr et
al.[5] ran a comprehensive simulation study. Two key findings were (1) that
ignoring CEs is more of a problem than falsely assuming constant hazards by
the use of the incidence density and (2) that the choice of the AE probability
estimator is crucial for the estimation of relative effects, i.e. comparison
between groups.

To illustrate and further solidify these simulation-based results with real
data the SAVVY consortium, a collaboration between nine pharmaceutical
companies and three academic institutions, meta-analyzed data from seventeen
randomized controlled trials (RCT).

In this review article we summarize the results of the empirical study,
reported in two separate publications: Stegherr et al.[6] was concerned with
estimation of AE risk in one treatment group and Rufibach et al.[7] with the
comparison of AE risks between two groups in an RCT. A cursory assessment how
relevant guidelines recommend to estimate AE risk is given with a call for
updates. We conclude with a discussion.

##  2 Scientific questions of the SAVVY project

The overarching scientific questions of SAVVY can be phrased as follows:

  1. 1)

For estimation of the probability of an AE, how biased are commonly used
estimators, especially the incidence proportion and one minus Kaplan-Meier, in
presence of censoring, varying follow-up between patients, CEs, and in the
case of incidence densities a restrictive parametric model?

  2. 2)

What is the bias of common estimators that quantify the relative risk of
experiencing an AE between two treatment arms in a RCT?

  3. 3)

Can trial characteristics be identified that help explain the bias in
estimators?

  4. 4)

How does the use of potentially biased estimators impact qualification of AE
probabilities and relative effects in regulatory settings?

Within the SAVVY project, these questions were approached in two ways: First,
in Stegherr et al.[5], via simulation of clinical trial data. This approach
has the advantage that the true underlying data generating mechanism is
specified by the researcher and therefore known. This allows to exactly
quantify the bias of a given estimator, i.e. to answer 1) and 2) above (it
would also allow to answer Question 4, but that was not addressed in Stegherr
et al.[5]). Second, in Stegherr et al.[6] and Rufibach et al.[7], biases of
commonly used estimators of absolute and relative AE risks were estimated by
comparing them to the best available estimator. Having real clinical trial
data available also allows to answer Question 3 above, through meta-analytic
methods.

###  2.1 Competing events and their connection to the ICH E9 estimands
addendum

In what follows we will use competing event (apart from direct quotes) and
consider it synonymous to competing risk.

An important but largely unrecognized aspect when quantifying AE risk is the
likely presence of CEs. Gooley et al.[8] define a CE as

> “We shall define a competing risk as an event whose occurrence either
> precludes the occurrence of another event under examination or fundamentally
> alters the probability of occurrence of this other event.”

whereas the ICH E9(R1) estimands addendum[9] defines an intercurrent event as

> “Events occurring after treatment initiation that affect either the
> interpretation or the existence of the measurements associated with the
> clinical question of interest.”

These two definitions appear to be, if not the same, then at least very
related. However, the ICH E9(R1) addendum does not discuss CEs, so it is not
entirely clear how to embed CEs into the addendum framework, i.e. whether and
if yes which of the proposed strategies of the estimand addendum applies to
the CE situation. More research and discussion is needed to align CEs (if
necessary at all), and the analysis of complex time-to-event data with the
addendum.

Stegherr et al.[6] and Rufibach et al.[7] took a pragmatic approach and
defined events as “competing” that preclude the occurrence or recording of the
AE under consideration in a time-to-first-event analysis. Specifically, one
important CE is death before AE. In addition, any event that would both be
viewed from a patient perspective as an event of his/her course of disease or
treatment and would stop the recording of the interesting AE was viewed as a
CE. Since all these CEs apart from death may be prone to some subjectivity in
the empirical analysis reported in Stegherr et al.[6] and Rufibach et al.[7] a
variant of the estimators with a CE of death only was also considered. Since
results were in line with the broader definition of CEs as given above we omit
the results for the death only variant here.

###  2.2 Estimation methods

A precise mathematical definition of all estimators of the probability of an
AE that were compared in SAVVY is provided in Stegherr et al.[10], a
prospectively published statistical analysis plan for the SAVVY project. A
short introduction in all estimators is also provided in Stegherr et al.[6].
In this overview article we only provide very brief descriptions of the
considered estimators. By far the most commonly used estimator, e.g. in
standard safety reporting that enters benefit-risk assessment for the approval
of new medicines, to estimate the risk of an adverse event up to a maximal
observation timepoint $\tau$ is the incidence proportion[4]. It simply divides
the number of patients with an observed AE on $[0,\tau]$ in group $A$ by the
number of patients in group $A$. The incidence proportion is an estimator of
the probability that an AE happens in the interval $[0,\tau]$, and that this
AE is observed, i.e. not censored. This illustrates that the incidence
proportion is not properly dealing with censored observations. However, it
correctly accounts for CEs, see Allignol et al.[4] for an exemplary
illustration of that feature. To account for the differing follow-up times
between patients researchers and guidelines (see Table 4) often advocate to
use the incidence density, where the number of AEs in the nominator is divided
by the total patient-time at risk instead of simply the number of patients. As
such, the incidence density does not directly estimate the probability of an
AE, but rather the AE hazard. As described in Stegherr et al.[10], this hazard
estimator can easily be transformed to indeed estimate the probability of an
AE. However, it only does so assuming that the AE hazard is constant, i.e. the
probability-transformed incidence density is a fully parametric estimator. In
addition, as such it does not correctly account for CEs, but it can be
modified to do so, leading to the probability transform incidence density
accounting for CEs. Researchers are often aware of the inability of the
incidence proportion to properly deal with varying follow-up times and
censoring. As a remedy they (and many guidelines, see Table 4) then advocate
to consider time to AE and estimate the probability of an AE by reading off
the one minus Kaplan-Meier estimator at a timepoint of interest, e.g. the end
of observation time $\tau$ or at an earlier timepoint. While indeed, this
estimator properly accounts for varying follow-up times and censoring the
question remains how to deal with CEs. Numerous papers have been written[8, 4]
and providing technical arguments explaining why the Kaplan–Meier estimator is
a biased estimator of the probability of an AE. Intuitively, one minus the
Kaplan–Meier estimator estimates the distribution function of the time of
interest, that is, it tends towards one as we move to the right on the time
axis. However, if we add up the probability of an AE and the probability of a
CE in a truly CE scenario that also must be equal to one, implying that the
probability of an AE is strictly smaller than one. As a consequence, the
Kaplan–Meier estimator to estimate the probability of an AE will be biased
upwards. Finally, there is an estimator that at the same time accounts for
(random or independent) censoring, respects varying follow-up times, accounts
for CEs in the right way, is fully nonparametric and therefore free of bias
introduced by any of these processes: the AJE [3]. It is therefore considered
the gold standard estimator. In the empirical analysis of the 17 RCTs in the
SAVVY project, it served as a benchmark against which all estimators were
measured against. The term bias was therefore used for deviations of the
estimators from this benchmark estimator, or gold standard. For an evaluation
of the true bias, i.e. the deviation of estimators to the true underlying
value, we refer to Stegherr et al.[5].

Table 1 in Stegherr et al.[6] concisely summarizes the properties of each
considered estimator with respect to whether it accounts for censoring, CEs,
and whether it makes a parametric assumption, and we therefore reproduce it
here in Table 1 for the estimators discussed here.

Table 1: Overview whether the estimators deal with the possible sources of bias, reproduced from Stegherr et al.[6].  | Accounts for | Makes no constant | Accounts for
---|---|---|---
| censoring | hazard assumption | CEs
Incidence proportion | No | Yes | Yes
1-Kaplan-Meier | Yes | Yes | No
gold-standard AJE | Yes | Yes | Yes

###  2.3 Quantification of bias - the SAVVY project

One of the goals of the SAVVY project is to quantify the bias of standard
estimators of the probability of an AE. Based on simulations, i.e. comparing
estimated to true underlying values from which the data was simulated, this
has been done in Stegherr et al.[5]. Key findings in this study were that
ignoring CEs is more of a problem than falsely assuming constant hazards. The
one minus Kaplan-Meier estimator may overestimate the true AE probability by a
factor of four at the largest observation time. Moreover, the choice of the AE
probability estimator is crucial for group comparisons.

To confirm these results in real clinical trial data three academic
institutions and nine pharmaceutical companies collaborated within the SAVVY
consortium. In order not to have to share the trial data macros to compute the
above estimators were developed centrally. Companies then ran these macros on
trial data of their choice and only returned high-level results (estimated
values of AE probabilities for each estimator and some basic trial
characteristics) to the central data analysis unit hosted at one of the
academic institutions. The central data analysis unit then meta-analyzed these
results to quantify biases for estimators and impact on regulatory decision-
making. A statistical analysis plan for these analyses was published[10].
Results for estimation of the probability of an AE in one arm are reported in
Stegherr et al.[6] and for the comparison between two groups in Rufibach et
al.[7].

In this overview paper we focus attention on summarizing results for the most
commonly used estimators of AE risks, namely the incidence proportion and one
minus Kaplan-Meier in comparison to the gold standard AJE. Furthermore, we
report results for the maximal evaluation time as defined in the above two
papers, i.e. the latest time at which a dataset has an observation, either
event or censored. Results for shorter evaluation times were in line with
those for the maximal evaluation time.

##  3 Results

Ten organisations provided 17 trials including 186 types of AEs (median $8$;
interquartile range $[3,9]$). Twelve (71.6% out of 17) trials were from
oncology, nine (52.9%) were actively controlled and eight (47.1%) were placebo
controlled. The trials included between 200 and 7171 patients (median: 443;
interquartile range $[411,1134]$).

For one of the 17 trials, details of the trial and the analysis of one of the
AEs by the different methods investigated in this paper are presented in
Stegherr et al.[5].

Note that for RRs only those AEs were considered where neither the estimated
AE probability in the experimental arm nor the estimated AE probability in the
control arm is 0 ($n=156$ for incidence proportion and $n=155$ for one minus
Kaplan-Meier). This also applies to the categories in Table 2.

###  3.1 Empirical bias of common estimators of absolute risk with respect to
the gold-standard AJE

For the comparison of the AE probabilities Stegherr et al.[6] focussed on the
experimental arm. Median follow-up was 927 days (interquartile range
$[449,1380]$). The median of the gold-standard AJE was 0.092 (minimum 0 and
maximum 0.961), i.e. the estimated probability of an AE was 9.2% on average
over all 186 AEs over all trials. The one minus Kaplan-Meier estimator was on
average about 1.2-fold larger than the AJE and the probability transform of
the incidence density ignoring CEs was even two-fold larger. The average
empirical bias (i.e. the difference between the considered estimator and the
gold-standard AJE) using the incidence proportion was less than 5%. Assuming
constant hazards using incidence densities was hardly an issue provided that
CEs were accounted for. However, beyond these average biases it was striking
how extreme the bias could become: For one minus Kaplan-Meier, in our
empirical analysis using real clinical trial data, we observed an
overestimation of the AE probability up to a factor of five, whereas for the
incidence proportion we observed underestimations of up to a factor of three.
This is in line with the findings of the simulation study in Stegherr et
al.[5] and illustrates that using too simplistic estimators for the
probability of an AE can be truly misleading. To evaluate what study
characteristics impact the bias a meta-regression was performed. For that, the
response variable was defined as the AE probability estimates obtained with
the different estimators divided by the AE probability obtained with the gold-
standard AJE, considered on the log-scale. Covariates were the proportion of
censoring, the evaluation time point $\tau$ (i.e., the maximal time to event
in years – AE, CE or censoring – observed), and the size of the AE probability
estimated by the gold-standard AJE. The meta-regression showed that the bias
depended primarily on the amount of censoring and on the amount of CEs.

Finally, according to the European Commission’s guideline on summary of
product characteristics (SmPC)[11] and based on the recommendations of the
Council for International Organizations of Medical Science (CIOMS) Working
Groups III and V[12], the AE risk is classified into frequency categories
which are defined by ‘very rare’, ‘rare’, ‘uncommon’, ‘common’, and ‘very
common’ when the risk is $<$ 0.01%,$<$ 0.1%,$<$ 1%,$<$ 10%,$\geq$ 10%,
respectively. Stegherr et al.[6] (Table 2) assigned these categories to AE
probability estimates from all estimators (and all 186 AEs) and compared them
to the categories resulting from the AJE. As an example, systematic
overestimation of the one minus Kaplan-Meier resulted in upgrading of 2/6 AEs
from ‘uncommon’ to ‘common’ and 14/86 from ‘common’ to ‘very common’.

###  3.2 Empirical bias of probability-based estimators of relative risk with
respect to the gold-standard AJE

Naively, one could ask the question whether when we have biased estimates of
the AE probability in two groups, that go in the same direction in both
groups, and want to compare them, the bias “cancels out”. To assess this
hypothesis Rufibach et al.[7] extended the work from Stegherr et al.[6] to
quantify the bias when comparing AE risks between groups in a randomized
trial. The focus of Rufibach et al.[7] is not to define what a fit-for purpose
estimand to quantify safety risk could be, but rather to evaluate statistical
properties of commonly used estimators in the presence of varying follow-up
and CEs. A thorough discussion of effect measures to quantify the relative
risk is given there as well. Without going into causal or estimand details,
the effects to be compared between groups are to be understood as total
effects, comparing patients’ AE risk in this world and in the presence of CEs.
The estimators that were finally assessed are

  * •

the risk ratio (RR) $\widehat{RR}=\hat{q}_{E}/\hat{q}_{C}$ for estimators
$\hat{q}_{E}$ and $\hat{q}_{C}$ of AE probabilities calculated at a specific
evaluation time within each treatment arm (E for experimental, C for control),

  * •

the hazard ratio (HR) for both, the AE of interest and the CE.

In the one-sample case, estimates of AE probabilities were benchmarked on the
gold-standard AJE. This, because the latter is a fully nonparametric estimator
that accounts for censoring, does not rely on a constant hazard assumption,
and accounts for CEs, as discussed above. So, as a straightforward extension
for the comparison of AE probabilities between two arms using the RR, we
benchmark the latter on the RR estimated using the AJE in each arm, with
variance derived using the delta rule. The gold-standard for estimates of the
HR will be the HR from Cox regression. This, because the latter is typically
used to quantify a treatment effect not only for efficacy, but also for time-
to-first-AE type endpoints. Variances of comparisons of different estimators
of the RR and HR will be received via bootstrapping, because of the dependency
of estimators computed on the same dataset, see Stegherr et al.[5].

For the one-sample case of estimation of absolute AE risk direction of biases
can be derived analytically: incidence proportion systematically under- and
one minus Kaplan-Meier systematically overestimates the true AE risk. However,
there is no such derivation possible for direction of the bias for the RR or
HR. Rufibach et al.[7] found the RR based on the incidence proportion
overestimates the RR computed using the gold-standard AJE by up to a factor of
three, and one minus Kaplan-Meier underestimates up to a factor of four, see
Figure 1. Interestingly, dividing the two biased estimates of the AE
probability based on the incidence proportion, which both tend to
underestimate the true AE probability, leads to an estimate of the RR that on
average performs comparably to the AJE. Apart from shedding light on
estimation quality of the incidence proportion and one minus Kaplan-Meier to
estimate the RR we conclude that different patterns of under- or
overestimation of absolute AE probabilities can lead to similar performance
for RR. This implies that in general, one cannot conclude how an estimator of
the relative AE risk performs based on looking how these same estimators
performs on estimation of arm-wise AE probabilities.

As discussed in Stegherr et al.[6] one reason for the good performance of the
incidence proportion might be a high amount of CEs before possible censoring.
However, not only the proportion of censoring but also the timing of the
censoring are relevant.

Meta-analyses confirmed the above impressions. Meta-regressions showed that
(1) the key difference between estimators lies in the value of the average RR
and (2) the impact of covariates is overall limited, compared to the average
RR the estimated coefficients are close to 1. This emphasizes that the choice
of the estimator is key, and that this holds true over a wide range of
possible data configurations quantified through the considered covariates.

Figure 1: Ratio of RRs estimated with estimator of interest and the gold-
standard AJE.

A key finding of Rufibach et al.[7] was that categorization of evidence based
on RR crucially depends on the estimator one uses to estimate the RR. We
therefore reproduce Table 1 in their paper here in a version trimmed down to
incidence proprtion and one minus Kaplan-Meier, see Table 2. Overall, we find
quite a number of switches to neighboring categories, more so than for
estimators of the absolute AE risk in one arm. Reasons for switches are wider
CIs of the AJE as well as RR estimates / CI bounds that are close to the
cutoffs between categories. As the incidence proportion on average estimates
the RR well, we see a similar number of switches to a higher ($n=8$, below the
diagonal in Table 2) and lower ($n=9$) evidence category. Not surprisingly,
for one minus Kaplan-Meier that underestimates the RR with respect to the
gold-standard AJE we see relevantly more switches to a lower than higher
evidence category, namely $n=41/n=16$, $32/8$, and $28/6$, respectively.

In summary, the choice of the estimator of the RR does have a relevant impact
on the conclusions. Note that there is no universally accepted standard how
one should combine a point estimate and its associated variability, in our
case RR, into evidence categories. As an example, Rufibach et al.[7] used a
categorization motivated by the methods put forward by the IQWiG[13] for
severe AEs (Table 14) to be used for the German benefit-risk assessment.

Table 2: The impact of the choice of relative effect estimators, incidence proprtion and one minus Kaplan-Meier, for AE probabilities on qualitative conclusions. Diagonal entries are set in bold face. Deviations from the gold-standard AJE are the off-diagonal entries. Off-diagonal zeros are omitted from the display. |  |  |  |  |  | gold-standard Aalen-Johansen
---|---|---|---|---|---|---
|  |  |  |  |  | (0) no effect | (a) minor | (b) considerable | (c) major
|  |  incidence |  proportion |  | (0) no effect | 84 | 5 | 0 | 0
|  |  | (a) minor | 3 | 10 | 2 | 0
|  |  | (b) considerable | 1 | 2 | 12 | 2
|  |  | (c) major | 1 | 0 | 1 | 33
one |  minus |  Kaplan- |  Meier |  | (0) no effect | 84 | 9 | 4 | 8
| (a) minor | 3 | 6 | 3 | 3
| (b) considerable | 2 | 1 | 7 | 5
| (c) major | 0 | 1 | 1 | 18

###  3.3 Empirical bias of hazard-based estimators of relative risk with
respect to the gold-standard AJE

As for hazard-based inference it is important to note that, even if the event
of interest is AE, effects on all other CEs is generally recommended for any
(hazard-based) analysis of CEs[14]. Rufibach et al.[7] found that the hazard
of AE is generally larger for the experimental compared to the control arm,
meaning that the instantaneous risk of AE is typically higher in the
experimental arm of an RCT, not unexpectedly. For the hazards of CEs, for both
types, what we find is that the hazard in the experimental arm is generally
lower than in the control arm, i.e. there is an effect of the experimental
treatment on the CE. If we simply censored at CEs we would thus introduce arm-
dependent censoring, a feature that may lead to biased effect estimates[15,
16]. The ratio of the incidence densities of the AE in the two arms
underestimates with respect to the Cox regression HR while for the other two
endpoints on the median they turn out to be approximately unbiased compared to
the Cox HR, with a tendency to overestimation when accounting for all CEs. To
appreciate the differences between the two estimators of the RR based on
hazards, i.e. the incidence density ratio and the gold-standard Cox regression
HR, recall the properties of the two methods: Both properly account for
censoring and they properly estimate event-specific hazards, or rather the
relative effect based on these. The only difference between the two methods is
what they assume about the shape of the underlying hazard: the incidence
density assumes them to be constant up to the considered follow-up time, which
also implies that they are proportional.

The impact of the use of the different estimators on the conclusions derived
from the comparison of treatment arms was again investigated by the use of
categories. These are typically derived from comparing the confidence interval
(CI) of the RR to thresholds. In contrast to the usual IQWiG procedure,
however, they did not only categorize the benefit of a therapy, but also the
harm, where they did not distinguish between a positive and a negative
treatment effect. Four categories were possible: (0) “no effect” if 1 is
included in the CI, (a) “minor” (“gering”) if the upper bound of the CI is in
the interval $[0.9,1)$ for a RR$<1$ or the lower bound in the interval
$(1,1.11]$ for a RR$>1$, (b) “considerable” (“beträchtlich”) if the upper
bound of the CI is in the interval $[0.75,0.9)$ for a RR$<1$ or the lower
bound in the interval $(1.11,1.33]$ for a RR$>1$, and (c) “major”
(“erheblich”) if the upper bound is smaller than 0.75 for a RR$<1$ or the
lower bound greater than 1.33 for a RR$>1$. The same categorization was used
for the HR instead of RR.

Rufibach et al.[7] have considered two effect measures to quantify the RR of
an AE in two arms: the RR based on AE probabilities and the HR computed from
Cox regression. Rufibach et al.[7]’s analyses revealed that all the considered
estimators are overall inferior to the two gold-standards that were
considered, either the RR based on the arm-wise AJE or the HR based on Cox
regression. A question that remains is whether the qualitative conclusions
drawn based on the two gold-standards are relevantly different when relying on
the criteria put forward by the IQWiG (Table 14 in their general methods
document[13]), see Table 9 in Rufibach et al.[7], which we reproduce here as
Table 3. Quite some different classifications based on the two estimators of
the RR were observed. However, this is not a surprise, as the estimand the two
methods look at is not the same (see Varadhan et al.[17] for an exposition of
this issue): Cox HR quantifies a relative effect based on an endpoint of AE
_hazard_ , whereas RR based on gold-standard Aalen-Johansen is based on a
comparison of _probabilities_ at a evaluation time, see Rufibach et al.[7] for
an extended discussion of this issue.

Table 3: Conclusions of the RR calculated with the gold-standard Aalen-Johansen estimator compared to the conclusions of the HR calculated with the Cox model. The table shows the analysis of those $n=94$ AEs that were observed with an absolute frequency of $\geq 10$ in each arm. Off-diagonal zeros are omitted from the display. |  | HR Cox for AE
---|---|---
|  | (0) no effect | (a) minor | (b) considerable | (c) major
RR gold-standard Aalen-Johansen | (0) no effect | 42 | 3 | 3 | 1
(a) minor | 9 | 2 | 1 | 0
(b) considerable | 4 | 1 | 3 | 2
(c) major | 2 | 0 | 4 | 17

###  3.4 Overall conclusions from the empirical analysis

To conclude this section, based on theoretical and empirical evidence,
Stegherr et al.[6] clearly recommend the AJE as the non-parametric, unbiased
estimator in the presence of both CEs and censoring. If a parametric analysis
based on incidence densities is considered, they strongly recommend to
incorporate incidence densities of CEs as they detail in their paper. This is
also preferable over the one minus Kaplan-Meier approach. A conclusion of the
empirical study of the SAVVY project not discussed here was also that ignoring
CEs appeared to be worse than assuming constant hazards. This illustrates the
importance of careful consideration of CEs when aiming at properly estimating
and comparing AE risks.

##  4 Analysis of safety in key guidelines

To understand the landscape of existing guidelines on safety reporting which,
ideally, at some point will be updated based on the conclusions from SAVVY, we
reviewed an opportunistic sample of these and collected results in Table 4.
Many of these guidelines mention at least that varying follow-up times may be
relevant to quantify AE risk. For example, the CIOMS Working Group VI
handbook[18], which forms the basis of several guidances, already admits that
“ICH Guideline E3 mentions survival analysis methods for analysing safety
data, but it appears that this has often not been followed in practice.”
Looking at Table 4 two aspects are remarkable: First, the heterogeneity in how
related guidelines treat the issues of varying follow-up, use of the incidence
density, constant hazard assumption for the latter, proposing life-table or
one minus Kaplan-Meier techniques. Second, the complete absence of any
consideration of CEs, although at least death seems to be quite an obvious CE
in estimation of AE risk, with many others potentially relevant. Taken
together with the frequent mentioning of life-table / one minus Kaplan-Meier
methods to account for varying follow-up time this is specifically concerning
given the findings in Stegherr et al.[6] and Rufibach et al.[7] about the bias
of one minus Kaplan-Meier in presence of CEs.

Overall, it appears somewhat suprising that all guidelines exhibit awareness
of the varying follow-up time issue and even discuss potential remedies, but
that in routine reporting of safety and quantification of AE risks the
incidence proportion appears still to be the overwhelmingly preferred
approach.

Table 4: Coverage of relevant time-to-event aspects for quantification of AE
risk for an opportunistic sample of safety guidelines. ”x” means that this
aspect is mentioned in the respective guideline.

Guideline |  Acknowledges |  |  |  |
---|---|---|---|---|---
varying FU |  Proposes |  |  |  |
incidence density |  Acknowledges constant |  |  |  |
hazard assumption |  Proposes life-table / |  |  |  |
one minus Kaplan-Meier |  Acknowledges |  |  |  |
CEs |  |  |  |  |
ICH E3[19] |  x |  |  |  x |
ICH E9[20] |  x |  x |  |  x |
SmPC[11] |  x |  x |  |  |
CIOMS[18] |  x |  x |  x |  x |
FDA premarketing[21] |  x |  x |  x |  x |
CONSORT Harm 2022 update[22] |  x |  |  |  |

##  5 Implementation

All methods introduced in the SAVVY project have been implemented in the R
package savvyr[23]. This package is available from CRAN.

##  6 Discussion

###  6.1 Main conclusion of the SAVVY project

The main conclusion from the SAVVY project is that commonly used methods such
as incidence proportions, incidence densities (with and without ignoring CEs),
and one minus Kaplan-Meier are all biased and therefore inappropriate to
quantify AE risk in the presence of varying follow-up times, CEs and
censoring. Estimators are biased for estimation for absolute as well as
relative AE risks. It is important to note that this bias is a statistical
property of any of these estimators and independent of the purpose we use any
of these estimators for, i.e., whether we quantify the risk for a prespecified
or emerging AE, or estimate AE risk in a given therapeutic area, or want to
detect a different AE signal between two treatment arms. Taking together,
Stegherr et al.[5], Stegherr et al.[6], and Rufibach et al.[7] provide
theoretical (i.e. based on simulated) as well as empirical (based on data from
17 RCTs) evidence for the bias of all estimators apart from the gold-standard
AJE , and also quantify this bias. This supports decade-long theoretical
investigations into the bias of, e.g, the one minus Kaplan-Meier estimator in
this field. We are of the strong opinion that all relevant stakeholders, among
them clinicians, statisticians, regulators, should collaborate to finally
implement fit-for-purpose methods to quantify AE risk, and update pertinent
guidelines. In our opinion, the implementation of the ICH E9(R1) estimands
addendum – so far primarily for efficacy – offers a window of opportunity to
push for a change also in reporting safety information. In drug development,
safety contributes as much to the benefit-risk assessment of a medicine as
efficacy, so the same estimand, estimation, and reporting standards should
apply to both.

###  6.2 SAVVY – template for sharing of summary data

A special feature of SAVVY was they way data from the 17 RCTs was shared and
analyzed: In a big collaborative effort data had been gathered within ten
sponsor organisations (nine pharmaceutical companies and one academic trial
center). In order to avoid challenges with data sharing SAVVY used an approach
familiar from Health Informatics, see e.g. Budin et al.[24]. A standardized
data structure was defined[10] based on which SAS and R macros were developed
by the academic project group members. These macros where then shared with all
participating sponsor organisations and run by them locally on their
individual trial data. Only aggregated data necessary for meta-analyses were
forwarded to the academic group members to centrally run meta-analyses.

This appears to be an approach that could also be applied in other
applications, as long as the analysis of interest can be done on summary data
provided by single organisations.

###  6.3 Roadmap for the SAVVY project and the analysis of safety data

Work within the SAVVY project continues. Concrete plans are an analysis
restricted to the oncology trials within the 17 RCTs, discussing in more
detail the issue of CEs in a typical oncology clinical trial. Collaborations
with clinicians in other therapeutic areas to define AEs of interest for which
”proper” estimation of their risk would be informative and what clinical
events to define as competing are envisaged.

SAVVY’s long-term vision is indeed to further familiarize trialists with the
AJE and have this method be recommended in future revisions of pertinent
guidelines. This, in connection with developing pragmatic approaches how to
properly identify and define CEs in therapeutic areas.

##  7 Abbreviations

AE: adverse event; CE: competing event; SAVVY: Survival analysis for AdVers
events with VarYing follow-up times; HR: hazard ratio, RR: relative risk; RCT:
randomized clinical trial

##  8 Declarations

##  9 Ethics approval and consent to participate

Not applicable.

##  10 Consent for publication

Not applicable.

##  11 Availability of data and materials

Individual trial data analyses were run within the sponsor organizations using
SAS and R software provided by the academic project group members. Only
aggregated data necessary for meta-analyses were shared and meta-analyses were
run centrally at the academic institutions. The SAVVY project has a webpage:
https://numbersman77.github.io/savvy. Methods are implemented in the R package
savvyr[23]. This package is available from CRAN.

##  12 Competing interests

KR is employee of F. Hoffmann-La Roche (Basel, Switzerland). JB has received
personal fees for consultancy from Pfizer and Roche, all outside the submitted
work. TF has received personal fees for consultancies (including data
monitoring committees) from Bayer, Boehringer Ingelheim, Janssen, Novartis and
Roche, all outside the submitted work. CS has received personal fees for
consultancies (including data monitoring committees) from Novartis and Roche,
all outside the submitted work. RS has declared no conflict of interest.

##  13 Funding

Not applicable.

##  14 Author’s contributions

KR, JB, TF, CS and RS conceived the idea for article and drafted it. All
authors critically reviewed the manuscript and approved its final version.

##  15 Acknowledgements

We thank Thomas Künzel and Daniel Sabanés Bové for implementing earlier code
into the R package savvyr[23].

## References

  * [1] O’Neill R. T.. Statistical analyses of adverse event data from clinical trials: Special emphasis on serious events. Drug Inf J.  1987;21:9–20.
  * [2] Proctor Tanja, Schumacher Martin. Analysing adverse events by time-to-event models: the CLEOPATRA study Pharmaceutical Statistics. ;15:306–314.
  * [3] Aalen Odd O, Johansen Søren. An empirical transition matrix for non-homogeneous Markov chains based on censored observations Scandinavian Journal of Statistics. :141–150.
  * [4] Allignol Arthur, Beyersmann Jan, Schmoor Claudia. Statistical issues in the analysis of adverse events in time-to-event data Pharm Stat.  2016;15:297–305.
  * [5] Stegherr Regina, Schmoor Claudia, Lübbert Michael, Friede Tim, Beyersmann Jan. Estimating and comparing adverse event probabilities in the presence of varying follow-up times and competing events Pharmaceutical Statistics.  2021c;20:1125–1146.
  * [6] Stegherr Regina, Schmoor Claudia, Beyersmann Jan, et al. Survival analysis for AdVerse events with VarYing follow-up times (SAVVY) – estimation of adverse event risks Trials.  2021b;22.
  * [7] Rufibach Kaspar, Stegherr Regina, Schmoor Claudia, et al. Comparison of Adverse Event Risks in Randomized Controlled Trials with Varying Follow-Up Times and Competing Events: Results from an Empirical Study Statistics in Biopharmaceutical Research. ;0:1-14.
  * [8] Gooley Ted A., Leisenring Wendy, Crowley John, Storer Barry E.. Estimation of Failure Probabilities in the Presence of Competing Risks: New Representations of Old Estimators Statistics in Medicine. ;18:695–706.
  * [9] ICH E9 working group . ICH E9 (R1): addendum on estimands and sensitivity analysis in clinical trials to the guideline on statistical principles for clinical trials 2019.
  * [10] Stegherr Regina, Beyersmann Jan, Jehl Valentine, et al. Survival analysis for AdVerse events with VarYing follow-up times (SAVVY): Rationale and statistical concept of a meta-analytic study. Biometrical Journal.  2021;63:650–670.
  * [11] EMA . A Guideline on Summary of Product Characteristics (SmPC) 2009.
  * [12] CIOMS Working Groups III and V . Guidelines for Preparing Core Clinical-Safety Information on Drugs Geneva:Council for International Organizations of Medical Sciences.
  * [13] IQWiG . General Methods, Version 5.0. Institute of Quality and Efficiency in Health Care .
  * [14] Latouche Aurelien, Allignol Arthur, Beyersmann Jan, Labopin Myriam, Fine Jason P. A competing risks analysis should report results on all cause-specific hazards and cumulative incidence functions Journal of Clinical Epidemiology. ;66:648–653.
  * [15] Schemper M, Smith T L. A note on quantifying follow-up in studies of failure time Controlled Clinical Trials. ;17:343–346.
  * [16] Clark Taane G, Altman Douglas G, De Stavola Bianca L. Quantification of the completeness of follow-up The Lancet. ;359:1309–1310.
  * [17] Varadhan Ravi, Weiss Carlos O, Segal Jodi B, Wu Albert W, Scharfstein Daniel, Boyd Cynthia. Evaluating health outcomes in the presence of competing risks: a review of statistical methods and clinical applications Medical care. ;48:S96–105.
  * [18] CIOMS Working Groups VI . Management of Safety Informationfrom Clinical Trials Geneva:Council for International Organizations of Medical Sciences 2005\.
  * [19] ICH Harmonised Tripartite Guideline . Structure and content of clinical study reports E3 1995.
  * [20] ICH Harmonised Tripartite Guideline . Statistical principals for clinical trials E9 1998.
  * [21] U.S. Food and Drug Administration . Premarketing Risk Assessment.
  * [22] Junqueira D. R., Zorzela L., Golder S., et al. CONSORT Harms 2022 statement, explanation, and elaboration: updated guideline for the reporting of harms in randomized trials J Clin Epidemiol. ;158:149–165.
  * [23] Kuenzel Thomas, Rufibach Kaspar, Stegherr Regina, Sabanés Bové Daniel. savvyr: Survival Analysis for AdVerse Events with VarYing Follow-Up Times.  R package version 0.1.0.
  * [24] Budin-Ljosne I., Burton P., Isaeva J., et al. DataSHIELD: an ethically robust solution to multiple-site individual-level data analysis Public Health Genomics. ;18:87–96.


# A Comparison Between Nonlinear Force-Free Field and Potential Field Models
Using Full-Disk SDO/HMI Magnetogram

Tilaye Tadesse1 1 Department of Physics, Drexel University, Philadelphia, PA
19104-2875, U.S.A, email: tasfaw@einstein.physics.drexel.edu,email:
<EMAIL_ADDRESS>
3 Max Planck Institut für Sonnensystemforschung, Max-Planck Str. 2, D–37191
Katlenburg-Lindau, Germany, email<EMAIL_ADDRESS>email:
<EMAIL_ADDRESS>
4 NASA, Goddard Space Flight Center, Code 674, Greenbelt, MD 20771, U.S.A.
email<EMAIL_ADDRESS>
5 National Solar Observatory, Sunspot, NM 88349, U.S.A. email:
<EMAIL_ADDRESS>
T. Wiegelmann3 1 Department of Physics, Drexel University, Philadelphia, PA
19104-2875, U.S.A, email: tasfaw@einstein.physics.drexel.edu,email:
<EMAIL_ADDRESS>
3 Max Planck Institut für Sonnensystemforschung, Max-Planck Str. 2, D–37191
Katlenburg-Lindau, Germany, email<EMAIL_ADDRESS>email:
<EMAIL_ADDRESS>
4 NASA, Goddard Space Flight Center, Code 674, Greenbelt, MD 20771, U.S.A.
email<EMAIL_ADDRESS>
5 National Solar Observatory, Sunspot, NM 88349, U.S.A. email:
<EMAIL_ADDRESS>
B. Inhester3 1 Department of Physics, Drexel University, Philadelphia, PA
19104-2875, U.S.A, email: tasfaw@einstein.physics.drexel.edu,email:
<EMAIL_ADDRESS>
3 Max Planck Institut für Sonnensystemforschung, Max-Planck Str. 2, D–37191
Katlenburg-Lindau, Germany, email<EMAIL_ADDRESS>email:
<EMAIL_ADDRESS>
4 NASA, Goddard Space Flight Center, Code 674, Greenbelt, MD 20771, U.S.A.
email<EMAIL_ADDRESS>
5 National Solar Observatory, Sunspot, NM 88349, U.S.A. email:
<EMAIL_ADDRESS>
K. Olson1 1 Department of Physics, Drexel University, Philadelphia, PA
19104-2875, U.S.A, email: tasfaw@einstein.physics.drexel.edu,email:
<EMAIL_ADDRESS>
3 Max Planck Institut für Sonnensystemforschung, Max-Planck Str. 2, D–37191
Katlenburg-Lindau, Germany, email<EMAIL_ADDRESS>email:
<EMAIL_ADDRESS>
4 NASA, Goddard Space Flight Center, Code 674, Greenbelt, MD 20771, U.S.A.
email<EMAIL_ADDRESS>
5 National Solar Observatory, Sunspot, NM 88349, U.S.A. email:
<EMAIL_ADDRESS>
P. J. MacNeice4 1 Department of Physics, Drexel University, Philadelphia, PA
19104-2875, U.S.A, email: tasfaw@einstein.physics.drexel.edu,email:
<EMAIL_ADDRESS>
3 Max Planck Institut für Sonnensystemforschung, Max-Planck Str. 2, D–37191
Katlenburg-Lindau, Germany, email<EMAIL_ADDRESS>email:
<EMAIL_ADDRESS>
4 NASA, Goddard Space Flight Center, Code 674, Greenbelt, MD 20771, U.S.A.
email<EMAIL_ADDRESS>
5 National Solar Observatory, Sunspot, NM 88349, U.S.A. email:
<EMAIL_ADDRESS>
A. Pevtsov5 1 Department of Physics, Drexel University, Philadelphia, PA
19104-2875, U.S.A, email: tasfaw@einstein.physics.drexel.edu,email:
<EMAIL_ADDRESS>
3 Max Planck Institut für Sonnensystemforschung, Max-Planck Str. 2, D–37191
Katlenburg-Lindau, Germany, email<EMAIL_ADDRESS>email:
<EMAIL_ADDRESS>
4 NASA, Goddard Space Flight Center, Code 674, Greenbelt, MD 20771, U.S.A.
email<EMAIL_ADDRESS>
5 National Solar Observatory, Sunspot, NM 88349, U.S.A. email:
<EMAIL_ADDRESS>

###### Abstract

Measurements of magnetic fields and electric currents in the pre-eruptive
corona are crucial to study solar eruptive phenomena, like flare and coronal
mass ejections(CMEs). However, spectro-polarimetric measurements of certain
photospheric lines permit a determination of the vector magnetic field at the
photosphere. Thus, substantial collection of magnetograms relate to the
photospheric surface field only. Therefore, there is considerable interest in
accurate modeling of the solar coronal magnetic field using photospheric
vector magnetograms as boundary data. This numerical modeling is carried out
by applying state-of-the-art nonlinear force-free field (NLFFF)
reconstruction. Cartesian nonlinear force-free field (NLFFF) codes are not
well suited for larger domains, since the spherical nature of the solar
surface cannot be neglected when the field of view is large. One of the most
significant results of Solar Dynamic Observatory (SDO) mission to date has
been repeated observations of large, almost global scale events in which large
scale connection between active regions may play fundamental role. Therefore,
it appears prudent to implement a NLFFF procedure in spherical geometry for
use when large scale boundary data are available, such as from the
Helioseismic and Magnetic Imager (HMI) on board SDO. In this work, we model
the coronal magnetic field above multiple active regions with the help of a
potential field and a NLFFF extrapolation codes in a full-disk using HMI data
as a boundary conditions. We compare projections of the resulting magnetic
field lines solutions with full-disk coronal images from the Atmospheric
Imaging Assembly (SDO/AIA) for both models. This study has found that the
NLFFF model reconstructs the magnetic configuration better than the potential
field model. We have concluded that much of trans-equatorial loops connecting
the two solar hemispheres are current-free.

###### keywords:

Active Regions, Magnetic Fields; Active Regions, Models; Magnetic fields,
Corona; Magnetic fields, Photosphere; Magnetic fields, Models

##  1 Introduction

Magnetic fields are believed to play a dominant role for active phenomena
carried out in the solar corona. The study of solar eruptive phenomena
requires that we understand how magnetic energy is stored in the pre-eruptive
corona. One of the key questions is: What is the three dimensional (3D)
structure of magnetic fields and electric currents in the pre-eruptive corona,
and how much free energy is stored in the field? [Schrijver and Title (2005),
Jiang and Feng (2012)]. To answer this question we must measure the magnetic
field in the coronal volume. However, spectro-polarimetric measurements of
certain photospheric lines permit a determination of the vector magnetic field
at the photosphere. Thus, substantial collection of magnetograms relate to the
photospheric surface field only. Although measurement of magnetic fields in
the chromosphere and the corona has considerably improved in recent decades
[Lin, Penn, and Tomczyk (2000), Liu and Lin (2008)], further developments are
needed before accurate data are routinely available. The problem of measuring
the coronal field and its embedded electrical currents thus leads us to use
numerical modeling to infer the field strength in the higher layers of the
solar atmosphere from the measured photospheric field. Due to the low value of
the plasma $\beta$ (the ratio of gas pressure to magnetic pressure), the solar
corona is magnetically dominated [Gary (2001)]. To describe the equilibrium
structure of the static coronal magnetic field when non-magnetic forces are
negligible, the force-free assumption is appropriate:

$(\nabla\times\textbf{B})\times\textbf{B}=0\ilabel{one}$ (1)
$\nabla\cdot\textbf{B}=0\ilabel{two}$ (2)

subject to the boundary condition

$\textbf{B}=\textbf{B}_{\textrm{obs}}\quad\mbox{on photosphere}\ilabel{three}$
(3)

where B is the magnetic field and $\textbf{B}_{\textrm{obs}}$ is measured
vector field on the photosphere. Equation (one) states that the Lorentz force
vanishes (as a consequence of $\textbf{J}\parallel\textbf{B}$, where J is the
electric current density) and Equation (two) describes the absence of magnetic
monopoles.

As an alternative to real measurements, force-free extrapolation of
photosphericmagnetic fields is currently being used as a the primary tool for
the modeling of coronal magnetic fields [Inhester and Wiegelmann (2006),
Valori, Kliem, and Keppens (2005), Wiegelmann (2004), Wheatland (2004),
Wheatland and Régnier (2009), Tadesse, Wiegelmann, and Inhester (2009),
Wheatland and Leka (2011), Amari and Aly (2010), Wiegelmann, T. and Thalmann,
J. K. and Inhester, B. and Tadesse, T. et al. (2012), Jiang and Feng (2012),
Jiang, Feng, and Xiang (2012), Aschwanden and Malanushenko (2012), Wiegelmann
and Sakurai (2012)].

Cartesian nonlinear force-free field (NLFFF) codes are not well suited for
larger domains, since the spherical nature of the solar surface cannot be
neglected when the field of view is large [Wiegelmann (2007), Tadesse,
Wiegelmann, and Inhester (2009), Tadesse et al. (2012a), Guo et al. (2012)].
One of the most significant results of Solar Dynamic Observatory (SDO) mission
to date has been repeated observations of large, almost global scale events in
which large scale connection between active regions may play fundamental role.
Therefore, it appears prudent to implement a NLFFF procedure in spherical
geometry for use when large scale boundary data are available, such as from
the Helioseismic and Magnetic Imager (HMI) on board SDO. DeRosa has studied
that different NLFFF models have markedly different field line configurations
and provided widely varying estimates of the magnetic free energy in the
coronal volume. The main reasons for that problem are (1) the forces acting on
the field within the photosphere, (2) the uncertainties on vector-field
measurements, particularly on the transverse component, and (3) the large
domain that needs to be modeled to capture the connections of an active region
to its surroundings [Wiegelmann (2007), Tadesse et al. (2011), Tadesse et al.
(2012b), Tadesse et al. (2012c)]. In this study, we have considered those
three points explicitly into account for modeling coronal magnetic field in a
full-disk.

In this work, we apply our spherical NLFFF procedure to a group of active
regions observed on November 09 2011 around 17:45UT by SDO/HMI. During this
observation, there were four active regions (ARs 11338, 11339, 11341 and
11342) along with other smaller sunspots spreading on the disk. We compare the
extrapolated (both potential and NLFFF) magnetic loops with extreme
ultraviolet (EUV) observations by the Atmospheric Imaging Assembly (AIA)
[Lemen, J. R. and Title, A. M. and Akin, D. J. et al. (2012)] on board SDO.
During comparison, we check whether the NLFFF model reconstructs the magnetic
configuration better than the potential field model. This comparison can be
used to evaluate how well the model field lines approximate the observed
coronal loops.

##  2 Instrumentation and data set

The Helioseismic and Magnetic Imager (HMI) [Schou, J. and Scherrer, P. H. and
Bush, R. I. et al. (2012)] is part of the Solar Dynamics Observatory (SDO) and
observes the full Sun at six wavelengths and full Stokes profile in the Fe
slowromancapi@ 617.3 nm spectral line. It consists of a refracting telescope,
a polarization selector, an image stabilization system, a narrow band tunable
filter and two 4096 pixel CCD cameras with mechanical shutters and control
electronics.

The transverse components of vector magnetic fields suffer from the so-called
$180^{\circ}$ ambiguity. The $180^{\circ}$ ambiguity for the HMI data in this
study has been resolved by an improved version of the minimum energy method
[Metcalf (1994), Metcalf, T. R. and Leka, K. D. and Barnes, G. et al. (2006),
Leka, Barnes, and Crouch (2009)]. As described in Leka:2009, in weak-field
areas, the minimization may not return a good solution due to large noise. The
noise level is $\approx$ 10G and $\approx$ 100G for the longitudinal and
transverse magnetic field, respectively. Therefore, in order to get a
spatially smooth solution in weak-field areas, we divide the magnetic field
into two regions, i.e., strong-field region and weak-field region, which is
defined to be where the field strength is below 200 G at the disk center, and
400 G on the limb. The values vary linearly with distance from the center to
the limb. The ambiguity solution in the strong-field area is derived by the
annealing in the minimization method and released in the data series
“hmi.B_720s_e15w1332_cutout“. Magnetic field azimuths in the weak-field area
are finally determined by the potential-field acute-angle method after the
annealing. This is different from the original minimum energy method, which
uses a neighboring-pixel acute-angle algorithm to revisit the weak field area.
Figure fig3a) shows the vector magnetic fields observed by HMI on November 09
2011 around 17:45UT.

##  3 Optimization principle and preprocessing of HMI data

Molodenskii69 and Aly89 pointed out that vector magnetic fields on a closed
surface that fully encloses any force-free domain have to satisfy the force-
free and torque-free conditions. To serve as suitable lower boundary condition
for a force-free modeling, on average a net tangential force acting on the
boundary and shear stresses along axes lying on the boundary have to reduce to
zero [Molodensky (1969), Aly (1989), Wiegelmann, Inhester, and Sakurai (2006),
Tadesse, Wiegelmann, and Inhester (2009)]. Preprocessing method as implemented
in Wiegelmann06sak has to be used to fulfill those conditions. The
preprocessing scheme of tilaye09 involves minimizing a two-dimensional
functional of quadratic form in spherical geometry similar to

$\textbf{B}=\emph{argmin}(L_{p}),$
$L_{p}=\mu_{1}L_{1}+\mu_{2}L_{2}+\mu_{3}L_{3}+\mu_{4}L_{4},\ilabel{3}$ (4)

where B is the preprocessed surface magnetic field from the input observed
field $\textbf{B}_{obs}$. Each of the constraints $L_{n}$ is weighted by an as
yet undetermined factor $\mu_{n}$. The first term $(n=1)$ corresponds to the
force-balance condition, the next $(n=2)$ to the torque-free condition, and
the last term $(n=4)$ controls the smoothing. The explicit form of $L_{1}$,
$L_{2}$, and $L_{4}$ can be found in tilaye09. The term $(n=3)$ controls the
difference between measured and preprocessed vector fields. After
preprocessing the HMI data, we solve the force-free equations (one) and (two)
using optimization principle [Wheatland, Sturrock, and Roumeliotis (2000),
Wiegelmann (2004)] in spherical geometry [Wiegelmann (2007), Tadesse,
Wiegelmann, and Inhester (2009), Tadesse et al. (2012c)].

##  4 Potential field model

The simplest way to model the coronal field is to assume that it is potential,
i.e. that it carries no electric current. Solutions for this model in plane
geometry have been obtained by schmidt, for the case where the vertical
component of the field is specified at the photospheric boundary. Potential
field model has now led to an almost routine type reconstruction, used for
observational purposes [Sakurai (1989)], but also for building initial
conditions for dynamical MHD numerical simulations [Amari et al. (1996)]. This
assumption has proven to be adequate for many quiescent, old active regions
and even for the non-eruptive global coronal-heliospheric interface [Wang and
Sheeley (1990), Schrijver and De Rosa (2003)].

Studies of the coupling of the coronal field into the heliosphere suggest that
the global coronal magnetic field is often largely potential. For the
practical calculation of the global field, the so-called source-surface model
has been introduced [Schatten, Wilcox, and Ness (1969)], in which the
influence of the solar wind is artificially taken account of by the
requirement that the field be radial at some exterior spherical (source)
surface typically at $R_{s}=2.5R_{\odot}$ from the sun’s center. The
potential-field source surface (PFSS) model, uses this concept to extrapolate
the line-of-sight surface magnetic field through the corona with the boundary
assumed to be at the source surface.

The potential field approximation has often been used to deduce the magnetic
structure of the solar corona from measurements of the photospheric field. It
obeys the equation $\nabla\times\textbf{B}=0$, so that the magnetic field can
be expressed as the gradient of a scalar potential, i.e.
$\textbf{B}=-\nabla\Phi$. Since the magnetic field is divergence-free
($\nabla\cdot\textbf{B}=0$), the scalar potential obeys the Laplace equation
$\nabla^{2}\Phi=0$. Together with the photospheric boundary condition, which
is traditionally provided by a map of the radial magnetic field, the Laplace
equation can be solved in the spherical coordinate system $(r,\theta,\phi)$,
where $\theta$ stands for colatitude. Using the radial magnetic field observed
at $R_{\odot}$ and the source surface assumption at $R_{s}$ as the boundary
condition, together with the spherical harmonic expansion in the domain
$R_{\odot}<r<R_{s}$, we finally obtain the potential field model.

##  5 Results

In this study, we have used potential field model and our spherical NLFFF
optimization procedure to a group of active regions observed on November 09
2011 around 17:45UT by SDO/HMI instrument. During this observation there were
four active regions (ARs 11338, 11339, 11341 and 11342, see Fig. fig2a) along
with other smaller sunspots spreading on the disk. In order to accommodate the
connectivity between those ARs and their surroundings, and to study trans-
equatorial loops (loops connecting north and south hemisphere), we adopt a non
uniform spherical grid $r$, $\theta$, $\phi$ with $n_{r}=225$,
$n_{\theta}=375$, $n_{\phi}=425$ grid points in the direction of radius,
latitude, and longitude, respectively, with the field of view of
$[r_{\mbox{min}}=1R_{\odot}:r_{\mbox{max}}=2R_{\odot}]\times[\theta_{\mbox{min}}=-50^{\circ}:\theta_{\mbox{max}}=50^{\circ}]\times[\phi_{\mbox{min}}=90^{\circ}:\phi_{\mbox{max}}=270^{\circ}]$.
The main purpose of this work is to study the magnetic and electric current
connectivities of the northern and southern solar hemispheres and to compare
which of the two models (potential or NLFFF model) is best suited to fully
describe those connectivities in full-disk environment.

We use our preprocessing routine in spherical geometry to derive suitable
boundary conditions for force-free modeling from the measured photospheric
data [Wiegelmann, Inhester, and Sakurai (2006), Tadesse, Wiegelmann, and
Inhester (2009)]. We used a standard preprocessing parameter set
$\mu_{1}=\mu_{2}=1$, $\mu_{3}=0.001$ and $\mu_{4}=0.05$, which are similar to
the values calculated from vector data used in previous studies [Wiegelmann,
T. and Thalmann, J. K. and Inhester, B. and Tadesse, T. et al. (2012), Tadesse
et al. (2012c)] for HMI data. The preprocessing influences the structure of
the magnetic vector data. It does not only smooths $\textbf{B}_{t}$
(transverse field) but also alters its values in order to reduce the net force
and torque. The change in $\textbf{B}_{t}$ is more pronounced than the radial
component $\textbf{B}_{r}$ (radial field) since $\textbf{B}_{t}$ is measured
with lower accuracy than the longitudinal magnetic field.

Figure 1.: Magnetic vector maps of HMI data on part of the lower boundary. The
color coding shows $B_{r}$ on the photosphere and the black arrow indicates
the transverse components of the field. The vertical and horizontal axes show
latitude, $\theta$ and longitude, $\phi$ on the photosphere respectively.fig1

For computing the potential field, we use the preprocessed radial component
$\textbf{B}_{r}$ of the HMI-data using spherical harmonic expansion. This
potential field solution has been used to initialize our NLFFF code. For
nonlinear force-free fields we minimize the functional Eq. (4). We implement
the new term $L_{photo}$ in Eq. (4) to work with boundary data of different
noise levels and qualities [Wiegelmann and Inhester (2010), Tadesse et al.
(2011)]. For those pixels, for which $\textbf{B}_{obs}$ was successfully
inverted, we allow deviations between the model field B and the input fields
observed $\textbf{B}_{obs}$ surface field using Eq. (4), so that the model
field can be iterated closer to a force-free solution even if the observations
are inconsistent. In order to control the speed with which the lower boundary
is injected during the NLFFF extrapolation, we have used the Langrangian
multiplier of $\nu=0.001$ as suggested by Tadesse:2013. For more details of
the method used in this work we direct the readers to the study by
Tilaye:2010.

We plot the surface potential field and NLFFF solutions in Figure fig1 for
selected region excluding the polar regions. In this figure, one can see that
there are substantial differences between the surface radial and transverse
field components of the two models. To quantify the degree of disagreement
between the respective vector components of potential and NLFFF model
solutions on the bottom surface, we calculate their pixel-wise correlations.
The correlation were calculated[Schrijver, C. J. and Derosa, M. L. and
Metcalf, T. R. et al. (2006), Metcalf, T. R. and De Rosa, M. L. and Schrijver,
C. J. et al. (2008)] from

$C_{\mathrm{vec}}=\frac{\sum_{i}\textbf{v}_{i}\cdot\textbf{u}_{i}}{\Big{(}\sum_{i}|\textbf{v}_{i}|^{2}\sum_{i}|\textbf{u}_{i}|^{2}\Big{)}^{1/2}},\ilabel{6}$
(5)

where $\textbf{v}_{i}$ and $\textbf{u}_{i}$ are the vectors at each grid point
$i$ on the bottom surface. If the vector fields are identical, then
$C_{vec}=1$; if $\textbf{v}_{i}\perp\textbf{u}_{i}$ , then $C_{vec}=0$. Table
table1 shows the correlation ($C_{vec}$) of the 2D surface magnetic field
vectors of potential and NLFFF models for the radial and transverse
components. The vector correlation between $\textbf{B}_{t}$ of the potential
and NLFFF surface vector maps is much more less than that of respective
$\textbf{B}_{r}$ components. This is due to the fact that the transverse
components ($\textbf{B}_{t}$) of the potential field model are computed using
spherical harmonics expansion from the radial components of the magnetic field
used as a boundary condition for NLFFF model. The average vector correlations
between the radial and transverse components of the two models in Table table1
indicates that the NLFFF solutions are somewhat far from potential as they
carry electric currents.

Table 1.: The correlations between the components of surface fields from potential and NLFFF models.table1 v | u | $C_{\mathrm{vec}}$
---|---|---
$(\textbf{B}_{pot})_{r}$ | $(\textbf{B}_{NLFFF})_{r}$ | $0.937$
$(\textbf{B}_{pot})_{t}$ | $(\textbf{B}_{NLFFF})_{t}$ | $0.813$
Figure 2.: a) Full-disk SDO/HMI magnetogram, the rectangles figure outline the
main active regions during this observation and the yellow line indicates the
solar equator ($\theta=0^{\circ}$), b) Full-disk AIA 171 Å image. Both data
sets were obtained on November 09, 2011 around 17:45UT, c) Field lines of the
potential field model at 17:45 UT overlaid on AIA 193Å image and d) Field
lines of the NLFFF model at 17:45 UT overlaid on AIA 193Å image. Green and red
lines represent open and closed magnetic field lines.fig2

In order to compare our reconstructions with observation, we plot the selected
fieldlines of the potential and NLFFF models in Figure fig2c and d. We overlay
the field lines with AIA 193 Å image. For observation time in Figure fig2c and
d, the field lines of the potential and NLFFF models are reconstructed from
the same footpoints. The potential field lines in Figure fig2(c) have an
obvious deviation from the observed EUV loops, since the projection of the
field lines at the bottom of AR 11339 leans toward the east, while the loops
toward the west. The loops connecting positive and negative polarities of AR
11341 are best overlaid by NLFFF lines than potential ones. The spatial
correspondence between the overall shape of the NLFFF field lines and the EUV
loops is much improved as shown in Figure fig2(d). Therefore, the qualitative
comparison between the model magnetic field lines and the observed EUV loops
indicates that the NLFFF model provides a more consistent field for full-disk
magnetic field reconstruction. In the absence of a more reliable quantitative
comparison, it remains the best option even if the qualitative nature is not
ideal.

In addition to the above comparison to quantify the degree of disagreement
between vector field solutions of the two models in the computational volume
that are specified on identical sets of grid points, we use five metrics that
compare either local characteristics (e.g., vector magnitudes and directions
at each point) or the global energy content in addition to the force and
divergence integrals as defined in Schrijver06. The vector correlation
($C_{vec}$) metric of Equation (6) is also used analogous to the standard
correlation coefficient for scalar functions. The second metric, $C_{CS}$ is
based on the Cauchy-Schwarz
inequality($|\textbf{a}\cdot\textbf{b}\leq|\textbf{a}||\textbf{b}|$ for any
vector a and b)

$C_{\rm CS}=\frac{1}{N}\sum_{i}\frac{{\bf B_{i}}\cdot{\bf b_{i}}}{|{\bf
B_{i}}||{\bf b_{i}}|},\ilabel{7}$ (6)

where $\textbf{B}=\textbf{B}_{\mbox{pot}}$,
$\textbf{b}=\textbf{B}_{\mbox{NLFFF}}$, and $N$ is the number of vectors in
the field. This metric is mostly a measure of the angular differences of the
vector fields: $C_{CS}=1$ when B and b are parallel and $C_{CS}=-1$ if they
are anti-parallel; $C_{CS}=0$ if $\textbf{B}_{i}\perp\textbf{b}_{i}$ at each
point. Next, we introduce two measures for the vector errors, one normalized
to the average vector norm, one averaging over relative differences. The
normalized vector error $E_{N}$ is defined as

$E_{\rm N}=\sum_{i}|{\bf b_{i}}-{\bf B_{i}}|/\sum_{i}|{\bf B_{i}}|,\ilabel{8}$
(7)

The mean vector error $E_{M}$ is defined as

$E_{\rm M}=\frac{1}{N}\sum_{i}\frac{|{\bf b_{i}}-{\bf B_{i}}|}{|{\bf
B_{i}}|}.\ilabel{9}$ (8)

Unlike the first two metrics, perfect disagreement of the two vector fields
results in $1-E_{M}=1-E_{N}=0$. As we are also interested in determining how
well the models estimate the energy contained in the field, we use the total
magnetic energy in NLFFF model field normalized to the total magnetic energy
in potential field as a global measure of the quality of the fit:

$\epsilon=\frac{\sum_{i}|{\bf b_{i}}|^{2}}{\sum_{i}|{\bf
B_{i}}|^{2}}.\ilabel{10}$ (9)

$\epsilon=1$, if there is no difference between the potential field and the
nonlinear force-free model solutions.

Table 2.: Comparision of figure of merits between potential field and NLFFF models. We have used pherical grids of $225\times 375\times 425$.table2 Model | $L_{f}$ | $L_{d}$ | $C_{\rm vec}$ | $C_{\rm CS}$ | $1-E_{N}$ | $1-E_{M}$ | $\epsilon$ | Time
---|---|---|---|---|---|---|---|---
Potential | $0.000$ | $0.002$ | $1$ | $1$ | $0$ | $0$ | $1$ | 30min
NLFFF | $0.591$ | $0.997$ | $0.878$ | $0.833$ | $0.735$ | $0.876$ | $1.220$ | 4h:39min

The degree of convergence towards a force-free and divergence-free model
solution can be quantified by the integral measures of the Lorentz force and
divergence terms in the minimization functional in Equation (4), computed over
the entire model volume V. $L_{f}$ and $L_{d}$ of Equation (4) measure how
well the force-free and divergence-free conditions are fulfilled,
respectively. In Table table2, we provide some quantitative measures to rate
the quality (figure of merit) of our reconstruction. Column $1$ names the
corresponding models. Columns $2-3$ show how well the force and solenoidal
condition are fulfilled for both models. The next five columns of Table table2
contain different measures which compare our NLFFF solution with potential
field. Those figures of merit indicate that there is a clear difference
between potential field and NLFFF models solutions for full-disk. The last
column shows the computing time on $1$ processor. The figures of merit show
that the potential field is far away from the true solutions and contains only
$81.97\%$ of the total magnetic energy.

We estimate the free magnetic energy to be the difference between the
extrapolated NLFFF and the potential field with the same normal boundary
conditions in the photosphere[Régnier and Priest (2007), Thalmann, Wiegelmann,
and Raouafi (2008)]. We therefore estimate the upper limit to the free
magnetic energy associated with coronal currents of the form

$E_{\mathrm{free}}=\frac{1}{8\pi}\int_{V}\Big{(}B_{nlff}^{2}-B_{pot}^{2}\Big{)}r^{2}sin\theta
drd\theta d\phi,\ilabel{ten}$ (10)

Our result for the estimation of free-magnetic energy in Table table3 shows
that the NLFFF model has $18.83\%$ more energy than the corresponding
potential field model. Therefore, whenever we try to analysis magnetic field
configuration and its magnetic field contents in full-disk, it is inaccurate
to use potential field model as it undermines the real features of the field
in the corona.

Table 3.: The magnetic energy associated with extrapolated potential and NLFFF field configurations from full disk SDO/HMI data.table3 Model | $E_{\mbox{total}}(10^{33}{\mbox{erg}})$ | $E_{\mbox{free}}(10^{33}{\mbox{erg}})$
---|---|---
Potential | $7.306$ | $0$
NLFFF | $8.913$ | $1.607$

In our previous work [Tadesse et al. (2012a)], we have studied the
connectivity between three neighbouring active regions. In this work, we study
the magnetic and electric connectivities between active regions in northern
and southern solar hemispheres using full-disk HMI data. During this
observation, there were three active regions in the northern hemisphere and
one active region with surrounding sunspots in the southern hemisphere. As a
result the total unsigned magnetic field flux in the northern hemisphere is
larger than that of the southern one. In order to quantify these
connectivities, we have calculated the magnetic flux and the electric currents
shared between those active regions. For the magnetic flux, e.g., we use

$\Phi_{\alpha\beta}=\sum_{i}|\textbf{B}_{i}\cdot\hat{r}|R^{2}_{\odot}\textrm{sin}(\theta_{i})\Delta\theta_{i}\Delta\phi_{i}\ilabel{ten}$
(11)

where the summation is over all pixels of $\mbox{AR}_{\alpha}$ from which the
field line ends in $\mbox{AR}_{\beta}$ or
$i\in\mbox{AR}_{\alpha}\|\,\mbox{conjugate footpoint}(i)\in\mbox{AR}_{\beta}$.
The indices $\alpha$ and $\beta$ denote the active region. For the electric
current we replace the magnetic field, $\textbf{B}_{i}\cdot\hat{r}$, by the
vertical current density $\textbf{J}_{i}\cdot\hat{r}$ in Equation (ten). Table
table4 shows the percentage of the total magnetic flux and electric current
shared between the ARs 11339 and 11338 and between AR 11341 and sunspots in
the southern hemisphere (see Fig. fig3). We have calculated total an unsigned
flux for each active regions 11339 and 11341 (both being in northern
hemisphere) and the flux due to those field lines ending in the AR 11338 and
magnetic patches in the southern hemisphere, respectively. The first column of
Table table4 shows that $7.31\%(6.17\%)$ of positive/negative polarity of AR
11339 (11341) in the northern hemisphere is connected to positive/negative
polarity of AR 11338 (magnetic patches) in southern hemisphere for potential
field configuration. The percentage share between those ARs are
$7.57\%(6.91\%)$ for NLFFF configuration. In Table table4, the percentages in
the bracket show the share in electric current between the corresponding
active regions and small magnetic patches. Figure fig3 display the trans-
equatorial loops connecting those active regions. In the figure, we rotate
those active regions to the limb to show loops connecting the two hemispheres.
The surface contour plots in Figures fig4 and fig5 show the magnetic flux and
electric current density flux of field lines crossing equatorial plane,
respectively. From those figures we can see that NLFFF model has more magnetic
and electric current density fluxes compared to the potential field. Figure
fig5 shows that there are more electric currents crossing the plane carried by
trans-equatorial loops connecting AR 11341 in northern hemisphere to sunspots
in the south than those currents connecting ARs 11339 and 11338. This indicate
that much of the trans-equatorial loops connecting ARs 11339 and 11338 are
potential. In general the two hemispheres are more magnetically connected than
electric current. However, one has to use NLFFF model over potential model to
study magnetic field structure in any phenomena involving the field.

Figure 3.: a) Selected potential field lines connecting ARs 11338 and 11339 rotated to limb. b) The corresponding field lines for NLFFF in a), c) Selected potential field lines connecting AR 11341 and Sun spots in southern hemisphere rotated to limb d) The corresponding field lines for NLFFF in c).fig3 Figure 4.: Surface contour plot of perpendicular component of magnetic fields crossing a plane of the solar equator ($\theta=0^{\circ}$). The color coding shows $B_{\perp}$ crossing the $\theta=0^{\circ}$ equatorial plane. The vertical and horizontal axes show height, $r$(in solar radius) and longitude, $\phi$(in degree) on the plane.fig4 Figure 5.: Surface contour plot of perpendicular component of electric current density crossing a plane of the solar equator ($\theta=0^{\circ}$). The color coding shows $J_{\perp}$ crossing the $\theta=0^{\circ}$ equatorial plane. The vertical and horizontal axes show height, $r$(in solar radius) and longitude, $\phi$(in degree) on the plane.fig5 Table 4.: The percentage of the total magnetic flux and electric current density shared between the ARs 11339 and 11338 and between AR 11341 and sunspots in the southern hemisphere ( see Fig. fig3). The percentage in brackets denote that of electric currents.table4 Model | $\mbox{AR 11339}\longmapsto\mbox{AR 11338}$ | $\mbox{AR 11341}\longmapsto\mbox{Patches in South}$
---|---|---
Potential | $7.31\,(0.0)$ | $6.17\,(0.01)$
NLFFF | $7.57\,(0.13)$ | $6.91\,(3.35)$

##  6 Conclusion and outlook

Both potential and nonlinear force-free field (NLFFF) codes in Cartesian
geometry are not well suited for larger domains, since the spherical nature of
the solar surface cannot be neglected when the field of view is large. One of
the most significant results of Solar Dynamic Obseratory (SDO) mission to date
has been repeated observations of large, almost global scale events in which
large scale connection between active regions may play fundamental role.
Therefore, it appears prudent to implement a NLFFF procedure in spherical
geometry for use when large scale boundary data are available, such as from
the Helioseismic and Magnetic Imager (HMI) on board SDO.

In this study, we have investigated the coronal magnetic field associated with
full solar disk on November 09 2011 by analysing SDO/HMI data using potential
and NLFFF models. During this particular observation, there were three active
regions in the northern hemisphere and one active region surrounded by
sunspots in the south. We have used our spherical NLFFF and potential codes to
compute the magnetic field solutions in full-disk. For computing the potential
field, we use the preprocessed radial component $\textbf{B}_{r}$ of the HMI-
data using spherical harmonic expansion. We implement our NLFFF code
initialized by the this potential field solution (except the lower observed
bottom boundary) during relaxation towards force-freeness state in the
computational volume.

We have compared the magnetic field solutions from both models. The
qualitative comparison between the model magnetic field lines and the observed
EUV loops indicates that the NLFFF model provides a more consistent field for
full-disk magnetic field reconstruction. In addition to this, the figures of
merits have been used to quantify the disagreements between the two models
field lines. The figures of merit indicate that the magnetic field lines
structures are not potential in full-disk field-of-views. However, much of the
trans-equatorial field lines connecting the active regions in the northern and
southern hemispheres are potential. This indicates that the two solar
hemispheres are more magnetically connected than electric current. The
magnetic field lines obtained from nonlinear force-free extrapolation bear
larger total magnetic energy that of the potential field model. Therefore, one
has to use NLFFF model over potential model to study magnetic field structure
in any phenomena involving the field in full-disk.

## Acknowledgements

Data are courtesy of NASA/SDO and the AIA and HMI science teams. SOLIS/VSM
vector magnetograms are produced cooperatively by NSF/NSO and NASA/LWS. The
National Solar Observatory (NSO) is operated by the Association of
Universities for Research in Astronomy, Inc., under cooperative agreement with
the National Science Foundation. This work was supported by NASA grant
NNX07AU64G and the work of T. Wiegelmann was supported by DLR-grant $50$ OC
$453$ $0501$.

## References

  * Aly (1989) Aly, J.J.: 1989,. Sol. Phys. 120, 19\.
  * Amari and Aly (2010) Amari, T., Aly, J.: 2010,. A&A 522, A52.
  * Amari et al. (1996) Amari, T., Luciani, J.F., Aly, J.J., Tagger, M.: 1996,. ApJL 466, L39.
  * Aschwanden and Malanushenko (2012) Aschwanden, M.J., Malanushenko, A.: 2012,. Sol. Phys., 182\. doi:10.1007/s11207-012-0070-1.
  * DeRosa, M. L. and Schrijver, C. J. and Barnes, G. et al. (2009) DeRosa, M. L. and Schrijver, C. J. and Barnes, G. et al.: 2009,. ApJ 696, 1780\.
  * Gary (2001) Gary, G.A.: 2001,. Sol. Phys. 203, 71\.
  * Guo et al. (2012) Guo, Y., Ding, M.D., Liu, Y., Sun, X.D., DeRosa, M.L., Wiegelmann, T.: 2012,. ApJ 760, 47\.
  * Inhester and Wiegelmann (2006) Inhester, B., Wiegelmann, T.: 2006,. Sol. Phys. 235, 201\.
  * Jiang and Feng (2012) Jiang, C., Feng, X.: 2012,. ApJ 749, 135\.
  * Jiang, Feng, and Xiang (2012) Jiang, C., Feng, X., Xiang, C.: 2012,. ApJ 755, 62\.
  * Leka, Barnes, and Crouch (2009) Leka, K.D., Barnes, G., Crouch, t.: 2009,. Sol. Phys. 260, 83\.
  * Lemen, J. R. and Title, A. M. and Akin, D. J. et al. (2012) Lemen, J. R. and Title, A. M. and Akin, D. J. et al.: 2012,. Sol. Phys. 275, 17\.
  * Lin, Penn, and Tomczyk (2000) Lin, H., Penn, M.J., Tomczyk, S.: 2000,. ApJL 541, L83.
  * Liu and Lin (2008) Liu, Y., Lin, H.: 2008,. ApJ 680, 1496\.
  * Metcalf (1994) Metcalf, T.R.: 1994,. Sol. Phys. 155, 235\.
  * Metcalf, T. R. and De Rosa, M. L. and Schrijver, C. J. et al. (2008) Metcalf, T. R. and De Rosa, M. L. and Schrijver, C. J. et al.: 2008,. Sol. Phys. 247, 269\.
  * Metcalf, T. R. and Leka, K. D. and Barnes, G. et al. (2006) Metcalf, T. R. and Leka, K. D. and Barnes, G. et al.: 2006,. Sol. Phys. 237, 267\.
  * Molodensky (1969) Molodensky, M.M.: 1969,. Soviet Astron.-AJ 12, 585\.
  * Régnier and Priest (2007) Régnier, S., Priest, E.R.: 2007,. ApJL 669, L53.
  * Sakurai (1989) Sakurai, T.: 1989,. Space Sci. Rev. 51, 11\.
  * Schatten, Wilcox, and Ness (1969) Schatten, K.H., Wilcox, J.M., Ness, N.F.: 1969,. Sol. Phys. 6, 442\.
  * Schmidt (1964) Schmidt, H.U.: 1964, In: The Physics of Solar Flares, 107\.
  * Schou, J. and Scherrer, P. H. and Bush, R. I. et al. (2012) Schou, J. and Scherrer, P. H. and Bush, R. I. et al.: 2012,. Sol. Phys. 275, 229\.
  * Schrijver and De Rosa (2003) Schrijver, C.J., De Rosa, M.L.: 2003,. Sol. Phys. 212, 165\.
  * Schrijver and Title (2005) Schrijver, C.J., Title, A.M.: 2005,. ApJ 619, 1077\.
  * Schrijver, C. J. and Derosa, M. L. and Metcalf, T. R. et al. (2006) Schrijver, C. J. and Derosa, M. L. and Metcalf, T. R. et al.: 2006,. Sol. Phys. 235, 161\.
  * Tadesse, Wiegelmann, and Inhester (2009) Tadesse, T., Wiegelmann, T., Inhester, B.: 2009,. A&A 508, 421\.
  * Tadesse et al. (2011) Tadesse, T., Wiegelmann, T., Inhester, B., Pevtsov, A.: 2011,. A&A 527, A30.
  * Tadesse et al. (2012a) Tadesse, T., Wiegelmann, T., Inhester, B., Pevtsov, A.: 2012a,. Sol. Phys. 277, 119\.
  * Tadesse et al. (2012b) Tadesse, T., Wiegelmann, T., Inhester, B., Pevtsov, A.: 2012b,. Sol. Phys. 281, 53\.
  * Tadesse et al. (2012c) Tadesse, T., Wiegelmann, T., Inhester, B., MacNeice, P., Pevtsov, A., Sun, X.: 2012c, Full-disk nonlinear force-free field extrapolation of SDO/HMI and SOLIS/VSM magnetograms. ArXiv e-prints.
  * Thalmann, Wiegelmann, and Raouafi (2008) Thalmann, J.K., Wiegelmann, T., Raouafi, N.E.: 2008,. A&A 488, L71.
  * Valori, Kliem, and Keppens (2005) Valori, G., Kliem, B., Keppens, R.: 2005,. A&A 433, 335\.
  * Wang and Sheeley (1990) Wang, Y.M., Sheeley, N.R. Jr.: 1990,. ApJ 365, 372\.
  * Wheatland (2004) Wheatland, M.S.: 2004,. Sol. Phys. 222, 247\.
  * Wheatland and Leka (2011) Wheatland, M.S., Leka, K.D.: 2011,. ApJ 728, 112\.
  * Wheatland and Régnier (2009) Wheatland, M.S., Régnier, S.: 2009,. ApJL 700, L88.
  * Wheatland, Sturrock, and Roumeliotis (2000) Wheatland, M.S., Sturrock, P.A., Roumeliotis, G.: 2000,. ApJ 540, 1150\.
  * Wiegelmann (2004) Wiegelmann, T.: 2004,. Sol. Phys. 219, 87 – 108.
  * Wiegelmann (2007) Wiegelmann, T.: 2007,. Sol. Phys. 240, 227 – 239.
  * Wiegelmann and Inhester (2010) Wiegelmann, T., Inhester, B.: 2010,. A&A 516, A107.
  * Wiegelmann and Sakurai (2012) Wiegelmann, T., Sakurai, T.: 2012,. Living Reviews in Sol. Phys. 9, 5\.
  * Wiegelmann, Inhester, and Sakurai (2006) Wiegelmann, T., Inhester, B., Sakurai, T.: 2006,. Sol. Phys. 233, 215\.
  * Wiegelmann, T. and Thalmann, J. K. and Inhester, B. and Tadesse, T. et al. (2012) Wiegelmann, T. and Thalmann, J. K. and Inhester, B. and Tadesse, T. et al.: 2012,. Sol. Phys. 281, 37\.