# A New Magnetic Parameter of Active Regions Distinguishing Large Eruptive and
Confined Solar Flares

Ting Li11affiliation: CAS Key Laboratory of Solar Activity, National
Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101,
China<EMAIL_ADDRESS>22affiliation: School of Astronomy and Space Science,
University of Chinese Academy of Sciences, Beijing 100049, China , Xudong
Sun33affiliation: Institute for Astronomy, University of Hawai‘i at Mānoa,
Pukalani, HI 96768, USA , Yijun Hou11affiliation: CAS Key Laboratory of Solar
Activity, National Astronomical Observatories, Chinese Academy of Sciences,
Beijing 100101, China<EMAIL_ADDRESS>22affiliation: School of Astronomy
and Space Science, University of Chinese Academy of Sciences, Beijing 100049,
China , Anqin Chen44affiliation: Key Laboratory of Space Weather, National
Center for Space Weather, China Meteorological Administration, Beijing 100081,
China , Shuhong Yang11affiliation: CAS Key Laboratory of Solar Activity,
National Astronomical Observatories, Chinese Academy of Sciences, Beijing
100101, China<EMAIL_ADDRESS>22affiliation: School of Astronomy and Space
Science, University of Chinese Academy of Sciences, Beijing 100049, China &
Jun Zhang55affiliation: School of Physics and Optoelectronics engineering,
Anhui University, Hefei 230601, China

###### Abstract

With the aim of investigating how the magnetic field in solar active regions
(ARs) controls flare activity, i.e., whether a confined or eruptive flare
occurs, we analyze 106 flares of Geostationary Operational Environmental
Satellite (GOES) class $\geq$M1.0 during 2010$-$2019\. We calculate mean
characteristic twist parameters $\alpha$FPIL within the “flaring polarity
inversion line” region and $\alpha$HFED within the area of high photospheric
magnetic free energy density, which both provide measures of the
nonpotentiality of AR core region. Magnetic twist is thought to be related to
the driving force of electric current-driven instabilities, such as the
helical kink instability. We also calculate total unsigned magnetic flux
($\Phi$AR) of ARs producing the flare, which describes the strength of the
background field confinement. By considering both the constraining effect of
background magnetic fields and the magnetic non-potentiality of ARs, we
propose a new parameter $\alpha$/$\Phi$AR to measure the probability for a
large flare to be associated with a coronal mass ejection (CME). We find that
in about 90% of eruptive flares, $\alpha$FPIL/$\Phi$AR and
$\alpha$HFED/$\Phi$AR are beyond critical values (2.2$\times$$10^{-24}$ and
3.2$\times$$10^{-24}$ Mm-1 Mx-1), whereas they are less than critical values
in $\sim$ 80% of confined flares. This indicates that the new parameter
$\alpha$/$\Phi$AR is well able to distinguish eruptive flares from confined
flares. Our investigation suggests that the relative measure of magnetic
nonpotentiality within the AR core over the restriction of the background
field largely controls the capability of ARs to produce eruptive flares.

Sun: activity—Sun: coronal mass ejections (CMEs)—Sun: flares

##  1 Introduction

Solar flares and coronal mass ejections (CMEs) are the rapid release of a huge
amount of magnetic energy accumulated in the solar corona through
magnetohydrodynamic instabilities and magnetic reconnection. Solar flares are
often, but not always, accompanied by CMEs. We refer to flares with a CME as
“eruptive flares” and flares not associated with a CME as “confined flares”.
It is revealed that flare-CME association rate increases with the flare
intensity (Andrews 2003; Yashiro et al. 2006). Recently, Li et al. (2020,
2021) found that flare-CME association rate decreases with total unsigned
magnetic flux ($\Phi$AR) of active regions (ARs) producing the flare, which
provides a global parameter relating to the strength of the background field
confinement.

Over the last two decades, significant progress has been made in understanding
the physical factors determining whether a flare event is associated with a
CME or not. It is suggested that eruptive flares tend to occur if the
overlying background magnetic fields are weaker or more quickly decay with
height (Török & Kliem 2005; Wang & Zhang 2007; Wang et al. 2017; Baumgartner
et al. 2018; Amari et al. 2018; Jing et al. 2018; Duan et al. 2019). Moreover,
the magnetic non-potentiality of ARs is thought to be another important factor
governing the eruptive character of solar flares (Nindos & Andrews 2004; Liu
et al. 2016; Cui et al. 2018; Vasantharaju et al. 2018; Thalmann et al. 2019;
Avallone & Sun 2020; Gupta et al. 2021), such as free magnetic energy,
relative helicity, magnetic twists, etc. Statistical studies have shown that
CME productivity is correlated with the twist parameter $\alpha$ (Falconer et
al. 2002, 2006), which characterizes the degree to which the photospheric
magnetic fields of an AR deviate from a potential field (Leka & Skumanich
1999; Yang et al. 2012). Bobra & Ilonidis (2016) found that the twist
parameter $\alpha$ and mean gradient of the horizontal field are two
relatively high-performing features in predicting CMEs based on machine-
learning algorithms. However, the measurement of previous known parameters
showed that there are high degrees of overlap between confined and eruptive
flares if only one factor (overlying confinement or magnetic non-potentiality)
was considered in selecting the parameters (Nindos & Andrews 2004; Wang et al.
2017; Vasantharaju et al. 2018).

Sun et al. (2015) suggested that AR eruptivity is related to the relative
value of magnetic nonpotentiality over the restriction of the background
field. In this Letter, following the idea of Sun et al. (2015), we consider
both the constraining effect of background magnetic fields and the magnetic
non-potentiality of ARs, and propose a new parameter $\alpha$/$\Phi$AR to
describe the eruptive character of a flare. By measuring the parameter of 106
large flares within the cores of 12 CME-active and 9 CME-quiet ARs, we find
that the new parameter $\alpha$/$\Phi$AR is well able to distinguish flares
associated with CMEs from flares that are not.

##  2 Database Selection and Parameter Calculations

We use a subset of 106 flare events111https://doi.org/10.12149/101087
$\geq$M1.0 (43 eruptive and 63 confined) from a large database of 322 M-class
flares222https://doi.org/10.12149/101030 during the period of June 2010 to
June 2019 (Li et al. 2020). The subset is selected based on the
characteristics of ARs, and the selected ARs must fulfill the following two
selection criteria. First, the ARs are flare-active and produced $\geq$3
M-class flares. Second, the ARs can be unambiguously classified into CME-
active and CME-quiet. We refer to an AR as a CME-active (CME-quiet) AR if the
flares from it are all eruptive (confined) or only one exceptional flare is
confined (eruptive). Finally, a total of 21 ARs fulfill these selection
criteria, including 12 CME-active ARs and 9 CME-quiet ARs.

For each event, based on the vector magnetograms from Space-Weather
Helioseismic and Magnetic Imager (HMI; Scherrer et al. 2012) AR Patches
(SHARP; Bobra et al. 2014) observed by _Solar Dynamics Observatory_ (_SDO_ ;
Pesnell et al. 2012), we calculate $\Phi$AR before the flare onset by summing
all pixels where vertical magnetic field $|$$B_{z}$$|$$>$100 G (Kazachenko et
al. 2017). The magnetograms are re-mapped using a cylindrical equal area (CEA)
projection with a pixel size of $\sim$0$\arcsec$.5 and presented as ($B_{r}$,
$B_{\theta}$, $B_{\phi}$) in heliocentric spherical coordinates corresponding
to ($B_{z}$, -$B_{y}$, $B_{x}$) in heliographic coordinates (Sun 2013). We
have identified a “flaring polarity inversion line” (FPIL) mask to demarcate
the core of an AR by using the method of Sun et al. (2015). We first find the
polarity inversion line (PIL) pixels from a smoothed vertical magnetic field
$B_{z}$, and dilate them with a circular kernel with a radius of 18 pixels
(about 6.5 Mm; other radii are also used and the results are not affected).
Then we isolate flare ribbons by using the 1600 Å image near the flare peak
from the Atmospheric Imaging Assembly (AIA; Lemen et al. 2012) on board the
_SDO_ (above 700 DN $s^{-1}$ which is about 10 times the standard deviation
above the mean of the quiet-Sun values), and dilate them with a large kernel
having a radius of 20 pixels (about 7.2 Mm). Finally the intersection of
dilated PIL and ribbon areas is considered as the FPIL mask. The FPIL mask is
similar to the strong-field, high-gradient PIL mask in Schrijver (2007),
however, it only involves part of the PIL mask.

We calculate distributions of vertical electric current density $J_{z}$ and
the mean characteristic twist parameter $\alpha$FPIL within the FPIL mask
region. We note that the FPIL mask region could be determined only after the
flare occurrence and thus $\alpha$FPIL can not be used for CME forecasting.
The energy release during a flare is generally believed to originate from
areas with high values of photospheric magnetic free energy density
$\rho_{\mathrm{free}}$ in ARs, and thus we select the region with
$\rho_{\mathrm{free}}$$>$4.0$\times$104 erg cm-3 (HFED region; Chen & Wang
2012) as a proxy for AR core region and calculate mean characteristic twist
parameter $\alpha$HFED and mean shear angle $\Psi$HFED within HFED region.
Compared with $\alpha$FPIL, $\alpha$HFED does not need ribbon information and
is more suitable for CME forecasting. Detailed formulas of the parameters are
listed in Table 1.

##  3 Statistical Results

Figure 1 shows four examples of two eruptive and two confined flares, which
includes 1600 Å images (left), photospheric magnetograms (middle) and derived
vertical electric current density $J_{z}$ maps (right). We can see that for
two eruptive flares (X2.1 flare in AR 11283 and M6.5 in AR 12371) the positive
and negative currents have a coherent structure around the PIL (Figures
1(a)-(b)), indicating the presence of “current ribbons” as in a coherent flux
rope. However, the two confined events (X3.1 in AR 12192 and M6.1 in AR 12222)
exhibit disordered current distributions and do not have any noticeable
structure (Figures 1(c)-(d)). It also can been seen that the FPIL mask regions
(orange and black contours) overlap the area of large currents, implying that
the FPIL mask corresponds to the AR core with the strongest magnetic non-
potentiality. We estimate the errors of $J_{z}$ to be about 10 mA m-2 based on
the noise level of transverse magnetic field ($\sim$ 100 G in Liu et al.
2012).

Based on the derived $J_{z}$ map, we then calculated the mean characteristic
twist parameter $\alpha$FPIL within the FPIL mask region for the 106 flares.
Figure 2(a) shows the scatter plot of $\alpha$FPIL versus $\Phi$AR. Blue (red)
circles are the eruptive (confined) flares. It needs to be noted that
$\alpha$FPIL is a signed parameter and in our study $\alpha$FPIL means its
absolute value. It can be seen that the events with $\alpha$FPIL$<$0.07 Mm-1
(black dotted line in Figure 2(a)) and $\Phi$AR$>$1.0$\times$$10^{23}$ Mx
(right green dashed line in Figure 2(a)) are all confined and those with
$\alpha$FPIL$\geq$0.07 Mm-1 and $\Phi$AR$<$3.5$\times$$10^{22}$ Mx (left green
dashed line in Figure 2(a)) are all eruptive. According to error propagation
theory, we estimate the errors of $\alpha$FPIL to be $10^{-4}$$-$$10^{-3}$
Mm-1 (by considering the error of $J_{z}$ $\sim$ 10 mA m-2). The errors are
small and they do not affect our results. The noise level of $B_{z}$ is on the
order of 10 G (Liu et al. 2012), and the errors in the calculation of $\Phi$AR
are estimated to be about $10^{18}$$-$$10^{19}$ Mx. The errors are much
smaller than $\Phi$AR and thus are not considered. Figure 2(b) shows the
scatter plot of flare peak X-ray flux versus $\alpha$FPIL/$\Phi$AR. About 93%
(40 of 43) of eruptive events have
$\alpha$FPIL/$\Phi$AR$\geq$2.2$\times$$10^{-24}$ Mm-1 Mx-1, and $\sim$ 83% (52
of 63) of confined flares have $\alpha$FPIL/$\Phi$AR$<$2.2$\times$$10^{-24}$
Mm-1 Mx-1 (black dash-dotted line in Figure 2(b)). This shows that the new
relative non-potential parameter $\alpha$FPIL/$\Phi$AR is well able to
distinguish the two populations of eruptive and confined flares. The
statistical comparison between the two distributions in Figure 2(b) was done
with a two-sample t-test, which shows a significant difference (P$<$0.001).
However, there is still a small overlap (14 of 106) by using our criterion and
the exceptional events are mainly from a CME-active AR 11302 and two CME-quiet
ARs 11476 and 12268.

