$\begin{split}F\left[x,\,y\right]\sim&\left\\{\begin{array}[]{lll}\displaystyle\frac{q-1}{x_{\text{IR}}^{1-q}}\,\frac{1}{x^{q}}&\text{for}&x_{0}\sqrt{\xi}=x_{\text{IR}}\leq
x\leq y\\\\[15.0pt] \hskip 10.11775pt0&&\text{otherwise}\leavevmode\nobreak\
.\end{array}\right.\end{split}$ (63)

In terms of arguments in Eq. (62), the finite support of the function $F$ is
defined by the interval $(x_{0}^{2}\,c_{1}\log
y^{\prime})^{1/2}\leq(y^{\prime}/y)^{1/2}x\leq y^{\prime}$ and, subsequently,
it defines the integration range $y_{1}\leq y^{\prime}\leq y_{2}$ in which $F$
in Eq. (62) is nonvanishing. There are two choices of $(y_{1},\,y_{2})$
depending on the value of $y^{-1/2}x=\frac{k}{\sqrt{m_{r}H}}$.

When $y^{-1/2}x$ is large (large momentum region), the integration is done
over $x^{2}/y=y_{1}\leq y^{\prime}\leq y_{2}=y$ where $x^{2}/y$ is the
intersection between two curves $y^{\prime}$ and $(y^{\prime}/y)^{1/2}x$. The
integration in Eq. (62) can be done straightforwardly:

$\begin{split}\frac{\partial\rho_{a}}{\partial k}(k,\,t)&=8\pi
f_{a}^{2}c_{1}m_{r}(x_{0}^{2}c_{1})^{\frac{q-1}{2}}y^{\frac{q-3}{2}}\frac{1}{x^{q}}\left(\frac{q-1}{2}\right)^{-\frac{q+3}{2}}\\\\[5.0pt]
&\quad\times\left[\Gamma\left(\frac{q+5}{2},\,\frac{q-1}{2}\log\frac{x^{2}}{y}\right)-\Gamma\left(\frac{q+5}{2},\,\frac{q-1}{2}\log
y\right)\right]\leavevmode\nobreak\ ,\end{split}$ (64)

where $\Gamma(s,x)=\int_{x}^{\infty}dt\,t^{s-1}e^{-t}$ is the incomplete Gamma
function. The variable $x^{2}/y=\frac{k^{2}}{m_{r}H}$ for the fixed momentum
$k$. When $\frac{m_{r}}{H}$ is taken to be large, the second argument of
$\Gamma$ function in Eq. (64) becomes large as well. Since
$\Gamma(s,\,x)\rightarrow x^{s-1}e^{-x}$ as $x\rightarrow\infty$, the result
in Eq. (64) can be approximated as, in the large $m_{r}/H$ limit,

$\begin{split}\frac{\partial\rho_{a}}{\partial
k}(k,\,t)\approx\frac{8H\mu_{\text{eff}}\sqrt{\xi}}{x_{0}}\left(\frac{Hx_{0}\sqrt{\xi}}{k}\right)^{q}f\left(\frac{m_{r}}{H},\,\frac{k}{m_{r}}\right)\leavevmode\nobreak\
,\end{split}$ (65)

where $\pi f_{a}^{2}\log y\sim\mu_{\text{eff}}$ and $c_{1}\log y\sim\xi$ and
the function $f(y,\,u)$ is defined as

$f(y,\,u)=\left(\frac{\log(yu^{2})}{u^{2}\log
y}\right)^{\frac{q+3}{2}}u^{4}-1\leavevmode\nobreak\ .$ (66)

Since $f(y,\,1)=0$ and $\partial f/\partial u<0$ for
$\exp\frac{q+3}{2(q-1)}<y^{1/2}u=\frac{k}{\sqrt{m_{r}H}}$, the function
$f(y,\,u)$ is positive decreasing function in $u$ for $u<1$. Therefore, we see
that $\partial\rho_{a}/\partial k$ rapidly decays faster than $\sim k^{-q}$
for the momentum $k\gtrsim e^{\frac{q+3}{2(q-1)}}\sqrt{m_{r}H}$, and its
contribution to the axion abundance will be accordingly suppressed.

When $y^{-1/2}x$ is small (low momentum region), the integration is done over
$-\frac{x_{0}^{2}c_{1}y}{x^{2}}W_{-1}(-\frac{x^{2}}{x_{0}^{2}c_{1}y})=y_{1}\leq
y^{\prime}\leq y_{2}=y$ where $y_{1}$ in terms of Lambert $W$ function is the
intersection between two curves $(y^{\prime}/y)^{1/2}x$ and
$(x_{0}^{2}c_{1}\log y^{\prime})^{1/2}$. The integration in Eq. (62) can also
be done straightforwardly and it is given by

$\begin{split}\frac{\partial\rho_{a}}{\partial k}(k,\,t)&=8\pi
f_{a}^{2}c_{1}m_{r}(x_{0}^{2}c_{1})^{\frac{q-1}{2}}y^{\frac{q-3}{2}}\frac{1}{x^{q}}\left(\frac{q-1}{2}\right)^{-\frac{q+3}{2}}\\\\[5.0pt]
&\times\left[\Gamma\left(\frac{q+5}{2},\,\frac{q-1}{2}\log\left(-\frac{x_{0}^{2}c_{1}y}{x^{2}}W_{-1}\left(-\frac{x^{2}}{x_{0}^{2}c_{1}y}\right)\right)\right)-\Gamma\left(\frac{q+5}{2},\,\frac{q-1}{2}\log
y\right)\right]\leavevmode\nobreak\ ,\\\\[5.0pt] &=8\pi
f_{a}^{2}c_{1}m_{r}(x_{0}^{2}c_{1})^{\frac{q-1}{2}}y^{\frac{q-3}{2}}\frac{1}{x^{q}}\left(\frac{q-1}{2}\right)^{-\frac{q+3}{2}}\\\\[5.0pt]
&\times\left[\Gamma\left(\frac{q+5}{2},\,-\frac{q-1}{2}W_{-1}\left(-\frac{x^{2}}{x_{0}^{2}c_{1}y}\right)\right)-\Gamma\left(\frac{q+5}{2},\,\frac{q-1}{2}\log
y\right)\right]\leavevmode\nobreak\ ,\end{split}$ (67)

where $\log(-wW_{k}(-w^{-1}))=-W_{k}(-w^{-1})$ was used in the second
relation. It can be easily shown that $-W_{-1}(-w^{-1})$ becomes large at
large $y$ for the low momentum. Using $\Gamma(s,\,x)\rightarrow x^{s-1}e^{-x}$
as $x\rightarrow\infty$ as before, the expression in Eq. (67) can similarly be
approximated as, in the large $\log\frac{m_{r}}{H}$ limit,

$\begin{split}\frac{\partial\rho_{a}}{\partial
k}(k,\,t)\approx\frac{8H^{2}\mu_{\text{eff}}\xi}{k}\left[\left(\frac{-W_{-1}\left(-\frac{H}{m_{r}}\log\frac{m_{r}}{H}\left(\frac{k}{x_{0}H\sqrt{\xi}}\right)^{2}\right)}{\log\frac{m_{r}}{H}}\right)^{2}-\left(\frac{k}{x_{0}H\sqrt{\xi}}\right)^{1-q}\right]\leavevmode\nobreak\
,\end{split}$ (68)

where the identity $\exp[-W_{k}(-w^{-1})]=-wW_{k}(-w^{-1})$ was used and $\pi
f_{a}^{2}\log y\sim\mu_{\text{eff}}$ and $c_{1}\log y\sim\xi$ were used for
the neat expression. The above expression can be further simplified by
approximating the Lambert function,

$\begin{split}\frac{\partial\rho_{a}}{\partial
k}(k,\,t)\approx\frac{8H^{2}\mu_{\text{eff}}\xi}{k}\left[\left(1-2\frac{\log\frac{k}{x_{0}H\sqrt{\xi}}}{\log\frac{m_{r}}{H}}\right)^{2}-\left(\frac{k}{x_{0}H\sqrt{\xi}}\right)^{1-q}\right]\leavevmode\nobreak\
.\end{split}$ (69)

From Eq. (69), we see that axion spectrum $\partial\rho_{a}/\partial k$ scales
as $\sim k^{-1}$ in the low momentum region $x_{0}H\sqrt{\xi}\lesssim
k\lesssim z_{k}\sqrt{m_{r}H}$
161616$z_{k}=(-\frac{x_{0}^{2}c_{1}}{2}W_{-1}(-\frac{2}{x_{0}^{2}c_{1}}))^{1/4}$
can be read off from the equality
$-\frac{x_{0}^{2}c_{1}y}{x^{2}}W_{-1}(-\frac{x^{2}}{x_{0}^{2}c_{1}y})=\frac{x^{2}}{y}$,
provided that $\frac{x^{2}}{x_{0}^{2}c_{1}y}<e^{-1}$, for our relevant choice
of parameters, and it is found to be order one constant. The large and low
momentum regions are separated below and above $\sim\sqrt{x_{0}m_{r}H}$ when
the IR cutoff of $H$ is used. . Comparing with [14], the net effect by the
modification of the IR cutoff from $x_{0}H$ to $x_{0}H\sqrt{\xi}$ is
equivalent to shifting $k^{0}=x_{0}H$ (in their notation) to
$k^{0}=x_{0}H\sqrt{\xi}$ without changing the overall factor.

Another place where the modified IR cutoff can affect is the average axion
field value, $\langle a^{2}\rangle^{1/2}/f_{a}$ where $\langle
a^{2}(t)\rangle=\int dkk^{-2}(\partial\rho_{a}(k,\,t)/\partial k)$. Using our
result in Eq. (69), we can evaluate $\langle a^{2}(t)\rangle$ (or by a simple
power counting to get the leading term) with the modified cutoff
$k_{\text{IR}}=x_{0}H\sqrt{\xi}$, and it gives rise to $\langle
a^{2}(t)\rangle\approx\frac{8\mu_{\text{eff}}}{x_{0}^{2}}(1/2-(1+q)^{-1}-\log^{-1}\frac{m_{r}}{H}+\log^{-2}\frac{m_{r}}{H})$.
That is, $\langle a^{2}(t)\rangle\sim 4\mu_{\text{eff}}$ at late times for
$q>1$, as opposed to $\langle a^{2}(t)\rangle\sim 4\mu_{\text{eff}}\xi$ when
the IR cutoff is $\sim H$. The ratio $\langle a^{2}\rangle^{1/2}/f_{a}$ is
reduced by the factor of $\sqrt{\xi}$. While $\langle
a^{2}(t_{\star})\rangle\sim 4\mu_{\text{eff}\star}\gg 1$ is still expected
around the time $t_{\star}$ of the QCD crossover based on our simulation, it
may be informative to consider the opposite case with $\langle
a^{2}(t)\rangle^{1/2}/f_{a}\ll 1$ and see how the modified IR cutoff by
$\sqrt{\xi}$ can change the parametric of the axion abundance. In this
hypothetical situation, the nonlinearities due to the axion potential will be
suppressed, and the axion number density at later time $t>t_{\star}$ will be
simply given by
$n_{a}^{\text{str}}(t)=(H(t)/H_{\star})^{3/2}\,n_{a}^{\text{str}}(t_{\star})$
where $n_{a}^{\text{str}}(t_{\star})=\int dkk^{-1}(\partial\rho_{a}/\partial
k)\approx 8H_{\star}\mu_{\text{eff}\star}\sqrt{\xi_{\star}}$, as opposed to
$8H_{\star}\mu_{\text{eff}\star}\xi_{\star}$ when $k_{\text{IR}}\sim H$ was
assumed. That is, the shifted cutoff by $\sqrt{\xi}$ suppresses the abundance
by the factor of $\sqrt{\xi}$.

Now getting back to the our situation where a large $\langle
a^{2}(t_{\star})\rangle^{1/2}/f_{a}$ is expected, and the nonlinearities
arising due to the axion potential can not be neglected. The axions after time
$t_{\star}$ still evolves as free relativistic fields, and the axion energy
density at time $t>t_{\star}$ is accordingly redshifted and it is given by, in
the range of $x_{0}\sqrt{\xi_{\star}H_{\star}H}<k<z_{k}\sqrt{m_{r}H}$,

$\begin{split}\frac{\partial\rho_{a}}{\partial
k}(k,\,t)\approx\frac{8H^{2}\mu_{\text{eff}\star}\xi_{\star}}{k}\left[\left(1-2\frac{\log\frac{k}{x_{0}\sqrt{\xi_{\star}H_{\star}H}}}{\log\frac{m_{r}}{H_{\star}}}\right)^{2}-\left(\frac{k}{x_{0}\sqrt{\xi_{\star}H_{\star}H}}\right)^{1-q}\right]\leavevmode\nobreak\
.\end{split}$ (70)

The evolution continues until the transition time, denoted by $t_{\ell}$,
where the axion energy density stored in the gradient terms becomes comparable
with the axion potential and during which the axion energy density is assumed
to promptly converted into non-relativistic ones, namely

$\rho_{\text{IR}}(t_{\ell})=\int_{x_{0}\sqrt{\xi_{\star}H_{\star}H}}^{c_{m}m_{a}(t_{\ell})}dk\frac{\partial\rho_{a}}{\partial
k}(k,\,t_{\ell})=c_{V}m_{a}^{2}(t_{\ell})f_{a}^{2}\leavevmode\nobreak\ ,$ (71)

where $c_{m}$ and $c_{V}$ are order one parameters that have to be determined
by numerical simulation. Axions with $k>c_{m}m_{a}$ will decay faster than
those contributing to the dominant abundance. The relation in Eq. (71) leads
to the condition (differs by the definition of $\kappa$ compared to [14]),

$8H^{2}_{\ell}\mu_{\text{eff}\star}\xi_{\star}\left[\log\kappa\left(1-2\frac{\log\kappa}{\log\frac{m_{r}}{H_{\star}}}+\frac{4}{3}\frac{\log^{2}\kappa}{\log^{2}\frac{m_{r}}{H_{\star}}}-\frac{1-\kappa^{1-q}}{q-1}\right)\right]=c_{V}m_{a}^{2}(t_{\ell})f_{a}^{2}\leavevmode\nobreak\
,$ (72)

where $\kappa=\frac{c_{m}m_{a}(t_{\ell})}{x_{0}\sqrt{\xi_{\star}H_{\star}H}}$
and $H_{\ell}=H(t_{\ell})$. Parametrizing the axion mass as
$m_{a}(t)=H_{\star}(H_{\star}/H)^{\alpha/4}$ with $m_{a}(t_{\star})=H_{\star}$
and introducing [14],

$z\equiv\left(\frac{m_{a}(t_{\ell})}{H_{\star}}\right)^{1+\frac{6}{\alpha}}\leavevmode\nobreak\
,$ (73)

we obtain, keeping only first term as the dominant one in Eq. (72) at late
times $\log\gg 1$,

$\frac{8\mu_{\text{eff}\star}\xi_{\star}}{c_{V}f_{a}^{2}}\log\left(\frac{c_{m}}{x_{0}\sqrt{\xi_{\star}}}z^{\frac{\alpha+2}{\alpha+6}}\right)=z^{\frac{2(\alpha+4)}{\alpha+6}}\leavevmode\nobreak\
.$ (74)

Comparing with the result in [14], the modified momentum IR cutoff amounts to
shifting $x_{0}$ to $x_{0}\sqrt{\xi_{\star}}$. The axion abundance at
$t=t_{\ell}$ is estimated as
$n_{a}^{\text{str}}(t_{\ell})=c_{n}\frac{\rho_{\text{IR}}(t_{\ell})}{m_{a}(t_{\ell})}=c_{n}c_{V}m_{a}(t_{\ell})f_{a}^{2}$
whereas the contribution from the misalignment with the order one angle is
$n_{a}^{\text{mis}}(t_{\ell})=c^{\prime}_{n}m_{a}(t_{\star})f_{a}^{2}(H_{\ell}/H_{\star})^{3/2}$.
Expressing $m_{a}(t_{\ell})$ in terms of the solution $z$ using Eq. (73) and
taking ratio of two different types of contributions, we obtain

$\frac{n_{a}^{\text{str}}(t_{\ell})}{n_{a}^{\text{mis}}(t_{\ell})}=\frac{c_{n}c_{V}}{c^{\prime}_{n}}z\leavevmode\nobreak\
,$ (75)

where the solution $z$ of Eq. (87) is given by

$z=\left[-\frac{\alpha+2}{2(\alpha+4)}\frac{8\mu_{\text{eff}\star}\xi_{\star}}{c_{V}f_{a}^{2}}W_{-1}\left(-\frac{2(\alpha+4)}{\alpha+2}\frac{c_{V}f_{a}^{2}}{8\mu_{\text{eff}\star}\xi_{\star}}\left(\frac{c_{m}}{x_{0}\sqrt{\xi_{\star}}}\right)^{-\frac{2(\alpha+4)}{\alpha+2}}\right)\right]^{\frac{\alpha+6}{2(\alpha+4)}}\leavevmode\nobreak\
,$ (76)

which is the same as that in [14] except for $x_{0}\rightarrow
x_{0}\sqrt{\xi}$. Using the relation,
$-W_{-1}(-w^{-1})=\log(w\log(w\log(\cdots)))$, the above ratio of two
different contributions is finally given by (in the similar form to [14])

$\frac{n_{a}^{\text{str},q>1}(t_{\ell})}{n_{a}^{\text{mis}}(t_{\ell})}=\frac{c_{n}c_{V}}{c^{\prime}_{n}}\left[\frac{4\mu_{\text{eff}\star}\xi_{\star}}{c_{V}f_{a}^{2}}\frac{\alpha+2}{\alpha+4}\log\left(\frac{4\mu_{\text{eff}\star}\xi_{\star}}{c_{V}f_{a}^{2}}\frac{\alpha+2}{\alpha+4}\left(\frac{c_{m}}{x_{0}\sqrt{\xi_{\star}}}\right)^{\frac{2(\alpha+4)}{\alpha+2}}\log(\cdots)\right)\right]^{\frac{\alpha+6}{2(\alpha+4)}}\leavevmode\nobreak\
.$ (77)

Comparing to [14], the final result differs by the replacement of
$x_{0}\rightarrow x_{0}\sqrt{\xi_{\star}}$ without changing the overall
factors. Upon replacements of $\mu_{\text{eff}\star}\sim\pi
f_{a}^{2}\log_{\star}$, $\xi_{\star}\sim c_{1}\log_{\star}$ inside $\log$ at
late times, or $\log_{\star}\gg 1$, the shifted momentum induces the
modification of the overall coefficient of the axion abundance (and sub-
leading terms) as

$\frac{n_{a}^{\text{str},q>1}(t_{\ell})}{n_{a}^{\text{mis}}(t_{\ell})}\approx\frac{c_{n}c_{V}}{c^{\prime}_{n}}\left[\frac{4\mu_{\text{eff}\star}\xi_{\star}}{c_{V}f_{a}^{2}}\frac{\alpha}{\alpha+4}\left(\log\log\frac{m_{r}}{H_{\star}}+\mathcal{O}(\log\log\log\frac{m_{r}}{H_{\star}})\right)\right]^{\frac{1}{2}\left(1+\frac{2}{\alpha+4}\right)}\leavevmode\nobreak\
,$ (78)

where the factor $\frac{\alpha}{\alpha+4}$ will change to
$\frac{\alpha+2}{\alpha+4}$ when using $k_{\text{IR}}=x_{0}H$ as in [14]. For
instance, it can cause roughly 20% reduction of the overall rate for the
typical choice $\alpha=8$.

###  E.2 Estimation for $q=1$

We extend the previous discussion to the case with $q=1$ for the sake of
completeness (and for the clear comparison with literature). We primarily
present the estimate with IR cutoff $x_{\text{IR}}=x_{0}\sqrt{\xi}$ while
commenting on the case with $x_{\text{IR}}=x_{0}$ when it is relevant.
$x_{\text{UV}}=1$ is chosen as before. The instantaneous emission for $q=1$ is
given by

$\begin{split}F\left[x,\,y\right]\sim&\left\\{\begin{array}[]{lll}\displaystyle\frac{1}{\log{\frac{y}{x_{\text{IR}}}}}\,\frac{1}{x}&\text{for}&x_{0}\sqrt{\xi}=x_{\text{IR}}\leq
x\leq y\\\\[15.0pt] \hskip 10.11775pt0&&\text{otherwise}\leavevmode\nobreak\
.\end{array}\right.\end{split}$ (79)

Since the finite support of $F$ is the same as in Section E.1, the integration
in Eq. (62) can be done similarly for large and small momenta $k$, or
$y^{-1/2}x$ over similar integration ranges to Section E.1.

When $y^{-1/2}x$ is large, the integration gives rise to

$\begin{split}\frac{\partial\rho_{a}}{\partial k}(k,\,t)&=\frac{2\pi
f_{a}^{2}c_{1}m_{r}}{yx}\left(\log^{2}{y}-\log^{2}{\frac{x^{2}}{y}}\right)\approx\frac{8H^{2}\mu_{\text{eff}}\xi}{k}\frac{\log\frac{k}{H}}{\log\frac{m_{r}}{H}}\left(1-\frac{\log\frac{k}{H}}{\log\frac{m_{r}}{H}}\right)\leavevmode\nobreak\
,\end{split}$ (80)

where $\pi f_{a}^{2}\log{y}\sim\mu_{\text{eff}}$ and $c_{1}\log{y}\sim\xi$
were used. While $\partial\rho_{a}/\partial k$ scales as $\sim k^{-1}$ as
expected, it is multiplied by the logarithmic suppression
$\frac{\log{(k/H)}}{\log{(m_{r}/H)}}$. Similarly, when $y^{-1/2}x$ is small,
the integration is given by

$\begin{split}\frac{\partial\rho_{a}}{\partial k}(k,\,t)&=\frac{2\pi
f_{a}^{2}c_{1}m_{r}}{yx}\left[\log^{2}{y}-\log^{2}\left(-\frac{{x_{0}}^{2}c_{1}y}{x^{2}}W_{-1}\left(-\frac{x^{2}}{{x_{0}}^{2}c_{1}y}\right)\right)\right]\\\\[5.0pt]
&=\frac{2\pi
f_{a}^{2}c_{1}m_{r}}{yx}\left[\log^{2}{y}-\left(-W_{-1}\left(-\frac{x^{2}}{{x_{0}}^{2}c_{1}y}\right)\right)^{2}\right]\leavevmode\nobreak\
,\end{split}$ (81)

where the relation $\log(-wW_{k}(-w^{-1}))=-W_{k}(-w^{-1})$ is used in the
second relation. Using the approximate relations $\pi
f_{a}^{2}\log{y}\sim\mu_{\text{eff}}$ and $c_{1}\log{y}\sim\xi$ and
approximating the Lambert function, the expression in Eq. (81) is further
simplified as

$\begin{split}\frac{\partial\rho_{a}}{\partial
k}(k,\,t)\approx\frac{8H^{2}\mu_{\text{eff}}\xi}{k}\frac{\log{\frac{k}{x_{0}H\sqrt{\xi}}}}{\log{\frac{m_{r}}{H}}}\left(1-\frac{\log{\frac{k}{x_{0}H\sqrt{\xi}}}}{\log{\frac{m_{r}}{H}}}\right)\leavevmode\nobreak\
.\end{split}$ (82)

The axion spectrum $\partial\rho_{a}/\partial k$ in the large momentum also
scales as $\sim k^{-1}$, as is expected, with the logarithmic suppression.
Unlike the situation with $q>1$, the current case with $q=1$ shows the power
law behavior $\sim k^{-1}$ over the entire momentum range as well as the
logarithmic suppression. These properties remain the same even if
$x_{\text{IR}}=x_{0}$ is taken. It does not change $\partial\rho_{a}/\partial
k$ except for replacing $x_{0}H\sqrt{\xi}$ with $x_{0}H$ and, subsequently,
the splitting of the momentum range into $x_{0}H<k\leq\sqrt{x_{0}m_{r}H}$ (for
low) and $\sqrt{x_{0}m_{r}H}<k$ (for high).

Switching from $q>1$ to $q=1$ case also affects the average axion field value.
Using the result in Eq. (82) with cutoff $k_{\text{IR}}=x_{0}H\sqrt{\xi}$, the
average axion field value is given by

$\langle
a^{2}(t)\rangle\approx\frac{2\mu_{\text{eff}}}{x_{0}^{2}\log{\frac{m_{r}}{H}}}\quad\text{for}\quad
x_{\text{IR}}=x_{0}\sqrt{\xi}\leavevmode\nobreak\ ,$ (83)

that is, $\langle a^{2}(t)\rangle\sim 2\pi f_{a}^{2}$ at late times with no
$\log{\frac{m_{r}}{H}}$ enhancement unlike to $q>1$ case. However, this
behavior is IR cutoff sensitive, and taking $x_{\text{IR}}=x_{0}$ instead
gives

$\langle
a^{2}(t)\rangle\approx\frac{2\mu_{\text{eff}}\xi}{x_{0}^{2}\log{\frac{m_{r}}{H}}}\quad\text{for}\quad
x_{\text{IR}}=x_{0}\leavevmode\nobreak\ ,$ (84)

which has a $\log$ enhancement. It implies that the size of $\langle
a^{2}(t_{\star})\rangle^{1/2}/f_{a}$ is not obvious in priori when $q=1$,
explicit numerical check may be necessary as in [15]. Here we discuss the
parametric behavior of axion number density for both cases, one with $\langle
a^{2}(t_{\star})\rangle^{1/2}/f_{a}\ll 1$ and the other one with $\langle
a^{2}(t_{\star})\rangle^{1/2}/f_{a}\gg 1$. In the former case, the scaling in
$\log{\frac{m_{r}}{H_{\star}}}$ is sensitive to the IR cutoff. As was
discussed in Section E.1, the axion number density at later time $t>t_{\star}$
is given by
$n_{a}^{\text{str}}(t)=(H(t)/H_{\star})^{3/2}\,n_{a}^{\text{str}}(t_{\star})$
where $n_{a}^{\text{str}}(t_{\star})=\int dkk^{-1}(\partial\rho_{a}/\partial
k)\approx\frac{8\pi f_{a}^{2}H_{\star}}{x_{0}}\sqrt{\xi_{\star}}$ for
$x_{\text{IR}}=x_{0}\sqrt{\xi}$. The abundance is enhanced by the factor of
$\sqrt{\xi_{\star}}$, when $x_{\text{IR}}=x_{0}$, or
$n_{a}^{\text{str}}(t_{\star})\approx\frac{8\pi
f_{a}^{2}H_{\star}\xi_{\star}}{x_{0}}$. However, this may not cause a large
deviation in the abundance $n_{a}^{\text{str}}(t_{\star})$ since $\xi_{\star}$
is limited to $\xi_{\star}\lesssim\frac{x_{0}^{2}}{2\pi}$ to satisfy $\langle
a^{2}(t_{\star})\rangle^{1/2}/f_{a}\lesssim 1$.

Now we discuss the latter case, where $\langle
a^{2}(t_{\star})\rangle^{1/2}/f_{a}$ is large. Similarly to the case for
$q>1$, the axion radiation around the transition time $t_{\ell}$ is estimated
by solving the relation in Eq. (71) with the axion energy density at time
$t>t_{\star}$,

$\begin{split}\frac{\partial\rho_{a}}{\partial
k}(k,\,t)\approx\frac{8H^{2}\mu_{\text{eff}\star}\xi_{\star}}{k}\frac{\log{\frac{k}{x_{0}\sqrt{\xi_{\star}H_{\star}H}}}}{\log{\frac{m_{r}}{H_{\star}}}}\left(1-\frac{\log{\frac{k}{x_{0}\sqrt{\xi_{\star}H_{\star}H}}}}{\log{\frac{m_{r}}{H_{\star}}}}\right)\leavevmode\nobreak\
,\end{split}$ (85)

in the range of $x_{0}\sqrt{\xi_{\star}H_{\star}H}<k<z_{k}\sqrt{m_{r}H}$.
Solving the relation in Eq. (71) and expressing in terms of
$\kappa=\frac{c_{m}m_{a}(t_{\ell})}{x_{0}\sqrt{\xi_{\star}H_{\star}H}}$ and
$H_{\ell}=H(t_{\ell})$ gives rise to

$8\pi
f_{a}^{2}H_{l}^{2}\xi_{\star}(\log\kappa)^{2}\left(\frac{1}{2}-\frac{1}{3}\frac{\log{\kappa}}{\log{\frac{m_{r}}{H_{\star}}}}\right)=c_{V}m_{a}^{2}(t_{\ell})f_{a}^{2}\leavevmode\nobreak\
.$ (86)

With the same parametrization of the axion mass $m_{a}(t)$ and introducing the
$z$ parameter as in Section E.1, we obtain, keeping only the first term in the
parenthesis on the left hand side in Eq. (86),

$\left(\frac{4\pi\xi_{\star}}{c_{V}}\right)^{\frac{1}{2}}\log\left(\frac{c_{m}}{x_{0}\sqrt{\xi_{\star}}}z^{\frac{\alpha+2}{\alpha+6}}\right)=z^{\frac{\alpha+4}{\alpha+6}}\leavevmode\nobreak\
,$ (87)

where choosing $x_{\text{IR}}=x_{0}$ amounts to the replacement of
$x_{0}\sqrt{\xi_{\star}}$ with $x_{0}$. As in Section E.1, the relative axion
abundance with respect to the misalignment can be estimated using the solution
$z$ to Eq. (87) which is given by

$z=\left[-\frac{\alpha+2}{\alpha+4}\left(\frac{4\pi\xi_{\star}}{c_{V}}\right)^{\frac{1}{2}}W_{-1}\left(-\frac{\alpha+4}{\alpha+2}\left(\frac{c_{V}}{4\pi\xi_{\star}}\right)^{\frac{1}{2}}\left(\frac{c_{m}}{x_{0}\sqrt{\xi}}\right)^{-\frac{\alpha+4}{\alpha+2}}\right)\right]^{\frac{\alpha+6}{\alpha+4}}\leavevmode\nobreak\
.$ (88)

Unlike the case for $q>1$, the above solution $z$ in Eq. (88) does not hold
for an arbitrarily large $\xi_{\star}$ for given fixed value of $c_{m}$. The
value of $-W_{-1}(-w^{-1})$ exists only for $w\geq e$, and the solution in Eq.
(88) makes sense only for $\xi_{\star}$ smaller than
$\left(\frac{1}{e}\frac{\alpha+2}{\alpha+4}\sqrt{\frac{4\pi}{c_{V}}}\right)^{\alpha+2}\left(\frac{c_{m}}{x_{0}}\right)^{\alpha+4}$,
denoted by $\xi_{\star\text{max}}$. The solution to Eq. (88) does not exist
for larger values of $\xi_{\star}>\xi_{\star\text{max}}$ since the low
momentum range $x_{0}\sqrt{\xi_{\star}H_{\star}H}<k\leq c_{m}m_{a}(t)$ (that
dominantly contributes to $\rho_{IR}$) becomes too narrow so that
$\rho_{IR}(t)$ can not be large enough to satisfy Eq. (88).

A large $\xi_{\star}$ may become compatible with the solution to Eq. (88) if
$c_{m}$ is allowed to depend on $\xi_{\star}$. While the actual values of
$c_{m}$ should be determined from numerical simulations, applying the property
$-W_{-1}(-w^{-1})\geq 1$ to the solution $z$ in Eq. (88), we obtain

$z\geq\left(\frac{\alpha+2}{\alpha+4}\left(\frac{4\pi\xi_{\star}}{c_{V}}\right)^{\frac{1}{2}}\right)^{\frac{\alpha+6}{\alpha+4}}\leavevmode\nobreak\
,$ (89)

where switching inequality to equality amounts to restoring the multiplicative
$-W_{-1}(-w^{-1})$ function which depends at most logarithmically on $w$.
Finally, the axion abundance from strings in the scaling regime is estimated
as

$\frac{n_{a}^{\text{str},q=1}(t_{\ell})}{n_{a}^{\text{mis}}(t_{\ell})}\approx\frac{c_{n}c_{V}}{c^{\prime}_{n}}\left[\frac{\alpha+2}{\alpha+4}\left(\frac{4\pi\xi_{\star}}{c_{V}}\right)^{\frac{1}{2}}\right]^{1+\frac{2}{\alpha+4}}\leavevmode\nobreak\
,$ (90)

with a possible multiplicative factor scaling at most as
$\log\log\frac{m_{r}}{H_{\star}}$.

## References

  * [1] R. D. Peccei and H. R. Quinn, CP Conservation in the Presence of Instantons, Phys. Rev. Lett. 38 (1977) 1440–1443.
  * [2] S. Weinberg, A New Light Boson?, Phys. Rev. Lett. 40 (1978) 223–226.
  * [3] F. Wilczek, Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons, Phys. Rev. Lett. 40 (1978) 279–282.
  * [4] M. Dine and W. Fischler, The Not So Harmless Axion, Phys. Lett. B 120 (1983) 137–141.
  * [5] J. Preskill, M. B. Wise, and F. Wilczek, Cosmology of the Invisible Axion, Phys. Lett. B 120 (1983) 127–132.
  * [6] L. F. Abbott and P. Sikivie, A Cosmological Bound on the Invisible Axion, Phys. Lett. B 120 (1983) 133–136.
  * [7] D. J. E. Marsh, Axion Cosmology, Phys. Rept. 643 (2016) 1–79, [arXiv:1510.07633].
  * [8] T. W. B. Kibble, Evolution of a system of cosmic strings, Nucl. Phys. B 252 (1985) 227. [Erratum: Nucl.Phys.B 261, 750 (1985)].
  * [9] M. Gorghetto, E. Hardy, and G. Villadoro, Axions from Strings: the Attractive Solution, JHEP 07 (2018) 151, [arXiv:1806.04677].
  * [10] M. Kawasaki, T. Sekiguchi, M. Yamaguchi, and J. Yokoyama, Long-term dynamics of cosmological axion strings, PTEP 2018 (2018), no. 9 091E01, [arXiv:1806.05566].
  * [11] A. Vaquero, J. Redondo, and J. Stadler, Early seeds of axion miniclusters, JCAP 04 (2019) 012, [arXiv:1809.09241].
  * [12] V. B. Klaer and G. D. Moore, Global cosmic string networks as a function of tension, JCAP 06 (2020) 021, [arXiv:1912.08058].
  * [13] M. Buschmann, J. W. Foster, and B. R. Safdi, Early-Universe Simulations of the Cosmological Axion, Phys. Rev. Lett. 124 (2020), no. 16 161103, [arXiv:1906.00967].
  * [14] M. Gorghetto, E. Hardy, and G. Villadoro, More axions from strings, SciPost Phys. 10 (2021), no. 2 050, [arXiv:2007.04990].
  * [15] M. Buschmann, J. W. Foster, A. Hook, A. Peterson, D. E. Willcox, W. Zhang, and B. R. Safdi, Dark matter from axion strings with adaptive mesh refinement, Nature Commun. 13 (2022), no. 1 1049, [arXiv:2108.05368].
  * [16] C. A. J. O’Hare, G. Pierobon, J. Redondo, and Y. Y. Y. Wong, Simulations of axionlike particles in the postinflationary scenario, Phys. Rev. D 105 (2022), no. 5 055025, [arXiv:2112.05117].
  * [17] M. Hindmarsh, J. Lizarraga, A. Lopez-Eiguren, and J. Urrestilla, Scaling Density of Axion Strings, Phys. Rev. Lett. 124 (2020), no. 2 021301, [arXiv:1908.03522].
  * [18] M. Yamaguchi, M. Kawasaki, and J. Yokoyama, Evolution of axionic strings and spectrum of axions radiated from them, Phys. Rev. Lett. 82 (1999) 4578–4581, [hep-ph/9811311].
  * [19] M. Yamaguchi, Scaling property of the global string in the radiation dominated universe, Phys. Rev. D 60 (1999) 103511, [hep-ph/9907506].
  * [20] M. Yamaguchi, J. Yokoyama, and M. Kawasaki, Evolution of a global string network in a matter dominated universe, Phys. Rev. D 61 (2000) 061301, [hep-ph/9910352].
  * [21] M. Yamaguchi and J. Yokoyama, Quantitative evolution of global strings from the Lagrangian view point, Phys. Rev. D 67 (2003) 103514, [hep-ph/0210343].
  * [22] T. Hiramatsu, M. Kawasaki, T. Sekiguchi, M. Yamaguchi, and J. Yokoyama, Improved estimation of radiated axions from cosmological axionic strings, Phys. Rev. D 83 (2011) 123531, [arXiv:1012.5502].
  * [23] M. Kawasaki, K. Saikawa, and T. Sekiguchi, Axion dark matter from topological defects, Phys. Rev. D 91 (2015), no. 6 065014, [arXiv:1412.0789].
  * [24] D. P. Bennett, The evolution of cosmic strings, Phys. Rev. D 33 (1986) 872. [Erratum: Phys.Rev.D 34, 3932 (1986)].
  * [25] D. P. Bennett, Evolution of cosmic strings. 2., Phys. Rev. D 34 (1986) 3592.
  * [26] R. A. Battye and E. P. S. Shellard, Global string radiation, Nucl. Phys. B 423 (1994) 260–304, [astro-ph/9311017].
  * [27] R. A. Battye and E. P. S. Shellard, Axion string constraints, Phys. Rev. Lett. 73 (1994) 2954–2957, [astro-ph/9403018]. [Erratum: Phys.Rev.Lett. 76, 2203–2204 (1996)].
  * [28] C. J. A. P. Martins and E. P. S. Shellard, String evolution with friction, Phys. Rev. D 53 (1996) 575–579, [hep-ph/9507335].
  * [29] W. Zhang, A. Myers, K. Gott, A. Almgren, and J. Bell, AMReX: Block-Structured Adaptive Mesh Refinement for Multiphysics Applications, arXiv:2009.12009.
  * [30] W. Zhang, A. Almgren, V. Beckner, J. Bell, J. Blaschke, C. Chan, M. Day, B. Friesen, K. Gott, D. Graves, M. Katz, A. Myers, T. Nguyen, A. Nonaka, M. Rosso, S. Williams, and M. Zingale, AMReX: a framework for block-structured adaptive mesh refinement, Journal of Open Source Software 4 (May, 2019) 1370.
  * [31] K. Clough, P. Figueras, H. Finkel, M. Kunesch, E. A. Lim, and S. Tunyasuvunakool, GRChombo : Numerical Relativity with Adaptive Mesh Refinement, Class. Quant. Grav. 32 (2015), no. 24 245011, [arXiv:1503.03436].
  * [32] A. Drew and E. P. S. Shellard, Radiation from global topological strings using adaptive mesh refinement: Methodology and massless modes, Phys. Rev. D 105 (2022), no. 6 063517, [arXiv:1910.01718].
  * [33] M. J. Berger and J. Oliger, Adaptive mesh refinement for hyperbolic partial differential equations, Journal of Computational Physics 53 (1984), no. 3 484–512.
  * [34] L. Fleury and G. D. Moore, Axion dark matter: strings and their cores, JCAP 01 (2016) 004, [arXiv:1509.00026].
  * [35] M. Yamaguchi and J. Yokoyama, Lagrangian evolution of global strings, Phys. Rev. D 66 (2002) 121303, [hep-ph/0205308].
  * [36] K. Saikawa, J. Redondo, A. Vaquero, and M. Kaltschmidt, Spectrum of global string networks and the axion dark matter mass, arXiv:2401.17253.
  * [37] T. W. B. Kibble, Topology of Cosmic Domains and Strings, J. Phys. A 9 (1976) 1387–1398.
  * [38] T. W. B. Kibble, Some Implications of a Cosmological Phase Transition, Phys. Rept. 67 (1980) 183.
  * [39] A. Vilenkin, Cosmic Strings, Phys. Rev. D 24 (1981) 2082–2089.
  * [40] A. Albrecht and N. Turok, Evolution of Cosmic Strings, Phys. Rev. Lett. 54 (1985) 1868–1871.
  * [41] D. P. Bennett and F. R. Bouchet, Evidence for a Scaling Solution in Cosmic String Evolution, Phys. Rev. Lett. 60 (1988) 257.
  * [42] B. Allen and E. P. S. Shellard, Cosmic string evolution: a numerical simulation, Phys. Rev. Lett. 64 (1990) 119–122.
  * [43] W. H. Press, B. S. Ryden, and D. N. Spergel, Dynamical Evolution of Domain Walls in an Expanding Universe, Astrophys. J. 347 (1989) 590–604.
  * [44] J. N. Moore, E. P. S. Shellard, and C. J. A. P. Martins, On the evolution of Abelian-Higgs string networks, Phys. Rev. D 65 (2002) 023503, [hep-ph/0107171].
  * [45] G. Vincent, N. D. Antunes, and M. Hindmarsh, Numerical simulations of string networks in the Abelian Higgs model, Phys. Rev. Lett. 80 (1998) 2277–2280, [hep-ph/9708427].
  * [46] A. Dabholkar and J. M. Quashnock, Pinning Down the Axion, Nucl. Phys. B 333 (1990) 815–832.
  * [47] R. L. Davis, Goldstone Bosons in String Models of Galaxy Formation, Phys. Rev. D 32 (1985) 3172.
  * [48] R. L. Davis, Cosmic Axions from Cosmic Strings, Phys. Lett. B 180 (1986) 225–230.
  * [49] D. Harari and P. Sikivie, On the Evolution of Global Strings in the Early Universe, Phys. Lett. B 195 (1987) 361–365.
  * [50] C. Hagmann, S. Chang, and P. Sikivie, Axions from string decay, Nucl. Phys. B Proc. Suppl. 72 (1999) 81–86, [hep-ph/9807428].
  * [51] D. H. Lyth, Axions and inflation: Vacuum fluctuations, Phys. Rev. D 45 (May, 1992) 3394–3404.
  * [52] P. Sikivie, Axion Cosmology, Lect. Notes Phys. 741 (2008) 19–50, [astro-ph/0610440].
  * [53] T. Hiramatsu, M. Kawasaki, K. Saikawa, and T. Sekiguchi, Axion cosmology with long-lived domain walls, JCAP 01 (2013) 001, [arXiv:1207.3166].
  * [54] T. Hiramatsu, M. Kawasaki, K. Saikawa, and T. Sekiguchi, Production of dark matter axions from collapse of string-wall systems, Phys. Rev. D 85 (2012) 105020, [arXiv:1202.5851]. [Erratum: Phys.Rev.D 86, 089902 (2012)].
  * [55] P. Fox, A. Pierce, and S. D. Thomas, Probing a QCD string axion with precision cosmological measurements, hep-th/0409059.
  * [56] G. Grilli di Cortona, E. Hardy, J. Pardo Vega, and G. Villadoro, The QCD axion, precisely, JHEP 01 (2016) 034, [arXiv:1511.02867].


#  UHH-LT at SemEval-2020 Task 12: Fine-Tuning of Pre-Trained Transformer
Networks for Offensive Language Detection

Gregor Wiedemann
Seid Muhie Yimam
Language Technology Group
Department of Informatics
University of Hamburg, Germany
{gwiedemann, yimam<EMAIL_ADDRESS>
Chris Biemann

###### Abstract

Fine-tuning of pre-trained transformer networks such as BERT yield state-of-
the-art results for text classification tasks. Typically, fine-tuning is
performed on task-specific training datasets in a supervised manner. One can
also fine-tune in unsupervised manner beforehand by further pre-training the
masked language modeling (MLM) task. Hereby, in-domain data for unsupervised
MLM resembling the actual classification target dataset allows for domain
adaptation of the model. In this paper, we compare current pre-trained
transformer networks with and without MLM fine-tuning on their performance for
offensive language detection. Our MLM fine-tuned RoBERTa-based classifier
officially ranks 1st in the SemEval 2020 Shared Task 12 for the English
language. Further experiments with the ALBERT model even surpass this result.

##  1 Offensive Language Detection

The automatic detection of hate-speech, cyber-bullying, or aggressive and
offensive language became a vividly studied task in natural language
processing (NLP) in recent years. The offensive language detection (OLD)
shared Task 6 of 2019’s International Workshop on Semantic Evaluation
(SemEval) [Zampieri et al., 2019b] attracted submissions from more than 100
teams. The Offensive Language Identification Dataset (OLID) used in this
shared task comprises three hierarchical classification sub-tasks: A)
offensive language detection, B) categorization of offensive language, and C)
offensive language target identification [Zampieri et al., 2019a]. For Task A,
14,100 Twitter tweets were manually labeled as either offensive (OFF) or not
offensive (NOT). Task B distinguishes the 4,640 offensive tweets from task A
into targeted (TIN) or general, untargeted (UNT) offensive language. Task C,
finally, separates targeted insults into the three categories: groups (GRP),
individuals (IND), and others (OTH).

