# Physical conditions that lead to the detection of the pair annihilation line
in the BOAT GRB221009A

Asaf Pe’er Bar-Ilan University, Ramat-Gan 5290002, Israel  Bing Zhang The
Nevada Center for Astrophysics, University of Nevada, Las Vegas, NV 89154, USA
Department of Physics and Astronomy, University of Nevada, Las Vegas, NV
89154, USA

###### Abstract

The brightest of all time (BOAT) GRB221009A show evidence for a narrow,
evolving MeV emission line. Here, we show that this line can be explained as
due to pair annihilation in the prompt emission region, and that its temporal
evolution is naturally explained as the high-latitude emission (emission from
higher angles from the line of sight) after prompt emission is over. We
consider both the high and low optical depth for pair production regimes, and
find acceptable solutions, with the GRB Lorentz factor $\Gamma\approx 600$ and
the emission radius $r\gtrsim 10^{16.5}$ cm. We discuss the conditions for the
appearance of such a line, and show that a unique combination of high
luminosity and Lorentz factor that is in a fairly narrow range are required
for the line detection. This explains why such an annihilation line is rarely
observed in GRBs.

Relativistic jets — High energy astrophysics — Gamma-ray bursts — Theoretical
models

##  1 Introduction

Gamma-ray bursts (GRBs) are the most extreme explosions in the universe,
releasing, in many cases more than $10^{52}-10^{53}$ ergs within a few
seconds. This energy release occurs in a very compact region, comparable to
the core of a collapsing, massive star, namely $10^{7}-10^{8}$ cm (see
Mészáros, 2006; Pe’er, 2015; Kumar & Zhang, 2015; Zhang, 2018; Mészáros, 2019;
Bošnjak et al., 2022, for a few recent reviews). Observations of very
energetic photons, in the MeV- GeV (and recently, TeV) range, thus imply that
the explosion energy must first be converted to kinetic energy in the form of
relativistic jet before dissipation causes the release of the gamma-ray
signal. Jet Lorentz factor of $\Gamma\geq 100$ is typically needed to solve
the compactness problem (Krolik & Pier, 1991; Woods & Loeb, 1995; Lithwick &
Sari, 2001); but see Dereli-Bégué et al. (2022). Explosion energy conversion
to kinetic energy can be mediated either by photons (the well known ”fireball”
model; Cavallo & Rees, 1978; Shemi & Piran, 1990; Rees & Meszaros, 1992;
Meszaros & Rees, 1993; Rees & Meszaros, 1994) or by strong magnetic fields
(Drenkhahn, 2002; Drenkhahn & Spruit, 2002; Mészáros & Rees, 1997; Zhang &
Pe’er, 2009; Zhang & Yan, 2011). In either scenario, a large number of
$e^{\pm}$-pairs is expected to be created, by interaction of energetic photons
having (comoving) energies exceeding the threshold energy of 0.511 MeV.

Although the creation of a large number of pairs was predicted a long time ago
(Pe’er & Waxman, 2004; Pe’er et al., 2006; Ioka et al., 2007), a firm
observational confirmation of the existence of an expected pair annihilation
line in the observed GRB spectra was so far lacking. This fact can partially
be explained as that the blue-shifted line is expected to be detected in the
hundred’s of MeVs (the Fermi-LAT) range, in which the detector’s sensitivity
is poorer than at lower energy bands. In addition, due to relativistic motion,
this line may be smeared if the Lorentz factor in different jet regions is
somewhat different. Nonetheless, the lack of detection of this line seem to
contradict with theoretical expectations.

This situation has changed recently with the observation of GRB221009A, which
became known as the ’brightest-of-all-times’ (BOAT) GRB. This extremely bright
burst had a reported fluence of $F\sim 0.2~{}{\rm erg/cm^{2}}$ (An et al.,
2023; Frederiks et al., 2023; Lesage et al., 2023; Burns et al., 2023). At
redshift $z=0.151$ this implies a record-breaking isotropic energy of
$E_{iso}\sim 1.5\times 10^{55}$ erg as detected by GECAM-C (An et al., 2023),
as well as a luminosity exceeding $10^{54}$ erg/s, which is also among the
highest known.

This GRB was so bright that the Fermi-GBM detector was saturated in between
219-277 s from the burst trigger. In addition to being so bright, this burst
had several other unique features. Starting at $\sim 230$ s from the trigger,
extensive TeV emission was detected, with more than 5000 photons detected
above TeV (LHAASO Collaboration et al., 2023) and the most energetic photon
having an energy as high as 13 TeV (Cao et al., 2023). The TeV emission first
rose until $\sim 242$ s, and then the flux decayed with time roughly as
$t^{-1}$. Furthermore, this burst showed a clear narrow emission line at
energy of $\sim 10$ MeV, separated, in both flux, energy as well as temporal
behaviour from the peak of the flux which was at the 100’s of keV regime. This
line was first identified immediately after the saturation end time (at $\sim
280$ s), and detected for tens of seconds (Edvige Ravasio et al., 2023; Zhang
et al., 2024b). Its energy was observed to decay in time, roughly as
$\nu_{line}(t)\propto(t-226~{}{\rm s})^{-1}$. Both Edvige Ravasio et al.
(2023) and recently Zhang et al. (2024c) have interpreted the line feature as
due to high-latitude pair annihilation line emission. However, they did not
study the detailed physical conditions to justify that a pair annihilation
line with such a luminosity can be realized in this burst.

Here, we explore the consequences of this observation. We focus on
understanding the physical conditions that enable the creation of this line,
as well as its uniqueness, namely why it has not been detected until recently.
As we show below, a relatively narrow range of parameters is required for this
line to be detected. We convincingly show that the observed signal decay is
indeed due to high latitude emission, and we explore the consequences of this
on the jet structure.

##  2 Key observations and basic model assumptions

The observed hard X-ray lightcurve of GRB22109A is very long, lasting for
several hundreds of seconds. It is variable, with variability time of a few
seconds, and show multiple peaks. The highest peak occurs at $\sim 244$ s
after the GRB trigger. In fact, this peak was so bright that the Fermi-GBM
detector was saturated in between the time interval 219-277 s (Edvige Ravasio
et al., 2023). The physical origin of this peak is uncertain. In particular,
the radius in which energy is dissipated is unknown. However, the fact that
the observed spectra is broad, and extends to high energies suggests that the
energy dissipation that produces this peak occurs above the photosphere. A
likely possibility is that a narrow Poynting flux dominated jet dissipate
energy at an emission radius beyond $10^{16}$ cm (see Zhang et al. (2024a) for
arguments).

After the saturation of the GBM detector, clear evidence for emission line
emerges. The line energy is in the range of $\sim 10$ MeV, and decays with
time very close to $h\nu_{line}(t)\propto(t-226~{}s)^{-1}$ (see Figure 1). The
decay of the flux varies considerably, but may be consistent with a $t^{-2}$
law (see Figure 2). These scaling laws are well consistent with the prediction
of high-latitude emission (e.g. Kumar & Panaitescu, 2000; Zhang et al., 2009;
Zhang, 2018).

As a figure of merit, the line luminosity is $L_{line}^{ob.}\approx 10^{50}$
erg/s. It is related to the comoving frame luminosity by
$L_{line}^{ob.}=D^{2}L^{\prime}$, where $D$ is the Doppler boost factor (see
below). Here and below, we use primes to describe quantities measured in the
comoving frame. An underlying assumption of this work is that this observed
line is due to pair annihilation. This assumption implies that the rate of
pair annihilation must be

$\dot{N}^{\prime\pm}={L^{\prime}\over 2m_{e}c^{2}}={6\times 10^{55}\over
D^{2}}~{}{\rm s^{-1}}.$ (1)

For simplicity, we assume that the energy dissipation that produces the flux
peak at $226$ s occurs at a single radius, $r$, which is the radius when the
prompt emission abruptly stops. To estimate this radius, we use the
characteristic observed lightcurve variability time scale, $\Delta t\approx
10$ s. This gives

$r\simeq\Gamma^{2}c\Delta t^{ob.}=1.08\times 10^{17}~{}\left({\Gamma\over
600}\right)^{2}\Delta t_{1}~{}{\rm cm}.$ (2)

Here and below, we estimate a fiducial value to the Lorentz factor of
$\Gamma=600$ (Zhang et al., 2024a) for the calculation. In section 4 below, we
compare the results with the data, to constrain the value of the Lorentz
factor. We further use the standard notation $Q_{x}\equiv Q/10^{x}$ in cgs
units. Note that this calculation takes into account both the radial and side
spreading of the expanding plasma shell.

As the observed spectrum is flat, the flux of energetic photons in the plasma
is high. As a result, a large number of pairs are produced. The produced pairs
rapidly cool, mainly radiatively, as well as adiabatically. The cold pairs
annihilate, producing the observed line. The observed line energy drops with
time due to the decrease of the Doppler factor of high latitude emission.

Using the observed (total) burst luminosity $\sim 10^{54}~{}{\rm
erg~{}s^{-1}}$, one can estimate the photospheric radius to be at

$r_{ph}={L\sigma_{T}\over 8\pi m_{p}c^{3}\Gamma^{3}}={5.8\times
10^{20}\over\Gamma^{3}}~{}L_{54}~{}{\rm cm}.$ (3)

Thus, for jet Lorentz factor of a few hundreds, the assumed radius of pair
emission line is much greater than the photosphere. However, this results
strongly depends on the value of the jet Lorentz factor, $\Gamma$.

Figure 1: Temporal evolution of the line energy. Left: the full data set;
Right: the zoom-in version. Figure made from data taken from the tables given
in Edvige Ravasio et al. (2023) and Zhang et al. (2024b).

Figure 2: Temporal evolution of the line luminosity. Left: the full data set;
Right: the zoom-in version. Figure made from data taken from the tables given
Edvige Ravasio et al. (2023) and Zhang et al. (2024b).

##  3 Theoretical expectation of the pair annihilation rate

The line luminosity is due to annihilation of pairs. In order to estimate the
expected luminosity, one needs to calculate the rate of pair annihilation.
Here we calculate the theoretical expectation of this rate under various
plausible plasma conditions.

###  3.1 Pair production rate

The pair annihilation rate depends on the density of pairs inside the plasma.
As pairs are both produced and annihilated, their density is set by the
balance between these two processes. Thus, we first calculate the rate of pair
production inside the plasma, from which one can deduce the pair density.

Calculation of the pair production rate is done by estimating the density of
photons inside the plasma that are energetic enough to produce pairs. The
comoving number density of photons having energies exceeding
$\epsilon^{\prime}\geq Am_{e}c^{2}$ is calculated assuming an approximately
flat spectra ($\nu F_{\nu}\propto\nu^{0}$), which gives an average photon
energy $\langle\epsilon^{\prime}\rangle\approx Am_{e}c^{2}$. Here, $A\geq 1$
is an uncertainty factor, which takes into account deviation from a flat
spectrum, variation of the cross section with energy, etc. Its value could be
much greater than unity but could be of order of unity, depending on the
physical parameters (see more discussion below). For a flat spectrum (which is
suggested by both Fermi and GECAM spectral data), the number density of
photons in the plasma is therefore

$n^{\prime}_{\gamma}(\epsilon^{\prime}=Am_{e}c^{2})={L\over 4\pi
Am_{e}c^{2}\Gamma^{2}cr^{2}}={7.7\times 10^{8}\over A}~{}L_{54}~{}\Delta
t_{1}^{-2}~{}\left({\Gamma\over 600}\right)^{-6}~{}{\rm cm^{-3}}.$ (4)