Figure 3 shows the maps of photospheric magnetic free energy density
$\rho_{\mathrm{free}}$ and magnetic shear angle $\Psi$ for the four examples
shown in Figure 1. Free magnetic energy is the amount of magnetic energy in
excess of the minimum energy attributed to the potential field, and magnetic
shear is defined as the angle between the horizontal components of the
observed magnetic field and a modeled potential magnetic field based on
photospheric $B_{z}$ map. They are commonly used parameters in describing the
magnetic complexity and non-potentiality (Wang et al. 1994; Su et al. 2014).
It can be seen that the maps of $\rho_{\mathrm{free}}$ and $\Psi$ exhibit
similar distributions, with their large values around the PILs of ARs. The
comparison of two eruptive flares (Figures 3(a)-(b)) with two confined ones
(Figures 3(c)-(d)) shows that the mean values of $\rho_{\mathrm{free}}$ and
$\Psi$ of eruptive flares are larger than those of confined events. Similar to
the appearance of $J_{z}$ maps (Figure 1), $\rho_{\mathrm{free}}$ and $\Psi$
of eruptive flares exhibit ribbon patterns, however, the $\rho_{free}$ and
$\Psi$ maps of confined flares show disordered distributions.

We determine high-$\rho_{\mathrm{free}}$ areas (HFED region) as a proxy for AR
core region and then calculate mean characteristic twist parameter
$\alpha$HFED and mean shear angle $\Psi$HFED within the HFED region. Figure 4
shows the calculation results for 43 eruptive and 63 confined flares. It can
be seen that the distributions of $\alpha$HFED (Figure 4(a)) are similar to
those of $\alpha$FPIL (Figure 2(a)). For $\alpha$HFED$<$0.1 Mm-1 (black dotted
line in Figure 4(a)), an overwhelming majority (about 93%, 26 out of 28) of
flares are confined. Almost all the eruptive flares (41 out of 43) have
$\alpha$HFED$\geq$0.1 Mm-1. If we consider the relative parameter
$\alpha$HFED/$\Phi$AR, the differences between confined and eruptive cases are
more evident (Figure 4(b)). About 91% (39 of 43) of eruptive flares have
$\alpha$HFED/$\Phi$AR$\geq$3.2$\times$$10^{-24}$ Mm-1 Mx-1, and $\sim$ 75% (47
of 63) of confined events have $\alpha$HFED/$\Phi$AR$<$3.2$\times$$10^{-24}$
Mm-1 Mx-1 (black dash-dotted line in Figure 4(b)). There is also a difference
of mean shear angle $\Psi$HFED within HFED region between eruptive and
confined flares (Figure 4(c)). The proportion of confined flares is $\sim$91%
(29 of 32) corresponding to $\Psi$HFED$<$60∘. Similarly, the relative
parameter $\Psi$HFED/$\Phi$AR can provide a good ability for distinguishing
the eruptive and confined events (Figure 4(d)). An overwhelming majority of
eruptive flares have $\Psi$HFED/$\Phi$AR$\geq$1.0$\times$$10^{-21}$ degree
Mx-1 and most of confined flares show
$\Psi$HFED/$\Phi$AR$<$1.0$\times$$10^{-21}$ degree Mx-1. The degree of overlap
between eruptive and confined events for $\Psi$HFED/$\Phi$AR is a little
higher than that for $\alpha$HFED/$\Phi$AR. The two-sample t-test was carried
out for the two parameters in Figures 4(b) and (d), and showed a significant
difference (P$<$0.001).

##  4 Summary and Discussion

In this study, we have analyzed the magnetic non-potentiality of ARs and the
constraining effect of background fields for 106 flares $\geq$M1.0-class from
12 CME-active and 9 CME-quiet ARs. We proposed a new parameter
$\alpha$FPIL/$\Phi$AR that is well able to distinguish ARs with the capability
of producing eruptive flares. About 93% of eruptive events have the
distributions of $\alpha$FPIL/$\Phi$AR$\geq$2.2$\times$$10^{-24}$ Mm-1 Mx-1,
and $\sim$ 83% of confined flares have
$\alpha$FPIL/$\Phi$AR$<$2.2$\times$$10^{-24}$ Mm-1 Mx-1. Moreover, we select
the areas of high photospheric magnetic free energy density before the flare
onset as a proxy for AR core region to calculate mean characteristic twist
parameter $\alpha$HFED and mean shear angle $\Psi$HFED. The statistical
results showed that the relative parameters $\alpha$HFED/$\Phi$AR and
$\Psi$HFED/$\Phi$AR can also provide a good ability for distinguishing the
eruptive and confined flares. Overall, about 81% (86 out of 106) of flare
events can be classified into the two populations of confined and eruptive
flares using the parameter $\alpha$HFED/$\Phi$AR (through determining a
critical value). The performance of $\Psi$HFED/$\Phi$AR is not as good as that
of $\alpha$HFED/$\Phi$AR, and about 74% (78 out of 106) of flare events can be
distinguished.

Parameter $\alpha$ is the average characteristic twist of the magnetic field
lines around the PILs of an AR, and provide measures of the nonpotentiality of
AR core region (Leka & Skumanich 1999; Benson et al. 2021). Previous studies
have shown that the parameter of magnetic twist plays an important role in
discriminating between confined and eruptive events and can be used to predict
whether an X- or M-class flaring AR would produce a CME (Bobra & Ilonidis
2016; Duan et al. 2019). Magnetic twist is thought to be related to the
driving force of electric current-driven instabilities, such as the helical
kink instability (Hood & Priest 1979). An AR containing a highly twisted
magnetic field tends to produce an eruption when the twist exceeds a certain
threshold. On the other hand, according to the study of Li et al. (2020),
$\Phi$AR has a high positive correlation with the critical decay index height
(related to the torus instability of a magnetic flux rope; Kliem & Török
2006), implying that $\Phi$AR describes the strength of the background field
confinement. Our statistical study reveals that the relative parameter
$\alpha$/$\Phi$AR has a better performance in distinguishing between the two
types of flares than only $\Phi$AR or $\alpha$ does. We suggest that the
relative parameter $\alpha$/$\Phi$AR indicates the balance between the upward
force that drives the eruptions and the downward force that suppresses the
eruptions.

Recently, several studies have shown the general trend that confined events
have a smaller reconnection flux fraction ($\Phi$ribbon/$\Phi$AR; defined by
the ratio between the flux swept by the flare ribbons to the total AR flux)
compared with eruptive flares (Toriumi et al. 2017; Li et al. 2020; Kazachenko
et al. 2021). However, there is a considerable overlap between the flux
fraction distributions of confined and eruptive flares. Comparably, our
parameters $\alpha$/$\Phi$AR and $\alpha$HFED/$\Phi$AR show more significant
differences in distributions between the two populations of flares. We
hypothesize that $\Phi$ribbon involves part of overlying background confining
fields of ARs participating the flare reconnection at the late stage of the
flare. In comparison, $\alpha$FPIL and $\alpha$HFED only correspond to the
pre-flare high-twist fields in an AR core region, which are strongly related
to the initial driving force of solar flares. Therefore, $\alpha$FPIL/$\Phi$AR
and $\alpha$HFED/$\Phi$AR can provide a better ability for distinguishing the
eruptive and confined events than $\Phi$ribbon/$\Phi$AR.

The idea of the new parameter $\alpha$/$\Phi$AR can be generalized to
“relative non-potentiality”, which refers to the ratio of magnetic flux (or
other physical quantities) in a flux rope to that in the surrounding magnetic
structures (Lin et al. 2021). A larger relative non-potentiality indicates a
higher probability for a flux rope to erupt (Toriumi et al. 2017; Thalmann et
al. 2019; Gupta et al. 2021). Recently, Lin et al. (2020) proposed a new
relative non-potentiality parameter of the magnetic flux in the highly twisted
region relative to its ambient background fields, which demonstrated a
moderate ability in discriminating between confined and ejective events. It
seems that our parameters $\alpha$/$\Phi$AR shows a better classification
performance than the parameter proposed in Lin et al. (2020). For instance,
two events of X1.6 and X3.1-class flares in AR 12192 failed to be correctly
classified in Lin et al. (2020, 2021), however, they can be correctly
classified in our study. However, due to the difference of selected samples, a
direct comparison can not be carried out.

In recent years, the importance of the confinement of flux ropes by overlying
loops (and thus the importance of structural relativity) has been actively
discussed in the theoretical and modeling studies. Leake et al. (2013)
simulated a magnetic flux emergence into pre-existing dipole coronal field and
found that it becomes a stable flux rope if the dipole coronal field is
orientated to minimize magnetic reconnection. Using the same simulation
results, Pariat et al. (2017) found that the ratio of the magnetic helicity of
the current-carrying magnetic field to the total relative helicity diagnoses
very clearly the eruptive potential of their parametric simulations. Toriumi &
Takasao (2017) conducted flux emergence models and showed that the confinement
of the flux rope or the access to the outer space depends on the large-scale
AR structures, which determines the CME eruption. The CME eruption is
predicted by the SHARP parameters that characterize the “relativity” between
the flaring zone and overall AR area. Amari et al. (2018) suggested that the
role of the magnetic cage is one deciding factor for the success or failure of
CMEs.

Our findings imply that the relative measure of magnetic nonpotentiality
within the AR core over the restriction of the background field largely
controls whether a flare is eruptive or confined. However, it needs to be
noted that there is still a small overlap (less than 20%) between the two
populations of confined and eruptive flares by using parameter
$\alpha$/$\Phi$AR. We suggest that other unknown mechanisms or intrinsic
stochasticity may also play a role in governing CME production.

We are grateful to Dr. Jiangtao Su for useful discussions. This work is
supported by the B-type Strategic Priority Program of the Chinese Academy of
Sciences (XDB41000000), the National Key R&D Program of China
(2019YFA0405000), the National Natural Science Foundations of China (11773039,
11903050, 12073001, 11790304, 11873059 and 11790300), Key Programs of the
Chinese Academy of Sciences (QYZDJ-SSW-SLH050), the Youth Innovation Promotion
Association of CAS (2014043), Yunnan Academician Workstation of Wang Jingxiu
(No. 202005AF150025) and NAOC Nebula Talents Program. X. D. Sun acknowledges
the support by NSF award 1848250. A. Q. Chen is supported by the Strategic
Priority Program on Space Science, Chinese Academy of Sciences, Grant No.
XDA15350203. _SDO_ is a mission of NASA’s Living With a Star Program.