OffensEval 2020, the SemEval 2020 Offensive Language Detection Shared Task,
does not provide an own manually labeled training set for the English language
[Zampieri et al., 2020]. Instead, a large ‘weakly labeled’ dataset was
published by the organizers, containing roughly nine million tweets. Each
tweet has been automatically classified by an ensemble of five different
supervised classifiers trained on the OLID dataset. The weakly labeled dataset
contains the raw tweet (with user mentions replaced by a special ‘USER’ token)
along with the average label probability and the variance of the five
classifier predictions. Since there is no way that such weak labels themselves
carry more useful information to a machine learning system than the original
dataset on which the five classifiers were trained, we decided not to use any
of the weakly labeled information. Instead, for our classification systems, we
rely on the 2019 OLID dataset only. However, the OffensEval 2020 dataset is an
ample source to build models using unsupervised learning, particularly for
domain-adaptation of a pre-trained language model such as BERT [Devlin et al.,
2019] or its successors which are based on the transformer neural network
architecture. Unfortunately, training a transformer-based language model in an
unsupervised manner is incredibly resource-consuming, making it impractical to
learn from large datasets without access to larger GPU clusters or TPU
hardware. Regarding this, the contribution of our paper is two-fold:

  1. 1.

We evaluate to what extent different pre-trained transformer-based neural
network models can be fine-tuned to detect offensive language and its sub-
categories. An ensemble based on the ALBERT [Lan et al., 2019] model achieves
the best overall performance.

  2. 2.