The energetic photons interact with low energy photons to produce electron-
positron pairs. For a photon with (comoving) energy $Am_{e}c^{2}$, the cross
section for pair production is maximal for producing a pair with a photon of
(comoving) energy $\epsilon^{\prime}\approx A^{-1}m_{e}c^{2}$, and is equal to
$\approx\sigma_{T}$. The typical (comoving) timescale for producing pairs is
therefore

$(t^{\prime})^{-1}=n^{\prime}_{\gamma}(\epsilon^{\prime})c\sigma_{T}={L\sigma_{T}\over
4\pi A^{-1}m_{e}c^{2}\Gamma^{2}r^{2}}=1.5\times 10^{-5}A~{}L_{54}~{}\Delta
t_{1}^{-2}~{}\left({\Gamma\over 600}\right)^{-6}~{}{\rm s^{-1}},$ (5)

or $t^{\prime}=6.5\times 10^{4}$ s. This value should be compared to the
comoving (dynamical) time, $t^{\prime}_{dyn}=r/\Gamma c=\Gamma\Delta t=6\times
10^{3}~{}(\Gamma/600)\Delta t_{1}$ s. The fact that $t^{\prime}\gg
t^{\prime}_{dyn}$ (for the fiducial parameter values chosen) implies that for
these parameter values, only a small fraction of the photons will produce
pairs within the available (dynamical) time.

The comoving width of the expanding plasma is $\Delta r^{\prime}=r/\Gamma$.
The optical depth for producing pairs is thus

$\tau_{\gamma\gamma}=\Delta
r^{\prime}n^{\prime}_{\gamma}(\epsilon^{\prime})\sigma_{T}={L\sigma_{T}\over
4\pi A^{-1}m_{e}c^{3}\Gamma^{3}r}=0.092A~{}L_{54}~{}\Delta
t_{1}^{-1}~{}\left({\Gamma\over 600}\right)^{-5}.$ (6)

Given the high uncertainty in the value of the jet Lorentz factor and the very
strong dependence of the optical depth on its unknown value, there is a large
uncertainty in the value of the optical depth, $\tau_{\gamma\gamma}$. Since a
qualitatively different result is expected for high ($\tau_{\gamma\gamma}>1$)
and low ($\tau_{\gamma\gamma}<1$) optical depths, we split the calculations
into these two regimes: low optical depth and high optical depth for pair
production.

###  3.2 Low optical depth for pair production

For optical depth $\tau_{\gamma\gamma}<1$, the intensity drops as
$I(\tau_{\gamma\gamma})\simeq I_{0}(1-\tau_{\gamma\gamma})$. This means that
only a fraction $\tau_{\gamma\gamma}$ of the photons interact to produce pairs
before escaping the plasma. In this regime, no sharp cutoff in the high energy
spectra that results from high energy photon annihilation is expected. Rather,
some gradual decay in the high energy spectra is predicted. We therefore
expect in this regime that the value of the uncertainty parameter $A$ will not
be too far from unity.

The rate of pair production by a single photon (number of pairs produced per
second) is given by $n^{\prime}_{\gamma}(\epsilon^{\prime})c\sigma_{T}$.
Therefore, the density of the pairs produced inside the plasma per unit time
is given by

$\begin{array}[]{lcl}\dot{n}^{\prime}_{e^{\pm}}&=&n^{\prime}_{\gamma}(\epsilon^{\prime})n^{\prime}_{\gamma}(\epsilon^{\prime})c\sigma_{T}\\\
&=&{L^{2}\over(4\pi
A^{-1}m_{e}c^{2}\Gamma^{2}cr^{2})^{2}}c\sigma_{T}=1.2\times
10^{4}A^{2}~{}L_{54}^{2}~{}\Delta t_{1}^{-4}~{}\left({\Gamma\over
600}\right)^{-12}~{}{\rm cm^{-3}~{}s^{-1}},\end{array}$ (7)

which is given in units of number per unit volume per unit time. Note the
extreme dependence on the values of the parameters, especially the Lorentz
factor, $\Gamma$.

Assuming that this rate is constant during the dynamical time,
$t^{\prime}_{dyn}=r/\Gamma c$, the density of pairs inside the plasma at the
end of the dynamical time is given by

$n^{\prime}_{e^{\pm}}={A^{2}L^{2}\sigma_{T}\over(4\pi)^{2}(m_{e}c^{2})^{2}\Gamma^{5}c^{2}r^{3}}.$
(8)

####  3.2.1 Density of pairs inside the plasma: low optical depth

We assume that the pairs that are produced rapidly cool, reaching typical
velocity $v<c$. In this case, the cross section for pair annihilation is
$v\sigma_{ann.}\simeq c\sigma_{T}$. The optical depth for pair annihilation is
therefore

$\tau_{e^{\pm}}=\Delta
r\sigma_{T}n^{\prime}_{e^{\pm}}={L^{2}\sigma_{T}^{2}\over(4\pi)^{2}(m_{e}c^{2})^{2}\Gamma^{6}c^{2}r^{2}}=8.5\times
10^{-3}A^{2}~{}L_{54}^{2}~{}\Delta t_{1}^{-2}~{}\left({\Gamma\over
600}\right)^{-10},$ (9)

where $\Delta r=r/\Gamma$ was used. This optical depth must be smaller than
unity, since $n^{\prime}_{e^{\pm}}<n^{\prime}_{\gamma}$ (this is indeed seen
in the equation). Therefore, the rate of annihilation is

$\dot{n}^{\prime}_{ann}=(n^{\prime}_{e^{\pm}})^{2}c\sigma_{T}.$ (10)

We may therefore write the rate equation for the time evolution of pair
density inside the plasma:

$\dot{n}^{\prime}_{e^{\pm}}=(n^{\prime}_{\gamma})^{2}c\sigma_{T}-(n^{\prime}_{e^{\pm}})^{2}c\sigma_{T}\equiv
C_{1}-C_{2}(n^{\prime}_{e^{\pm}})^{2}$ (11)

(where $C_{1}$ and $C_{2}$ are constants). The solution to this equation is

$n^{\prime}_{e^{\pm}}=n^{\prime}_{\gamma}\tanh\left[(n^{\prime}_{\gamma}c\sigma_{T})t\right].$
(12)

Here,
$n^{\prime}_{\gamma}=\sqrt{n^{\prime}_{\gamma}(\epsilon)n^{\prime}_{\gamma}(\epsilon^{\prime})}$.

This means that in a steady state,

$n^{\prime}_{e^{\pm}}\underset{t\rightarrow\infty}{\rightarrow}n^{\prime}_{\gamma}=7.7\times
10^{8}A~{}L_{54}~{}\Delta t_{1}^{-2}~{}\left({\Gamma\over
600}\right)^{-6}~{}{\rm cm^{-3}},$ (13)

and that the characteristic timescale for approaching this solution is

$t^{\prime}_{ss}={1\over n^{\prime}_{\gamma}c\sigma_{T}}={4\pi
A^{-1}m_{e}c^{2}\Gamma^{2}r^{2}\over L\sigma_{T}}=6.5\times
10^{4}A^{-1}~{}L_{54}^{-1}~{}\Delta t_{1}^{2}~{}\left({\Gamma\over
600}\right)^{6}~{}{\rm s}.$ (14)

(Note again that all times are measured in the comoving frame). This time
should be compared to the (comoving) dynamical time,
$t^{\prime}_{dyn}=r/\Gamma c=6\times 10^{3}~{}(\Gamma/600)\Delta t_{1}$ s. One
thus concludes: (i) that the time to reach steady state is about an order of
magnitude longer than the available time; and (ii) this time depends on
$\Gamma$. The condition for $t^{\prime}_{ss}>t^{\prime}_{dyn}$ is translated
into

$\Gamma>370A^{1/5}~{}L_{54}^{1/5}~{}\Delta t_{1}^{-1/5}.$ (15)

For lower values of $\Gamma$, the density of pairs inside the plasma
approaches the asymptotic value given in Equation 13 (with an inverse
dependence on $\Gamma$!). For a larger Lorentz factor, a steady state is not
reached.

We proceed by assuming that $t^{\prime}_{ss}\gg t^{\prime}_{dyn}$. In this
case, one can approximate $\tanh(x)\approx x$, which is valid for $x\ll 1$.
Equation 12 then gives the comoving pair density inside the plasma,

$\begin{array}[]{lcl}n^{\prime}_{e^{\pm}}(t^{\prime})&\approx&n^{\prime}_{\gamma}(n^{\prime}_{\gamma}c\sigma_{T})t^{\prime}=(n^{\prime}_{\gamma})^{2}c\sigma_{T}t^{\prime}\\\
&=&{L^{2}\sigma_{T}\over(4\pi)^{2}A^{-1}(m_{e}c^{2})^{2}\Gamma^{4}cr^{4}}t^{\prime}=1.18\times
10^{4}t^{\prime}~{}AL_{54}^{2}~{}\Delta t_{1}^{-4}~{}\left({\Gamma\over
600}\right)^{-12}~{}{\rm cm^{-3}}.\end{array}$ (16)

Note the dependence on the uncertain factor of $A$ is linear.

We may now calculate the rate of annihilation, namely the total number of
photons in the emission line that are emitted per unit comoving time (second)
from the entire comoving volume,

$\dot{N}^{\prime}_{ann}=\dot{n}^{\prime}_{ann}V^{\prime}=(n^{\prime}_{e^{\pm}})^{2}c\sigma_{T}V^{\prime}$
(17)

(see Equation 10). Here, $V^{\prime}=4\pi r^{2}ct_{dyn}=2.6\times
10^{49}~{}(\Gamma/600)^{5}\Delta t_{1}^{3}~{}{\rm cm^{3}}$ is the comoving
volume111We note that we assume here a constant comoving volume. The error
introduced by this assumption is of the order unity.. Putting numbers, one
finds:

$\dot{N}^{\prime}_{ann}={L^{4}\sigma_{T}^{3}\over(4\pi)^{3}A^{-2}(m_{e}c^{2})^{4}\Gamma^{9}cr^{5}}t^{\prime
2}=7.4\times 10^{43}t^{\prime 2}~{}A^{2}L_{54}^{4}~{}\Delta
t_{1}^{-5}~{}\left({\Gamma\over 600}\right)^{-19}~{}{\rm s^{-1}}.$ (18)

At $t^{\prime}=t^{\prime}_{dyn}$, this reads

$\dot{N}^{\prime}_{ann}(t^{\prime}_{dyn})={L^{4}\sigma_{T}^{3}\over(4\pi)^{3}A^{-2}(m_{e}c^{2})^{4}\Gamma^{11}c^{3}r^{3}}=2.7\times
10^{51}A^{2}~{}L_{54}^{4}~{}\Delta t_{1}^{-3}~{}\left({\Gamma\over
600}\right)^{-17}~{}{\rm s^{-1}}.$ (19)

One therefore finds a huge dependence of the pair annihilation rate on the
value of the Lorentz factor, $\Gamma$. This dependence enables a very wide
range of values to be achieved. However, this result holds only as long as the
condition in Equation 15 is met, otherwise the optical depth becomes larger
than unity, and the parametric dependence becomes different.