## References

  * (1) Amari, T., Canou, A., Aly, J.-J., et al. 2018, Nature, 554, 211.
  * (2) Andrews, M. D. 2003, Sol. Phys., 218, 261
  * (3) Avallone, E. A. & Sun, X. 2020, ApJ, 893, 123
  * (4) Baumgartner, C., Thalmann, J. K., & Veronig, A. M. 2018, ApJ, 853, 105
  * (5) Benson, B., Pan, W. D., Prasad, A., et al. 2021, Sol. Phys., 296, 163
  * (6) Bobra, M. G. & Ilonidis, S. 2016, ApJ, 821, 127.
  * (7) Bobra, M. G., Sun, X., Hoeksema, J. T., et al. 2014, Sol. Phys., 289, 3549
  * (8) Chen, A. Q. & Wang, J. X. 2012, A&A, 543, A49
  * (9) Cui, Y., Wang, H., Xu, Y., et al. 2018, Journal of Geophysical Research (Space Physics), 123, 1704
  * (10) Duan, A., Jiang, C., He, W., et al. 2019, ApJ, 884, 73
  * (11) Falconer, D. A., Moore, R. L., & Gary, G. A. 2002, ApJ, 569, 1016
  * (12) Falconer, D. A., Moore, R. L., & Gary, G. A. 2006, ApJ, 644, 1258
  * (13) Gupta, M., Thalmann, J. K., & Veronig, A. M. 2021, A&A, 653, A69
  * (14) Hood, A. W. & Priest, E. R. 1979, Sol. Phys., 64, 303
  * (15) Jing, J., Liu, C., Lee, J., et al. 2018, ApJ, 864, 138
  * (16) Kazachenko, M. D., Lynch, B. J., Savcheva, A., et al. 2021, arXiv:2111.06048
  * (17) Kazachenko, M. D., Lynch, B. J., Welsch, B. T., et al. 2017, ApJ, 845, 49
  * (18) Kliem, B., & Török, T. 2006, Phys. Rev. Lett., 96, 255002
  * (19) Leake, J. E., Linton, M. G., & Török, T. 2013, ApJ, 778, 99
  * (20) Leka, K. D. & Skumanich, A. 1999, Sol. Phys., 188, 3
  * (21) Lemen, J. R., Title, A. M., Akin, D. J., et al. 2012, Sol. Phys., 275, 17
  * (22) Li, T., Hou, Y., Yang, S., et al. 2020, ApJ, 900, 128
  * (23) Li, T., Chen, A., Hou, Y., et al. 2021, ApJ, 917, L29
  * (24) Lin, P. H., Kusano, K., & Leka, K. D. 2021, ApJ, 913, 124
  * (25) Lin, P. H., Kusano, K., Shiota, D., et al. 2020, ApJ, 894, 20
  * (26) Liu, L., Wang, Y., Wang, J., et al. 2016, ApJ, 826, 119
  * (27) Liu, Y., Hoeksema, J. T., Scherrer, P. H., et al. 2012, Sol. Phys., 279, 295
  * (28) Nindos, A., & Andrews, M. D. 2004, ApJ, 616, L175
  * (29) Pariat, E., Leake, J. E., Valori, G., et al. 2017, A&A, 601, A125
  * (30) Pesnell, W. D., Thompson, B. J., & Chamberlin, P. C. 2012, Sol. Phys., 275, 3
  * (31) Scherrer, P. H., Schou, J., Bush, R. I., et al. 2012, Sol. Phys., 275, 207
  * (32) Schrijver, C. J. 2007, ApJ, 655, L117
  * (33) Su, J. T., Jing, J., Wang, S., et al. 2014, ApJ, 788, 150
  * (34) Sun, X. 2013, arXiv e-prints, arXiv:1309.2392
  * (35) Sun, X., Bobra, M. G., Hoeksema, J. T., et al. 2015, ApJ, 804, L28
  * (36) Thalmann, J. K., Moraitis, K., Linan, L., et al. 2019, ApJ, 887, 64
  * (37) Toriumi, S., Schrijver, C. J., Harra, L. K., et al. 2017, ApJ, 834, 56
  * (38) Toriumi, S. & Takasao, S. 2017, ApJ, 850, 39
  * (39) Török, T. & Kliem, B. 2005, ApJ, 630, L97
  * (40) Vasantharaju, N., Vemareddy, P., Ravindra, B., et al. 2018, ApJ, 860, 58
  * (41) Wang, D., Liu, R., Wang, Y., et al. 2017, ApJ, 843, L9
  * (42) Wang, H., Ewell, M. W., Zirin, H., et al. 1994, ApJ, 424, 436
  * (43) Wang, Y., & Zhang, J. 2007, ApJ, 665, 1428
  * (44) Yang, X., Zhang, H., Gao, Y., et al. 2012, Sol. Phys., 280, 165
  * (45) Yashiro, S., Akiyama, S., Gopalswamy, N., et al. 2006, ApJ, 650, L143

Table 1: Parameters Used to Distinguish the Eruptive and Confined Flares11Adapted from Chen & Wang (2012) and Bobra et al. (2014). Parameters | Description | Unit | Formula
---|---|---|---
$\Phi$AR | Total unsigned flux | Mx |  $\Phi$AR=$\Sigma$$|$$B_{z}$$|$dA
Jz | Mean vertical electric current density | mA m-2 | Jz=$\frac{1}{N\mu}$$\Sigma$($\nabla$$\times$B)z 22$\mu$ is the magnetic permeability in vacuum (4$\pi$$\times$$10^{-3}$ G m $A^{-1}$).
$\alpha$ | Mean characteristic twist parameter | Mm-1 |  $\alpha$=$\frac{\mu\Sigma J_{z}B_{z}}{\Sigma B_{z}^{2}}$
$\rho_{free}$ | Magnetic free energy density | erg cm-3 |  $\rho_{free}$=$\frac{1}{8\pi}$$|$Bobs-Bpot$|$2
$\Psi$ | Mean shear angle | degree |  $\Psi$=$\arccos$$\frac{\textbf{B}_{obs}\cdot\textbf{B}_{pot}}{|B_{obs}B_{pot}|}$
Figure 1: Four examples of two eruptive and two confined flares showing
SDO/AIA 1600 Å images (left), SDO/HMI photospheric magnetograms $B_{z}$
(middle) and derived vertical electric current density $J_{z}$ maps (right).
From top to bottom: eruptive X2.1-class flare in AR 11283, eruptive M6.5-class
flare in AR 12371, confined X3.1-class flare in AR 12192 and confined
M6.1-class flare in AR 12222. AIA 1600 Å images were remapped with CEA
projection. Orange and black contours outline the FPIL mask regions within
which the mean characteristic twist parameter ($\alpha$FPIL) in Figure 2 was
calculated.

Figure 2: Scatter plots of mean characteristic twist parameter $\alpha$FPIL
versus total unsigned magnetic flux $\Phi$AR and flare peak X-ray flux FSXR
versus $\alpha$FPIL/$\Phi$AR. Blue (red) circles are eruptive (confined)
flares. Two vertical green lines in panel (a) correspond to $\Phi$AR of
3.5$\times$$10^{22}$ Mx and 1.0$\times$$10^{23}$ Mx, respectively. The
horizontal black line in panel (a) refers to $\alpha$FPIL of 0.07 Mm-1. The
vertical black line in panel (b) denotes $\alpha$FPIL/$\Phi$AR of
2.2$\times$$10^{-24}$ Mm-1 Mx-1.

Figure 3: Maps of photospheric free magnetic energy density
$\rho_{\mathrm{free}}$ and magnetic shear angle $\Psi$ of four examples (from
top to bottom: eruptive X2.1-class flare in AR 11283, eruptive M6.5-class
flare in AR 12371, confined X3.1-class flare in AR 12192 and confined
M6.1-class flare in AR 12222). The white and green contours are the magnetic
fields $B_{z}$ at $\pm$800 G levels. Mean characteristic twist parameter
$\alpha$HFED and mean shear angle $\Psi$HFED in Figure 4 are calculated within
the areas of $\rho_{\mathrm{free}}$$>$4.0$\times$104 erg cm-3.

Figure 4:  Similar to Figure 2, but for mean characteristic twist parameter
$\alpha$HFED and mean shear angle $\Psi$HFED within HFED region. Blue (red)
circles are eruptive (confined) flares. The horizontal dotted line in panel
(a) corresponds to $\alpha$HFED of 0.1 Mm-1. The vertical dash-dotted line in
panel (b) refers to $\alpha$HFED/$\Phi$AR of 3.2$\times$$10^{-24}$ Mm-1 Mx-1.
The horizontal dotted line in panel (c) denotes $\Psi$HFED of 60∘. The
vertical dash-dotted line in panel (d) corresponds to $\Psi$HFED/$\Phi$AR of
1.0$\times$$10^{-21}$ degree Mx-1.


together with the Dirichlet boundary conditions $J_{m}(q_{mn}R)=0$. Comparing
Eqs. (180e), (180a) and (180b), it is easy to see that
${\overline{\mathcal{E}}}=\mathcal{P}/3$ with
$\mathcal{P}=\frac{2}{3}\mathcal{P}_{R}+\frac{1}{3}\mathcal{P}_{z}$ being the
isotropic pressure. Moreover, the relations in Eq. (180) are compatible with
the Euler relation (76).

The relations in Eq. (180) are valid for both real and imaginary rotation. In
the latter case, ${\overline{\mathcal{M}}}$ becomes

${\overline{\mathcal{M}}}=-\frac{8}{L^{2}\beta^{2}}\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}\int_{0}^{\infty}\frac{dk\,m\,e^{\beta\omega}\sin(\beta\Omega_{I}m)}{e^{2\beta\omega}-2e^{\beta\omega}\cos(\beta\Omega_{I}m)+1}.$
(182)

Thus, ${\overline{\mathcal{M}}}$ vanishes in the imaginary rotation case when
$\nu=1/2$, as was the case also in the unbounded system [see Eq. (166c)].

The results for $\overline{\Phi^{2}}$ and ${\overline{\mathcal{E}}}$ are shown
in the top and bottom panels of Fig. 13. As before, we set $\nu=p^{\prime}/10$
with $0\leq p^{\prime}\leq 5$, leading to irreducible fractions $p/q$ with
$q\in\\{1,2,5,10\\}$. The horizontal axis shows $L=2\pi R/\beta$, where $R$
represents the radius of the bounding cylinder. The dotted black lines
represent same quantities computed for the unbounded system using Eqs. (160a)
and (160b) for the same value of $L$. For $L\lesssim 10$, the boundary effects
lead to strong quenching of both $\phi^{2}$ and $T^{tt}$, such that both
$\overline{\Phi^{2}}$ and ${\overline{\mathcal{E}}}$ tend to $0$ as
$L\rightarrow 0$ in the bounded case. This is contrary to the unbounded case,
where the $L\rightarrow 0$ limit is finite in both cases. As already seen in
Fig. 10, with increasing $L$, the boundary effects become localized around a
small vicinity of the boundary and the bounded and unbounded results approach
each other. While for $\overline{\Phi^{2}}$, visible discrepancies remain even
for $L\gtrsim 10^{3}$, in the case of ${\overline{\mathcal{E}}}$, the results
obtained in the bounded case start following the ones corresponding to the
unbounded case already when $L\gtrsim 10$. As in the previous sections, the
fractal structure reveals itself at large values of $L$.

###  VI.3 Slow rotation: moment of inertia and shape

Figure 15: Ratios ${\overline{\mathcal{F}}}(0)/F_{0}$, $K_{2}/2!$ and
$K_{4}/4!$ computed at various values of the normalized transverse size of the
system, $L=2\pi R/\beta$ for the bounded system discussed in Sec. VI.

Finally, we discuss the expansion in Eq. (82) of the free energy density
${\overline{\mathcal{F}}}$ in the case of slow rotation. In particular, we
focus on the free energy in the absence of rotation,
${\overline{\mathcal{F}}}(0)$, as well as the first two coefficients, the
dimensionless moment of inertia, $K_{2}$, and the dimensionless shape
coefficient, $K_{4}$, which we evaluate using:

$\displaystyle{\overline{\mathcal{F}}}(0)$
$\displaystyle=-\frac{4}{L^{2}\beta^{2}}\sum_{m=-\infty}^{\infty}\sum_{n=1}^{\infty}\int_{0}^{\infty}\frac{dk\,k^{2}/\omega_{mn}}{e^{\beta\omega_{mn}}-1},$
(183a) $\displaystyle\overline{K}_{2}$
$\displaystyle=-\frac{16\pi^{2}}{\beta^{3}L^{4}}\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}\int_{0}^{\infty}\frac{m^{2}\,dk}{\sinh^{2}(\beta\omega_{m,n}/2)},$
(183b) $\displaystyle\overline{K}_{4}$
$\displaystyle=-\frac{16\pi^{4}}{L^{6}\beta^{3}}\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}\int_{0}^{\infty}\frac{dk\,m^{4}[2+\cosh(\beta\omega_{mn})]}{\sinh^{4}(\beta\omega_{mn}/2)},$
(183c)

where the notation $\overline{K}_{2n}=K_{2n}{\overline{\mathcal{F}}}(0)$ was
employed for brevity. The coefficients ${\overline{\mathcal{F}}}(0)$, $K_{2}$
and $K_{4}$ are studied as functions of the transverse size of the system $L$
with respect to their unbounded counterparts.

The ratios ${\overline{\mathcal{F}}}(0)/F_{0}$, $K_{2}/2!$ and $K_{4}/4!$ are
represented in Fig. 15, where the denominators of these expressions are the
classical expectations given in Eq. (83). As we already noted in Subsect.
IV.3, a strong quenching due to the boundary can be seen at small values of
$L$, which is however less pronounced for the shape coefficient $K_{4}$
compared to the dimensionless moment of inertia $K_{2}$ and the free energy
${\overline{\mathcal{F}}}_{0}$ itself. At $L=100$, these coefficients already
reach their expected asymptotic values $K_{2n}=(2n)!$, Eq. (83).

In a cylinder of a finite radius $R$, a large-size $L\to\infty$ limit
corresponds also to the high-temperature limit, $T\to\infty$ since $L=2\pi
R/\beta\equiv 2\pi RT$. Figure 15 shows that the dimensionless moment of
inertia $K_{2}$ approaches the asymptotic value $K_{2}=2$ from below,
indicating that the moment of inertia should decrease as temperature
decreases. This effect is related to the presence of an effective energy gap
between the states with zero, $m=0$, and non-zero, $m\neq 0$, orbital momenta
due to the finite size of the system. Therefore, at lower temperatures, the
system mostly resides in the $m=0$ state and the rotational modes, which
contribute to the moment of inertia (183b), are not excited. Since the latter
modes do not participate in rotation at low $T$, the moment of inertia of the
system decreases as the system gets colder. This effect should evidently also
occur for $K_{4}$ and higher coefficients. Interestingly, the same qualitative
behaviour for the moment of inertia $K_{2}$ is also observed in the first-
principle simulations of gluon plasma in the high-temperature phase of Yang-
Mills theory: as the temperature increases, the moment of inertia approaches
the high-temperature value $K_{2}=2$, Eq. (83), from below Braguta  _et al._
(2023).

##  VII Conclusions

In the present work, we studied the thermodynamic properties of massless
scalar fields subjected to rigid rotation and an inter-relation of the real
rotation with its imaginary analogue. The latter concept – a rigid rotation
with an imaginary frequency Chernodub (2021a); Chen  _et al._ (2022) – has a
practical interest since rotating systems cannot be implemented in Euclidean
path-integral formalism, suitable, for example, for numerical first-principle
calculations on the lattice Yamamoto and Hirono (2013); Braguta  _et al._
(2020, 2021, 2023); Chernodub  _et al._ (2022). In this sense, rotation shares
the deficiency suffered by finite-density systems, namely the sign problem
Yamamoto and Hirono (2013), and needs to be implemented in Euclidean spacetime
via its imaginary version supplemented with the subsequent analytical
continuation to real angular frequencies Braguta  _et al._ (2020, 2021, 2023).

Using the $1+1$-dimensional toy model of a scalar field under rigid rotation
on a ring, we explicitly demonstrated that the analytical no-go theorem
Chernodub (2022b), which describes the impossibility of continuation of the
imaginary-angular-frequency thermodynamics to the real angular frequencies is
related to the development of the fractality of thermodynamics for the former.
The result applies in the thermodynamic limit. Within this model,
thermodynamic functions such as pressure and energy density can be expressed
analytically via the Dedekind $\eta$ function (31). The latter function tends
to the fractal, non-analytical Thomae function (17) as its argument approaches
the real axis Nechaev and Polovnikov (2017), which corresponds to the
infinite-volume (thermodynamic) limit.

In the case of the $3+1$-dimensional Minkowski space, we first considered a
classical description of scalar particles under rotation using relativistic
kinetic theory. In the absence of boundaries, rigid rotation with a real
rotation parameter leads to a violation of causality and subsequent divergence
of all observables on the light cylinder. Imaginary rotation can be described
by a real distribution function only as an average of clockwise and
counterclockwise rotations, leading a seemingly non-equilibrium state. As
expected, observables such as the energy-momentum tensor $T^{\mu\nu}$ decrease
as the distance to the rotation axis is increased, at a faster rate for faster
rotation.

Under the quantum-field theoretical treatment, rigid rotation with a real
rotation parameter leads to the divergence of $T^{\mu\nu}$ at each point of
the space-time, also inside the light cylinder. Under imaginary rotation,
$T^{\mu\nu}$ evaluated on the rotation axis can be expressed via Bernoulli
polynomials (137). Away from the rotation axis, we were able to demonstrate
the fractalization of both the field fluctuations $\phi^{2}$ and of
$T^{\mu\nu}$ in the case of imaginary rotation, in complete analogy to the
$1+1$-dimensional toy model discussed above. In all cases, our results exhibit
a periodicity with respect to the imaginary rotation parameter which is not
present in the classical, kinetic analysis. For this reason, we found
agreement with kinetic theory only in the limit of slow rotation and only in
the vicinity of the rotation axis (before fractalization sets in).

For comparison with the case of real rotation, we took results obtained using
a perturbative calculation for slow rotation (with respect to a stationary
background) and found on the rotation axis an analytical result (145) which
can be related to the one obtained for imaginary rotation in an essentially
non-analytical way (146). We conclude that even on the axis of rotation –
which appears to be static in the rotating system – the non-analyticity is
strong and unavoidable.

We demonstrated the same analytic–non-analytical transition using numerical
calculations for the thermal scalar field system enclosed in a cylinder,
undergoing rotation with imaginary angular frequencies. As the radius of the
cylinder grows, the pressure becomes a non-analytical function of temperature
expressed, again, via the Thomae function (demonstrated in Figs. 8–14). In
this limit, the boundary effects become less important and the results
obtained in the unbounded case provide a good approximation for our
observables inside the cylindrical boundary. For values of the rotation
parameter respecting the causality constraints (i.e., when the light cylinder
is outside the boundary), the fractalization features do not appear in the
case of imaginary rotation. For sufficiently slow rotation and high
temperature, our numerically-obtained results are compatible with the flip
$\Omega_{I}^{2}\rightarrow-\Omega^{2}$ from imaginary to real rotation,
signaling the restoration of analytical continuation in this limit.

The exotic fractalization properties discussed above are related to a ninionic
deformation (12) of statistical distributions at imaginary angular frequencies
Chernodub (2022b). For this reason, we conclude that the results obtained in
the infinite-volume system subjected to imaginary rotation cannot be
analytically related to the properties of the physically rotating system (with
a real angular frequency). However, we explicitly demonstrated both for the
analytically treatable case on the ring and numerically accessible case of
rotating cylinder that the imaginary rotation in a spatially bounded system in
Euclidean space can be continued to the real-frequency domain in Minkowski
spacetime provided that the causality is respected for the latter spacetime.

We have also shown that for real-frequency rotation, the dimensionless moment
of inertia $K_{2}$, normalized per one degree of freedom, is equal to two,
$K_{2}=2$, in the thermodynamic limit of large radius $R$ of the cylinder. The
quantity $K_{2}$ determines the correction to the free energy (82), or,
equivalently, to the pressure of a non-interacting gas of scalar bosons,

$\displaystyle
P(\Omega)=P(0)\biggl{(}1+\frac{1}{2}K_{2}R^{2}\Omega^{2}+\dots\biggr{)}\,,$
(184)

due to small nonzero angular frequency, $\Omega\to 0$. This result matches
well the first-principle result of Ref. Braguta  _et al._ (2023) on the
behavior of gluons in high-temperature limit of Yang-Mills theory.

Below we summarize our main findings:

  1. 1.

We demonstrated the no-go theorem Chernodub (2022b) regarding the
impossibility of continuation of the imaginary-angular-frequency
thermodynamics to the real angular frequencies using an analytical
$1+1$-dimensional toy model, revealing the development of the fractality of
thermodynamics for the former.

  2. 2.

Since fractalization does not show up in the classical, kinetic theory
treatment of imaginary rotation, we conclude that this is a purely quantum
effect.

  3. 3.

We found similar fractalization on the unbounded Minkowski space for imaginary
rotation, as well as in the thermodynamic (infinite volume) limit of the
bounded system.

  4. 4.

For the case of a causal boundary that excludes the light cylinder, we were
able to restore the analytical continuation from imaginary to real rotation in
the limit of slow rotation and large temperatures.

  5. 5.

We attributed the exotic fractalization properties to a ninionic deformation
(12) of statistical distributions at imaginary angular frequencies Chernodub
(2022b).

The results obtained in this paper shed light on the implications of the
effects of imaginary rotation obtained in the context of first-principle
lattice simulations and constitute the basis for the analysis of more
complicated systems, e.g., free fermions or the chiral phase transition in the
effective QCD models such as the Nambu–Jona-Lasinio model or nonlinear sigma
models.

###### Acknowledgements.

V.E.A. gratefully acknowledges the support through a grant of the Ministry of
Research, Innovation and Digitization, CNCS - UEFISCDI, project number PN-
III-P1-1.1-TE-2021-1707, within PNCDI III.