We study how an additional fine-tuning step with masked language modeling
(MLM) of the best individual model RoBERTa [Liu et al., 2019b] conducted on
in-domain data affects the model performance. An ensemble of models trained
with this strategy was submitted as our official contribution to the
OffensEval 2020 shared task for the English language and achieved first place
in the competition.

##  2 Related Work

#### Offensive language detection:

Nowadays, a number of public datasets are available to train machine
classifiers for detecting English offensive language. Unfortunately,
underlying data sources, category definitions, data sampling strategies, and
annotation guidelines differ to a large extent between these datasets. Hence,
results of different datasets are hardly comparable, and training sets usually
cannot be combined to obtain a more robust classification system. ?), and ?)
conducted insightful surveys on this rapidly growing field. ?), ?), and ?)
recently organized shared tasks on the topic. Although winning systems can
achieve striking prediction accuracy, OLD is far from being a solved problem.
Prediction performance usually drops severely if the target data comprises
different characteristics than the training data. ?), for instance, show that
many machine learning architectures can be fooled easily just by adding the
word “love” to an offensive tweet to make it appear as non-offensive. ?)
highlight that linguistic information alone is not enough in many cases to
decide whether a tweet is hateful or not. Also context information, e.g. about
tweeting users themselves [Ribeiro et al., 2018], or mentioned users in tweets
[Wiedemann et al., 2018] can be a useful feature for automatic OLD.