It is easy to show that the condition $t^{\prime}_{ss}>t^{\prime}_{dyn}$ is
equivalent to the condition $\tau_{\gamma\gamma}<1$. Thus, if the conditions
are such that $t^{\prime}_{ss}<t^{\prime}_{dyn}$, then
$\tau_{\gamma\gamma}>1$, and one needs to consider the regime of high optical
depth, discussed below.

###  3.3 High optical depth for pair production

We here explore the other regime, in which the optical depth for pair
production, $\tau_{\gamma\gamma}>1$ (see Equation 6). The condition
$\tau_{\gamma\gamma}>1$ is equivalent to the condition that the comoving time
for producing pairs (Equation 5) is smaller than the dynamical time, i.e.,

$\tau_{\gamma\gamma}>1~{}\leftrightarrow t^{\prime}<t^{\prime}_{dyn}.$

Therefore, in this regime, production of pairs occurs on a timescale shorter
than the dynamical time. The pair production rate can therefore be estimated
as the rate in which energy is produced in photons having energies
$\epsilon^{\prime}_{ph}\geq Am_{e}c^{2}$. This rate (per unit volume) is
simply

$\dot{n}^{\prime}_{e^{\pm}}\approx{L\over(4\pi r^{2}\Gamma^{2}c)}{1\over
Am_{e}c^{2}t^{\prime}_{dyn}}.$ (20)

The produced pairs rapidly cool. At low energies, the cross section for pair
annihilation is also $\sim\sigma_{T}$. Therefore, the rate of pair
annihilation is
$\dot{n}^{\prime}_{ann}\sim{n^{\prime}_{e^{\pm}}}^{2}c\sigma_{T}$. For large
optical depth $\tau_{\gamma\gamma}\gg 1$, there is a balance between pair
production and pair annihilation rates, namely
$\dot{n}^{\prime}_{e^{\pm}}=\dot{n}^{\prime}_{ann}$, giving

${n^{\prime}_{e^{\pm}}}^{2}={L\over(4\pi
r^{2}\Gamma^{2}Am_{e}c^{3}t^{\prime}_{dyn})}{1\over
c\sigma_{T}}={L\sigma_{T}t_{dyn}\over(4\pi
r^{2}\Gamma^{2}Am_{e}c^{2})}\left({1\over
c\sigma_{T}t^{\prime}_{dyn}}\right)^{2}={\tau_{\gamma\gamma}A^{2}\over(c\sigma_{T}t^{\prime}_{dyn})^{2}},$
(21)

where we made use of $r\sim\Gamma ct^{\prime}_{dyn}$.

One therefore finds that in this regime of high optical depth to pair
production, the number density of pairs inside the plasma is expected to be

$n^{\prime}_{e^{\pm}}\sim{\tau_{\gamma\gamma}^{1/2}A\over
c\sigma_{T}t^{\prime}_{dyn}}=2.5\times
10^{9}~{}A^{3/2}~{}L_{54}^{1/2}~{}\Delta t_{1}^{-3/2}~{}\left({\Gamma\over
600}\right)^{-7/2}~{}{\rm cm^{-3}}.$ (22)

This result further implies that the optical depth to pair annihilation is

$\tau_{e^{\pm}}=\Delta
rn^{\prime}_{e^{\pm}}\sigma_{T}=\tau_{\gamma\gamma}^{1/2}A,$ (23)

namely, at any given radius, $\tau_{\gamma\gamma}$ and $\tau_{e^{\pm}}$ are
either both greater or smaller than unity.

The rate of annihilation, namely the total number of photons in the emission
line that are emitted per unit comoving time from the entire comoving volume,
can now be calculated using equation (20),

$\begin{array}[]{lcl}\dot{N}^{\prime}_{ann}=\dot{n}^{\prime}_{ann}V^{\prime}=\dot{n}^{\prime}_{e^{\pm}}V^{\prime}&=&{L\over(4\pi
r^{2}\Gamma^{2}c)}{1\over Am_{e}c^{2}t^{\prime}_{dyn}}\times 4\pi
r^{2}ct^{\prime}_{dyn}={L\over\Gamma^{2}Am_{e}c^{2}}\\\ &=&3.4\times
10^{54}~{}A^{-1}~{}L_{54}~{}\left({\Gamma\over 600}\right)^{-2}~{}{\rm
s^{-1}}.\end{array}$ (24)

Note that essentially, this result means that the brightness of the line is
equal to the rate of energy loss from high energy photons by annihilation.

The rate of pair annihilation, as calculated in Equations (19) and (24) for
the low and high optical depth regimes is plotted in Figure 3, for different
luminosities, together with the inferred value from the data. In the plot,
observed variability time $\Delta t=10$ s was assumed. Values of $A=1$ (solid
lines) and $A=8$ (dashed lines) were taken for demonstration. This helps
explaining the rareness of this line; see further discussion below.

Figure 3: Rate of pair annihilation as a function of the Lorentz factor,
$\Gamma$, from Equations 19 and 24

At low $\Gamma$ the system is in the high optical depth to pair production
regime, and the rate of annihilation is proportional to $\Gamma^{-2}$, while
at higher values of $\Gamma$ (the low optical depth regime), it sharply decays
as $\Gamma^{-17}$. Solid lines are for $A=1$, while dashed lines are for $A=8$
(note that $A=10$ would result in lines overlap).

##  4 The Doppler boost: high latitude emission

The observed line show a clear temporal decay, with peak energy that is well
fitted as $\epsilon_{\rm line}=h\nu_{\rm line}\propto(t-226~{}s)^{-1}$. This
temporal behaviour is typical for emission originating at high-lattitude, and
is due to the evolving Doppler boost.

The comoving pair annihilation line energy ($\epsilon^{\prime}=0.5$ MeV) is
Doppler boosted towards the observer: $\epsilon^{ob.}=D\epsilon^{\prime}$.
Here,

$D={1\over\Gamma(1-\beta\cos\theta)},$

where $\theta$ is the angle to the line of sight. Under the simplified
hypothesis we use here, all the photons are emitted simultaneously at the same
radius, $r$, but at different angles to the line of sight. Since photons
emitted at angle $\theta$ are delayed with respect to photons emitted on the
line of sight by an observed time

$\delta t^{ob.}={r\over c}(1-\cos\theta)\approx{r\over c}{\theta^{2}\over 2},$

one finds the Doppler boost to be

$\begin{array}[]{l}\theta<{1\over\Gamma}\rightarrow D\simeq
2\Gamma\propto{\delta t^{ob.}}^{0},\\\ \theta>{1\over\Gamma}\rightarrow
D\simeq{2\over\Gamma\theta^{2}}={r\over\Gamma c\delta t^{ob.}}.\end{array}$
(25)

Note that here we use the symbol $\delta t^{ob.}$ to describe the difference
between the observed time of the first photon (emitted on the line of sight)
and a photon emitted at angle $\theta$ to the line of sight, while $\Delta t$
is used to describe the characteristic observed variability time, from which
the emission radius is estimated.

One may now compare this theoretical expectation to the data. About $\delta
t^{ob.}\simeq 60$ s after the peak of the flux, the pair annihilation line is
at $\sim 15$ MeV, namely $D\simeq 30$. Using $r=\Gamma^{2}c\Delta t$, one
finds $r/(\Gamma c)=\Gamma\Delta t=1800$, or $(\Gamma/600)\Delta t_{1}=0.3$.
This result can be combined with the result in Equation (25) (bottom) inside
Equation (1), using observed time $\delta t^{ob.}=60$ s, to write the
observational constraint in Equation (1) as a constraint on the rate of pair
annihilation at the end of the dynamical time,

$\dot{N}^{\pm}={L^{\prime}\over 2m_{e}c^{2}}={L^{ob.}\over
2m_{e}c^{2}}\left({\delta t^{ob.}\over\Gamma\Delta t}\right)^{2}=6.8\times
10^{52}~{}s^{-1}.$ (26)

where we took $L^{ob.}=10^{50}$ erg/s as the luminosity of the line at $\delta
t^{ob.}=60$ from the peak time.

###  4.1 Low optical depth for pair production

The pair annihilation rate calculated in Equation (26) is to be compared to
the theoretical expectation given in Equation (19). Basic algebra reveals that
the two results match for $\Gamma=620~{}L_{54}^{2/7}$ and $\Delta t=3.0$ s.
Assuming these value of the Lorentz factor, the dissipation radius is
$r=3.5\times 10^{16}$ cm, which is much higher than the photospheric radius,
at $r_{ph}=2.4\times 10^{12}$ cm. The optical depth for pair production is
$\tau_{\gamma\gamma}=0.26~{}A$, i.e., close to unity, and for pair
annihilation it is $\tau_{e^{\pm}}=0.068$. For these parameters, The time to
reach steady state, $t^{\prime}_{ss}=7.1\times 10^{3}$ s is indeed longer than
the dynamical time, $t^{\prime}_{dyn}=1.86\times 10^{3}$ s (although by a
factor of a few only), implying that the approximation we used above is valid.

###  4.2 High optical depth for pair production

An alternative scenario may be that the emission occurs in region of high
optical depth to pair production. In this scenario, the constraint set by the
data (Equation 26) is to be compared to the theoretical expectation as given
in Equation (24).

The two results match, provided that $(\Gamma/600)=7~{}L_{54}^{1/2}A^{-1/2}$.
However, too high a value of $\Gamma$ is unacceptable, as in this case the
optical depth for pair production (Equation 6) will be smaller than unity.
Using the constraint $(\Gamma/600)\Delta t_{1}=0.3$ inside Equation (6) leads
to the constraint

$0.3AL_{54}>\left({\Gamma\over 600}\right)^{4},$

which, combined with the constraint on the line luminosity above, results in
the condition

$A>20~{}L_{54}.$ (27)

This is translated to an expected spectral cutoff at energy

$\epsilon^{ob.}>20Dm_{e}c^{2}=20\left({1800~{}s\over
t^{ob.}}\right)m_{e}c^{2}=10\left({1800~{}s\over t^{ob.}}\right)~{}{\rm MeV}.$
(28)

this regime is not excluded by existing data, and is therefore valid as well.

##  5 Conclusions and Discussion

###  5.1 Pair annihilation as a source of the observed line

In this paper, we considered the narrow emission line seen in GRB221009A. We
argue that this line originates from $e^{\pm}$ pair annihilation. We
calculated the expected rate of pair annihilation, as a function of the
observed GRB parameters, namely the luminosity, a spectral cutoff ($A$) and
variability time, and the unknown jet Lorentz factor. We calculated this rate
in the two regimes: low and high optical depths for pair production. The
results are presented in Equations (19) and (24), respectively.

When comparing the theoretical expectations to the observed line luminosity,
we found that both scenarios are plausible, and both lead to high Lorentz
factor, of $\sim 600$. The reason why such two different scenarios can provide
acceptable explanations is the uncertainty in the cutoff energy,
$Am_{e}c^{2}$. Therefore, observational discrimination between the two
scenarios is possible, by detecting a spectral cutoff ($A$). For the available
data of GRB221009A, this cutoff is expected in the sub-GeV range (Equation
27). This calculation can be trivially generalized to the data of any other
source.

###  5.2 High latitude emission

We explain the observed line energy, as well as the temporal decay of this
line as due to emission from high-latitude after the line-of-sight prompt
emission is over, i.e., at high angles to the line of sight. This
interpretation naturally explains why this line is detected only in the 10 MeV
range, rather than the naive expectation of $\Gamma m_{e}c^{2}\sim 300$ MeV
range, as would be expected had the emission originated from the line of
sight. Indeed, according to our model, at very early times, comparable to
$\Delta t\approx 3$ s after the peak (at 226 s), the Doppler boost is expected
to be $2\Gamma=1200$, and the emission line is expected at the sub-GeV range.
However, during this time, the detectors are saturated.