## References

  * Beams (1968) J. W. Beams, Physical Review Letters 21, 1093 (1968).
  * Barnett (1915) S. J. Barnett, Physical Review 6, 239 (1915).
  * Einstein and De Haas (1915) A. Einstein and W. De Haas, in  _Proc. KNAW_ , Vol. 181 (1915) p. 696.
  * Arabgol and Sleator (2019) M. Arabgol and T. Sleator, Physical Review Letters 122 (2019), 10.1103/physrevlett.122.177202.
  * Adamczyk  _et al._ (2017) L. Adamczyk _et al._ (STAR), Nature 548, 62 (2017), arXiv:1701.06657 [nucl-ex] .
  * Deng and Huang (2016) W.-T. Deng and X.-G. Huang, Phys. Rev. C 93, 064907 (2016), arXiv:1603.06117 [nucl-th] .
  * Jiang  _et al._ (2016) Y. Jiang, Z.-W. Lin, and J. Liao, Phys. Rev. C 94, 044910 (2016), [Erratum: Phys.Rev.C 95, 049904 (2017)], arXiv:1602.06580 [hep-ph] .
  * Becattini and Lisa (2020) F. Becattini and M. A. Lisa, Ann. Rev. Nucl. Part. Sci. 70, 395 (2020), arXiv:2003.03640 [nucl-ex] .
  * Huang  _et al._ (2020) X.-G. Huang, J. Liao, Q. Wang, and X.-L. Xia,  (2020), 10.1007/978-3-030-71427-7_9, arXiv:2010.08937 [nucl-th] .
  * Chen  _et al._ (2016) H.-L. Chen, K. Fukushima, X.-G. Huang, and K. Mameda, Phys. Rev. D 93, 104052 (2016), arXiv:1512.08974 [hep-ph] .
  * Jiang and Liao (2016) Y. Jiang and J. Liao, Phys. Rev. Lett. 117, 192302 (2016), arXiv:1606.03808 [hep-ph] .
  * Chernodub and Gongyo (2017a) M. N. Chernodub and S. Gongyo, JHEP 01, 136 (2017a), arXiv:1611.02598 [hep-th] .
  * Chernodub and Gongyo (2017b) M. N. Chernodub and S. Gongyo, Phys. Rev. D 95, 096006 (2017b), arXiv:1702.08266 [hep-th] .
  * Wang  _et al._ (2019) X. Wang, M. Wei, Z. Li, and M. Huang, Phys. Rev. D 99, 016018 (2019), arXiv:1808.01931 [hep-ph] .
  * Zhang  _et al._ (2020) Z. Zhang, C. Shi, X.-T. He, X. Luo, and H.-S. Zong, Phys. Rev. D 102, 114023 (2020), arXiv:2012.01017 [hep-ph] .
  * Sadooghi  _et al._ (2021) N. Sadooghi, S. M. A. Tabatabaee Mehr, and F. Taghinavaz, Phys. Rev. D 104, 116022 (2021), arXiv:2108.12760 [hep-ph] .
  * Braguta  _et al._ (2020) V. V. Braguta, A. Y. Kotov, D. D. Kuznedelev, and A. A. Roenko, Pisma Zh. Eksp. Teor. Fiz. 112, 9 (2020).
  * Chen  _et al._ (2021) X. Chen, L. Zhang, D. Li, D. Hou, and M. Huang, JHEP 07, 132 (2021), arXiv:2010.14478 [hep-ph] .
  * Chernodub (2021a) M. N. Chernodub, Phys. Rev. D 103, 054027 (2021a), arXiv:2012.04924 [hep-ph] .
  * Fujimoto  _et al._ (2021) Y. Fujimoto, K. Fukushima, and Y. Hidaka, Phys. Lett. B 816, 136184 (2021), arXiv:2101.09173 [hep-ph] .
  * Braguta  _et al._ (2021) V. V. Braguta, A. Y. Kotov, D. D. Kuznedelev, and A. A. Roenko, Phys. Rev. D 103, 094515 (2021), arXiv:2102.05084 [hep-lat] .
  * Golubtsova  _et al._ (2022) A. A. Golubtsova, E. Gourgoulhon, and M. K. Usova, Nucl. Phys. B 979, 115786 (2022), arXiv:2107.11672 [hep-th] .
  * Chen  _et al._ (2022) S. Chen, K. Fukushima, and Y. Shimada, Phys. Rev. Lett. 129, 242002 (2022), arXiv:2207.12665 [hep-ph] .
  * Golubtsova and Tsegelnik (2022) A. A. Golubtsova and N. S. Tsegelnik, (2022), arXiv:2211.11722 [hep-th] .
  * Zhao  _et al._ (2022) Y.-Q. Zhao, S. He, D. Hou, L. Li, and Z. Li, (2022), arXiv:2212.14662 [hep-ph] .
  * Chernodub  _et al._ (2022) M. N. Chernodub, V. A. Goy, and A. V. Molochkov, (2022), arXiv:2209.15534 [hep-lat] .
  * Braguta  _et al._ (2023) V. V. Braguta, M. N. Chernodub, A. A. Roenko, and D. A. Sychev, (2023), arXiv:2303.03147 [hep-lat] .
  * Ambruş and Winstanley (2014) V. E. Ambruş and E. Winstanley, Phys. Lett. B 734, 296 (2014), arXiv:1401.6388 [hep-th] .
  * Ambrus and Winstanley (2016) V. E. Ambrus and E. Winstanley, Phys. Rev. D 93, 104014 (2016), arXiv:1512.05239 [hep-th] .
  * Yamamoto and Hirono (2013) A. Yamamoto and Y. Hirono, Phys. Rev. Lett. 111, 081601 (2013), arXiv:1303.6292 [hep-lat] .
  * Chernodub (2022a) M. N. Chernodub, (2022a), arXiv:2208.04808 [hep-th] .
  * Chernodub (2022b) M. N. Chernodub, (2022b), arXiv:2210.05651 [quant-ph] .
  * Kapusta and Gale (2006) J. I. Kapusta and C. Gale,  _Finite-temperature field theory: Principles and applications_ (Cambridge University Press, 2006).
  * Rothe (2012) H. J. Rothe,  _Lattice gauge theories: an introduction_ (World Scientific Publishing Company, 2012).
  * de Forcrand (2009) P. de Forcrand, PoS LAT2009, 010 (2009), arXiv:1005.0539 [hep-lat] .
  * Chernodub (2021b) M. N. Chernodub, Symmetry 13, 1569 (2021b), arXiv:2104.05032 [cond-mat.mes-hall] .
  * Devaney (1999) R. L. Devaney, The American Mathematical Monthly 106, 289 (1999).
  * Nechaev and Polovnikov (2017) S. Nechaev and K. Polovnikov, arXiv preprint arXiv:1702.06757 (2017).
  * Milton (2001) K. A. Milton,  _The Casimir Effect_ (WORLD SCIENTIFIC, 2001).
  * Tolman and Ehrenfest (1930) R. Tolman and P. Ehrenfest, Phys. Rev. 36, 1791 (1930).
  * Tolman (1930) R. C. Tolman, Phys. Rev. 35, 904 (1930).
  * Duffy and Ottewill (2003) G. Duffy and A. C. Ottewill, Phys. Rev. D 67, 044002 (2003), arXiv:hep-th/0211096 .
  * Christensen (1976) S. M. Christensen, Phys. Rev. D 14, 2490 (1976).
  * Davies  _et al._ (1996) P. C. W. Davies, T. Dray, and C. A. Manogue, Phys. Rev. D 53, 4382 (1996), arXiv:gr-qc/9601034 .
  * Ambrus (2014) V. E. Ambrus,  _Dirac fermions on rotating space-times_ , Ph.D. thesis, Sheffield U. (2014).
  * Becattini and Grossi (2015) F. Becattini and E. Grossi, Phys. Rev. D 92, 045037 (2015), arXiv:1505.07760 [gr-qc] .
  * Ambrus (2017) V. E. Ambrus, Phys. Lett. B 771, 151 (2017), arXiv:1704.02933 [hep-th] .
  * Becattini  _et al._ (2021) F. Becattini, M. Buzzegoli, and A. Palermo, JHEP 02, 101 (2021), arXiv:2007.08249 [hep-th] .


EPJ Web of Conferences CONF12  11institutetext: Instituto de Física de
Cantabria / Instituto de Física Teórica (IFCA/IFT-CSIC) 22institutetext:
Instituto de Física Teórica UAM/CSIC, Universidad Autónoma de Madrid,
C/ Nicolás Cabrera 13-15, Cantoblanco, Madrid 28049, Spain 33institutetext:
Universidad Autónoma de Madrid, Cantoblanco, Madrid 28049, Spain
44institutetext: Theoretical Physics Department, CERN, 1211 Geneva 23,
Switzerland 55institutetext: INFN, Sezione di Tor Vergata, c/o Dipartimento di
Fisica, UniversitÃ di Roma Tor Vergata,
Via della Ricerca Scientifica 1, 00133 Rome, Italy

# Controlling quark mass determinations non-perturbatively in three-flavour
QCD

Isabel Campos 11 Patrick Fritzsch 2244<EMAIL_ADDRESS>Carlos Pena 2233
David Preti 22 Alberto Ramos 44 Anastassios Vladikas 55

###### Abstract

The determination of quark masses from lattice QCD simulations requires a non-
perturbative renormalization procedure and subsequent scale evolution to high
energies, where a conversion to the commonly used $\kern 1.0pt\overline{\kern
0.0pt\rm MS\kern-1.0pt}\kern 1.0pt$ scheme can be safely established. We
present our results for the non-perturbative running of renormalized quark
masses in $N_{\rm f}=3$ QCD between the electroweak and a hadronic energy
scale, where lattice simulations are at our disposal. Recent theoretical
advances in combination with well-established techniques allows to follow the
scale evolution to very high statistical accuracy, and full control of
systematic effects.

##  1 Quantum Chromodynamics at all scales

For numerical simulations one starts from one of many valid discretisations of
the (bare) Euclidean action density of Quantum Chromodynamics (QCD),

$\displaystyle\mathcal{L}_{\tiny QCD}$
$\displaystyle=\frac{1}{2g_{0}^{2}}\,\hbox{tr}\,\\{F_{\mu\nu}F_{\mu\nu}\\}+\sum_{i=1}^{N_{\rm
f}}\overline{\psi}_{i}\big{(}\gamma_{\mu}D_{\mu}+m_{0,i}\big{)}\psi_{i}\;,$
$\displaystyle D_{\mu}$ $\displaystyle=\partial_{\mu}+A_{\mu}\;,$ (1)

with a fixed number of dynamical quark flavours $N_{\rm f}$. With $N_{\rm
f}+1$ input parameters the action is supposed to describe the strong
interaction at all energy scales. Lattice QCD is a gauge-invariant
regularisation and thus does not require a gauge-fixing and Faddeev–Popov
term. It is the only regularisation allowing to determine colorless bound
states $(m_{\tiny\pi^{0}},m_{\tiny\rm K^{0}},m_{\rm p},\ldots)$ between quarks
from first principle calculations, which due to computational and theoretical
advances can be systematically improved over time. However, incorporating all
physical mass-scales that are relevant for describing ($N_{\rm f}=6$) QCD—as
it seems to be realised in nature—is an impossible task to present state-of-
the-art simulations due to the immense costs to accommodate many orders of
magnitude of energy scales in a single simulation, cf. Figure 1. Being mainly
interested in non-perturbative, long-distance effects of QCD, the
phenomenological applicability of LQCD is limited to a fixed window111for
convenience, we shall refer to all physical scales well separable from the
cutoffs as ’low-energy’ window in the following, i.e., $\mu\ll 1/a\sim
2-4\,{\rm GeV}$

$\displaystyle 40\,{\rm MeV}\;\sim\;\Lambda_{\tiny\rm IR}\;\ll\;\mu_{\rm
had},\,m_{\tiny\pi^{0}},\,m_{\rm p},\,m_{\tiny\rm
K^{0}},\,\ldots\;\ll\;\Lambda_{\tiny\rm UV}\;\sim\;3\,{\rm GeV}\;.$ (2)

Any physically interesting scale $\mu$ has to be well separated from the
imposed infrared and ultraviolet cutoffs, $\Lambda_{\tiny\rm
IR}\ll\mu\ll\Lambda_{\tiny\rm UV}$, to allow for a controlled assessment of
the respective cutoff effects, see caption of Figure 1. Although, the lattice
community can afford to simulate with three or four dynamical flavours
nowadays, one still needs to balance the $N_{\rm f}+1$ chosen (bare) input
parameters, given in terms of $N_{\rm f}+1$ physical scales from the real
spectrum, against the intrinsic cutoff values and algorithmic costs. As no
practical perfect lattice setup exists, different views on how to best achieve
this balance lead to controversies among lattice practitioners. Another of
those disputable topics is the application of perturbation theory at energies
as low as those given by the low-energy window of lattice QCD as specified in
(2), including $\mu=\Lambda_{\tiny\rm UV}\sim\,3\,{\rm GeV}$, see also
Brida:2016flw .