#### Pre-trained language models for text classification:

Transfer learning with deep neural networks, in general, has proven to be
superior over supervised learning for text classification, especially for
small training data situations. This is illustrated exemplarily in our last
year’s approach [Wiedemann et al., 2019] to the OLD SemEval shared task which
employed unsupervised pre-training of a recurrent neural network architecture
with a topic cluster prediction task. Practically all winners of the
aforementioned shared task competitions employ some form of a fine-tuned
bidirectional transformer-based language model, a neural network architecture
for which ?) published with BERT the seminal work. This architecture has been
proven highly successful for transfer learning. A base model is pre-trained
with a MLM task and a next-sentence prediction (NSP) task in an unsupervised
manner on very large datasets. The knowledge about language regularities and
semantic coherence encoded in the network during this step can then be
employed successfully in later training steps of fine-tuning the network
weights to the actual classification task. For instance, ?) fine-tuned the
pre-trained BERT model winning the 2019 SemEval OLD shared task. Also ?), and
?) used it successfully for offensive language and hate speech detection. ?)
test a wide range of BERT fine-tuning methods for text classification and
develop best practice recommendations. Since BERT, a number of successor
models improving the network architecture, the pre-training strategy, or the
pre-training dataset have been published. A selection of these models will be
evaluated in Section 3.

##  3 Fine-tuning Transformer Networks

We investigate two questions regarding the fine-tuning of pre-trained
transformer networks for OLD. First, which pre-trained model performs best on
the 2020 OLD shared task? Second, we investigate how much language model fine-
tuning on in-domain data prior to classification fine-tuning improves the
performance of the best model.

###  3.1 Model Selection of Transformer Networks

As we have indicated in Section 2, transformer networks have been successfully
employed for several text classification tasks. We test the following
transformer-based pre-trained models for the OffensEval 2020 OLD shared task.

#### BERT – Bidirectional Encoder Representations from Transformers:

this seminal transformer- based language model employs an attention mechanism
that enables to learn contextual relations between (sub-)words in a text
sequence [Devlin et al., 2019]. BERT uses two training strategies: 1) MLM
where 15 % of the tokens in a sequence are replaced (masked) for which the
model learns to predict the original tokens, and 2) NSP where the model
receives pairs of sentences as input and learns to predict whether or not the
second sentence is a successor of the first one in their original document
context.

#### RoBERTa – A Robustly Optimized BERT Pretraining Approach:

this is a replication of BERT developed by Facebook [Liu et al., 2019b] with
the following modifications 1) training the model longer with bigger batches
as well as more and cleaner data, 2) discard the NSP objective, 3) training on
longer sequences, and 4) dynamically change the masking patterns, e.g. taking
care of masking complete multi-word units. RoBERTa outperformed BERT on most
tasks of the GLUE NLP benchmark (ibid.).

#### XLM-RoBERTa – XLM-R:

this is a cross-lingual version of RoBERTa which is trained on several
languages at once [Conneau et al., 2019]. The model itself is equivalent to
RoBERTa, but the training data consists of texts from more than 100 languages
filtered from the CommonCrawl111https://commoncrawl.org dataset.

#### ALBERT – A Lite BERT for Self-supervised Learning of Language
Representations:

this is a modification on BERT especially to mitigate memory limitations and
training time issues [Lan et al., 2019]. The main contributions that ALBERT
makes over the design choices of BERT are 1) decomposing the embedding
parameters into smaller matrices that will be projected to the hidden space
separately, 2) share parameters across layers to improve or stabilize the
learned parameters, and 3) inter-sentence coherence loss, which is based on
sentence order prediction (SOP), in contrast to BERT’s simpler NSP objective.

###  3.2 Masked Language Model Fine-tuning

?) showed that further pre-training of BERT with the masked language model
task can improve later results of supervised task-specific fine-tuning. The
authors tested within-task, in-domain and cross-domain further pre-training.
Evaluations show that the first strategy is susceptible to overfit the
training set and, thus, may harm classification performance. The last strategy
does not help since BERT is already trained on general-domain data. In-domain
further pre-training, however, helps to improve later classification
performance if there is a substantial overlap in language characteristics
between further pre-training data and supervised training data.

The ‘weakly labeled’ dataset of the 2020 OLD shared task most certainly is a
valuable dataset for further in-domain pre-training. However, with ca. 9
million tweets it is also rather large. Pre-training on the complete dataset
is not possible regarding our hardware limitations.222MLM of the RoBERTa-large
model with the full dataset on a single GPU with 12 GB RAM would take
estimated 40 days. However, due to increasing memory consumption of the Adam
optimizer during training, the process will stop unfinished way earlier due to
a memory exception. Therefore, we conduct MLM pre-training only on a small
sample of the original data. We strip URLs and user mentions from tweets,
remove duplicates and, finally, randomly sample 5 % of the original dataset
size, i.e. 436.123 tweets for further pre-training. We further pre-train the
presumably best model RoBERTa-large [Liu et al., 2019b] (cp. Section 4) for
one epoch (batch size 4, and learning rate 2e-5).

###  3.3 Ensembling

For our official OffensEval 2020 test set submission as team UHH-LT, we
aggregated predictions from classifiers with different ensemble approaches.

#### Ensemble of model variants:

We fine-tuned different transformer models with the OffensEval 2019 training
data using the corresponding test data for validation. The following models
were tested: BERT-base and BERT-large (uncased), RoBERTa-base and RoBERTa-
large, XLM-RoBERTa, and four different ALBERT models (large-v1, large-v2,
xxlarge-v1, and xxlarge-v2). Each model was fine-tuned for 6 epochs with a
learning rate of 5e-6, maximum sequence length of 128, and batch size 4. After
each epoch, the model was evaluated on the validation set. The best performing
epoch was saved for the ensembling. We tested two ensemble approaches: 1)
majority vote from all models, and 2) majority vote from one model type but
with different parameter sizes such as BERT-base and BERT-large.

#### MLM RoBERTa ensemble:

To be able to learn from the entire 2019 OLID dataset (training and test set),
as well as to smooth instabilities of predictions originating from random
effects during model training, we also aggregated predictions using 10-fold
cross-validation. For this, the further MLM pre-trained RoBERTa-large model is
fine-tuned 10 times, each time with 90 % of the OLID data for training and the
remaining 10 % as validation set. The best model after 6 epochs of training
with learning rate 5e-6 and batch size 8 is used to predict the OLD 2020 test
data. The final predictions for submission were obtained via majority vote on
the 10 predictions per test data instance.

##  4 Results

| NOT | OFF |  |
---|---|---|---|---
Model | P | R | F1 | P | R | F1 | Macro F1 | Acc.
Baselines
All NOT | 72.21 | 100.00 | 41.93 | - | 0.00 | 0.00 | 41.93 | 72.21
All OFF | - | 0.00 | 0.00 | 27.78 | 100.00 | 43.49 | 21.74 | 27.79
Single pre-trained transformer models
BERT-base | 99.06 | 90.2 | 94.42 | 79.34 | 97.78 | 87.60 | 91.01 | 92.31
BERT-large | 99.65 | 90.35 | 94.77 | 79.81 | 99.17 | 88.44 | 91.60 | 92.80
RoBERTa-base | 99.45 | 90.70 | 94.88 | 80.33 | 98.70 | 88.57 | 91.73 | 92.93
RoBERTa-large | 99.53 | 90.92 | 95.03 | 80.73 | 98.89 | 88.89 | 91.96 | 93.13
XLM-RoBERTa | 99.03 | 91.31 | 95.01 | 81.22 | 97.69 | 88.69 | 91.85 | 93.08
ALBERT-large-v1 | 98.87 | 90.24 | 94.36 | 79.32 | 97.31 | 87.40 | 90.88 | 92.20
ALBERT-large-v2 | 98.87 | 90.20 | 94.34 | 79.26 | 97.31 | 87.36 | 90.85 | 92.18
ALBERT-xxlarge-v1 | 98.35 | 91.09 | 94.58 | 80.57 | 96.02 | 87.62 | 91.10 | 92.46
ALBERT-xxlarge-v2 | 98.47 | 91.73 | 94.98 | 81.76 | 96.30 | 88.44 | 91.71 | 93.00
Ensembles of pre-trained transformer models
All | 99.65 | 90.95 | 95.10 | 80.83 | 99.17 | 89.06 | 92.08 | 93.23
BERT | 99.42 | 91.16 | 95.11 | 81.11 | 98.61 | 89.01 | 92.06 | 93.23
RoBERTa | 99.57 | 90.84 | 95.01 | 80.62 | 98.98 | 88.86 | 91.93 | 93.11
ALBERT-all | 98.23 | 92.66 | 95.36 | 83.37 | 95.65 | 89.00 | 92.23 | 93.49
ALBERT-xxlarge | 98.70 | 92.16 | 95.32 | 82.62 | 96.85 | 89.17 | 92.25 | 93.47
Table 1: Performance (in %) of baselines, single models, and ensemble models
on the OLID test set.

Table 1 shows results of binary offensive language detection for a naive
baseline (assuming all tweets as either offensive or not), as well as for the
individual fine-tuned transformer models and their corresponding ensembles.
All transformer models largely outperform the naïve baseline, some of them
(e.g. XLM-RoBERTa) even outperform most of the other system submissions in the
competition.333https://sites.google.com/site/offensevalsharedtask/results-and-
paper-submission

Our best individual model is RoBERTa-large with an F1-score of 91.96 %. Hence,
we select this model as the basis for further MLM pre-training. From Table 2,
we can see that the MLM fine-tuned RoBERTa model achieved consistently better
results in terms of Macro F1 than the single pre-trained transformer models
(to lower random effects of neural network training, the table shows average
values of 10 runs).

Regarding the ensembles of model variants, we see in Table 1 that all
approaches consistently perform better than the individual models. Here, the
ensemble averaging the predictions from the two ALBERT-xxlarge models
performed best with an F1-score of 92.25 %.

For the OffensEval 2020 Shared Task, we decided to submit the results form the
MLM pre-trained RoBERTa ensemble.444Of course, during the submission phase of
the shared task, the test set labels were not available. We, thus, based our
decision for this specific model on its performance on last year’s OLID test
set. Table 3 presents the official results of our system in the sub-tasks A,
B, and C together with their ranks achieved in the competition. While our
ensemble reaches the top rank for task A, there are a handful of competing
approaches achieving better results for offensive language categorization (B)
and target identification (C). The post-submission experiments on the official
test set as presented in this paper (cp. Table 1) show that the ALBERT-based
ensembles would have even beat this first ranked submission. However, the MLM
fine-tuned RoBERTa model was considerably more successful on task C than the
best ALBERT-xxlarge ensemble, especially to detect the “OTH” class.

| NOT | OFF |  |
---|---|---|---|---
Model | P | R | F1 | P | R | F1 | Macro F1 | Acc.
RoBERTa-large | 98.96 | 91.49 | 95.07 | 81.54 | 97.49 | 88.76 | 91.93 | 93.15
RoBERTa-large MLM-ft | 99.15 | 91.53 | 95.18 | 81.66 | 97.96 | 89.06 | 92.12 | 93.31
Table 2: Performance (in %) of MLM fine-tuned models on the OLID test set (average of 10 runs).  Team | UHH-LT
---|---
| Macro F1 (%) | Rank
Task A | 92.04 | 1 out of 84
Task B | 65.98 | 6 out of 42
Task C | 66.83 | 3 out of 38
Table 3: Official results of our OffenseEval 2020 test set submissions for
tasks A, B, and C (RoBERTa-large MLM test ensemble, cp. Table 2).

Figure 2 shows the corresponding confusion matrices for the submitted
predictions. One observation that can be revealed from the matrices is that
the predictions for Tasks B and C are considerably biased towards the majority
class. For Task A, however, we see more false positive cases for the offensive
class which is underrepresented in the training data. A qualitative look into
a sample of false predictions (cp. Fig. 1) reveals tweets wrongly predicted as
offensive (false positives), some of which seem inconsistently annotated in
the gold standard. Parts of expressions in Examples 1, and 2 qualified as
offensive in many other training examples. Examples 3, and 4 contain some
negative words that may have triggered a higher offensiveness prediction
score. For the false negative samples, it is not really obvious why the models
missed the correct class, since except for example 6, they actually contain
strong offensive vocabulary.

1. FP: _@USER the gov and nyc mayor are the biggest joke except maybe for the idiots who elected them_ 2. FP: _@USER What the fuck_ 3. FP: _@USER I know it was bad but I used to love it_ 4. FP: _men who voted yes are childish, sorry are you 17??Men, would you have a problem if a girl said if she’s not receiving head she’s not giving head?_ 5. FN: @USER It’s as if every single one of her supporters are just as stupid as her. Wehdone.. 6. FN: I’m gonna say he doesn’t really. You should check the zip code demographics of his various properties Only liars could deny that Al Sharpton hates white people 7. FN: @USER Fuck the chargers 8. FN: Last night I watched the Democrats throwing shit up against the wall, and none of it stuck. Figure 1: False positive (FP) and false negative (FN) examples of our UHH-LT Task A submission.

Figure 2: Confusion matrices from the submitted UHH-LT ensemble predictions.

##  5 Conclusion

After last year’s SemEval shared task on offensive language detection was
already dominated by the then newly published BERT model, for the year 2020
competition we were successful in fine-tuning BERT’s successor models to
create the best performing system. The predictions obtained from fine-tuning
of the ALBERT model on the OLID dataset achieved 92.25 % macro F1-score as the
best overall result (including our post-submission experiments) on the
official test set of the Shared Task A for the English language. We also found
the ‘weak labels’ distributed along with 9 million tweets by the shared task
organizers not useful for training our classifiers. However, the tweets
themselves provided useful in-domain data for unsupervised pre-training. With
92.04 % macro-F1, our predictions based on further language model pre-training
of the RoBERTa model on ca. 440.000 tweets, before fine-tuning on last year’s
OLID dataset achieved the first rank in the official SemEval competition for
sub-task A, and also high ranks (6, and 3) for the other two sub-tasks. We
conclude that further pre-training of a transformer model with in-domain data
is useful for offensive language detection. However, for tasks B and C, our
models are clearly biased towards the majority class resulting in somewhat
lower ranks. Hence, taking the high class imbalance of the OLID dataset better
into account could further improve our results.

## Acknowledgements

The paper was supported by BWFG Hamburg (Germany) within the “Forum 4.0”
project as part of the ahoi.digital funding line, and the DFG project “FAME”
(WI 4949/2-1, grant no. 406289255).