Furthermore, this interpretation provides a natural explanation to the
temporal decay of the line frequency as $(t-226~{}s)^{-1}$, and its luminosity
decay as $(t-226~{}s)^{-2}$, although this is less pronounced by the data.
This high latitude interpretation further indicates that the dissipation
episode that leads to the emission of this line is instantaneous – it occurs
over a short observed period of a few seconds, rather than being continuous.
Furthermore, the data clearly indicates that the energy dissipation episode
occurs far above the photosphere, and therefore is not directly related to the
on-going events that take place in the central engine during the prompt phase.

Indeed, the high latitude emission model has a firm prediction. Assuming a
$\delta$-function emission in time and radius (but that extends over a wide
angular range), the observed signal should be constant (time-independent) for
a duration of $\Delta t^{ob.}\approx R/(2\Gamma^{2}c)$, and after that the
flux should drop as $\Delta t^{-2}$. This is a universal prediction, which is
independent on the origin of the emission, provided it occurs at a constant
radius.

This prediction aligns very well with the observed properties of the narrow
emission line in GRB221009A. However, this only gives an upper limit on the
decay rate of a given signal. Indeed, other parts of the spectra do not follow
this behaviour. For example, the peak of the flux seen by Fermi and GECAM is
much more variable, and show high peaks at later times, up to 600 s. This
indicates that multiple energy dissipation episodes occurred, though none of
them provides sufficient energy to produce an observed pair signal.

The TeV flux shows a very different temporal behaviour. While it starts its
rise at the same time, 226 s from the trigger, it rises for $\approx 20$ s,
before decaying as $F_{\rm TeV}\propto t^{-1.1}$, which is much shallower than
the limit set by the high latitude emission. The shape of lightcurve is
consistent with synchrotron self-Compton scattering from the external shock
(LHAASO Collaboration et al., 2023). Since the TeV photons are observed, they
are not attenuated to produce pairs, likely because the emission radius is
larger than the prompt emission radius where $\tau_{\gamma\gamma}$ is much
smaller.

###  5.3 Rarity of this line

The analysis carried here can be used to explain the lack of detection of pair
annihilation line in GRBs so far, despite the clear theoretical prediction.
One characteristic of this particular GRB (221009A) is its extreme brightness,
$L\simeq 10^{54}$ erg/s (Burns et al., 2023). This is about two - three orders
of magnitude brighter than the average GRB.

Assume that for less bright GRBs the detectors would not saturate, continuous
observations would occur during the dynamical time. Thus, one could
approximate the Doppler boost as $D\simeq 2\Gamma$. The observed line
luminosity is expected to be

$L_{line}^{ob.}=D^{2}L^{\prime}=4\Gamma^{2}\times 2m_{e}c^{2}\dot{N}^{\pm}.$

The results of Equations (19) and (24) provide the theoretical expectation,

$L_{line}^{ob.}=\left\\{\begin{array}[]{ll}6.2\times
10^{51}~{}~{}L_{54}^{4}~{}\Delta t_{1}^{-3}~{}\left({\Gamma\over
600}\right)^{-15}~{}~{}{\rm erg/s}&({\rm low~{}optical~{}depth})\\\ 7.8\times
10^{54}~{}A^{-1}~{}L_{54}~{}~{}{\rm erg/s}&({\rm
high~{}optical~{}depth}).\end{array}\right.$ (29)

For typical GRB with luminosity $L\sim 10^{51}$ erg/s, the transition between
high and low optical depth occurs at $\Gamma=90~{}L_{51}^{1/5}~{}\Delta
t_{1}^{-1/5}$ (see Equation 15). Thus, for typical values of $\Gamma\geq 100$,
the valid regime is always the low optical depth regime, where the line
luminosity depends on the GRB luminosity to the power 4. In this regime, the
theoretical calculations show a huge dependence on the value of $\Gamma$,
$L_{line}^{ob.}\propto\Gamma^{-15}$. This dependence implies that $\Gamma$
must have a narrow range of values in order for the line to be detected. E.g.,
for $L=10^{51}$, $\Gamma>150$ will produce a too dim emission line.

For brighter GRBs, this constraint is somewhat releases. For example, for
$L=10^{52}$ erg/s, annihilation line from GRBs with Lorentz factor of 200-250
should be within the detected range, at least for a typical time comparable to
the variability time, $\Delta t$. In such cases, this line is expected in the
same energy range as the line seen in GRB221009A. However, for such a moderate
luminosity GRB located at a typical cosmological distance, the line would have
too few photons to be detected.

###  5.4 The dissipation radius

One of the most interesting results of our analysis is that the radius in
which energy is dissipated to produce the observed line is large, $r_{dis}\sim
10^{16.5}$ cm. This radius is far above the photospheric radius. On the other
hand, this radius cannot be associated with the onset of the afterglow phase,
as the data indicates that the prompt emission phase lasts much longer.

The origin of the dissipation process at this radius is therefore not entirely
clear. Due to the very large amount of energy released in this radius, We find
internal shocks to be a less appealing scenario, as only the differential
energy can be released, and as a result the efficiency of energy conversion is
very low, typically a few %. The two alternative options are therefore either
dissipation of a large amount of magnetic energy (e.g., Zhang & Yan, 2011;
Zhang et al., 2024a), or alternatively, an abrupt change in the external
conditions. As was shown recently (Pe’er & Ryde, 2024), this can naturally
occur by interaction of the jet with the contact discontinuity that is
naturally expected close to the edge of the wind bubble produced by the
progenitor stellar wind.

## Acknowledgements

AP acknowledge the support from the European Union (EU) via ERC consolidator
grant 773062 (O.M.J.).

## References

  * An et al. (2023) An, Z.-H., Antier, S., Bi, X.-Z., et al. 2023, arXiv e-prints, arXiv:2303.01203, doi: 10.48550/arXiv.2303.01203
  * Bošnjak et al. (2022) Bošnjak, Ž., Barniol Duran, R., & Pe’er, A. 2022, Galaxies, 10, 38, doi: 10.3390/galaxies10020038
  * Burns et al. (2023) Burns, E., Svinkin, D., Fenimore, E., et al. 2023, ApJ, 946, L31, doi: 10.3847/2041-8213/acc39c
  * Cao et al. (2023) Cao, Z., Aharonian, F., An, Q., et al. 2023, Science Advances, 9, eadj2778, doi: 10.1126/sciadv.adj2778
  * Cavallo & Rees (1978) Cavallo, G., & Rees, M. J. 1978, MNRAS, 183, 359
  * Dereli-Bégué et al. (2022) Dereli-Bégué, H., Pe’er, A., Ryde, F., et al. 2022, Nature Communications, 13, 5611, doi: 10.1038/s41467-022-32881-1
  * Drenkhahn (2002) Drenkhahn, G. 2002, A&A, 387, 714, doi: 10.1051/0004-6361:20020390
  * Drenkhahn & Spruit (2002) Drenkhahn, G., & Spruit, H. C. 2002, A&A, 391, 1141, doi: 10.1051/0004-6361:20020839
  * Edvige Ravasio et al. (2023) Edvige Ravasio, M., Sharan Salafia, O., Oganesyan, G., et al. 2023, arXiv e-prints, arXiv:2303.16223, doi: 10.48550/arXiv.2303.16223
  * Frederiks et al. (2023) Frederiks, D., Svinkin, D., Lysenko, A. L., et al. 2023, ApJ, 949, L7, doi: 10.3847/2041-8213/acd1eb
  * Ioka et al. (2007) Ioka, K., Murase, K., Toma, K., Nagataki, S., & Nakamura, T. 2007, ApJ, 670, L77, doi: 10.1086/524405
  * Krolik & Pier (1991) Krolik, J. H., & Pier, E. A. 1991, ApJ, 373, 277, doi: 10.1086/170048
  * Kumar & Panaitescu (2000) Kumar, P., & Panaitescu, A. 2000, ApJ, 541, L51, doi: 10.1086/312905
  * Kumar & Zhang (2015) Kumar, P., & Zhang, B. 2015, Phys. Rep., 561, 1, doi: 10.1016/j.physrep.2014.09.008
  * Lesage et al. (2023) Lesage, S., Veres, P., Briggs, M. S., et al. 2023, ApJ, 952, L42, doi: 10.3847/2041-8213/ace5b4
  * LHAASO Collaboration et al. (2023) LHAASO Collaboration, Cao, Z., Aharonian, F., et al. 2023, Science, 380, 1390, doi: 10.1126/science.adg9328
  * Lithwick & Sari (2001) Lithwick, Y., & Sari, R. 2001, ApJ, 555, 540, doi: 10.1086/321455
  * Mészáros (2006) Mészáros, P. 2006, Reports on Progress in Physics, 69, 2259, doi: 10.1088/0034-4885/69/8/R01
  * Mészáros (2019) —. 2019, Mem. Soc. Astron. Italiana, 90, 57.  https://arxiv.org/abs/1904.10488
  * Meszaros & Rees (1993) Meszaros, P., & Rees, M. J. 1993, ApJ, 405, 278, doi: 10.1086/172360
  * Mészáros & Rees (1997) Mészáros, P., & Rees, M. J. 1997, ApJ, 482, L29, doi: 10.1086/310692
  * Pe’er (2015) Pe’er, A. 2015, Advances in Astronomy, 2015, 907321, doi: 10.1155/2015/907321
  * Pe’er et al. (2006) Pe’er, A., Mészáros, P., & Rees, M. J. 2006, ApJ, 642, 995, doi: 10.1086/501424
  * Pe’er & Ryde (2024) Pe’er, A., & Ryde, F. 2024, arXiv e-prints, arXiv:2406.03841, doi: 10.48550/arXiv.2406.03841
  * Pe’er & Waxman (2004) Pe’er, A., & Waxman, E. 2004, ApJ, 613, 448, doi: 10.1086/422989
  * Rees & Meszaros (1992) Rees, M. J., & Meszaros, P. 1992, MNRAS, 258, 41P
  * Rees & Meszaros (1994) —. 1994, ApJ, 430, L93, doi: 10.1086/187446
  * Shemi & Piran (1990) Shemi, A., & Piran, T. 1990, ApJ, 365, L55, doi: 10.1086/185887
  * Woods & Loeb (1995) Woods, E., & Loeb, A. 1995, ApJ, 453, 583, doi: 10.1086/176421
  * Zhang (2018) Zhang, B. 2018, The Physics of Gamma-Ray Bursts, doi: 10.1017/9781139226530
  * Zhang & Pe’er (2009) Zhang, B., & Pe’er, A. 2009, ApJ, 700, L65, doi: 10.1088/0004-637X/700/2/L65
  * Zhang et al. (2024a) Zhang, B., Wang, X.-Y., & Zheng, J.-H. 2024a, Journal of High Energy Astrophysics, 41, 42, doi: 10.1016/j.jheap.2024.01.002
  * Zhang & Yan (2011) Zhang, B., & Yan, H. 2011, ApJ, 726, 90, doi: 10.1088/0004-637X/726/2/90
  * Zhang et al. (2009) Zhang, B.-B., Zhang, B., Liang, E.-W., & Wang, X.-Y. 2009, ApJ, 690, L10, doi: 10.1088/0004-637X/690/1/L10
  * Zhang et al. (2024b) Zhang, Y.-Q., Xiong, S.-L., Mao, J.-R., et al. 2024b, Science China Physics, Mechanics, and Astronomy, 67, 289511, doi: 10.1007/s11433-023-2381-0
  * Zhang et al. (2024c) Zhang, Z., Lin, H., Li, Z., et al. 2024c, arXiv e-prints, arXiv:2405.12977, doi: 10.48550/arXiv.2405.12977


# Cross-Domain Semi-Supervised Audio Event Classification using Contrastive
Regularization

###### Abstract

In this study, we proposed a novel semi-supervised training method that uses
unlabeled data with a class distribution that is completely different from the
target data or data without a target label. To this end, we introduce a
contrastive regularization that is designed to be target task-oriented and
trained simultaneously. In addition, we propose an audio mixing based simple
augmentation strategy that performed in batch samples. Experimental results
validate that the proposed method successfully contributed to the performance
improvement, and particularly showed that it has advantages in stable training
and generalization.

Index Terms—  Audio event classification, semi-supervised learning,
contrastive learning

##  1 Introduction

Sound contains a lot of information, but unfortunately, most of the audio-
based studies have focused on speech. One of the barriers is the lack of a
high-quality dataset. Since labeling audio is time-consuming and expensive, it
is difficult to establish a large-scale labeled audio dataset. Datasets such
as AudioSet [1] acquire a large amount of data, but the quality of the labels
is not guaranteed, which causes an additional problem, weakly supervised
learning.

Self-supervised learning is one approach to use datasets with such large scale
and unreliable labels. This can be done in a label-free manner and the
training is based on the assumption that the meaning of the data must remain
in the presence of noise or transformations. It aims to learn general-purpose
expressions that can be used anywhere without a target task and can be used
for various downstream tasks through transfer learning or fine-tuning.
Especially, the recent contrastive learning-based approach is promising as it
shows similar or superior performance to supervised training in the image
domain [2, 3]. In the case of the audio domain, various self-learning methods
for speech have been proposed [4, 5], and recently, it has been expanded to
various attempts targeting general audio [6, 7]. However, the downside is that
self-supervised learning itself is so general that the performance of certain
downstream tasks is not always guaranteed [8].

Semi-supervised learning [9] also utilizes unlabeled data. It differs from
self-supervised learning in that there is a clear target task in the training
procedure. In other words, the goal of semi-supervised learning is to perform
the target task with the help of unlabeled data. In general, most studies
consider unlabeled data sets to contain a large number of target classes [10,
11]. Therefore, the focus is on finding available data in unlabeled datasets.
In the implementation, the performance validation is done using labels for
parts of large datasets. However, this assumption may not always be correct.
Unlabeled data is literally unknown data, so there is no guarantee that the
target data will be included.

In this study, we propose a cross-domain semi-supervised learning method that
can be used when unlabeled data have a different class distribution than
labeled data. The proposed method combines the concepts of semi-supervised and
self-supervised training. More specifically, the network is trained to perform
target tasks while regularized by unlabeled data using a contrastive learning
approach. One advantage of our approach is that any unlabeled data can be
used. It is also useful as it is applied in the form of additional losses, so
it can be applied to almost all common networks. To the best of our knowledge,
this is the first study to perform neural network-based semi-supervised
training using unlabeled data with different label distributions.

Our contributions are as follows:

  * •

We propose a cross-domain semi-supervised training that can also be applied to
data of completely different classes from the target data.

  * •

We present a simple but efficient mixing strategy for applying contrastive
learning to the audio domain, which is batch-split mixing.

Figure 1: The overview of the proposed cross-domain semi-supervised training.
$f,g_{c},g_{p}$ are the shared feature extractor, classification layer, and
L2-normalized projection layer, respectively. The contrastive regularization
is applied as an additional loss function.

##  2 contrastive regularization for semi-supervised learning

Figure 2: The proposed audio augmentation strategy. The original audio sample
is augmented with a) random offset, b) mixing with another source, and c)
adding artificial noise.