In course of our determination of quark masses, we have to fix a
renormalization prescription at scale $\mu_{\rm had}$ that lies well inside
the window $[\Lambda_{\tiny\rm IR},\Lambda_{\tiny\rm UV}]$. Then the natural
question arises how to compare or connect our result to other determinations
as summarised for instance by the Particle Data Group PDG:2014 . Both
naturally differ in the renormalization scheme and scale. After removing the
UV cutoff dependence, we thus have to evolve our (continuum) result, $\kern
1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern 1.0pt_{i}(\mu_{\rm had})$, using
renormalization group (RG) transformations to the same physical scale and
scheme, typically $\kern 1.0pt\overline{\kern 0.0pt\rm MS\kern-1.0pt}\kern
1.0pt$.

Figure 1: Sketch of energy scales relevant to QCD phenomenology. The masses
for the six quark flavours found in nature so far (u, d, s, c, b, t) are taken
as quoted in PDG:2014 . Also the masses of the neutral pseudoscalar meson
bound states ($m_{\tiny\pi^{0}}$, $m_{\tiny\rm K^{0}}$, $m_{\tiny\rm D^{0}}$,
$m_{\tiny\rm B^{0}}$) as well as the depicted scale dependence of the strong
coupling $\alpha_{\rm s}$ are taken from that reference. Although,
$\alpha_{\rm s}(\mu)$ can be computed in any sensible renormalization scheme
one conveniently quotes it for the 5-flavour theory at the electroweak scale
$\mu=m_{\tiny\rm Z}$ in the $\kern 1.0pt\overline{\kern 0.0pt\rm
MS\kern-1.0pt}\kern 1.0pt$ scheme.
The lattice regularisation of QCD provides the only known (non-perturbative
and gauge-invariant) framework to compute bound states of the strong
interaction. For phenomenological applications one is restricted to a window
($\Lambda_{\rm IR}\leq\mu\leq\Lambda_{\rm UV}$) in the low-energy regime of
QCD in order to incorporate the important long-distance effects. They are of
the order of the Compton wavelength of the lightest particle,
$\lambda_{\pi^{0}}=1/m_{\tiny\pi^{0}}$. It has to be well separated from the
infrared cutoff $\Lambda_{\rm IR}^{-1}\sim L\gg\lambda_{\pi^{0}}$ in order to
not be distorted by finite volume effects. Typical simulations have a physical
extent of $L\gtrsim 3\,{\rm fm}$ and a number of lattice points of
$N=(L/a)^{4}\approx(50-100)^{4}$, thus restricting the lattice spacing to
$a\sim\Lambda_{\rm UV}^{-1}\gtrsim 0.045\,{\rm fm}$. To extract physical
quantities one also has to stay away from the ultraviolet cutoff, i.e.
$\mu\ll\Lambda_{\rm UV}$, otherwise control over finite lattice spacing
effects is lost, and the continuum limit $a\to 0$ cannot be taken. Although,
various advances in LQCD steadily increase the quality and reliability of
numerical simulations with different number of flavours $N_{\rm f}$, the cost
of simulations and algorithmic difficulties still prevent us from reaching
even smaller lattice spacings.
For any non-perturbative renormalization problem posed at fixed
renormalization scale $\mu$ in the low-energy regime of LQCD the continuum
limit has to be taken in a controlled way. Then the remaining question is how
to relate the obtained results to experiments at much higher energies,
typically two orders of magnitude. For high precision physics the supposedly
most convenient renormalization scheme, $\kern 1.0pt\overline{\kern 0.0pt\rm
MS\kern-1.0pt}\kern 1.0pt$, is of little help at scales much below $10\,{\rm
GeV}$ (or $\alpha_{\rm s}\gtrsim 0.2$). Its intrinsic perturbative nature, the
behaviour of asymptotic series expansions in general, and missing non-
perturbative contributions prevent us from quantifying its (non-)reliability
at low energies from within the scheme itself. As detailed in the main text
the renormalization group running of any operator can be determined purely
non-perturbative, thus connecting the low-energy regime of LQCD in a
controlled and quantifiable way to scales well above $10\,{\rm GeV}$.
Recently, it has become possible to non-perturbatively test the accuracy of
continuum perturbation theory for $\alpha_{\rm s}\lesssim 0.2$ Brida:2016flw ,
stressing the relevance and necessity of a careful assessment beyond
perturbation theory in present high precision physics. Non-perturbatively it
is even possible nowadays to follow the renormalization group in QCD down to
scales of about $200\,{\rm MeV}$, cf. Ref. DallaBrida:2016kgh . Both results
have been reported at this conference Sint:12Conf2016 ; DallaBrida:12Conf2016
.

In a mass-independent renormalization scheme for QCD, the renormalization
group equations (RGE) for the running coupling and running quark mass(es) read

$\displaystyle\mu\dfrac{\partial}{\partial\mu}\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt(\mu)$
$\displaystyle=\beta(\kern 1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern
1.0pt)\;,$ $\displaystyle\dfrac{\mu}{\kern
1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern
1.0pt}\dfrac{\partial}{\partial\mu}\kern
1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern 1.0pt(\mu)$
$\displaystyle=\tau(\kern 1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern
1.0pt)\;,$ (3)

and can be formally integrated. The resulting, exact solutions are known as
renormalization group invariants (RGI), conveniently written as

$\displaystyle{\Lambda}$ $\displaystyle\equiv\mu\big{[}{b_{0}\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern
1.0pt^{2}(\mu)}\big{]}^{-{b_{1}}/({2b_{0}^{2}})}\,{\rm e}^{-1/({2b_{0}\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern
1.0pt^{2}(\mu)})}\exp\bigg{(}\\!-\\!\\!\int_{0}^{\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern
1.0pt^{\vphantom{A}}(\mu)}\\!\\!\\!\\!{\rm
d}{g}\,\bigg{[}\dfrac{1}{\beta(g)}+\dfrac{1}{b_{0}g^{3}}-\dfrac{b_{1}}{b_{0}^{2}g}\bigg{]}\bigg{)}\;,$
(4) $\displaystyle{M_{i}}$ $\displaystyle\equiv\kern
1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern 1.0pt_{i}(\mu)\big{[}{2b_{0}\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern
1.0pt^{2}(\mu)}\big{]}^{-{d_{0}}/({2b_{0}})}\,\exp\bigg{(}\\!-\\!\\!\int_{0}^{\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern
1.0pt^{\vphantom{A}}(\mu)}\\!\\!\\!\\!{\rm
d}{g}\,\bigg{[}\dfrac{\tau(g)}{\beta(g)}-\dfrac{d_{0}}{b_{0}g}\bigg{]}\bigg{)}\;,\quad
i\in\\{\rm u,d,s,c,b,t\\}\,.$ (5)

Only the leading, universal (scheme-independent) coefficients
$b_{0},b_{1},d_{0}$ of the perturbative expansions of the beta-function and
mass-anomalous dimension,

$\displaystyle\beta(g)$ $\displaystyle\mathrel{\mathop{\sim}^{g\to
0}}-g^{3}({b_{0}}+{b_{1}}g^{2}+{b_{2}^{s}}g^{4}+\ldots)\;,$
$\displaystyle\tau(g)$ $\displaystyle\mathrel{\mathop{\sim}^{g\to
0}}-g^{2}({d_{0}}+{d_{1}^{s}}g^{2}+\ldots)\;,$ (6)

appear such that the individual integrands are finite when $g\to 0$. Higher
order coefficients are scheme-dependent—indicated by the superscript $s$—and
so are $\beta$ and $\tau$. Thanks to recently completed efforts these
approximations are now known to 5-loop order in the minimal subtraction
scheme(s), cf. Gross:1973id ; Politzer:1973fx ; Caswell:1974gg ;
Tarasov:1980au ; vanRitbergen:1997va ; Czakon:2004bu ; Baikov:2016tgj and
Tarrach:1980up ; Tarasov:1982gk ; Chetyrkin:1997dh ; Vermaseren:1997fq ;
Baikov:2014qja . However, it is well-known that by using an asymptotic series
expansions one is missing non-perturbative contributions222typically
associated with instantons, renormalons, … which become increasingly more
relevant towards low energies, i.e., when the strong coupling becomes large.
One has to appreciate that these issues can be overcome by determining
$\beta(g)$ and $\tau(g)$ non-perturbatively. In that case, eqs. (4) and (5)
determine the fundamental parameters $\\{\Lambda,M_{i}\\}_{N_{\rm f}}$
uniquely for any input scale $\mu$.

##  2 An intermediate non-perturbative renormalization scheme

Figure 2: Non-perturbative $\beta$-function for the SF ($\alpha_{\rm
SF}\lesssim 0.2$) and GF ($\alpha_{\rm GF}\gtrsim 0.17$) running coupling
scheme; picture from DallaBrida:2016kgh . There also the two schemes have been
matched non-perturbatively at fixed physical volume corresponding to about
$4\,{\rm GeV}$.

For massless renormalization schemes (RS) it is most natural to first solve
the RG equation for the coupling and then for the mass, cf. (3). In LQCD it
has become customary to determine the scale evolution of coupling and mass
using a _finite-volume renormalization scheme_ in combination with recursive
finite-size scaling. For that purpose the physical size $L$ of the simulated
volume $L^{4}$ is identified with the inverse renormalization scale for the
sole purpose of solving the RG equations by rescaling the box size, $L\to sL$.
The constant factor $s$ is usually chosen to be $2$. In that way one is able
to cleanly disentangle the renormalization procedure from the low-energy
window of LQCD and to connect the hadronic renormalization scale $\mu_{\rm
had}=1/L_{\rm had}$ to the electroweak scale $\mu=m_{\tiny\rm Z}$, where for
instance a change of renormalization schemes towards the commonly quoted
$\kern 1.0pt\overline{\kern 0.0pt\rm MS\kern-1.0pt}\kern 1.0pt$ scheme can
typically be achieved in a controlled way,333for the time being we stay at
fixed $N_{\rm f}$ and neglect changes from $N_{\rm f}=3\to 5$ cf. Figure 1.
For that purpose one should ideally cover a range of scales compatible with

$\displaystyle L_{\rm max}$ $\displaystyle\geq L_{\rm had}\;,$ $\displaystyle
L_{\rm min}$ $\displaystyle\leq 1/m_{\tiny\rm Z}\;,$ $\displaystyle L_{\rm
max}/L_{\rm min}$ $\displaystyle=s^{N_{s}}\;,$ $\displaystyle s>1\;.$ (7)

For example, to connect a scale $\mu_{\rm had}$ as low as $200\,{\rm MeV}$ to
the electroweak scale of about $100\,{\rm GeV}$ using $s=2$, at least
$N_{s}=9$ individual steps are necessary. Each of the $N_{s}$ steps requires a
series of lattice simulations at different resolutions $L/a$ but fixed
renormalized input parameters, $\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt^{2}(L)=u$ and $L\kern
1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern 1.0pt_{i}(L)=0$. In this way
$\mu=1/L$ is kept fixed while only the lattice spacing $a$ is being varied, or
equivalently, the bare parameters $\\{g_{0}^{2},m_{0,i}\\}$ as a function of
$a$. A second set of simulations at the same bare parameters but $sL/a$ number
of points in each direction permits to trace the change of renormalization
scale $L\to sL$ at finite $a$. With adequate renormalization conditions for
the coupling and quark mass, this allows to take the continuum limit $a\to 0$
of their lattice approximants in order to determine the response w.r.t. a
change of scale $s$ in the continuum, equivalent to

$\displaystyle-\ln(s)$ $\displaystyle=\int^{\,\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt(sL)}_{\,\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt(L)}\\!\\!\frac{{\rm
d}g}{\beta(g)}\;,$ $\displaystyle\ln(\sigma_{\rm P})$
$\displaystyle=-\int^{\,\kern 1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern
1.0pt(sL)}_{\,\kern 1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern
1.0pt(L)}\\!\\!\\!{\rm d}g\frac{\tau(g)}{\beta(g)}\;,$ for fixed
$\displaystyle\quad\begin{cases}\;\;\;\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt^{2}(L)=u\\\ \,L\kern
1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern 1.0pt_{i}(L)=0\end{cases}\;.$ (8)