## References

  * [Basile et al., 2019]  Valerio Basile, Cristina Bosco, Elisabetta Fersini, Debora Nozza, Viviana Patti, Francisco Manuel Rangel Pardo, Paolo Rosso, and Manuela Sanguinetti.  2019\.  SemEval-2019 Task 5: Multilingual Detection of Hate Speech Against Immigrants and Women in Twitter.  In Proceedings of the 13th International Workshop on Semantic Evaluation, pages 54–63, MN, USA. ACL.
  * [Conneau et al., 2019]  Alexis Conneau, Kartikay Khandelwal, Naman Goyal, Vishrav Chaudhary, Guillaume Wenzek, Francisco Guzmán, Edouard Grave, Myle Ott, Luke Zettlemoyer, and Veselin Stoyanov.  2019\.  Unsupervised cross-lingual representation learning at scale.  CoRR, abs/1911.02116.
  * [Devlin et al., 2019]  Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  2019\.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers), pages 4171–4186, MN, USA. ACL.
  * [Fortuna and Nunes, 2018]  Paula Fortuna and Sérgio Nunes.  2018\.  A survey on automatic detection of hate speech in text.  ACM Computing Surveys (CSUR), 51(4):85.
  * [Gröndahl et al., 2018]  Tommi Gröndahl, Luca Pajola, Mika Juuti, Mauro Conti, and N. Asokan.  2018\.  All you need is “love”: Evading hate speech detection.  In Proceedings of the 11th ACM Workshop on Artificial Intelligence and Security, page 2–12, NY, USA. ACM.
  * [Lan et al., 2019]  Zhenzhong Lan, Mingda Chen, Sebastian Goodman, Kevin Gimpel, Piyush Sharma, and Radu Soricut.  2019\.  ALBERT: A lite BERT for self-supervised learning of language representations.  CoRR, abs/1909.11942.
  * [Liu et al., 2019a]  Ping Liu, Wen Li, and Liang Zou.  2019a.  NULI at SemEval-2019 task 6: Transfer learning for offensive language detection using bidirectional transformers.  In Proceedings of the 13th International Workshop on Semantic Evaluation, pages 87–91, MN, USA, June. ACL.
  * [Liu et al., 2019b]  Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov.  2019b.  RoBERTa: A robustly optimized BERT pretraining approach.  CoRR, abs/1907.11692.
  * [Mandl et al., 2019]  Thomas Mandl, Sandip Modha, Prasenjit Majumder, Daksh Patel, Mohana Dave, Chintak Mandlia, and Aditya Patel.  2019\.  Overview of the HASOC track at FIRE 2019: Hate speech and offensive content identification in Indo-European languages.  In Proceedings of the 11th Forum for Information Retrieval Evaluation, pages 14–17, Kolkata, India.
  * [Mozafari et al., 2019]  Marzieh Mozafari, Reza Farahbakhsh, and Noël Crespi.  2019\.  A BERT-based transfer learning approach for hate speech detection in online social media.  In Hocine Cherifi, Sabrina Gaito, José Fernendo Mendes, Esteban Moro, and Luis Mateus Rocha, editors, Proceedings of the 8th International Conference on Complex Networks and their Applications, pages 928–940, Lisbon, Portugal.
  * [Ribeiro et al., 2018]  Manoel Horta Ribeiro, Pedro H. Calais, Yuri A. Santos, Virgílio A. F. Almeida, and Wagner Meira Jr.  2018\.  Characterizing and detecting hateful users on Twitter.  In Proceedings of the 12th International Conference on Web and Social Media, ICWSM 2018, pages 676–679, CA, USA. AAAI Press.
  * [Risch et al., 2019]  Julian Risch, Anke Stoll, Marc Ziegele, and Ralf Krestel.  2019\.  hpiDEDIS at GermEval 2019: Offensive language identification using a German BERT model.  In Proceedings of the 15th Conference on Natural Language Processing (KONVENS 2019), pages 405–410, Erlangen, Germany.
  * [Schmidt and Wiegand, 2017]  Anna Schmidt and Michael Wiegand.  2017\.  A survey on hate speech detection using natural language processing.  In Proceedings of the 5th International Workshop on Natural Language Processing for Social Media, pages 1–10, Valencia, Spain. ACL.
  * [Struß et al., 2019]  Julia Maria Struß, Melanie Siegel, Josep Ruppenhofer, Michael Wiegand, and Manfred Klenner.  2019\.  Overview of GermEval task 2, 2019 shared task on the identification of offensive language.  In Proceedings of the 15th Conference on Natural Language Processing (KONVENS 2019), pages 354–365, Erlangen, Germany.
  * [Sun et al., 2019]  Chi Sun, Xipeng Qiu, Yige Xu, and Xuanjing Huang.  2019\.  How to fine-tune BERT for text classification?  In Maosong Sun, Xuanjing Huang, Heng Ji, Zhiyuan Liu, and Yang Liu, editors, Chinese Computational Linguistics, pages 194–206, Cham. Springer.
  * [Wiedemann et al., 2018]  Gregor Wiedemann, Eugen Ruppert, Raghav Jindal, and Chris Biemann.  2018\.  Transfer Learning from LDA to BiLSTM-CNN for Offensive Language Detection in Twitter.  In Proceedings of GermEval Task 2018, 14th Conference on Natural Language Processing (KONVENS), pages 85–94, Vienna, Austria.
  * [Wiedemann et al., 2019]  Gregor Wiedemann, Eugen Ruppert, and Chris Biemann.  2019\.  UHH-LT at SemEval-2019 task 6: Supervised vs. unsupervised transfer learning for offensive language detection.  In Proceedings of the 13th International Workshop on Semantic Evaluation, pages 782–787, MN, USA. ACL.
  * [Zampieri et al., 2019a]  Marcos Zampieri, Shervin Malmasi, Preslav Nakov, Sara Rosenthal, Noura Farra, and Ritesh Kumar.  2019a.  Predicting the type and target of offensive posts in social media.  In Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics (NAACL), pages 1415–1420, MN, USA. ACL.
  * [Zampieri et al., 2019b]  Marcos Zampieri, Shervin Malmasi, Preslav Nakov, Sara Rosenthal, Noura Farra, and Ritesh Kumar.  2019b.  SemEval-2019 Task 6: Identifying and Categorizing Offensive Language in Social Media (OffensEval).  In Proceedings of The 13th International Workshop on Semantic Evaluation (SemEval), MN, USA. ACL.
  * [Zampieri et al., 2020]  Marcos Zampieri, Preslav Nakov, Sara Rosenthal, Pepa Atanasova, Georgi Karadzhov, Hamdy Mubarak, Leon Derczynski, Zeses Pitenis, and Çağrı Çöltekin.  2020\.  SemEval-2020 Task 12: Multilingual Offensive Language Identification in Social Media (OffensEval 2020).  In Proceedings of The 14th International Workshop on Semantic Evaluation (SemEval), Barcelona, Spain. ACL.


# Emerging Dimming as Coronal Heating Episodes

Anna V. Payne NASA Graduate Fellow  Institute for Astronomy, University of
Hawai‘i at Mānoa, Honolulu, HI 96822, USA  Xudong Sun (孙旭东) Institute for
Astronomy, University of Hawai‘i at Mānoa, Pukalani, HI 96768, USA;
<EMAIL_ADDRESS>
(Received January 25, 2021; Revised March 11, 2021; Accepted March 13, 2021)

###### Abstract

Emerging dimming occurs in isolated solar active regions (ARs) during the
early stages of magnetic flux emergence. Observed by the Atmospheric Imaging
Assembly, it features a rapid decrease in extreme-ultraviolet (EUV) emission
in the 171 Å channel images, and a simultaneous increase in the 211 Å images.
Here, we analyze the coronal thermodynamic and magnetic properties to probe
its physical origin. We calculate the time-dependent differential emission
measures for a sample of 18 events between 2010 and 2012. The emission measure
(EM) decrease in the temperature range $5.7\leq\log_{10}T\leq 5.9$ is well
correlated with the EM increase in $6.2\leq\log_{10}T\leq 6.4$ over eight
orders of magnitude. This suggests that the coronal plasma is being heated
from the quiet-Sun, sub-MK temperature to 1–2 MK, more typical for ARs.
Potential field extrapolation indicates significant change in the local
magnetic connectivity: the dimming region is now linked to the newly emerged
flux via longer loops. We conclude that emerging dimming is likely caused by
coronal heating episodes, powered by reconnection between the emerging and the
ambient magnetic fields.

Solar extreme ultraviolet emission (1493); Solar coronal heating (1989); Solar
magnetic flux emergence (2000)

††facilities: SDO

##  1 Introduction

Our nearest star, the Sun, serves as a unique laboratory for stellar
processes. In particular, strong magnetic fields in active regions (ARs)
interact with the coronal plasma, changing its thermal structure as they
emerge. Studying the formation, evolution, and decay of ARs allows for a
deeper understanding of solar variability, and how it impacts us on Earth in
the form of space weather. The relation between the magnetic field and the
plasma thermodynamic properties may hold the key to solving the coronal
heating problem.

One interesting feature recently discovered is the “emerging dimming” (Zhang
et al., 2012) in isolated ARs. This type of ARs emerge with no other pre-
existing AR in the vicinity, so the local magnetic flux is likely balanced
with little open flux. As initially reported, 24 isolated ARs between 2010
June to 2011 May exhibited decreases of emission in an extreme-ultraviolet
(EUV) channel dominated by the lower-temperature Fe IX 171 Å line (0.6 MK)
during the early emerging stages. For the higher-temperature lines, in
particular Fe XIV 211 Å (2 MK), the emission increased continuously (Figure
1). No dimming was observed in other channels. The emerging dimming regions
are situated next to or around the emerging flux, with a fan or halo shape,
respectively. Zhang et al. (2012) speculated that coronal magnetic
reconnection between the emerging and background fields heats up the coronal
plasma, thus causing the cooler channel to dim and the warmer channels to
brighten.

Here we analyze the coronal thermodynamic and magnetic properties of emerging
dimming. We specifically test the hypothesis that the observed dimming is
caused by plasma heating, rather than other mechanisms such as mass loss.
Using multi-wavelength EUV data, we calculate the time-dependent differential
emission measure (DEM) and the emission measure (EM) in selective temperature
ranges. Further, we study the change of the magnetic field connectivity for a
representative case using a potential field extrapolation model. These
calculations yield insights on how the density of the optically-thin plasma at
different temperatures changes temporally and in response to the surface
magnetic field, thus allowing us to probe the underlying physical processes.

Figure 1: An example of emerging dimming around AR 11570, which is used as a
representative example throughout the paper. Left and right columns show SDO
observations taken at the start and the maximum of dimming, about 9 hr apart.
Top and middle rows show AIA 171 and 211 Å channels; bottom row shows HMI LOS
magnetograms. The red contours outline the emerging dimming region. The cyan
cross denotes the location of the pixel used in Figure 2 as an example.

##  2 Data & Methods

###  2.1 Observation

With the advent of the Solar Dynamics Observatory (SDO; Pesnell et al., 2012),
launched in 2010, the solar corona is continuously monitored over 10 UV/EUV
wavelengths with a cadence of 12 s and spatial sampling of 0$\farcs$6 with the
Atmospheric Imaging Assembly (AIA; Lemen et al., 2012). In addition, the
Helioseismic and Magnetic Imager (HMI; Schou et al., 2012), also on board SDO,
measures both line-of-sight (LOS) and vector magnetic fields in the
photosphere. The LOS and vector magnetograms are taken at a cadence of 45 s
and 12 minute respectively, with a spatial sampling of 0$\farcs$5.

With high cadence and a continuous, voluminous dataset, SDO is a rich resource
for studying how emerging dimming forms, and how they change through time. By
observing the Sun in different, optically-thin EUV channels with AIA,
information about the coronal electron density and temperature is recorded
through time. Contemporaneous EUV and magnetic field observations allows for
direct association between the coronal dynamics and their photospheric
sources.

We used the list of emerging ARs reported in Schunker et al. (2016) to search
for cases of emerging dimming. The two main properties provided by the list
are (1) their visibility in the continuum, and (2) the degree of spatial
isolation, i.e., whether the emergence occurs in relatively quiet areas of the
Sun that do not contain, and are not in the vicinity of, pre-existing ARs. The
second criterion allows us to select only isolated ARs where emerging dimming
was first identified (Zhang et al., 2012). It also avoids the possible
influence of AR-AR interaction so the interpretation of the results is more
straightforward.

###  2.2 Data Acquisition & Processing