The proposed semi-supervised method uses unlabeled data to learn differences
between samples. As shown in Fig. 1, proposed regularization is performed in
the form of multitask learning. However, it is different from the original
multitask learning setup in that the regularization task uses a different
unlabeled dataset than the target dataset. Rather, it is similar to the
constraint method using other pre-trained networks [12], but the networks used
for constraints are also trained at the same time. Such simultaneous training
is expected to prevent common transfer learning problems such as catastrophic
forgetting.

###  2.1 Contrastive regularization

Contrastive learning refers to a training method that recognizes identities
and differences between each sample rather than classification through class
labels [13, 14, 15]. This can also be applied to unlabeled data and is used to
allow the network to learn semantics that is robust to data transformations.
In particular, introducing NTXent loss [2] to solve contrastive learning as an
in-batch classification is a promising method in recent self-supervised
studies.

$NTXent=-log\frac{exp(x\cdot x^{+}/\tau)}{exp(x\cdot
x^{+}/\tau)+\sum_{x^{-}}exp(x\cdot x^{-}/\tau)}$ (1)

The proposed method introduces NTXent loss in semi-supervised learning. It
differs from the original method [2] in that the input of the NTXent loss is
one original sample and one augmented sample. We postulate that it will be
more useful to the target task as it directly includes the target dataset as
opposed to the original method of assuming a more general-purpose task. This
is a similar approach to the case where contrastive learning is used for a
specific purpose. In [16], contrastive learning is applied to extract audio
fingerprint. The learning process of contrastive learning can be viewed as a
fingerprint task performed on a subset of the database. Therefore, in this
case, the successful learning of contrastive learning guarantees the
performance of the target task.

1

Input: labeled dataset $(X,Y)$, unlabeled dataset $U$

Nets: shared feature extractor $f(.)$, classificaton layer $g_{c}(.)$, and L2
projection layer $g_{p}(.)$

2

3 $f,g_{c}\leftarrow$ WarmUp$(X,Y)$

// standard training for $f$ and $g_{c}$

4

5while  _epoch $<$ MaxEpoch_ do

6 while  _iter $<$ NumIters_ do

7 from $(X,Y)$, draw a mini-batch $\\{x_{b},y_{b}\\}$

8 from $U$, draw a mini-batch $\\{u_{b}\\}$

9

10 $s_{1},s_{2}$ = Split(Shuffle(Concat($x_{b},u_{b}$)))

// all batch data divided in half

11 $\overline{s}_{1}=$ Mixing($s_{1},s_{2},noise,SNR_{s1}$)

12 $\overline{s}_{2}=$ Mixing($s_{2},s_{1},noise,SNR_{s2}$)

// augmented samples

13

14 $L_{clf}=CrossEntropy(y_{b},g_{c}\cdot f(x_{b}))$

// classification loss

15

16 $L_{reg}^{(1)}$ =$NTXentLoss(g_{p}\cdot f(s_{1}),g_{p}\cdot
f(\overline{s}_{1}))$

17 $L_{reg}^{(2)}$ =$NTXentLoss(g_{p}\cdot f(s_{2}),g_{p}\cdot
f(\overline{s}_{2}))$

// regularization loss

18 $Loss=L_{clf}+\lambda_{reg}(L_{reg}^{(1)}+L_{reg}^{(2)})$

// total loss

19 $f,g_{c},g_{p}\leftarrow$ Optimize$(Loss)$

20

21

Algorithm 1 Training procedure of the cross-domain semi-supervised learning.

###  2.2 Audio augmentation strategy

We only used audio mixing as an augmentation method based on the unique
properties of audio. One attribute of audio that differs from an image is that
individual objects could be recognized even if multiple sources occur at the
same time. Unlike an image where one masks the other when another source is
added, the audio can be distinguished even if the two sources are mixed. We
thought that this property alone could make audio augmentation difficult
enough and could be applied to any kind of audio data. Audio mixing is often
used in contrastive learning [6, 7], but mainly to add noise or background.
The same idea as ours have been applied to triplet-based learning [17], but
not to contrast learning, to the best of our knowledge.

To do this in a batch, we split the batch in two and mix it with another
split. The mixing process involves mixing with random offsets, random
amplitudes, and additional artificial noise, as shown in Fig. 2. Audio events
can be shorter than the length of the network input, so after the audio mixing
process, they can exist individually on the time axis rather than as a mixture
at the same time. This can be a relatively easy problem that only requires
ensuring time-domain invariance, so additional noise is added to make it
difficult. The reason for using artificial noise is that the traditional
approach of using environmental sounds as noise, rather than the traditional
target class [18], is limited in use by what the target sound is.

Table 1: The results for semi-supervised training using the data with mismatched label on the ESC10. Performance refers to the averaged test accuracy (%) and standard deviation for predefined 5-folds. best means the highest performance in the entire training epochs and last means the averaged performance of the last 20 epochs. Note that $\dagger$ is evaluated on different dataset distributions. Data fraction | 10% | 25% | 50% | 75% | 100%
---|---|---|---|---|---
Supervised | best | 59.0 ± 5.0 | 67.2 ± 2.2 | 75.0 ± 2.9 | 76.2 ± 3.7 | 81.2 ± 1.6
| last | 53.0 ± 5.8 | 60.2 ± 2.4 | 65.9 ± 1.6 | 66.8 ± 5.0 | 72.3 ± 1.3
Supervised | best | 54.0 ± 4.8 | 76.2 ± 7.1 | 86.5 ± 1.5 | 89.5 ± 1.7 | 93.0 ± 3.6
(ImageNet init.) | last | 49.9 ± 4.8 | 73.0 ± 7.7 | 82.5 ± 1.7 | 84.9 ± 1.9 | 89.1 ± 3.6
Proposed | best | 65.2 ± 1.7 | 78.5 ± 3.7 | 90.0 ± 1.4 | 93.0 ± 1.3 | 95.2 ± 4.0
| last | 61.8 ± 0.4 | 75.6 ± 3.2 | 85.4 ± 1.9 | 87.6 ± 3.2 | 92.4 ± 2.5
Supervised [19] |  | 67.8 ± 4.0 | 82.5 ± 5.7 | 88.3 ± 3.0 | - | 91.7 ± 2.0
Deep-Co-Training$\dagger$ [19] |  | 75.7 ± 5.3 | 89.2 ± 4.0 | 91.7 ± 5.1 | - | -

###  2.3 Training procedure

The actual training procedure is listed in Algorithm 1. At the first stage,
the target classification network, which includes a shared feature extractor
and classification layer, is initialized through warm-up training. We found
that training the network from scratch eventually yields similar performance,
but takes more time to converge. Therefore, warm-up training is performed with
a relatively high learning rate to shorten the training time.

In the semi-supervised learning stage, the labeled data is used for typical
classification as well as warm-up training. The difference is that
regularization that uses both labeled and unlabeled at the same time has been
added. This results in NTXent loss after batch split mixing, as mentioned
earlier. The total loss is expressed as the sum of the classification loss and
regularization loss.

##  3 Experiments

###  3.1 Dataset

  * •

ESC10 and ESC50: The dataset for environmental sound classification [20] has a
total of 2,000 audio samples for various audio events. ESC10 and ESC50 have 40
excerpts of 5 seconds for each class. It has the predefined 5-fold validation
split.

  * •

US8K: The urban sound 8K [21] contains 8,732 labeled sound excerpts less than
4 seconds of 10 environmental sounds.
The audio length varies from 50 milliseconds to 4 seconds and it has a
predefined 10-fold validation split.