In this way, the relevant information is encoded in the _step-scaling
functions_ (SSFs) Luscher:1993gh ; Capitani:1998mq ,

$\displaystyle\sigma(s,u)$ $\displaystyle=\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt^{2}(sL)\;,$
$\displaystyle\sigma_{\rm P}(s,u)$ $\displaystyle=\kern
1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern 1.0pt_{i}(L)\big{/}\kern
1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern 1.0pt_{i}(sL)\;.$ (9)

While various ways exist to determine the SSFs numerically via lattice
simulations of QCD, the so-called _Schrödinger functional_ (SF) setup
Luscher:1992an ; Sint:1993un ; Sint:1995rb has been develop for exactly that
purpose. For instance, in contrast to standard (anti-)periodic boundary
conditions in all space time directions, Dirichlet boundary conditions are
imposed in time direction. This allows for a natural non-perturbative
definition of a strong coupling via non-vanishing boundary gauge fields as
well as to simulate at vanishing quark mass. For additional details we point
to the review article Sommer:2015kza and references therein.

The scale evolution of the strong coupling has been determined in refs.
Brida:2016flw ; DallaBrida:2016kgh , and a short account of the whole
procedure can be found in Bruno:2016gvs . However, a technical complication is
inherited from that determination: instead of a single RS, two different
schemes are being used and matched non-perturbatively at a scale of approx.
$4\,{\rm GeV}$. At present, its physical motivation is two-fold: 1) control
the perturbative matching at the electroweak scale to high precision
Brida:2016flw using the Schrödinger functional coupling, $\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt_{\rm SF}$, and 2) reach
good statistical accuracy down to energies of about $200\,{\rm MeV}$
DallaBrida:2016kgh by applying the Gradient flow coupling scheme, $\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt_{\rm GF}$ Fritzsch:2013je ;
Ramos:2015baa . Compared to previous estimates of $\alpha_{\rm s}(m_{\tiny\rm
Z})$ this represents a new quality of rigor from lattice QCD determinations
and for the first time allowed to accurately determine the corresponding non-
perturbative QCD $\beta$-functions, see Figure 2.

For our determination of the scale evolution of renormalized quark masses in
three-flavour QCD, this does not pose any additional obstacles but only
complicates the overall presentation. In fact, assuming a diagonal quark mass
matrix of rank $N_{\rm f}$ in two massless schemes, one can write Sint:1998iq

$\displaystyle\mu^{\prime}$ $\displaystyle=c\mu\,,\;c>0,$ $\displaystyle\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt^{\prime 2}$
$\displaystyle=\kern 1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern
1.0pt^{2}{\cal X}_{\rm g}(\kern 1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern
1.0pt)\,,$ $\displaystyle\kern 1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern
1.0pt_{j}^{\prime}$ $\displaystyle=\kern
1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern
1.0pt_{j}^{\vphantom{{}^{\prime}}}\,{\cal X}_{\rm
m}^{\vphantom{{}^{\prime}}}(\kern 1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern
1.0pt)\,,$ $\displaystyle j=1,\ldots,N_{\rm f}\,.$ (10)

Invariance of a physical observable $P$ under this change of variables implies

$\displaystyle P^{\prime}\big{[}\,\mu^{\prime}(\mu),\,\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt^{\prime}(\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt),\,\\{\kern
1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern 1.0pt_{j}^{\prime}(\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt,\kern
1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern
1.0pt_{j}^{\vphantom{{}^{\prime}}})\\}\,\big{]}$
$\displaystyle=P\,\big{[}\,\mu,\,\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt,\,\\{\kern
1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern 1.0pt_{j}\\}\big{]}\;,$ (11)

where $P^{\prime}$ satisfies the Callan–Symanzik equation in the primed scheme
w.r.t.

$\displaystyle\beta^{\prime}(\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt^{\prime})$
$\displaystyle=\left\\{\beta(\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt)\frac{\partial\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt^{\prime}}{\partial\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt}\right\\}_{\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt=\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt(\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt^{\prime})}\;,$
$\displaystyle\tau^{\prime}(\kern 1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern
1.0pt^{\prime})$ $\displaystyle=\left\\{\tau(\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt)+\beta(\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern
1.0pt)\frac{\partial}{\partial\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt}\ln{\cal X}_{\rm m}(\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt)\right\\}_{\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt=\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt(\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt^{\prime})}\;.$ (12)

While this is the general connection between two massless schemes at the level
of RG functions, every step in our determination is done non-perturbatively.
In practise, when switching between the two schemes, we do not change the
imposed renormalization condition—which is relevant for determining
$\sigma_{\rm P}(2,u)$—but merely the implicit, parametric dependence on the
renormalized coupling $u\equiv\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt^{2}(L)$. The connection
between SF and GF schemes has been established non-perturbatively at a well-
chosen fixed physical box size, $L_{0}\approx(4\,{\rm GeV})^{-1}$, and reads
DallaBrida:2016kgh

$\displaystyle\kern 1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt_{\rm
SF}^{2}(L_{0})$ $\displaystyle=2.012\;,$ $\displaystyle\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt_{\rm GF}^{2}(2L_{0})$
$\displaystyle=2.6723(64)\;.$ (13)

##  3 The running quark mass

Figure 3: _Left:_ Representative continuum extrapolations ($u$-by-$u$) for the
step-scaling function $\sigma_{\rm P}(u)\equiv\lim_{a/L\to
0}\Sigma^{(1)}(u,a/L)$ in the SF scheme. For further details please consult
Campos:2015fka ; Lattice16 . _Right:_ Non-perturbatively determined step-
scaling functions $\sigma_{\rm P}(u)$ in the SF ($u\lesssim 2.0$) and GF
($u\gtrsim 2.1$) running coupling scheme, together with different
interpolating fits. We note that for $u_{\rm SF}\lesssim 1.5$ the data agrees
with NLO perturbation theory such that the matching with SF-PT can be safely
established.

Given the discussions in the previous sections, our strategy of a non-
perturbatively controlled determination of quark masses in the three-flavour
theory with two renormalization schemes proceeds as follows. We determine the
flavour-independent RG factor ${M}/{\kern
1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern 1.0pt(\mu_{\rm had})}$, cf. eq.
(5), that connects the RGI mass to the quark mass in the hadronic, low-energy
regime

$\displaystyle\frac{M}{\kern 1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern
1.0pt(\mu_{\rm had})}$ $\displaystyle=\left[\frac{\kern
1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern 1.0pt(\mu_{\rm swi})}{\kern
1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern 1.0pt(\mu_{\rm had})}\right]_{\rm
GF-NP}\\!\\!\times\;\left[\frac{\kern
1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern 1.0pt(\mu_{\rm PT})}{\kern
1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern 1.0pt(\mu_{\rm swi})}\right]_{\rm
SF-NP}\\!\\!\times\;\left[\frac{M}{\kern
1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern 1.0pt(\mu_{\rm PT})}\right]_{\rm
SF-PT}\;.$ (14)

Note that we intend to make use of perturbation theory at a scale $\mu_{\rm
PT}$ close to $m_{\tiny\rm Z}$ where truncation errors from using the known
3-/2-loop orders in the perturbative $\beta$-/$\tau$-functions in the SF
scheme are sufficiently suppressed,444at high energies the (non-perturbative)
SF scheme is parametrically close to the (perturbative) $\kern
1.0pt\overline{\kern 0.0pt\rm MS\kern-1.0pt}\kern 1.0pt$ scheme  see Figure 3.
The true challenge is to determine the two non-perturbative factors to a
sufficient precision while controlling systematic effects at the sub-percent
level for each individual contribution. When this is achieved, the scale- and
scheme-independent RGI mass of flavour $j$ can be determined from

$\displaystyle M_{j}$ $\displaystyle=\frac{M}{\kern
1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern 1.0pt(\mu_{\rm had})}\times\kern
1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern 1.0pt_{j}(\mu_{\rm had})\;,$ (15)

where $\kern 1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern 1.0pt_{j}(\mu_{\rm
had})$ constitutes a renormalized quark mass computed in the hadronic regime.
Ideally, $M_{j}$ is the fundamental parameter to be compared to other
determinations. However, it has become customary to quote $\kern
1.0pt\overline{\kern 0.0pt\rm MS\kern-1.0pt}\kern 1.0pt$ masses. In order to
do so one subsequently has to employ perturbation theory and, depending on the
individual flavour and renormalization scale, appropriately match at the charm
and bottom thresholds.

Figure 4: Non-perturbative mass anomalous dimension $\tau(g^{2})$ in the SF
($\kern 1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt^{2}_{\rm
SF}\lesssim 2.45$) and GF ($\kern 1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern
1.0pt^{2}_{\rm GF}\gtrsim 2.10$) running coupling scheme. The vertical lines
in the inner plot represent the respective couplings at the switching scale.
Note that for $u=\kern 1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt_{\rm
GF}^{2}$ our data exactly covers the plotted fit range.

An important fact worth emphasising is that the continuum running factor in
eq. (14) can be used by other lattice practitioners which may employ a
different lattice discretization in their computations. They just have to
recompute $\kern 1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern
1.0pt_{j}(\mu_{\rm had})$ using the same—but so far
unmentioned—renormalization prescription, such that eq. (15) remains well-
defined. For the employed Wilson fermion action, the quark mass $\kern
1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern 1.0pt_{j}$ can be obtained, on
the lattice up to a finite renormalization, from the non-singlet axial Ward
identity, or current quark mass

$\displaystyle\kern 1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern
1.0pt_{j}^{\rm\tiny awi}(\mu)$ $\displaystyle=\lim_{a/L\to
0}\left[\frac{{Z_{\rm A}}(g_{0}^{2})}{{Z_{\rm
P}}(g_{0}^{2},a\mu)}m_{j}^{\rm\tiny awi}\right]_{\mu=L^{-1}}\;,$
$\displaystyle\partial_{\mu}\left[\overline{\psi}_{j}\gamma_{\mu}\gamma_{5}\psi_{j^{\prime}}(x)\right]$
$\displaystyle=(m_{j}^{\rm\tiny awi}+m_{j^{\prime}}^{\rm\tiny
awi})\,\overline{\psi}_{j}(x)\gamma_{5}\psi_{j^{\prime}}(x)\;.$ (16)

The bare mass $m_{j}^{\rm\tiny awi}\equiv m_{j}^{\rm\tiny
awi}(g_{0}^{2},\\{am_{0,i}\\})$ still depends on the bare parameters of the
action, cf. eq. (1). From the previous discussions the connection to the mass-
SSF and notation should become evident

$\displaystyle\sigma_{\rm P}(s,u)\equiv\frac{\kern
1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern 1.0pt_{j}^{\rm\tiny
awi}(L)}{\kern 1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern
1.0pt_{j}^{\rm\tiny awi}(sL)}$ $\displaystyle=\lim_{a/L\to
0}\left[\frac{{Z_{\rm P}}(g_{0}^{2},sL/a)}{{Z_{\rm
P}}(g_{0}^{2},L/a)}\right]_{\kern 1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern
1.0pt^{2}(L)=u,m=0}\;,$ (17)

cf. eqs. (8) and (9). ${Z_{\rm P}}$ is the multiplicative renormalization
factor of the non-singlet pseudoscalar current, which itself is determined
from a SF renormalization condition as in Capitani:1998mq .