We search for emerging dimming events over a two-year period from 2010 to 2012
based on visual inspection. Using the starting times of flux emergence given
in Schunker et al. (2016), we check the AIA 171 Å images using JHelioviewer
(Müller et al. 2017), which provides an easy way to query and visualize SDO
data. If we observe significant dimming next to or around the emerging flux,
we manually define (1) a time range (typically a few days) to encapsulate the
full history of the AR, including times before and during full emergence, and
(2) a fixed field of view (FOV) in the native Helioprojective coordinate. We
subsequently extract SDO data cubes corresponding to this time range and FOV
for each event at 30 minute cadence. The data cubes include six AIA EUV
channels (94, 131, 171, 193, 211, and 335 Å) and HMI LOS magnetograms. Our
final sample consists of 18 events (Table 1): nine were included in Zhang et
al. (2012); nine are new from 2012. We discuss the efficacy of this sample in
Section 4.

We perform all data reductions using the SolarSoftWare (SSW) IDL package. SSW
query of the SDO database returns level-1 AIA images, in which the raw images
have been de-spiked, had bad pixels removed, and flat-fielded. We use the
module aia_prep.pro to process the images into level 1.5, adjusting images to
a common plate scale and ensuring all images are centered at the same pixel.
The images are also de-rotated and co-aligned so as to enable a direct
comparison between different times.

Our analysis focuses on an “emerging dimming region” defined as a contiguous,
enclosed area whose 171 Å brightness declined post-emergence of the central
AR. The contours for each AR were selected based on the 171 Å image pixel
values during times when the AR’s surrounding area visibly darkened. This
subregion is fixed for each AR through its entire time series. This is
illustrated by the red contours in Figure 1 for a sample event in AR 11570.
The dimming regions were extracted based on a contour applied to the 171 Å
images after the time of emergence given in Schunker et al. (2016). The
contour threshold was set to be below half of the average surrounding medium
around the AR in the 171 Å data to form a contiguous entity surrounding the
AR. The same contour was then applied to the 211 Å images and magnetograms. We
consider the start of the dimming time, $t_{0}$, around the time of emergence
given in Schunker et al. (2016). The time of maximum dimming, $t_{*}$, was
then determined as the minimum of the spatially integrated EMlow curve,
described in the following section.

Table 1: Thermodynamic properties of all 18 emerging dimming events analyzed
in this work.

AR | $t_{*}$ | $t_{*}-t_{0}$ (hr) | $I(t_{*})/I_{\mathrm{QS}}$ |  $\sum\mathrm{EM}(t_{*})-\sum\mathrm{EM}(t_{0})$ (cm-3) | $\sum\mathrm{EM}(t_{*})/\sum\mathrm{EM}(t_{0})$ | Morphology
---|---|---|---|---|---|---
171 Å | 211 Å |  low |  high |  low |  high |
11122 | 2010-11-06T05:00 | 17.0 | 0.50 | 2.19 |  $-7.14\times 10^{55}$ |  $9.58\times 10^{56}$ |  0.06 |  1.66 | Fan
11179 | 2011-03-21T12:30 | 16.5 | 0.50 | 2.61 |  $-2.56\times 10^{51}$ |  $5.15\times 10^{52}$ |  0.02 |  7.55 | Fan
11194 | 2011-04-13T04:30 | 6.5 | 0.75 | 3.23 |  $-2.45\times 10^{49}$ |  $2.20\times 10^{51}$ |  0.06 |  2.21 | Fan
11198 | 2011-04-22T06:30 | 7.0 | 0.68 | 3.90 |  $-1.54\times 10^{51}$ |  $3.33\times 10^{52}$ |  0.09 |  1.34 | Halo
11211 | 2011-05-08T03:30 | 5.5 | 0.74 | 0.68 |  $-7.73\times 10^{51}$ |  $9.11\times 10^{52}$ |  0.24 |  1.20 | Fan
11214 | 2011-05-13T23:12 | 10.5 | 0.55 | 1.15 |  $-2.50\times 10^{51}$ |  $9.75\times 10^{51}$ |  0.28 |  1.01 | Fan
11215 | 2011-05-12T04:30 | 13.5 | 0.53 | 1.33 |  $-1.07\times 10^{55}$ |  $5.09\times 10^{55}$ |  0.19 |  1.12 | Fan
11220 | 2011-05-22T05:00 | 11.0 | 1.16 | 2.32 |  $-4.82\times 10^{53}$ |  $1.57\times 10^{55}$ |  0.11 |  1.87 | Fan
11221 | 2011-05-22T03:30 | 3.5 | 0.80 | 1.46 |  $-1.06\times 10^{50}$ |  $1.66\times 10^{52}$ |  0.71 |  1.08 | Fan
11400 | 2012-01-13T22:48 | 14.0 | 0.47 | 1.25 |  $-9.94\times 10^{54}$ |  $1.68\times 10^{56}$ |  0.16 |  1.22 | Fan
11414 | 2012-02-04T11:24 | 10.0 | 0.58 | 1.33 |  $-3.29\times 10^{54}$ |  $1.30\times 10^{56}$ |  0.12 |  1.47 | Fan
11431 | 2012-03-04T15:48 | 15.0 | 0.73 | 2.62 |  $-8.00\times 10^{53}$ |  $1.32\times 10^{55}$ |  0.02 |  4.43 | Fan
11437 | 2012-03-16T14:48 | 7.0 | 0.63 | 1.17 |  $-5.69\times 10^{55}$ |  $1.17\times 10^{57}$ |  0.06 |  1.65 | Halo
11446 | 2012-03-22T22:24 | 10.0 | 0.53 | 1.74 |  $-2.80\times 10^{54}$ |  $3.04\times 10^{56}$ |  0.25 |  1.25 | Fan
11570 | 2012-09-12T04:00 | 9.0 | 0.73 | 1.94 |  $-4.27\times 10^{54}$ |  $8.55\times 10^{56}$ |  0.27 |  1.38 | Fan
11603 | 2012-10-31T00:48 | 12.0 | 0.64 | 2.96 |  $-2.13\times 10^{54}$ |  $1.11\times 10^{56}$ |  0.16 |  1.35 | Fan
11607 | 2012-11-04T22:18 | 5.5 | 0.59 | 1.92 |  $-3.03\times 10^{56}$ |  $2.33\times 10^{57}$ |  0.11 |  1.23 | Fan
11624 | 2012-11-28T04:00 | 16.0 | 0.74 | 1.42 |  $-4.88\times 10^{55}$ |  $1.38\times 10^{57}$ |  0.25 |  1.36 | Fan

Notes. The first nine events were included in the original sample in Zhang et
al. (2012); the last nine are new cases. $t_{*}$ and $t_{0}$ indicate the time
of maximum dimming and the beginning of flux emergence, respectively.
$I(t_{*})/I_{\mathrm{QS}}$ is the average brightness within the dimming
contour at $t_{*}$ divided by the average brightness of the nearby quiet Sun
(QS). $\sum\mathrm{EM}(t_{*})-\sum\mathrm{EM}(t_{0})$ is the change of the EM
integrated over the emerging dimming region;
$\sum\mathrm{EM}(t_{*})/\sum\mathrm{EM}(t_{0})$ shows the relative change. The
“low” and “high” columns refer to EM calculated for the temperature range
$5.7\leq\log_{10}T\leq 5.9$ and $6.1\leq\log_{10}T\leq 6.3$, respectively.
Random errors of all variables are negligible. Morphology refers to whether
the emerging dimming region sits next to the emerging magnetic flux on one
side (fan), or fully surround the flux (halo) as classified by Zhang et al.
(2012).

Figure 2: Example DEM solution before (light blue, solid line) and during
emerging dimming (dark blue, dashed line) for AR11570. This DEM solution is
obtained from one pixel within the map shown in Figure 1. Error bars are
caused by random noise, determined by the method of Monte-Carlo resampling.
The two transparent orange regions denote the temperature bins used for the EM
calculations, where the low temperature bin ranges from $\log_{10}T=5.7-5.9$
and the high temperature bin ranges from $\log_{10}T=6.2-6.4$.

###  2.3 Differential Emission Measure (DEM) Analysis

To understand the physical origin of emerging dimming, we analyze the coronal
DEM:

$\mathrm{DEM}(T)\,dT=\int^{\infty}_{0}n_{e}^{2}(T,z)\,dz.$ (1)

Here, DEM is a function of the temperature $T$; $n_{e}(T,z)$ is the electron
number density as function of $T$ and the spatial coordinate $z$. By
convention, $z$ is $0$ at the coronal base and increases toward the observer
along the LOS. The integral of DEM over a finite temperature range $T_{0}\leq
T\leq T_{1}$ is called the emission measure (EM):

$\mathrm{EM}(T)=\int^{T_{1}}_{T_{0}}\mathrm{DEM}(T)\,dT.$ (2)

The DEM is related to the narrow-band, optically-thin EUV observations by an
integral in the temperature space:

$y_{i}=\int^{\infty}_{0}K_{i}(T)\,\mathrm{DEM}(T)\,dT,$ (3)

where $y_{i}$ is the exposure-normalized pixel value (i.e., count rate) in the
i-th EUV channel, and $K_{i}(T)$ is the known instrument- and channel-
dependent temperature response function. Given a set of AIA observations
$y_{i}$, our goal is to solve for DEM from Equation 3 so as to learn about the
thermodynamic parameters $n_{e}$ and $T$. This process is called DEM
inversion.

We employ the DEM algorithm described in Cheung et al. (2015), which computes
DEM solutions using a linear programming (also called linear optimization)
approach based on the concept of sparsity. This differs from other commonly
used procedures, which are mostly based on $\chi^{2}$-minimization. The
sparsity constraint reduces the risk of overfitting for underdetermined
systems.

After creating data cubes for each emerging dimming event in the six AIA
channels, we use the SSW module aia_sparse_em_solve.pro to calculate the DEM
solutions over time. In the end, for each event we have a cube of a DEM
solution for each AIA pixel, aligned over time to cover the whole duration.
The calculations are typically performed in the logarithm temperature space
with a $\Delta\log_{10}T=0.1$ bin size.

Our analysis is mostly focused on the total EM integrated over the emerging
dimming region, denoted as $\sum\mathrm{EM}$. We integrate for following two
temperature ranges: $\sum\mathrm{EM_{low}}$ for $5.7\leq\log_{10}T\leq 5.9$,
and $\sum\mathrm{EM_{high}}$ for $6.1\leq\log_{10}T\leq 6.3$. These ranges are
chosen to encapsulate the characteristic temperatures, around the peak of the
response function of the 171 and 211 Å channels.

Figure 3: Spatially integrated EM for AR 11570 over the emerging dimming
region over time, elapsed from 2012-09-11T19:00 UT which corresponds to the
emergence time of the AR. Top: $\mathrm{EM_{low}}$; bottom:
$\mathrm{EM_{high}}$. The transparent vertical line denotes the time of
maximum dimming. The gap in the curves is due to missing AIA data. Error bars
are too small to be visible.

The algorithm from Cheung et al. (2015) does not include a method to estimate
the random error associated with the noise in EUV images. In order to estimate
the error bars, we devise a method of resampling in a Monte-Carlo-like manner.
We utilize the SSWIDL module aia_bp_estimate_error.pro for an estimate of the
uncertainty $\sigma(y_{i})$ for each nominal AIA pixel value $y_{i}$ in the
$i$-th channel. Assuming the noise is Gaussian like, we create for each pixel
a sample following the normal distribution
$\mathcal{N}(y_{i},\sigma^{2}(y_{i}))$, and repeat for all pixels
independently. For a sample size of $N$, this effectively creates $N$
realizations of the same AIA image; this is done for each channel at each time
step. After calculating $N$ unique DEM solutions, we estimate the error as the
68% range centered at the medium of the returned DEM and EM values. We find a
sample size $N=100$ returns sufficiently stable quotes of the errors. Due to a
large number of contributing pixels, the random error in $\sum\mathrm{EM}$ is
negligible. The overall uncertainty is expected to be dominated by the poorly
understood systematics, which is on the order of $20\%$ (Judge, 2010).

###  2.4 Potential Field Extrapolation

Zhang et al. (2012) hypothesized that coronal magnetic reconnection occurs
between the emerging and background fields during emerging dimming events. If
so, magnetic connectivity should change, that is, field lines originate from
inside emerging dimming regions should have different end points after the AR
emergence. To quantitatively assess the changes, we use a local potential
field extrapolation algorithm to model the coronal magnetic field. The
algorithm is based on a Green’s function method (Sakurai, 1989) and is
implemented by Wiegelmann (2004).