###  3.2 Network Architecture

A network takes 5-seconds of audio as input. It then converted to dB-scale log
mel spectrogram, a common feature used in audio analysis. For the STFT, the
nfft size was 1,024 and the hop size was 230, and the frequency information
from 50Hz to 10,000Hz was converted to a 128-bin mel frequency. Thereafter,
the mono channel mel spectrogram is stacked three times on the channel axis
with the same value, to use the pre-trained weights [22, 23], resulting in the
shape of (128, 480, 3).

ResNet50 [24], which is often used in self-supervised studies [7, 25] with
comparable performance in audio event classification, and is known to be
effective for similar temporal unit analysis [26], was used as a feature
extractor to output a 2,048-dimensional representation. The classification
task branch consists of a set of 128-dimensional fully connected layer, batch
normalization, and ReLU activation, followed by the final output layer. The
contrastive regularization branch is similarly connected through
128-dimensional fully connected layer, batch normalization, and ReLU
activation, and as an output, an L2-normalized 64-dimensional embedding
vector. The final output of this branch means an embedding for each data,
comparisons with another embedding can be done through a dot product, the
larger it means the two data are more similar.

###  3.3 Experiment configuration

Two experiments were conducted with data in and out of the dataset.

  * •

Semi-supervised training using the data with mismatched label: ESC10 was used
as labeled data. We defined ESC40 as ESC50 excluding ESC10 and used as
unlabeled data in the proposed method. As with self-supervised studies [2],
all unlabeled labels were always used, and parts of labels were used.
Deep-Co-Training$\dagger$ [19] is a method of pseudo-labeling unlabeled data
using two networks that are trained to provide different predictive values and
complement each other. A portion of ESC10 is used as labeled and the rest is
used as unlabeled.

  * •

Cross-dataset semi-supervised training: Both ESC50 and US8K were used in the
experiment, one for labeled data and the other for unlabeled data. That is,
when classifying ESC50, US8K was used for contrastive regularization and vice
versa.

To reduce the influence of other variables in the experiment, no changes such
as augmentation were made to the data pipe in target classification. The
networks are initialized with the weights learned from the ImageNet. This is
because the kernel trained to analyze the image is also known effectively for
the classification of environmental sounds [22, 23].

###  3.4 Implementation details

All audio data is resampled to 22,500 Hz and the scale is normalized based on
energy. When each audio is used for network input, the short audio is zero-
padded for up to 5 seconds. Audio mixing was done in the SNR range of [-5, 20]
dB. A total of 6 artificial noises including blue, brown, grey, pink, violet,
and white were used and added to less than 6 dB.

Warm-up was performed until the training accuracy was high enough, ESC10 and
ESC50 performed 20 epochs, and US8K performed 10 epochs. Adam optimizer [27]
with learning rate 1e-4 is used for warm-up training and with learning rate
1e-5 used for semi-supervised training. For the proposed method, we trained
the network 100 epochs with 8 labeled and 32 unlabeled samples, and for
supervised learning, 200 epochs were used because the test accuracy was
relatively slow to converge.

Weight for regularization loss ($\lambda_{reg}$) is set to 0.05 taking into
account the scale of the loss value. The temperature of the regularization
loss is set to 0.01. We experimentally confirmed that training was successful
in the range of [0.001, 0.1]. All experiments were implemented using
Tensorflow [28], and the network structure and weights of ResNet50V2 included
in the library were used.

##  4 Results and Discussion

###  4.1 Cross-domain semi-supervised learning on ESC10

Table 2: Experimental results for cross-dataset semi-supervised training. Performance refers to the average accuracy and standard deviation for predefined k-fold validation. Target data | Unlabeld data | Accuracy (%)
---|---|---
ESC10 | - | 89.1 ± 3.6
ESC10 | ESC40 | 92.4 ± 2.5
ESC50 | - | 73.6 ± 1.0
ESC50 | US8K | 83.5 ± 3.0
US8K | - | 78.6 ± 4.3
US8K | ESC50 | 79.7 ± 4.3

The results for semi-supervised training using the data with mismatched labels
are listed in Table 1. best is the average of the highest accuracy observed
within the entire epochs at each fold and can be seen as the maximum
performance that can be achieved in this experimental setup. last is the
average accuracy of the last 20 epochs and is generally expected to be
achieved. Except for 10%, an almost constant performance difference is
observed between supervised learning using different initialization methods.
Initializing the weights from the image domain had a huge impact on the
performance improvement, reconfirming that the kernel trained in the image
domain still helped the spectrogram domain. However, the reversed performance
at 10% shows that these weighted initialization methods may harm the network
when data is scarce. In this case, since the data is very small, 3 samples per
class, it is thought that the network focuses on other elements learned from
images rather than semantic elements, but more research is needed.

The proposed method shows the highest accuracy in all experimental settings.
Since the proposed method does not affect the classification layer
structurally, it can be seen that a better semantic representation is induced
by adding contrastive regularization. Multitask learning can improve target
performance depending on the composition of the task [29, 30], and in this
case, it can be seen that it behaved as we intended. It also shows that even
if each branch uses data from a different label distribution, or a completely
different dataset (Table 2), the multitask learning assumptions still work.
The performance improvement of US8K is less than that of ESC50. Unlike ESC50,
US8K consists of audio clips shorter than 5 seconds, and we assume that
energy-based mixing is expected to adversely affect.

The proposed method has relatively low performance and performance improvement
compared to regular semi-supervised training [19] that uses the rest of
labeled data as unlabeled data. Excluding the 10% condition, the standard
semi-supervised training shows a higher performance improvement and converges
on performance with the full data at about 50%. However, at 100%, the proposed
method showed higher performance even though there are no various performance
improvement techniques. This shows the difference between the standard method
and the proposed method. The standard method always uses the same class of
data under strong assumptions, so it can reach the desired performance faster
with less data, but there is a limit to the performance gain.

###  4.2 Generalization effect

Figure 3: The learning curve for each learning method under the condition of
using 100% label data of ESC10. Each accuracy represents the average accuracy
of each epoch in a 5-fold experiment.

We assume that the main factors in improving performance are stable training
and generalization. The difference between last and best in the proposed
method is always smaller than that of supervised, and smaller or similar to
that of the weight initialized network. It implies that contrastive
regularization covers the target task as we assumed and guides the network in
a direction that helps generalization even when the data is scarce. These
trends also appear in Fig 3. The supervised training shows a typical learning
curve. On the other hand, in the proposed method, the training accuracy
increases rapidly, and both training and test accuracy remain constant. The
rapid training seems to be due to contrastive regularization as well as warm-
up training. In addition to fast training set convergence, preventing
overfitting also seems to be an effect of the proposed method. High training
accuracy was quickly achieved through warm-up training, but test accuracy was
also consistently high through contrastive regularization. Through this, it is
expected that much effort to find the optimal point could be reduced.

###  4.3 Ablation study

Table 3: Ablation study on ESC50 dataset. Proposed | -Mixing | -Unlabeled | Supervised
---|---|---|---
83.5 ± 3.0 | 82.0 ± 3.0 | 79.5 ± 1.6 | 73.6 ± 1.0

Table 3 shows the results of excluding each element from the overall proposed
method. In -Unlabeled condition, contrast regularization is added only within
a labeled data set, with a performance improvement of 5.9% points over the
normal supervised condition. It shows that the proposed method is also helpful
in a typical supervised setting. Without the data mixing strategy of the
proposed method, -Mixing, which applies only time offset augmentation, there
is a performance improvement of 8.4% points over the supervised condition.
This implies that the increased diversity of tasks that external data can
provide is a major factor in contrastive regularization.

##  5 Conclusion

In this study, we proposed cross-domain semi-supervised training that works
with completely different data distributions. We introduced contrastive
learning to the semi-supervised framework and proposed a novel augmentation
method considering the audio characteristics. Experimental results have proven
the effectiveness of the proposed method for audio event classification.
Considering that the proposed method was implemented in a simplified form and
the characteristics of the proposed method that it can be applied to any
network, there is a lot of room for performance improvement. In particular, we
expect that various augmentation methods and large-batch, which are the main
factors in contrastive learning, may improve the performance further.

## References

  * [1] J. F. Gemmeke _et al._ , “Audio set: An ontology and human-labeled dataset for audio events,” in _2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2017, pp. 776–780.
  * [2] T. Chen, S. Kornblith, M. Norouzi, and G. Hinton, “A simple framework for contrastive learning of visual representations,” in _International conference on machine learning_. PMLR, 2020, pp. 1597–1607.
  * [3] T. Chen, S. Kornblith, K. Swersky, M. Norouzi, and G. Hinton, “Big self-supervised models are strong semi-supervised learners,” _arXiv preprint arXiv:2006.10029_ , 2020.
  * [4] S. Pascual, M. Ravanelli, J. Serrà, A. Bonafonte, and Y. Bengio, “Learning problem-agnostic speech representations from multiple self-supervised tasks,” _Proc. Interspeech 2019_ , pp. 161–165, 2019.
  * [5] S. Schneider, A. Baevski, R. Collobert, and M. Auli, “wav2vec: Unsupervised pre-training for speech recognition,” _Proc. Interspeech 2019_ , pp. 3465–3469, 2019.
  * [6] D. Niizumi, D. Takeuchi, Y. Ohishi, N. Harada, and K. Kashino, “Byol for audio: Self-supervised learning for general-purpose audio representation,” _arXiv preprint arXiv:2103.06695_ , 2021.
  * [7] L. Wang and A. v. d. Oord, “Multi-format contrastive learning of audio representations,” _arXiv preprint arXiv:2103.06508_ , 2021.
  * [8] T. Xiao, X. Wang, A. A. Efros, and T. Darrell, “What should not be contrastive in contrastive learning,” _arXiv preprint arXiv:2008.05659_ , 2020.
  * [9] X. Zhou and M. Belkin, “Semi-supervised learning,” in _Academic Press Library in Signal Processing_. Elsevier, 2014, vol. 1, pp. 1239–1269.
  * [10] D. Berthelot, N. Carlini, I. Goodfellow, N. Papernot, A. Oliver, and C. Raffel, “Mixmatch: A holistic approach to semi-supervised learning,” _arXiv preprint arXiv:1905.02249_ , 2019.
  * [11] K. Sohn _et al._ , “Fixmatch: Simplifying semi-supervised learning with consistency and confidence,” _arXiv preprint arXiv:2001.07685_ , 2020.
  * [12] J. Johnson, A. Alahi, and L. Fei-Fei, “Perceptual losses for real-time style transfer and super-resolution,” in _European conference on computer vision_. Springer, 2016, pp. 694–711.
  * [13] J. Bromley, J. W. Bentz, L. Bottou, I. Guyon, Y. LeCun, C. Moore, E. Säckinger, and R. Shah, “Signature verification using a “siamese” time delay neural network,” _International Journal of Pattern Recognition and Artificial Intelligence_ , vol. 7, no. 04, pp. 669–688, 1993.
  * [14] R. Hadsell, S. Chopra, and Y. LeCun, “Dimensionality reduction by learning an invariant mapping,” in _2006 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’06)_ , vol. 2. IEEE, 2006, pp. 1735–1742.
  * [15] F. Schroff, D. Kalenichenko, and J. Philbin, “Facenet: A unified embedding for face recognition and clustering,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2015, pp. 815–823.
  * [16] S. Chang, D. Lee, J. Park, H. Lim, K. Lee, K. Ko, and Y. Han, “Neural audio fingerprint for high-specific audio retrieval based on contrastive learning,” in _ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2021, pp. 3025–3029.
  * [17] A. Jansen, M. Plakal, R. Pandya, D. P. Ellis, S. Hershey, J. Liu, R. C. Moore, and R. A. Saurous, “Unsupervised learning of semantic audio representations,” in _2018 IEEE international conference on acoustics, speech and signal processing (ICASSP)_. IEEE, 2018, pp. 126–130.
  * [18] J. Thiemann, N. Ito, and E. Vincent, “The diverse environments multi-channel acoustic noise database (demand): A database of multichannel environmental noise recordings,” in _Proceedings of Meetings on Acoustics ICA2013_ , vol. 19, no. 1. Acoustical Society of America, 2013, p. 035081.
  * [19] L. Cances and T. Pellegrini, “Comparison of deep co-training and mean-teacher approaches for semi-supervised audio tagging,” in _IEEE 46th International Conference on Acoustics, Speech and Signal Processing (ICASSP 2021)_ , 2021.
  * [20] K. J. Piczak, “ESC: Dataset for Environmental Sound Classification,” in _Proceedings of the 23rd Annual ACM Conference on Multimedia_. ACM Press, 2015, pp. 1015–1018. [Online]. Available: http://dl.acm.org/citation.cfm?doid=2733373.2806390
  * [21] J. Salamon, C. Jacoby, and J. P. Bello, “A dataset and taxonomy for urban sound research,” in _22nd ACM International Conference on Multimedia (ACM-MM’14)_ , Orlando, FL, USA, Nov. 2014, pp. 1041–1044.
  * [22] K. Palanisamy, D. Singhania, and A. Yao, “Rethinking cnn models for audio classification,” _arXiv preprint arXiv:2007.11154_ , 2020.
  * [23] A. Guzhov, F. Raue, J. Hees, and A. Dengel, “Esresnet: Environmental sound classification based on visual domain models,” _arXiv preprint arXiv:2004.07301_ , 2020.
  * [24] K. He, X. Zhang, S. Ren, and J. Sun, “Identity mappings in deep residual networks,” in _European conference on computer vision_. Springer, 2016, pp. 630–645.
  * [25] J.-B. Grill, F. Strub, F. Altché, C. Tallec, P. H. Richemond, E. Buchatskaya, C. Doersch, B. A. Pires, Z. D. Guo, M. G. Azar, _et al._ , “Bootstrap your own latent: A new approach to self-supervised learning,” _arXiv preprint arXiv:2006.07733_ , 2020.
  * [26] M. Won, A. Ferraro, D. Bogdanov, and X. Serra, “Evaluation of cnn-based automatic music tagging models,” _arXiv preprint arXiv:2006.00751_ , 2020\.
  * [27] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [28] M. Abadi _et al._ , “TensorFlow: Large-scale machine learning on heterogeneous systems,” 2015, software available from tensorflow.org. [Online]. Available: https://www.tensorflow.org/
  * [29] R. Caruana, “Multitask learning,” _Machine learning_ , vol. 28, no. 1, pp. 41–75, 1997.
  * [30] D. Lee, J. Lee, J. Park, and K. Lee, “Enhancing music features by knowledge transfer from user-item log data,” in _ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2019, pp. 386–390.