In Figure 3 (left panel) we show preliminary results for some continuum
extrapolations as in eq. (17) with $s=2$. The corresponding continuum SSF data
points are then plotted against the initially fixed target couplings $u$ for
both the SF- and GF-coupling schemes (right panel). Note that these values of
$u$ can be chosen at will but should cover the range of interest well enough
in order to reach the desired precision. Typical fit ansaetze as shown in Fig.
3 are polynomials like $\sigma_{\rm P}(u)=1+p_{\rm LO}u+p_{2}u^{2}+\cdots$ or
Padé approximants with the correct, leading (universal) asymptotics fixed.
With such a given smooth interpolation formula one is then able to recursively
construct the running factor and its uncertainty in steps of $s=2$ as
summarised in (7), i.e., the two ratios in (14) become

$\displaystyle\frac{\kern 1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern
1.0pt(\mu_{a})}{\kern 1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern
1.0pt(\mu_{b})}$ $\displaystyle=\prod_{k=1}^{N_{s}}\sigma_{\rm P}(u_{k})\;,$
$\displaystyle\begin{cases}u_{k+1}=\sigma(u_{k})\\\ u_{1}\quad\\!=\kern
1.0pt\overline{\kern-1.0ptg\kern-1.0pt}\kern 1.0pt^{2}(\mu_{a})\\\
\end{cases}\\!\\!\\!\\!,$ $\displaystyle s_{ab}$
$\displaystyle\equiv\frac{\mu_{a}}{\mu_{b}}=\frac{L_{b}}{L_{a}}=2^{N_{s}}\;.$
(18)

Accordingly, we can identify the switching scale $L_{\rm swi}$ appearing in
eq. (14) with the scale $2L_{0}$ of (13), which we inherit together with the
coupling-SSF $\sigma(u)$ from the complementary determination of the running
coupling Brida:2016flw ; DallaBrida:2016kgh ; Bruno:2016gvs ;
DallaBrida:12Conf2016 ; Sint:12Conf2016 . This concludes the computation of
(14) and we return to the determination of eq. (15), which we can now rewrite
as

$\displaystyle M_{j}$ $\displaystyle=\frac{M}{\kern
1.0pt\overline{\kern-1.0ptm\kern-1.0pt}\kern 1.0pt(\mu_{\rm
had})}\times\lim_{a/L\to 0}\left[\frac{{Z_{\rm A}}(g_{0}^{2})}{{Z_{\rm
P}}(g_{0}^{2},a\mu_{\rm had})}m_{j}^{\rm\tiny awi}\right]\equiv\lim_{a/L\to
0}\left[{Z_{\rm M}}(g_{0}^{2})\,m_{j}^{\rm\tiny
awi}(g_{0}^{2},\\{am_{0,i}\\})\right]\;.$ (19)

The disappearance of the hadronic scale $\mu_{\rm had}$ reflects the scale-
independence of $M_{j}$, which itself is connected to the current quark mass
$m_{j}^{\rm\tiny awi}$ only by a scale-independent renormalization ${Z_{\rm
M}}$. Hence, by determining ${Z_{\rm M}}$ the renormalization problem is fully
solved non-perturbatively using the massless (intermediate) SF scheme. While
the costs of simulating the Schrödinger functional setup are good to control
in general, towards larger physical volumes the simulations get more and more
involved and thus set a natural limit on the values the hadronic scale can
take. At the moment we are exploring different ways to fix the value of
$a\mu_{\rm had}$ in order to determine the ${Z_{\rm P}}(g_{0}^{2},a\mu_{\rm
had})$, and thus ${Z_{\rm M}}$, with high precision.
As the reader may have noticed, forcing $s_{ab}$ in eq. (18) to be a multiple
of $2$ is a quite stringent condition, especially as we already had to fix the
scheme-switching scale in (14). Thanks to the high statistical accuracy and
the achieved control of a few systematic effects which are propagated into the
final error, we are able to lift this restriction for the first time. By
reassessing eqs. (8) and (9), one notices that with the data points at hand
and a given parameterisation of the $\beta$-function, as shown in Figure 2 and
taken from Brida:2016flw ; DallaBrida:2016kgh , only a well motivated ansatz
for $\tau(g)$ is required for a sensible least square minimisation. We present
a preliminary but already very encouraging result in Figure 4, again for both
the SF- and GF-running coupling schemes. It should be noted that this result
is actually independent of the originally chosen scale-change factor $s=2$. As
a non-trivial cross-check we can use the resultant to numerically reconstruct
$\sigma_{\rm P}(u)$ for arbitrary $u$ from eq. (8). This is shown in Figure 3
and coincides very well with the original estimate. At the time of this
conference we were still accumulating statistics at the two strongest
renormalized couplings in use. Now we are finalising our analysis by also
employing global fit ansaetze in various ways and include some yet missing
correlations between our data, which is going to be published soon.
Accordingly, we here refrain from quoting any quantitative numbers.

Compared to the aforementioned simulations, the $m_{j}^{\rm\tiny awi}$
estimates in (19) are obtained from large-scale lattice simulations that
define the accessible low-energy window in the first place. We will employ
$N_{\rm f}=2+1$ ensembles jointly produced by the Coordinated Lattice
Simulations (CLS) effort Bruno:2014jqa , with which we share the lattice
discretization in use. Although there have been many advances in lattice QCD,
the accessible parameters $am_{0,i}$, in general, do not correspond directly
to the set of parameters at which $m_{\pi}$ and $m_{\rm K}$ take their
physical values for different values of the cutoffs (even after correcting for
isospin splitting and QED effects). Beside setting the physical scale
Bruno:2016plf , one thus still has to apply chiral perturbation theory. In
that respect, eq. (19) is a minor simplification in our discussion, but
nevertheless a step that has to be carefully evaluated in the near future.

##  4 Conclusions

We have presented a full strategy to determine renormalized quark masses from
first principle lattice calculations in three-flavour QCD. Compared to other
approaches in the community, we separate and thus disentangle the problem of
renormalization from large-scale lattice QCD simulations. In this way the
massless renormalization group running for the coupling and quark masses could
be mapped out non-perturbatively to high accuracy in the employed Schrödinger
functional scheme, and a connection to other schemes can be easily established
via renormalization group invariant quantities $\\{\Lambda,M_{i}\\}_{N_{\rm
f}}$. This strategy fully circumvents the unanswerable question regarding the
applicability of perturbation theory at parametrically large values of the
strong coupling. Furthermore, it provides a quantifiable and systematically
improvable uncertainty at each stage of the calculation, something
perturbation theory cannot really provide.

For the first time two different non-perturbative schemes, based on the
Schrödinger functional coupling and the younger Gradient flow coupling, have
been used and matched non-perturbatively. This specific combination has been
chosen for the avail of improving on statistical and systematic uncertainties.
Also the non-perturbative determination of the continuum quark mass anomalous
dimension is a unique achievement for QCD.

The prospects for determining fundamental parameters of QCD valuable for high
precision physics are excellent. To obtain an even more realistic
approximation of nature this study can be extended to cross the charm flavour
threshold $(N_{\rm f}=3\to N_{\rm f}=4)$ non-perturbatively at some point in
the future. We finally remark that also the renormalization and scale
dependence of the tensor current—relevant in rare meson decays (especially B
decays), and studies of the neutron electric dipole moment—is being computed
along the same lines Fritzsch:2015lvo .

## Acknowledgments

The simulations were jointly performed on the Altamira HPC facility, on
Finisterrae-2 at CESGA, and on the GALILEO supercomputer at CINECA (INFN
agreement). We thankfully acknowledge the computer resources and technical
support provided by the University of Cantabria (IFCA), CESGA and CINECA.

P.F. acknowledges financial support from the Spanish MINECO’s “Centro de
Excelencia Severo Ochoa” Programme under grant SEV-2012-0249, as well as from
the grant FPA2015-68541-P (MINECO/FEDER).

## References

  * (1) M. Dalla Brida, P. Fritzsch, T. Korzec, A. Ramos, S. Sint, R. Sommer (ALPHA), Phys. Rev. Lett. 117, 182001 (2016), 1604.06193
  * (2) K. Olive et al. (Particle Data Group), Chinese Physics C 38, 090001 (2014)
  * (3) M. Dalla Brida, P. Fritzsch, T. Korzec, A. Ramos, S. Sint, R. Sommer (ALPHA) (2016), 1607.06423
  * (4) S. Sint (ALPHA), EPJ Web of Conferences CONF12 (2016), this conference
  * (5) M. Dalla Brida (ALPHA), EPJ Web of Conferences CONF12 (2016), this conference
  * (6) D.J. Gross, F. Wilczek, Phys.Rev.Lett. 30, 1343 (1973)
  * (7) H.D. Politzer, Phys.Rev.Lett. 30, 1346 (1973)
  * (8) W.E. Caswell, Phys.Rev.Lett. 33, 244 (1974)
  * (9) O. Tarasov, A. Vladimirov, A.Y. Zharkov, Phys.Lett. B93, 429 (1980)
  * (10) T. van Ritbergen, J. Vermaseren, S. Larin, Phys.Lett. B400, 379 (1997), hep-ph/9701390
  * (11) M. Czakon, Nucl.Phys. B710, 485 (2005), hep-ph/0411261
  * (12) P.A. Baikov, K.G. Chetyrkin, J.H. KÃhn (2016), 1606.08659
  * (13) R. Tarrach, Nucl.Phys. B183, 384 (1981)
  * (14) O.V. Tarasov (1982), in Russian, preprint JINR P2-82-900
  * (15) K. Chetyrkin, Phys.Lett. B404, 161 (1997), hep-ph/9703278
  * (16) J. Vermaseren, S. Larin, T. van Ritbergen, Phys.Lett. B405, 327 (1997), hep-ph/9703284
  * (17) P.A. Baikov, K.G. Chetyrkin, J.H. KÃhn, JHEP 10, 076 (2014), 1402.6611
  * (18) M. Lüscher, R. Sommer, P. Weisz, U. Wolff, Nucl.Phys. B413, 481 (1994), hep-lat/9309005
  * (19) S. Capitani, M. Lüscher, R. Sommer, H. Wittig (ALPHA), Nucl.Phys. B544, 669 (1999), hep-lat/9810063
  * (20) M. Lüscher, R. Narayanan, P. Weisz, U. Wolff, Nucl.Phys. B384, 168 (1992), hep-lat/9207009
  * (21) S. Sint, Nucl.Phys. B421, 135 (1994), hep-lat/9312079
  * (22) S. Sint, Nucl.Phys. B451, 416 (1995), hep-lat/9504005
  * (23) R. Sommer, U. Wolff (ALPHA), Nucl.Part.Phys.Proc. 261-262, 155 (2015), 1501.01861
  * (24) M. Bruno et al. (ALPHA), _The determination of $\alpha_{s}$ by the ALPHA collaboration_, in _Proceedings, 6th Workshop on Theory, Phenomenology and Experiments in Flavour Physics : Interplay of Flavour Physics with electroweak symmetry breaking. (Capri 2016): Anacapri, Capri, Italy, June 11-13, 2016_ (2017), Vol. 285-286, pp. 132–138, 1611.05750
  * (25) P. Fritzsch, A. Ramos, JHEP 1310, 008 (2013), 1301.4388
  * (26) A. Ramos, S. Sint, Eur. Phys. J. C76, 15 (2016), 1508.05552
  * (27) S. Sint, P. Weisz (ALPHA), Nucl.Phys. B545, 529 (1999), hep-lat/9808013
  * (28) I. Campos et al., PoS LATTICE2015, 249 (2016), 1508.06939
  * (29) I. Campos et al., PoS LATTICE2016, 210 (2016)
  * (30) M. Bruno et al., JHEP 1502, 043 (2015), 1411.3982
  * (31) M. Bruno, T. Korzec, S. Schaefer (2016), 1608.08900
  * (32) P. Fritzsch, C. Pena, D. Preti, PoS LATTICE2015, 250 (2016), 1511.05024