We use the photospheric magnetic field in the radial direction ($B_{r}$)
derived from HMI vector magnetograms (Hoeksema et al., 2014) as the lower
boundary condition. We ignore the local curvature and perform the
extrapolation in a Cartesian coordinate with a 364 km resolution. We
subsequently trace field lines from the emerging dimming regions identified in
co-aligned AIA images. Field lines with starting points close to the side
boundary are not included as to minimize the boundary effect.

Figure 4: Spatially integrated EM for all ARs shown in Table 1. For each AR,
the top panel shows the normalized $\mathrm{EM_{low}}$ evolution and the
bottom panel shows the normalized $\mathrm{EM_{high}}$ evolution. The
transparent vertical line denotes the time of maximum dimming. The gap in the
curves is due to missing AIA data. Error bars are too small to be visible.

Figure 5: Scatter plot of the absolute change of the spatially integrated EM,
$\left|\sum\mathrm{EM}(t_{*})-\sum\mathrm{EM}(t_{0})\right|$, in the low and
high temperature bins for all events. Line shows the result of linear
regression in the logarithm space. Error bars are too small to be visible. The
best fit linear regression line, shown as the red dotted line, returns a slope
of $0.95\pm 0.12$ in log-log space corresponding to $r^{2}=0.95$. The $95\%$
confidence limits and prediction limits are shown by the transparent region
and dashed gray lines, respectively.

##  3 Results

The SDO observations shown in Figure 1 reveal several interesting features.
The emerging dimming region resides to the west of the emerging flux with a
fan-shaped boundary. In the AIA 171 Å channel, the mean intensity in the
emerging dimming region at maximum dimming $I(t_{*})$ is not only darker than
the pre-dimming counterpart, but also much darker than the surrounding quiet
Sun (QS) $I_{\mathrm{QS}}$, with $I(t_{*})/I_{\mathrm{QS}}=0.73$. In the AIA
211 Å channel, the same region becomes brighter with
$I(t_{*})/I_{\mathrm{QS}}=1.94$. Diffuse loop-like structure appear to connect
the emerging dimming region to the emerging flux.

An example of DEM solution is shown in Figure 2. For this pixel, the DEM has a
single peak at $\log_{10}T=6.2$ prior to the start of the dimming. At maximum
dimming, however, the DEM peak shifts to $\log_{10}T=6.3$, and the peak value
significantly increases. The DEM values in the range $5.7\leq\log_{10}T\leq
5.9$, i.e., sub-MK coronal plasma typical for QS, decreases drastically during
this period to near depletion. At the same time, the DEM values in the range
$6.1\leq\log_{10}T\leq 6.3$ increases by over 100$\%$. The amount of the
plasma over 1 MK has significantly increased. Such behavior is expected from
the AIA observations.

Figure 3 shows the evolution of the spatially integrated EM in AR 11570.
$\sum\mathrm{EM_{low}}$ increases first, then decreases until the maximum
dimming is reached. $\sum\mathrm{EM_{high}}$, on the other hand, continues to
increase throughout. The ratio between the total EM at the maximum dimming and
pre-dimming times $\sum\mathrm{EM}(t_{*})/\sum\mathrm{EM}(t_{0})$ are 0.27 and
1.38 for low- and high-temperature bins, respectively. The maximum dimming
occurs 9 hr after the initial flux emergence.

Figure 6: Top view of selective field lines from the potential field model for
AR 11570. The two panels correspond to the two columns in Figure 1. Field
lines in two panels have the identical starting foot points inside the dimming
contours (shaded yellow). Colors indicate their lengths, with blue (red) being
shorter (longer). The background image shows $B_{r}$ saturated at $\pm$500 G.

In our sample, the EM variation in the low-temperature channel exhibits
significant variations in time, whereas the high-temperature channel increases
more smoothly (Figure 4). We summarize in Table 1 selective variables to
quantify these evolution. Several trends are obvious.

  1. 1.

The duration of emerging dimming, defined as the time difference between the
maximum and the start of dimming $t_{*}-t_{0}$, is, on average, 10.9 hours.

  2. 2.

At the maximum dimming, the AIA 171 Å emission in the emerging dimming region
is lower than the nearby QS, whereas in 211 Å it is higher.

  3. 3.

The change of total EM, $\sum\mathrm{EM}(t_{*})-\sum\mathrm{EM}(t_{0})$, is
always negative (positive) for the low-$T$ (high-$T$) bin, as expected from
our case selection criteria. The relative change,
$1-\sum\mathrm{EM}(t_{*})/\sum\mathrm{EM}(t_{0})$, can be over $95\%$ in the
low-$T$ channel in the most extreme cases (AR 11179, AR 11431).

  4. 4.

The total EM and its change are typically 1 to 2 orders of magnitudes greater
in the high-$T$ channel. It may be partially due to the larger integration
range in the real $T$ space. A bin size of 0.2 in $\log_{10}T$ translates to
$\Delta T=0.3$ and $0.9$ MK for the low- and high-$T$ channels, respectively.

  5. 5.

Similar to Zhang et al. (2012), most events here display a fan-shaped
morphology. Only 1 of the 9 events has a halo morphology.

We find a strong correlation between the changes of total EM in the two
temperature bins, as illustrated in Figure 5. The EM values cover a wide range
of 8 orders of magnitude. A best-fit linear regression in logarithm space
reveals a slope of $0.95\pm 0.12$.

Figure 7: Top: histograms of the vertical magnetic field $B_{r}$ inside the
dimming region in AR11570, at two times shown in Figure 6. Bottom: histograms
of the length of the field lines originating from the dimming region.

Figure 6 shows selective potential field lines originated within the dimming
contour for AR 11570. Figure 7 shows the histograms of $B_{r}$ and the field
line lengths before and during emerging dimming. While the magnetic field
within the contour does not appear to change much (Figure 1 and Figure 7), the
local magnetic connectivity does change significantly over the emerging
dimming period. Many field lines that originally close locally now appear to
connect to the newly emerged flux via longer, higher-arching loops. This is
consistent with the loop-like features observed in AIA 211 Å.

##  4 Discussion & Conclusion

We have analyzed a total of 18 emerging dimming events in a quantitative
manner to understand their physical origin. Our DEM analysis shows a marked
change of the total EM in the emerging dimming regions. The simultaneous
decrease in the range $5.7\leq\log_{10}T\leq 5.9$ and increase in
$6.2\leq\log_{10}T\leq 6.4$ well explain the behaviors observed in AIA 171 and
211 Å channels.

We find two lines of evidence in support of the hypothesis that emerging
dimming is caused by coronal heating. First, the changes of $\sum\mathrm{EM}$
in the two temperature bins are well correlated. The fact that such
correlation holds over 8 orders of magnitude strongly suggest a common
physical cause. As the heating acts on plasma of all temperatures, a
proportionality between different bins are perhaps not surprising. Another
common mechanism for EUV dimming, mass loss via plasma ejection or expansion
(e.g., Harra & Sterling, 2001), will result in a decrease of EM in all
temperature bins and can thus be ruled out. Second, the change of magnetic
connectivity suggests ongoing reconnection, which provides a viable energy
source for coronal heating. Such heating episodes due to new flux emergence
have been quantitatively studied before (Tarr et al., 2014). Additionally, the
emerging dimming region magnetic field remains QS-like. We find that the net
magnetic flux is well balanced. Most loops are closed, with very few open
field lines despite the small computation domain. A mass-loss scenario is thus
unlikely.

Our sample size is relatively small, and the selection criteria are somewhat
arbitrary (for example, no event is included between June 2011 and January
2012). Nevertheless, the fact that all events studied evolve in a similar
fashion suggests that our conclusion is not biased by the sample selection. We
further note that the two parent samples our study is based on, i.e., the
emerging dimming sample (Zhang et al., 2012) and the flux emergence sample
(Schunker et al., 2016), are both complete.

We discuss two interesting aspects mentioned in Zhang et al. (2012) and
reproduced in this study. First, there is a delay of hours between the start
of flux emergence and the maximum dimming. Both the flux emergence rate and
the AR total flux were shown to negatively correlate with the delay, and
positively correlate with the emerging dimming duration. Such is consistent
with a magnetic reconnection origin. Second, most emerging dimming events have
a fan-like morphology rather than a halo. Such is likely determined by the
property of the ambient field. If it is mostly unipolar (as in a coronal
hole), the minority polarity of the emerging bipole will reconnect in all
directions. This naturally leads to a dome-shaped separatrix (e.g., Tarr et
al., 2014), whose quasi-circular footprints map to the halo dimming region.

A recent study found reduced 171 Å emission in moat-like regions around seven
isolated, well-formed ARs (Singh et al., 2021). These extended dark moats are
also visible in the 304, 131, 94, and 335 Å channels, and less so in the 193
and 211 Å channels. There are no signs of brightening. A DEM analysis
indicates a reduction of plasma in the entire $5.7\leq\log_{10}T\leq 6.2$
range. These observations point to a physical origin different from the
emerging dimming. Following Antiochos & Noci (1986), the authors propose that
the magnetic loops from the moats are pushed to low altitudes by the strong AR
fields, and are thus restricted to lower-than-coronal temperatures.

Similar dimming is reported for Sun-as-a-star synthetic observations when a
large, isolated AR is transiting the disc (Toriumi et al., 2020). The EUV
irradiance in all AIA channels are positively correlated with the sunspot
magnetic flux except for 171 Å. The AR studied is relatively stable with no
significant flux emergence. Hence, the cause may be more in line with that of
Singh et al. (2021). If it is related to emerging dimming, it is possible that
the depleted low-temperature plasma is never replenished. For more active
stars with larger stellar spots, the effect may be more significant.

We note that magnetic flux emergence simulations can now reproduce
realistically many observed structures on AR scales (e.g., Toriumi & Hotta,
2019). When extended to include a coronal domain (e.g., Cheung et al., 2019)
or coupled with a coronal model, they will provide a self-consistent “digital
laboratory” for testing the relevant physics for emerging dimming.

We conclude that emerging dimming is likely related to coronal heating
episodes powered by reconnection between the emerging and the ambient magnetic
fields. An assessment of the mass and energy budget will be deferred to future
investigation.

This work is supported by the state of Hawai‘i and NSF award #1848250. We
thank M. Cheung and W. Liu for discussion and help on the DEM error analysis.
The SDO data are courtesy NASA, the HMI, and the AIA science teams.

## References

  * Antiochos & Noci (1986) Antiochos, S. K. & Noci, G. 1986, ApJ, 301, 440
  * Cheung et al. (2015) Cheung, M. C. M., Boerner, P., Schrijver, C. J., et al. 2015, ApJ, 807, 143
  * Cheung et al. (2019) Cheung, M. C. M., Rempel, M., Chintzoglou, G., et al. 2019, NatAstron, 3, 160
  * Harra & Sterling (2001) Harra, L. K. & Sterling, A. C. 2001, ApJ, 561, L215
  * Hoeksema et al. (2014) Hoeksema, J. T., Liu, Y., Hayashi, K., et al. 2014, Sol. Phys., 289, 3483
  * Judge (2010) Judge, P. G. 2010, ApJ, 708, 1238
  * Lemen et al. (2012) Lemen, J. R., Title, A. M., Akin, D. J., et al. 2012, Sol. Phys., 275, 17
  * Müller et al. (2017) Müller, D., Nicula, B., Felix, S., et al. 2017, A&A, 606, A10
  * Pesnell et al. (2012) Pesnell, W. D., Thompson, B. J., & Chamberlin, P. C. 2012, Sol. Phys., 275, 3
  * Sakurai (1989) Sakurai, T. 1989, Space Sci. Rev., 51, 11
  * Schou et al. (2012) Schou, J., Scherrer, P. H., Bush, R. I., et al. 2012, Sol. Phys., 275, 229
  * Schunker et al. (2016) Schunker, H., Braun, D. C., Birch, A. C., Burston, R. B., & Gizon, L. 2016, A&A, 595, A107
  * Singh et al. (2021) Singh, T., Sterling, A. C., & Moore, R. L. 2021, ApJ, 909, 57
  * Tarr et al. (2014) Tarr, L. A., Longcope, D. W., McKenzie, D. E., & Yoshimura, K. 2014, Sol. Phys., 289, 3331
  * Toriumi et al. (2020) Toriumi, S., Airapetian, V. S., Hudson, H. S., et al. 2020, ApJ, 902, 36
  * Toriumi & Hotta (2019) Toriumi, S. & Hotta, H. 2019, ApJ, 886, L21
  * Wiegelmann (2004) Wiegelmann, T. 2004, Sol. Phys., 219, 87
  * Zhang et al. (2012) Zhang, J., Yang, S., Liu, Y., & Sun, X. 2012, ApJ, 760, L29