# Spontaneous symmetry breaking in finite systems and anomalous order-
parameter correlations

F. K. Diakonos<EMAIL_ADDRESS>Faculty of Physics, University of Athens,
GR-15784 Athens, Greece  Y. F. Contoyiannis<EMAIL_ADDRESS>Department of
Electrical and Electronics Engineering, University of West Attica, GR-12244,
Egaleo, Greece  Faculty of Physics, University of Athens, GR-15784 Athens,
Greece  S. M. Potirakis<EMAIL_ADDRESS>Department of Electrical and
Electronics Engineering, University of West Attica, GR-12244, Egaleo, Greece

###### Abstract

We show that the spontaneous symmetry breaking can be defined also for finite
systems based on the properly defined jump probability between the ground
states in the 2d and 3d Ising models on a square and a cubic lattice
respectively. Our analysis reveals the existence of an interval in the
temperature (control parameter) space within which the spontaneous symmetry
breaking takes place. The upper limit of this region is the pseudocritical
point where the symmetric vacuum bifurcates in energetically degenerate non-
symmetric vacua, initiating the spontaneous symmetry breaking process. The
lower limit, identified as the temperature value at which the spontaneous
symmetry breaking is completed, is characterized by maximal characteristic
time for the decay of magnetization (order parameter) auto-correlations. We
argue that this anomalous enhancement of auto-correlations is attributed to
the transition from type I to on-off intermittency in the order parameter
dynamics. Possible phenomenological implications of this behaviour are briefly
discussed.

Critical phenomena are the origin of the physical world’s structure at all
scales. Their realization in physical systems in usually related to the
spontaneous breaking of symmetries implied by the parametric dependence of the
vacuum state of a system on control parameters. According to this scenario a
symmetry respecting ground state of a physical system, becoming dynamically
unstable, bifurcates in two (or more) dynamically stable ground states
whenever a critical value of the control parameter (critical point) is
crossed. The emerging ground states form a degenerate space and the original
symmetry transform, leaving the dynamics of the system invariant, generates a
mapping between the different ground states. Then, the underlying symmetry is
spontaneously broken since the ground states do not remain invariant under the
corresponding transformation. Although this scenario is rigorously defined in
infinite systems, it does not apply to real physical systems which posses
finite size. In fact, in a finite system a symmetry can never break
spontaneously Binder1985  in the strict mathematical sense. Obviously, the
spontaneous breaking of a symmetry in a finite system cannot take place at the
pseudocritical point signalling the occurrence of the bifurcation of the
vacuum state, which is a necessary but not sufficient condition. In fact, well
after the crossing of the bifurcation point in the control parameter space,
the system behaves as having a ground state which is still symmetric, since
the new vacua, which emerge, communicate through the underlying dynamics. In
this sense the physical system is in a symmetric superposition of the
corresponding vacua, a term borrowed by quantum mechanics. Effectively, this
quantum mechanical picture is always valid also for a classical field (like
the order parameter density) defined in a finite spatial region. This is due
to the fact that there is a finite jump probability from the one ground state
to the other when the system possesses a fixed size.

The present letter aims to present a definition of spontaneous symmetry
breaking (SSB) in a finite system. To achieve this task, we study the jump
probabilities between the new vacua occurring after crossing the
pseudocritical control parameter value. We show that these probabilities, as a
function of the control parameter, posses an exponential suppression term. The
latter defines a characteristic control parameter value, after which, the jump
probability becomes negligibly small and therefore the symmetry can be
considered as broken. In practice, this situation is realized in nature since
the time interval during which an observer probes a finite system is also
finite. We will call this characteristic control parameter value as the SSB
completion point. Thus, the pseudocritical point and the SSB completion point
define an interval in the control parameter space within which the symmetry,
despite the emergence of the new vacua, remains unbroken. Departing from this
interval in the control parameter space, the communication between the new
vacua, through the order parameter dynamics, is practically lost for finite
observation time. The width of the control parameter region, within which the
completion of the SSB takes place, depends on the size of the system and
increases as the latter decreases. We will in the following refer to this
region as the completion region of SSB. Analyzing the order parameter dynamics
in the SSB completion region we reveal an anomalous behaviour of the
corresponding auto-correlations expressed by an enhancement in the border of
this region. Furthermore, employing the Method of Critical Fluctuations (MCF)
Contoyiannis2002  to analyse the order parameter time-series at different
control parameter values within the SSB completion region, we observe a
transition of the order parameter dynamics from intermittency type I
Manneville1980 ; Schuster1995  to the on-off intermittency regime Platt1993 .
Finally, as we discuss at the end of this letter, despite its theoretical
character, our proposal can explain eventual discrepancy between experimental
and simulated results concerning the location of the critical point in a
physical system Gavai2016 ; Ratti2018 .

To further illustrate our claims it is necessary to employ a concrete system.
Here we will consider the magnetization dynamics of the ferromagnetic Ising
model with next neighbour interactions in 2 and 3 dimensions, defined in a
square and cubic lattice respectively. It is given by the Hamiltonian:

$H=-J\displaystyle{\sum_{\langle
i,j\rangle}}s_{i}s_{j}~{}~{}~{}~{}~{};~{}~{}~{}~{}~{}J>0$ (1)

where $\langle i,j\rangle$ denotes summation over next neighbours. As usual,
the coupling $J$ can be absorbed into the inverse temperature $\beta=T^{-1}$.
The system obeys periodic boundary conditions. The simulations are performed
in lattices of relatively small size to reveal the phenomenon addressed
previously in the clearest way. The Metropolis algorithm is used in all
simulations, taking into account autocorrelations in order to produce
statistically independent spin configurations. Our calculations are performed
using $\mathcal{N}=2\cdot 10^{6}$ sweeps/spin while we have checked the
convergence of our results for specific values of the size $L$ and the
temperature $T$, in both the 2d as well as the 3d model, using several
$\mathcal{N}$ values ranging from $10^{5}$ to $10^{7}$ sweeps/spin.

We focus on the properties of the average magnetization
$M=\displaystyle{{\sum_{i}s_{i}\over\Omega}}$ at equilibrium. With $\Omega$ we
denote the volume of the considered lattice. Firstly, we investigate the
distribution $P(M)$ for various temperature values. In Fig. 1a we show $P(M)$
at $\beta^{-1}=2.3$ for a square lattice of $128\times 128$ size. This
temperature value corresponds to the pseudo-critical one, notated as $T_{pc}$,
at which a plateau region around $M=0$ is formed in $P(M)$. This property
reflects the fact that the fixed point $M=0$ becomes marginally stable as an
intermediate step in its transition from a stable to an unstable one with
decreasing temperature. In Fig. 1b we show $P(M)$ for $\beta^{-1}=2.27$ where
the formation of the two new stable minima in the magnetization effective
potential $V_{eff}(M)$ is clearly seen. Notice that the minima of $V_{eff}(M)$
appear as maxima in the distribution $P(M)$ since $P(M)\propto e^{-\Omega
V_{eff}(M)}$. However, despite the presence of these two minima in
$V_{eff}(M)$, the dynamics in the magnetization space are still respecting the
$Z(2)$ symmetry, since there is finite, symmetric in the direction,
probability to jump from the one minimum to the other. For a clearer display
of this behaviour we include in this plot an inset, presenting the
magnetization time-series after equilibrium is reached. The jumps from the
left stable minimum to the right and vice versa are clearly observable in this
inset.

Figure 1: The magnetization distribution in 2d Ising square lattice with
$L=128$ and periodic boundary conditions for: (a) $T\approx T_{pc}=2.3$ and
(b) $T\approx T_{SSB}=2.27$. In the inset in (b) is shown a part of the
corresponding magnetization time-series.

In the following we define as left-right and right-left jump probabilities
$p_{LR}$ and $p_{RL}$ as:

$p_{LR}=\displaystyle{\frac{n_{-+}}{\mathcal{N}}}~{}~{}~{}~{}~{}~{};~{}~{}~{}~{}~{}~{}p_{RL}=\displaystyle{\frac{n_{+-}}{\mathcal{N}}}$
(2)

where $n_{-+}$ ($n_{+-}$) is the number of sign changes of the magnetization
$M$ from negative (positive) to positive (negative) values respectively, in a
given number of time steps $\mathcal{N}$ measured in sweeps. Since in the
considered temperature range, due to $Z(2)$ symmetry, it holds
$p_{LR}=p_{RL}$, we restrict our study on the jump probability $p_{RL}$
without loss of generality. Our goal is to calculate $p_{RL}$ at different
subcritical temperatures close to the pseudocritical $T_{pc}$ and for
different lattice sizes $L$ to obtain $p_{RL}(L,T)$. The results of our
calculations are presented in Figs. 2(a,b). In Fig. 2a we show the function
$Q(T,L)=-{1\over
b(L)}\ln\displaystyle{\left(\frac{p_{RL}(T,L)}{p_{RL}(T_{pc}(L),L)}\right)}$
for the 2d Ising model and various sizes of the square lattice $L=64$ (black
open circles), $L=80$ (red crosses), $L=96$ (green up triangles), $L=112$
(blue stars) and $L=128$ (violet triangles down). The scaling factor $b(L)$ is
used to reveal the universal behaviour obeyed by the calculated jump
probabilities. On the x-axis we set the temperature difference $T_{pc}(L)-T$
with $T_{pc}(L)$ the pseudocritical temperature corresponding to each
$L$-value, calculated with separate simulations and cross-checked with the
literature values Malakis2004 ; Kaupuzs2011 . We observe that the results for
the different lattice sizes are concentrated on a single curve providing a
universal law for the decay of the jump probability with temperature in the
immediate neighbourhood of the pseudocritical point within the subcritical
region. The dashed black line is the function $(T_{pc}-T)^{{3\over 2}}$. As
shown in Fig. 2b, similar results for $Q(T,L)$ are also obtained for the 3d
Ising case using cubic lattices with $L=16$ (black open circles), $L=20$ (red
crosses), $L=24$ (green up triangles), $L=28$ (blue stars) and $L=32$ (violet
triangles down).

Figure 2: The quantity $Q(T,L)=-{1\over
b(L)}\ln\displaystyle{\left(\frac{p_{RL}(T,L)}{p_{RL}(T_{pc}(L),L)}\right)}$
as a function of $T_{pc}-T$ for: (a) the 2d Ising model using lattices with
size $L=64,~{}80,~{}96,~{}112,~{}128$ and (b) the 3d Ising model using
lattices with size $L=16,~{}20,~{}24,~{}28,~{}32$. The scaling factors $b(L)$
in each case are obtained though fitting.

The universal behaviour in Figs. 2a and 2b leads to a modified exponential
decay law given by the function:

$p_{RL}(T,L)=p_{RL}(T_{pc}(L),L)e^{-b(L)(T_{pc}-T)^{{3\over 2}}}$ (3)

The characteristic exponent $b(L)$ contains a temperature scale which
determines the characteristic temperature $T_{SSB}$ below which the jump
probability $p_{RL}$ becomes negligibly small in a lattice of size $L$. In
fact we can write $b(L)=\frac{\sigma(L)}{T_{pc}^{{3\over 2}}}$ with
$\sigma(L)$ dimensionless. The jump probability becomes negligible small when
$\sigma(L)\left(\frac{T_{pc}-T}{T_{pc}}\right)^{{3\over
2}}\approx\left(\frac{8\Gamma[{2\over 3}]}{3}\right)^{3\over 2}$ 111It is a
linear approximation based on equality of areas.. In this limit $T\approx
T_{SSB}$ leading to:

$T_{SSB}\approx T_{pc}\left(1-\frac{8\Gamma[{2\over
3}]}{3}\left(\frac{1}{\sigma(L)}\right)^{{2\over 3}}\right)$ (4)

As shown in Fig. 3 the coefficient $\sigma(L)$ has a power-law dependence on
$L$ in the form $\sigma(L)\sim L^{2}$ for the 2d Ising and $\sigma(L)\sim
L^{\frac{7}{3}}$ in the 3d Ising case.

Figure 3: The coefficient $\sigma(L)$ used in Eq. (3)
($b(L)=\frac{\sigma(L)}{T_{pc}^{3/2}}$) as a function of $L$ for the 2d Ising
(olive circles) and the 3d Ising (blue stars) case. The red lines,
representing the scaling laws $0.09L^{2}$ (2d Ising) and $L^{{7\over 3}}$ (3d
Ising), are displayed to guide the eye.

Combining Eq. (4) with the results shown in Fig. 3 we obtain that the size
$\Delta T=T_{pc}-T_{SSB}$ of the SSB completion region in temperature space
depends on the lattice size $L$, shrinking with increasing $L$ according to
the power-laws:

$\Delta T_{2d}\sim L^{-4/3}~{}~{}~{}~{};~{}~{}~{}~{}\Delta T_{3d}\sim
L^{-14/9}$ (5)

Our next goal is to explore in more detail the magnetization dynamics within
the SSB completion region. In a previous work we have shown, in the context of
the 3d Ising model, that the magnetization dynamics at the pseudocritical
point possesses characteristics of type I intermittency, expressed by a power-
law form of the distribution
$N(\lambda)\sim\lambda^{-\frac{\delta+1}{\delta}}$ ($\delta$ is the isothermal
critical exponent) of the waiting times $\lambda$ (laminar lengths) within the
neighbourhood of the bifurcating fixed point $M=0$ Contoyiannis2002 . Here, we
argue that varying the temperature in the SSB completion region, departing
from $T_{pc}$ and approaching $T_{SSB}$, a gradual transition from type I
intermittent dynamics to on-off intermittency Platt1993  takes place. The
benchmark of this transition is the collective spin dynamics around the fixed
point at $M=0$. For $T\approx T_{pc}$ this dynamics is dominated by the
marginal instability of the $M=0$ fixed point leading to a typical type I
intermittent behaviour ala Manneville and Pomeau Manneville1980 . As the
temperature approaches $T_{SSB}$ the intermittent behaviour is expressed by
long oscillations around the new stable minima interrupted by burst dynamics
describing the transition from the neighbourhood of the one stable minimum to
the other. The stable minima in $V_{eff}(M)$ represent the stable manifolds in
the collective spin dynamics. The communication between them is achieved
through the $M=0$ channel formed by spin configurations in which the clusters
with positively aligned spins ($s_{i}=1$) are equal (in total size) with the
clusters with negatively aligned spins ($s_{j}=-1$). This channel acts as the
unstable manifold and becomes more and more narrow (statistically suppressed)
as $T\to T_{SSB}$ due to the coalescence of clusters of the same orientation.
As a consequence the on-off intermittency becomes established, a behaviour
which is expected to be observable in the modification of the laminar length
distribution $N(\lambda)$ approaching $T_{SSB}$. As shown in Cenys1997  for
noisy on-off intermittency the laminar length distribution becomes a power-law
with exponent $\frac{3}{2}$. We have calculated the laminar length
distribution around $M=0$ for the 3d Ising system at the temperatures
$T_{pc}=4.515$ and $T=4.45$ which lies very close to $T_{SSB}$ according to
the relation in Eq. (4). We pick out in Fig. 4 the bulk of the distribution
where the power-law behaviour is clearly seen. The reason for choosing this
region of $\lambda$-values is that we want to avoid the saturation region
close to $\lambda=1$ (notice that the time is measured in sweeps) as well as
the exponential tail of the distribution, induced by the finite size $L$ of
the considered system. The increase of the slope of the power-law decay from
$\approx 1.2$ to $\approx 1.5$ as $T_{SSB}$ is approached, can be clearly
seen.

Figure 4: The distribution $N(\lambda)$ (in log-log scale) of the waiting
times $\lambda$ around the marginally unstable fixed point $M=0$ for the 3d
Ising model at $T=4.515$, which is approximately $T_{pc}$ (blue stars) and at
$T=4.45$, which is approximately $T_{SSB}$ (black circles). The red dashed
lines indicate the power-laws with exponents $1.2$ and $1.58$ (fit results)
respectively and they are displayed to guide the eye.

Related to this transition from type I to on-off intermittency in the
magnetization dynamics as the temperature is varied from $T_{pc}$ to $T_{SSB}$
is an anomalous behaviour of the magnetization autocorrelation function
displayed in Figs. 5(a,b).

Figure 5: In (a) we show the autocorrelation function of the magnetization
times-series for the 2d Ising model at three different temperatures $T=2.25$
($T$ outside ot the SSB completion region), $T_{SSB}=2.27$ and $T_{pc}=2.31$.
In (b) it is shown the corresponding quantity for the 3d Ising model now at
temperatures $T=4.3$ ($T$ outside ot the SSB completion region), $T_{SSB}=4.4$
and $T_{pc}=4.5$.

We clearly observe in both cases an increase of the autocorrelations as $T$
approaches $T_{SSB}$. This anomalous enhancement of autocorrelations in the
on-off intermittency is due to the fact that, although the mean magnetization
in this case remains close to zero, the magnetization time-series contains
mainly small oscillations around the minima of the effective potential
$V_{eff}(M)$ (maxima of the distribution $P(M)$) which in turn may lead to
large waiting times in this magnetization region.

The analysis presented in this letter can have significant impact on our
understanding of critical phenomena as they are realized in true physical
systems. In particular, it clarifies the fact that the establishment of a non-
vanishing mean value of the order parameter, taking place below $T_{SSB}$,
does not coincide with the pseudocritical point $T_{pc}$ of a system where the
stable fixed point bifurcates becoming unstable. Only when the, stable in the
symmetric phase and unstable in the spontaneously broken phase, fixed point
disconnects from the other (stable) fixed points in the order parameter
dynamics the spontaneous breaking of the underlying symmetry is completed. A
characteristic property of the SSB completion region is that the distance
$T_{pc}-T_{SSB}$ increases with decreasing system’s size. As a consequence in
a microscopic system this distance can be large. For example, such a case
could appear in the Lattice QCD simulations at zero chemical potential. There,
the critical point is estimated through the formation of a non-zero
$\langle\bar{q}q\rangle$ condensate, which in fact takes place at $T_{SSB}$.
However, since the calculations are performed in quite small lattices
(systems) the corresponding pseudocritical temperature could be significantly
higher. Such a scenario could explain inconsistencies occurring in the attempt
to apply Lattice QCD results to the phenomenology of ultra-relativistic ion
collisions, related to the experimental search for the QCD critical endpoint
Gavai2016 ; Ratti2018 .

## References

  * (1) K. Binder, J. Comp. Phys. 59, 1 (1985).
  * (2) Y. F. Contoyiannis, F. K. Diakonos and A. Malakis, Phys.Rev. Lett. 89, 035701 (2002).
  * (3) ”Deterministic Chaos: An Introduction”, H. G. Schuster, Wiley-VCH, 1995.
  * (4) Y. Pomeau and P. Manneville, Comm. in Math. Phys. 74, 189 (1980).
  * (5) N. Platt, E. A. Spiegel and C. Tresser, Phys. Rev. Lett. 70, 279 (1993).
  * (6) R. V. Gavai, Cont. Phys. 57, 350 (2016).
  * (7) C. Ratti, Rep. Prog. Phys. 81, 084301 (2018).
  * (8) A. Malakis, A. Peratzakis and N. G. Fytas, Phys. Rev. E 70, 066128 (2004).
  * (9) J. Kaupužs, J. Rimšāns and R. V. N. Melnik, Ukr. J. Phys. 56, 845 (2011).
  * (10) A. Čenys, A. N. Anagnostopoulos and G. L. Bleris, Phys. Lett. A 224, 346 (1997).