# Energy distribution of small-scale flares derived using genetic algorithm

Toshiki Kawai Institute for Space-Earth Environmental Research, Nagoya
University,
Furo-cho, Chikusa-ku, Nagoya, Japan  Shinsuke Imada  Institute for Space-Earth
Environmental Research, Nagoya University,
Furo-cho, Chikusa-ku, Nagoya, Japan

(Received May 1, 2020; Revised November 9, 2020; Accepted November 10, 2020)

###### Abstract

To understand the mechanism of coronal heating, it is crucial to derive the
contribution of small-scale flares, the so-called nanoflares, to the heating
up of the solar corona. To date, several studies have tried to derive the
occurrence frequency distribution of flares as a function of energy to reveal
the contribution of small-scale flares. However, there are no studies that
derive the distribution with considering the following conditions: (1)
evolution of the coronal loop plasma heated by small-scale flares, (2) loops
smaller than the spatial resolution of the observed image, and (3)
multiwavelength observation. To take into account these conditions, we
introduce a new method to analyze small-scale flares statistically based on a
one-dimensional loop simulation and a machine learning technique, that is,
genetic algorithm. First, we obtain six channels of SDO/AIA light curves of
the active-region coronal loops. Second, we carry out many coronal loop
simulations and obtain the SDO/AIA light curves for each simulation in a
pseudo-manner. Third, using the genetic algorithm, we estimate the best
combination of simulated light curves that reproduce the observation.
Consequently, the observed coronal loops are heated by small-scale flares with
energy flux larger than that typically required to heat up an active region
intermittently. Moreover, we derive the occurrence frequency distribution
which have various power-law indices in the range from 1 to 3, which partially
supports the nanoflare heating model. In contrast, we find that $90\%$ of the
coronal heating is done by flares that have energy larger than
$10^{25}~{}\mathrm{erg}$.

Sun: corona — Sun: flares — hydrodynamics

††journal: ApJ

##  1 Introduction

Understanding the mechanism of the heating up of the solar corona to a few
million Kelvin or above is one of the long-standing problems in solar physics.
To date, two primary mechanisms have been proposed to explain how the corona
is heated, small-scale magnetic reconnection and wave dissipation. In the
former model, the corona heats up due to small-scale impulsive heating events,
so-called nanoflares ($E\simeq 10^{24}~{}\mathrm{erg}$), which are related to
coronal magnetic reconnections (e.g., Parker, 1988). In the wave dissipation
model, the corona heats up due to Alfvén waves propagated from the surface,
which are excited due to convection (e.g., Antolin & Shibata, 2010). In recent
years, much progress has been made in regard to theories and observations for
both models; however, a definitive solution to this problem is yet to be
achieved (e.g., Klimchuk, 2006, 2015). From the perspective of the nanoflare
model, quantifying the contribution of small-scale flares in the heating up of
the corona is crucial to understand their contribution to coronal heating.

Parker (1983, 1988) proposed that small-scale flares such as microflares and
nanoflares are magnetic reconnections between the coronal magnetic fields
tangled by the random motion of their foot points due to the convection of the
surface. Cargill (1994) succeeded in reproducing the emission measure
observation of skylab and pre-skylab by modelling the active-region closed
coronal loops as a bundle of numerous elemental loops that are randomly heated
by nanoflares. From the spectroscopic observations, Schmelz et al. (2001) and
Warren et al. (2008) reported that a coronal loop consists of plasma that has
a wide range of temperatures, which suggest that coronal loop is composed of
finer tubes. The volumetric filling factor of the loop is estimated as
approximately $10\%$ by Warren et al. (2008). Vekstein & Katsukawa (2000)
proposed the assumption that the cross-sectional area of each elemental loop
can be derived from the balance between the plasma pressure inside the loop
and the outer magnetic pressure. Based on this assumption, Sakamoto et al.
(2009) revealed that the major differences between the soft X-ray
($>2~{}\mathrm{MK}$) and EUV ($\simeq 1~{}\mathrm{MK}$) nanoflare-heated
coronal loops are their magnetic field strength ($40~{}\mathrm{G}$ for soft
X-ray and $8~{}\mathrm{G}$ for EUV) and released energies of nanoflares
($10^{24}~{}\textendash~{}10^{25}~{}\mathrm{erg}$ for soft X-ray and
$10^{23}~{}\mathrm{erg}$ for EUV). They also estimated the volumetric filling
factors of the loops as approximately $10\%$ and $70\%$ for soft X-ray and
EUV, respectively.

A microflare was observed by the hard X-ray balloon observation in 1980 for
the first time (Lin et al., 1984). After that, from the Soft X-ray Telescope
(SXT: Tsuneta et al., 1991) onboard the Yohkoh satellite (Ogawara et al.,
1991), it was revealed that many small explosive events occur in active
regions (Shimizu et al., 1992, 1994). To validate the possibility of the
coronal heating due to small-scale flares, it is necessary to derive the
occurrence frequency distribution of flares as a function of the energy. It is
known that the occurrence frequency of the flares is distributed as a power
law with the following equation:

$\frac{dN}{dE}=AE^{-\alpha}$ (1)

where $N$, $E$, and $\alpha$ are the event number, energy of each event, and
the power law constant, respectively (Hudson, 1991). From this equation, the
total energy released by all detected flares $P$ can be calculated as

$P=\int_{E_{\mathrm{min}}}^{E_{\mathrm{max}}}\frac{dN}{dE}EdE=\frac{A}{-\alpha+2}\left(E_{\mathrm{max}}^{-\alpha+2}-E_{\mathrm{min}}^{-\alpha+2}\right)$
(2)

where, $A$ is a constant. Therefore, small-scale flares significantly heat up
the corona when $E_{\mathrm{min}}$ is small enough and $\alpha$ is greater
than 2. However, in the case of solar and stellar flares ($10^{27}\lesssim
E\lesssim 10^{35}~{}\mathrm{erg}$), $\alpha$ is estimated as approximately
$1.8$ (Drake, 1971; Datlowe et al., 1974; Lin et al., 1984; Dennis, 1985;
Collura et al., 1988; Shakhovskaya, 1989), which does not support the
nanoflare heating model. On the other hand, the energy distributions of
smaller energy flares were analyzed using various observational equipments and
methods. In some studies, $\alpha$ was found to be greater than 2 (e.g.,
Parnell & Jupp, 2000; Benz & Krucker, 2002), which suggests that the frequency
distribution at smaller energies is steeper and therefore may heat the corona
sufficiently, while in other studies this is not the case (e.g., Shimizu,
1995; Aschwanden et al., 2000; Tajfirouze et al., 2016; Jess et al., 2019).

The smallest event detected by satellite observation is
$10^{25}~{}\mathrm{erg}$ (Testa et al., 2014). Current satellite observations
cannot detect extremely small-scale flare ($E\simeq 10^{23}~{}\mathrm{erg}$)
due to lack of sensitivity. On the other hand, a sounding rocket experiment,
the High-resolution Coronal Imager (Hi-C: Kobayashi et al., 2014), detected a
nanoflare that has $E\simeq 10^{23}~{}\mathrm{erg}$ at the foot point of a
coronal loop in an active region by EUV observation (Testa et al., 2013).
Moreover, another sounding rocket, the second flight of the Focusing Optics
Solar X-ray Imager (FOXSI-2: Christe et al., 2016) detected the plasma of the
upper 10 MK without any evident X-ray flare emissions (Ishikawa et al., 2017).
This observation indicates that the plasma is heated by nanoflares.

By using the 0-dimensional hydrodynamic coronal loop model EBTEL (Enthalpy
Based Thermal Evolution of Loops Klimchuk et al., 2008; Klimchuk, 2015;
Cargill et al., 2012) and probabilistic neural network, Tajfirouze et al.
(2016) estimated the power-law index $\alpha$ and the number of elemental
loops in one pixel of the EUV channel of Atmospheric Imaging Assembly (AIA:
Lemen et al., 2012) onboard the Solar Dynamics Observatory (SDO: Pesnell et
al., 2012). They concluded that $\alpha$ is approximately 1.5 and that 1000
loops are included in one pixel of an active region. However, they analyzed
only the case of $\alpha=1.5~{}\mathrm{and}~{}2.5$, and the range of flare
heating rate was very narrow (approximately from $0.01$ to
$0.5~{}\mathrm{erg~{}cm^{-3}~{}s^{-1}}$). Moreover, the location where
nanoflare occurs is neglected because they used the 0-dimensional model.
According to the model of Parker (1988), nanoflares may occur anywhere along
the coronal loop. The plasma evolution depends on the nanoflare location even
in one-dimensional coronal loop model (e.g.,  Reale & Orlando, 2008).
Therefore, to improve on this point, at least, a one-dimensional model must be
used for analysis. The modelling of a coronal loop has been done for a few
decades (e.g., Priest, 1978). Generally, it is assumed that the temperature
and density are almost the same as the surrounding corona, the magnetic
pressure is dominant ($\beta\ll 1$), the confined plasma is a compressible
fluid, and the energy is transported only along the magnetic field lines
(e.g., Rosner et al., 1978; Vesecky et al., 1979). Therefore, the evolution of
the coronal loop plasma can be described using a one-dimensional model only in
case neglecting the curvature, twisting, currents, waves, and non-uniform
cross-section (see review by Reale, 2014, for details).

In this paper, we derive the occurrence frequency distribution of small-scale
flares as a function of energy while considering the following:

  1. 1.

Evolution of coronal loop plasma heated by small-scale flares in the one-
dimensional model

  2. 2.

Elemental loops smaller than the spatial resolution of the observed image

  3. 3.

Multiwavelength observation

As far as we know, there are few studies that derive the frequency
distributions of small-scale flares with taking into account all the above
conditions. We succeeded in carrying out such an analysis by using one of the
machine learning techniques called genetic algorithm (GA), which is effective
at optimizing parameter combinations. This paper is organized as follows: The
active region observation, which we applied in our new method, is shown in
section 2. We describe the set of small-scale flare heating coronal loop
simulations in section 3. In section 4, we present how we derive the flare
parameter and energy distribution by using GA. Section 5 gives the obtained
results. Finally, we summarize this paper and discuss the validity of our
method and the nanoflare heating model in section 6.

##  2 Data and Observation

Figure 1: A snapshot of NOAA active region 12712 between 03:00 – 04:00 UT on
30 May 2018 obtained by the SDO/AIA 171 Å. The red crosses represent the
locations of loop foot points. The red line represents the vertical bisector
of both foot points.  Figure 2: Intensity distribution along the vertical
bisector of loop foot points (Red line in figure 1). Red and blue crosses
represent peaks and bottoms of large enhancements, respectively.

In this study, we use the coronal EUV images of NOAA active region 12712
between 03:00 and 04:00 UT on 30 May 2018, using six filters of the SDO/AIA,
which consists of 94 ($\log{T}\simeq 6.8$), 131 ($\log{T}\simeq 7.0$), 171
($\log{T}\simeq 5.8$), 193 ($\log{T}\simeq 6.0$ for quiet region and
$\log{T}\simeq 7.3$ for flaring region), 211 ($\log{T}\simeq 6.3$), and 335
($\log{T}\simeq 6.4$) Å (Lemen et al., 2012). This active region is suitable
for small-scale flare analysis because it is relatively calm and produces no
eruptions greater than a GOES B-class flare. The spatial and time resolution
of each filter is 0.6” and 12 s, respectively. The SDO/AIA data are calibrated
by the `aia_prep` routine in SolarSoftWare (SSW: Freeland & Handy, 1998). We
obtain light curves at the loop tops for our analysis. Figure 1 shows the
observed active region obtained from the AIA 171 Å filter. The red crosses and
red line indicates the foot points of coronal loops and their vertical
bisector, respectively. This red line corresponds to the observation
locations. The length and location of the red line is defined to avoid the
active region core and is extended to where the existence of the loop can be
visually confirmed.

Figure 2 shows the intensity distribution along the red line in figure 1. We
regard spikes in this profile as coronal loops. Each loop we detect is
enhancements above $3\sigma_{\mathrm{ph}}$, where $\sigma_{\mathrm{ph}}$
represents photon noise. In this study, we estimate $\sigma_{\mathrm{ph}}$ as
the square root of mean intensity along the red line. The red and blue crosses
in figure 2 represents the peak and bottom of each large spikes. We estimate
the minimum coronal loop radius as the minimum FWHM of these which is 5.2
pixels along the vertical bisector. Therefore, light curves are summed in each
$4\times 4$ pixels along the red line ($5.2/\sqrt{2}\simeq 4$).

Figure 3 presents a cropped map of figure 1. Each red square indicates pixels
used in the light curve. The size of each square is $4\times 4$ pixels as
described above. The number of the observation locations is 40.

The latitude of this active region is approximately 15 deg north. To track the
movement of this active region due to the rotation, we calculate the apparent
velocity as follows:

$v_{\mathrm{app}}=\frac{2\pi
R_{\odot}\cos{\phi}}{360}v_{\mathrm{rot}}~{}\mathrm{km~{}day^{-1}}$ (3)

where, $R_{\odot}$, $\phi$, and $v_{\mathrm{rot}}$ represent the solar radius,
latitude, and the rotation rate. According to the result of Snodgrass & Ulrich
(1990), the rotation rate of this active region is estimated as approximately
$14.54~{}\mathrm{deg~{}day^{-1}}$. Therefore, $v_{\mathrm{app}}$ can be
estimated as $1.7\times 10^{5}~{}\mathrm{km~{}day^{-1}}$ We assume that 1
arcsec along the solar-X at this active region is equivalent to 730 km. This
assumption is reasonable because the longitude of this active region is almost
zero. We calculate the apparent velocity of this active region as
$9.8~{}\mathrm{arcsec~{}h^{-1}}$. The apparent motion is corrected for by
moving the pixels in the light curve accordingly.

We obtain the light curves of the six filters in each red box. For example,
Figure 4 represents the observed light curves in one red box. Some
enhancements can be seen in the cooler channels (171, 193, 211 Å), while the
intensities of hotter channels (335, 94, 131 Å) are very low and noisy. The
loop lengths can be roughly estimated as
$200~{}\textendash~{}400~{}\mathrm{Mm}$.

Figure 3: Trimmed image of figure 1. Forty red squares indicate the area where
the light curves are obtained. Figure 4: An example of light curves of one
observational region (top-right red square in Figure 3) obtained from six
filters of SDO/AIA.

##  3 Numerical Simulation

We calculate the evolution of a coronal loop plasma heated by flares with a
one-dimensional model. In this study, we use CANS (coordinated astronomical
numerical software 111Details are available at the website
http://www.astro.phys.s.chiba-u.ac.jp/netlab/astro/index2-e.html. We used the
version of 21 September 2019 distributed at http://www-
space.pdf.s.u-tokyo.ac.jp/~yokoyama/etc/cans/index-e.html.) 1D solar flare
package. The fundamental equations are as follows:

$\displaystyle\frac{\partial}{\partial t}(\rho S)+\frac{\partial}{\partial
x}\left(\rho V_{x}S\right)=0$ (4) $\displaystyle\frac{\partial}{\partial
t}\left(\rho V_{x}S\right)+\frac{\partial}{\partial x}\left[\left(\rho
V_{x}^{2}+p\right)S\right]=\rho gS$ (5) $\displaystyle\frac{\partial}{\partial
t}\left[\left(\frac{p}{\gamma-1}+\frac{1}{2}\rho V_{x}^{2}\right)S\right]+$
$\displaystyle\frac{\partial}{\partial
x}\left[\left(\frac{\gamma}{\gamma-1}p+\frac{1}{2}\rho
V_{x}^{2}\right)V_{x}S-\kappa\frac{\partial T}{\partial x}S\right]=$
$\displaystyle\left(\rho gV_{x}+H-R+H_{f}\right)S$ (6) $\displaystyle
p=\frac{k_{B}}{m}\rho T$ (7)

where, $p$, $T$, $v_{x}$, $\rho$, $\gamma=5/3$, $S$, $g$, $H$, $R$, $H_{f}$,
$k_{B}$, $\kappa$, and $m$ represent the pressure, temperature, plasma
velocity along the loop, density, heat capacity ratio, cross-sectional area,
gravitational acceleration, static heating, radiative cooling, flare heating,
Boltzmann constant, thermal conductivity, and mean particle mass,
respectively. The simulation assumes that the length and cross section of the
loop do not change with time; the cross section is uniform along the loop; the
flow is inviscid and compressible; and the location where the flare occurs is
fixed at the loop top. The Spitzer thermal conductivity (Spitzer, 1956) and
gravity are taken into account as follows.

$\displaystyle\kappa=\kappa_{0}T^{5/2}$ (8) $\displaystyle
g=g_{0}\cos\left(\frac{\pi}{2L}x\right)$ (9)

where, $\kappa_{0}=10^{-6}~{}\mathrm{cgs}$, $g_{0}=270~{}\mathrm{m/s^{2}}$ and
$L$ represent thermal conductivity strength, a gravitational acceleration at
the photosphere, and the half loop length.

An approximation to correct for the effects of high-density plasma is included
in the radiative cooling model as following equations:

$\displaystyle R=\rho^{2}\Lambda_{\rho}(\rho)\Lambda(T)$ (10)
$\displaystyle\Lambda_{\rho}(\rho)=\frac{\rho_{\mathrm{cl}}}{\rho}\tanh\left(\frac{\rho}{\rho_{\mathrm{cl}}}\right)$
(11)

where, $\rho_{\mathrm{cl}}=10^{12}~{}\mathrm{cm^{-3}}$. The radiative cooling
function $\Lambda(T)$ is approximated by an algebraic function as follow:

$\displaystyle\Lambda(T)\approx\Lambda_{0}10^{\Theta(T)}$ (12)
$\displaystyle\theta\equiv\log_{10}\left(\frac{T}{T_{\mathrm{cl}}}\right)$
(13)
$\displaystyle\Theta(\theta)=0.4\theta-3+3\times\frac{2}{e^{1.5(\theta+0.08)}+e^{-2(\theta+0.08)}}$
(14)

where, $\Lambda_{0}=8\times 10^{-22}~{}\mathrm{cgs}$ and
$T_{\mathrm{cl}}=2\times 10^{5}~{}\mathrm{K}$.

This simulation includes not only the corona but also the transition region
and chromosphere. The flare energy input is represented by the following
equations:

$\displaystyle H_{f}=H_{f0}\cdot q(t)\cdot f(x)\cdot g_{1}(x)\cdot g_{2}(x)$
(15) $\displaystyle
q(t)=\frac{1}{4}\left\\{1+\tanh{\frac{t-t_{i0}}{0.1\tau_{0}}}\right\\}\left\\{1-\tanh{\frac{t-t_{i1}}{0.1\tau_{0}}}\right\\}$
(16) $\displaystyle
f(x)=\frac{1}{\sqrt{2\pi}}\exp{\left[-\frac{(x-x_{i})^{2}}{2w_{f}^{2}}\right]}$
(17) $\displaystyle
g_{1}(x)=\frac{1}{2}\left\\{1+\tanh{\frac{x-20\mathcal{H}_{0}}{3\mathcal{H}_{0}}}\right\\}$
(18) $\displaystyle
g_{2}(x)=\frac{1}{2}\left\\{1-\tanh{\frac{x-(2L-20\mathcal{H}_{0})}{3\mathcal{H}_{0}}}\right\\}$
(19)

where, $w_{f}=6000~{}\mathrm{km}$, $\mathcal{H}_{0}=200~{}\mathrm{km}$, and
$\tau_{0}=20~{}\mathrm{s}$ represent the width along the loop of the flare,
the scale height, and sound wave traveling time at the surface ($x=0$),
respectively. $t_{i0}$, $t_{i1}$, and $x_{i}$ represent the beginning and
finishing time and location of heating by the $i$th flare, respectively.
$q(t)$ is a function of time to make the heating impulsive. The role of
$g_{1}(x)$ and $g_{2}(x)$ is to prevent the heat pulse from entering directly
into the chromosphere. This simulation uses the modified Lax–Wendroff scheme
which is second-order accurate in both space and time. The boundary conditions
at $x=0$ and $x=2L$ are as follows:

$\displaystyle\frac{\partial\rho}{\partial x}=0$ (20)
$\displaystyle\frac{\partial p}{\partial x}=0$ (21) $\displaystyle V_{x}=0$
(22)

This simulation setup is mostly the same as that of Hori et al. (1997); Imada
& Zweibel (2012); Kawai et al. (2020).

Table 1: Free parameters in the simulation and their ranges Parameter | Minimum | Maximum
---|---|---
Heating rate $[\mathrm{erg/s/cm^{3}}]$ | $10^{-2}$ | $10^{2}$
Flare duration [$\mathrm{s}$] | 1 | 300
Flare occurrence time [$\mathrm{s}$] | $t=0$ | $t=5600$
Flare location | $x=0$ | $x=2L$
Loop length [$\mathrm{Mm}$] | 100 | 500

Figure 5: Results of the hydrodynamic simulation. Each panel presents the time
evolution of the temperature distribution (left-top), density (right-top),
pressure (left-bottom), and plasma velocity along the loop (right-bottom) with
30-s cadence. Each horizontal axis represents the loop coordinates from the
surface to the loop top. Line color indicates the progress of time (from blue
to red).

To demonstrate the simulation, we show a simple example in figure 5. Each
panel presents the temporal variation and spatial distribution of temperature,
density, pressure, and plasma velocity along the loop, respectively. As each
horizontal axis indicates coordinates along the loop, the left-hand edge
($x=0$) is the surface and the right-hand edge is the loop top. The region
where the temperature and pressure change rapidly ($x\simeq 0.3$ Mm) is the
transition region, and the chromosphere is in the left of this region. Line
color indicates the progress of time in the simulation, from blue to red. This
result is for a single flare occurrence at the loop top at the beginning of
the simulation. When a flare occurs, the temperature of the loop top is
increased and is transported to the foot points of the loop by thermal
conduction. Then, the temperature and pressure of the chromosphere are rapidly
increased by the incoming high-temperature plasma. As a result, high-density
plasma in the chromosphere is ejected into the corona by the pressure gradient
force; this is referred to as chromospheric evaporation. Consequently, the
coronal loop is filled with high-density plasma, and emits soft X-ray and EUV
irradiance.

To use various results in GA later, we carry out the simulation with random
flare heating rates, heating duration, occurrence time, occurrence location,
number, and loop length. Table 1 represents the parameters and their possible
ranges. According to the study of Testa et al. (2014), an energy flux of a
nanoflare ($10^{24}\lesssim E\lesssim 10^{25}~{}\mathrm{erg}$) is
approximately $7\times 10^{7}$ – $2\times 10^{9}~{}\mathrm{erg/s/cm^{2}}$.
Therefore, the range of flare heating rate in the simulation is roughly
$10^{-0.9}$ – $10^{0.5}~{}\mathrm{erg/s/cm^{3}}$ because the width along the
loop of flare is $6000~{}\mathrm{km}$. However, to include weaker and stronger
heating, we defined the range of heating rate as $10^{-2}$ –
$10^{2}~{}\mathrm{erg/s/cm^{3}}$. The heating duration of a nanoflare in Testa
et al. (2014) is from 10 to 30 s. We define the range of heating duration in
our simulation from 1 (10 times shorter than the minimum) to 300 (10 times
longer than the maximum) s. Testa et al. (2014) also suggested that the
duration of heating event is similar to that of transition region
brightenings. The heating duration range in our simulation almost covers the
typical lifetimes of compact transition region brightenings, that is from a
few tens of seconds to several minutes (Tiwari et al., 2019). The probability
distributions of heating rate and flare duration are uniform in a logarithmic
scale while those of flare occurrence time, location, and loop length are
uniform in their ranges. This maximum number of flares is roughly determined
to distinguish each heating event in each run even when flares have the
longest duration. The wider the parameter range, the more various the
simulation results, but the more calculations are required to reproduce the
observed light curves well. These parameter ranges are wider than those of
Tajfirouze et al. (2016) and enough to reproduce the observations because we
focus on relatively calm coronal loops. The cadence of the simulation output
is 2 s.

Figure 6: A schematic architecture of a gene used in this paper. Each gene
represents a synthetic light curve $I_{\mathrm{synth}}(t,\lambda)$, which is
described as a linear superposition of $l$ of simulated light curves.

The original CANS solar flare package calculates only half of the loop (foot
point to top) under the symmetry assumption as shown in figure 5. his
assumption should work well when a flare occurs at the top like in figure 5.
However, according to the Parker’s model, nanoflares can occur everywhere
along the loop and this symmetric assumption might be violated. Therefore, in
this study, we extend it to calculate the whole loop (one foot point to
another foot point) to consider various heating locations. Moreover, to
simplify the discussion, we remove static heating ($H=0$) and assume the loop
is heated by only flares. There are two reasons for this assumption. One is
that the mechanism and intensity of background heating is also unclear,
therefore, another assumption is needed anyway. The second reason is that the
contribution of background heating must be relatively small even when the loop
is heated by small-scale flares.

In some conditions, the calculation is broken due to an overheating or
overcooling of the loop. We do not use simulation results, which has density
of over $10^{12}~{}\mathrm{cm^{-3}}$ around the loop top. This high density
loop is sometimes caused by overcooling due to the imbalance between the flare
heating and radiative and conductive cooling. This imbalance mostly occurs
when the flare heating is too large.

Because the observation duration is one hour, as described in Section 2, the
simulation also needs to be calculated for one hour. However, to make various
initial conditions, we calculate time in advance of observation duration. As
can be seen in the simulation results, the variety of light curve becomes the
maximum around $t=2000~{}\mathrm{s}$ even for the longest loops. Therefore, we
calculate for 5600 s and remove the first 2000 s for the optimization. As
describe above, the number of flares are 15 in each run, however, some or all
of them can occur during the first 2000 s. Therefore, the range of number of
flares is from 0 to 15 during the latter one hour calculation. In this study,
we ran approximately 5,000 random simulations.

As described in section 2, we obtain a light curve of the SDO/AIA $4\times 4$
pixels ($\simeq 1.8\times 1.8~{}\mathrm{Mm^{2}}$) at the loop top by the
observation. Therefore, we derive the SDO/AIA light curve in the
$1.8~{}\mathrm{Mm}$ range at the loop top of each calculation as well. Using
the temperature response function of each SDO/AIA channel (Lemen et al.,
2012), we obtain the intensities in each time and grid from the simulation
result. The response functions are available from the SolarSoftWare (SSW:
Freeland & Handy, 1998) procedure `aia_get_response.pro`. We assume the
emission measure to be $\rho^{2}$ and derive the light curve per cross-
sectional area ($\mathrm{DN~{}s^{-1}~{}cm^{-2}}$) for each calculation. We
assign the identification to each run to use them for optimization, which is
described in the next section.

##  4 Genetic Algorithm

We estimate the best combinations of the results of simulations that reproduce
actually observed coronal loop light curves by using GA. GA is a machine
learning technique that is effective at solving optimization and searching
problems and is based on the concept of Darwin’s theory of evolution. “Genes”
and an “environment” are necessary for GA. GA can search the gene that is
optimized best for the environment through natural selection, crossover, and
mutation. In this study, each gene is a list of the simulation IDs that are
assigned in section 3. Each gene expresses the synthetic light curve by linear
superposition of its simulation results. Therefore, a synthetic light curve
$I_{\mathrm{synth}}^{X}$ of a gene $X$ can be described as follows:

$I_{\mathrm{synth}}^{X}(t,\lambda)=\sum_{i}^{l}I_{\mathrm{sim}}^{k_{i}}(t,\lambda)$
(23)

where, $l$ and $I_{\mathrm{sim}}^{k}$ represent the length of gene (the number
of simulated loops to combine) and simulated light curve of simulation ID $k$,
Figure 6 represents the schematic organization of the genes. In this paper,
gene length $l$ and the number of gene $n$ are 100 and $500$, respectively.
$l$ is defined based on the analysis of Tajfirouze et al. (2016), which
suggests that the number of loops in each SDO/AIA pixel is approximately 1000.
However, only one flare occurs on each coronal loop in their study while at
most 15 flares occur on each loop in our simulation. Therefore, we roughly
estimate the number of elemental loops in each red square in figure 3 as 100.
The accuracy of reproduction of light curves increases as the number of genes
$n$ increases. We define $n=500$ based on available computational resources.
However, even if $n$ is smaller, the reproducibility can be improved by
increasing the number of iteration.

Figure 7: A schematic flow of GA. Firstly, initial genes generated randomly
are evaluated based on the correlations between observed
$I_{\mathrm{obs}}(t,\lambda)$ and synthetic $I_{\mathrm{synth}}(t,\lambda)$
light curves. Second, based on the evaluation, better genes are left or
duplicated probabilistically, while worse genes are eliminated. Third, the
remaining genes swap their information partially with each other. Sometimes a
mutation occurs to avoid a local minimum. The best gene that reproduces the
observed light curves can be estimated by continuing this iteration.

The main components of GA are a natural selection, crossover, and mutation.
Figure 7 schematically shows the procedures. Genes have a random list of IDs
at the beginning. First, in natural selection, while better genes are left or
multiplied, worse genes are eliminated stochastically. Each gene is evaluated
by the averaged correlation between observed and synthetic light curves in six
AIA filters. Therefore, the selection probability of gene $i$ can be described
as follows:

$p_{i}=\frac{c_{i}}{\sum_{k=1}^{n}c_{k}}$ (24)

where, $c_{i}$ represents the averaged correlation of gene $i$. Genes in the
next generation are decided by this probability. However, in this case, even
the best gene among the generation has a possibility to be eliminated.
Therefore, we leave at least one of the genes which have the highest
correlation to the next generation, which is typical for GA. This can
accelerate the learning while it may increase the likelihood of getting a
local solution due to the lack of variety of genes.

Second, in crossover, we search for new better solutions by swapping
information partially between the pair of the genes. In this study, we choose
two-point crossover for this procedure. We randomly choose pairs of genes and
two crossover points (black vertical lines in figure 7) from within the genes.
The IDs between the crossover points are exchanged with each other. These new
genes are populations for the next generation. $90~{}\%$ of all genes are
crossed over in each generation.

Third, some genes are mutated to maintain diversity and avoid local minima. In
a mutation, a few IDs in the genes are randomly replaced to the another ID
with a low probability. In this paper, $1~{}\%$ of all IDs in the genes are
mutated in each generation.

We estimate the gene that can best reproduce the observed light curves by
iterating these procedures 1,000 times for each observational area. The gene
is a list of simulation results and we know the inputs of flare parameters of
each simulation. Therefore, we can estimate the combination of flare
parameters that reproduce observational results best by GA.

When synthetic light curves $I_{\mathrm{synth}}$ are satisfactorily correlated
to observational results $I_{\mathrm{obs}}$, the relationship between them can
be written as follows:

$I_{\mathrm{obs}}(t,\lambda)\simeq A(\lambda)\times
I_{\mathrm{synth}}(t,\lambda)+B(\lambda)$ (25)

where, $A$ and $B$ are coefficients derived by least squares. $A$ and $B$
indicate the ratio of amplitudes of fluctuations between $I_{\mathrm{obs}}$
and $I_{\mathrm{synth}}$ and the background component of $I_{\mathrm{obs}}$,
respectively. As mentioned above, $I_{\mathrm{sim}}$ is the time series of
emissions from the coronal loop whose cross-sectional area is
$1~{}\mathrm{cm^{2}}$. Therefore, the cross-sectional area of elemental
coronal loops in each observation region can be estimated as $\Delta
S_{\mathrm{loop}}(\lambda)=A(\lambda)~{}\mathrm{[cm^{2}]}$.

##  5 Result

At first, to verify the robustness of our procedure, we test whether the GA
can estimate the correct combination of simulations from light curves
generated by random set of genes. We create SDO/AIA synthetic light curves
from 100 sets of simulation results as a target. The created light curves have
a random gaussian noise whose standard deviation is $10\%$ of the intensity.
The estimated genes are not exactly the same with the target gene, however,
the derived distribution of flares are very similar to the input of that.
Figure 8 represents an example of occurrence frequency distribution of flares
in unit of cross sectional area. Red and black histograms represent
distributions of GA reproduction and input of target simulations. In this
case, the mean correlation between target and reproduced light curves is
approximately 0.97. We continue this test 25 times to define the uncertainty
of this method. We define the uncertainty of the occurrence frequency
distribution, $\sigma_{f}$, as a standard deviation of the difference between
target and reproduced one in each energy bin.

Figure 8: An example of frequency distribution of flares as a function of
energy. Black and red lines indicate the distribution of input of the target
simulations and of estimation by the GA, respectively. Figure 9: Averaged
correlations between observed and the best synthetic light curves in each
observation region. The horizontal axis indicates the observation location
(from outside to inside). Figure 10: Sample of observed SDO/AIA light curves
(solid lines) and synthetic ones which are optimized to reproduce the
observation the most by GA (red dashed lines), which has the best average
correlation ($\simeq 0.9$).  Figure 11: Filling factors of heated coronal loop
in each observation location and filter. Horizontal axis represents the peak
temperature of response function of each SDO/AIA channel. Colder (Warmer)
color indicates farther distance from (closer to) the core.

Figure 9 represents the average correlations between observed and optimized
synthetic light curves. The horizontal axis indicates observation locations
(red squares in figure 1), from outside to inside. The correlations are
approximately $0.4$ – $0.9$, though they depend on the regions.

Figure 10 shows an example of observed light curves obtained from SDO/AIA six
channels (solid lines) and synthetic ones estimated by GA, which reproduces
the observation best (dashed lines). These light curves are obtained from the
region where the average correlation is the best ($\simeq 0.9$). Large
fluctuations are reproduced accurately, however, smaller fluctuations other
than photon noise are not reproduced well.

Figure 12: Time series of energy flux of flares in each observation region
derived by GA. Colder (Warmer) color indicates farther distance from (closer
to) the core. The panel without lines indicates that the GA determines that
the calculations without flares from $t=2000$ to $t=5600$ reproduce the
observed light curves best.

We define a cross sectional area along the loop of each observation region
(red square in figure 1) as $S_{\mathrm{obs}}=w^{2}$, where $w$ represents the
width of each red square. In this study, the size of each red square is
$4\times 4$ pixels, and hence, $w\simeq 1.75~{}\mathrm{Mm}$ and
$S_{\mathrm{obs}}\simeq 3.1~{}\mathrm{Mm^{2}}$. We define the filling factor
of the coronal loop $\phi$ as follows:

$\phi(\lambda)=\frac{l\times\Delta
S_{\mathrm{loop}}(\lambda)}{S_{\mathrm{obs}}}$ (26)

where, $l$ represents the length of each gene. Figure 11 presents the
estimated filling factors of all observation areas for each AIA filter. The
warmer color indicates observation locations which are closer to the core. The
horizontal axis represents the peak temperature of the response function of
each AIA channel. As a result, the higher the plasma temperature, the lower
the filling factor. In addition, the location closer to the core tends to have
larger filling factor.

The energy flux by flares can be described as:

$F(t)=\frac{1}{l}\sum_{i}^{l}H_{f}^{k_{i}}(t)~{}\mathrm{[erg/s/cm^{2}]}$ (27)

where $H_{f}^{k}$ is the flux caused by flares in simulation ID $k$. Each
panel of figure 12 represents the estimated time series of energy flux by
flares in each observation area. The color indicates the observation location
as same as figure 11. The panel without lines shows that the GA determines
that the calculations without flares from $t=2000$ to $t=5600$ are the best to
reproduce the observed light curves. Almost all observed regions are heated by
flares intermittently, which have energy flux above the typical requirement
($10^{7}~{}\mathrm{erg/s/cm^{2}}$: Withbroe & Noyes, 1977).

Figure 13 presents the frequency distribution of detected flares as a function
of energy. Vertical lines indicate uncertainties defined by the test described
in the beginning of this section as follow:

$\sigma_{f}^{\prime}=N_{\mathrm{obs}}\times\Delta
S_{\mathrm{loop}}(94)\times\sigma_{f}.$ (28)

where, $N_{\mathrm{obs}}$ represents the number of observation regions (40 in
this study). The volume of the heated loop is necessary for the energy
estimation; however, it depends on the filter. Therefore, in this study, we
estimate the energy with $\Delta S_{\mathrm{loop}}(94)$, which is the most
similar to the X-ray observational study of Shimizu (1995). The number of all
detected flares is approximately $37,000$. Generally, the occurrence frequency
distribution of flares as a function of energy show the power-law. However,
the result does not have a single power-law index through the energy range.
Therefore, we calculate the power-law index of the distribution in each energy
bin as shown in Figure 14. The size of each bin is defined as the width by
which the flare energy on a logarithmic scale changes by one. Horizontal axis
indicates median of flare energy of each bin. In most cases, the power-law
index increases with flare energy. The power-law index is greater than 2 in
the energy range from $10^{26}$ to $10^{27}~{}\mathrm{erg}$ and above. On the
other hand, the index is approximately 1 in the energy range smaller than
$10^{25.5}~{}\mathrm{erg}$. This energy range is probably lower than the
detection limit of our method because the index of simulation input is also
approximately $1$.

The blue line (left axis) in figure 15 represents the contribution of flares
to heat the corona in each energy range. The red line (right axis) represents
the cumulative contribution from high energy. Vertical lines indicate the
errors calculated from $\sigma_{f}$ defined in the beginning of this section.
Flares that contribute the most to the heating of the corona are those in the
energy range of $10^{26}\lesssim E\lesssim 10^{27}~{}\mathrm{erg}$. Moreover,
the $90\%$ of the coronal heating is done by flares that have energies greater
than $10^{25}~{}\mathrm{erg}$.

Figure 13: Frequency distribution of flares as a function of energy estimated
by GA. vertical lines indicate errors which defined by the test described in
the beginning of Section 5. The number of all detected flares is 37,458.
Figure 14: Series of power-law indices of the occurrence frequency
distribution shown in Figure 13 in each energy bin. The size of each bin is
the width that the flare energy on a logarithmic scale changes one. Dashed
horizontal line indicates where power-law index becomes 2 which implies
whether smaller flares are dominant in the coronal heating or not. Figure 15:
Contribution of flares to the coronal heating in each energy range (left axis,
blue) and the cumulative contribution from high energy (right axis, red).
Vertical lines indicate the errors calculated from the definition of the
standard deviation $\sigma_{f}$ described in the beginning of Section 5.
Figure 16: Schematic picture of coronal loops in the observation region.
$\Delta S_{\mathrm{loop}}(T_{H})$ and $\Delta S_{\mathrm{loop}}(T_{L})$
represents the cross sectional areas of loops composed of hotter and cooler
plasmas, respectively. Red square represents the observation region depicted
in Figure 1.

##  6 Discussion and Summary

In this paper, we introduced a new method for detection and energy estimation
of small-scale flares by using one-dimensional simulation and GA, which is a
machine learning technique. We applied our method to the active-region coronal
loop observation by SDO/AIA and obtained the time series of energy flux of
flares, occurrence frequency distribution as a function of energy, and filling
factors.

The occurrence frequency distribution generally can be fitted by the power-law
distribution, however, the power-law index of the derived distribution depends
on the fitting range. As a result, the power-law index is found to be greater
than 2 in flare energy range of $10^{26}<E<10^{28}~{}\mathrm{erg}$ (Figure
14). In this energy range, smaller flares are dominant in heating the corona,
which is in line with the nanoflare heating model. On the other hand, the
power-law index is approximately 1 in the energy range of smaller than
$10^{25.5}~{}\mathrm{erg}$. This is probably caused by the detection
limitation of our method because the distribution is similar to that of
simulation inputs. Moreover, we found that the coronal loops are heated by
flares that have enough energy flux to heat the corona intermittently (Figure
12). This is mainly because we calculate the fine structure of the loop that
cannot be resolved. From figure 11, we found that the volume that is heated by
flares is much smaller than the observational resolution though it depends on
the temperature and region. On the other hand, we found that the flares in
energy range of $10^{26}\lesssim E\lesssim 10^{27}~{}\mathrm{erg}$ contribute
to the coronal heating the most. In addition, $90\%$ of the energy flux comes
from flares that release energies greater than $10^{25}~{}\mathrm{erg}$. It
has been thought that smaller flares
($10^{23}~{}\textendash~{}10^{24}~{}\mathrm{erg}$) contribute the most to heat
the corona because the energy flux of detected flares are not enough. There
are some reasons why our result is incompatible with previous studies. We
calculate the plasma evolution of coronal loops taking into account smaller
loops than the observational resolution unlike previous studies. It is
possible that the flare energy and occurrence frequency could be derived more
accurately than previous studies by considering such a process. This is
because our energy estimation is based on physics unlike the method that
defines a flare energy as a difference of thermal energy of a loop (e.g.,
Shimizu, 1995). However, our method has some shortcomings. First, our results
do not explain the balance of energy on the whole active region because we
focus only on the coronal loop and not on the core region. Second, as we
mentioned above, our method might be unable to detect small-scale flares which
have an energy less than approximately $10^{24}~{}\mathrm{erg}$. Third, our
combined synthetic light curves reproduce large fluctuations in the
observation while smaller fluctuations are not reproduced well.

We found that the filling factor increases as we get closer to the core, and
the variation between the minimum and maximum is approximately three orders of
magnitude. This is probably because the loops closer to the core can
experience more magnetic pressure from the outer loops. This tendency is
acceptable because the loops closer to the core are brighter than those in the
outer region. As shown in figure 11, the higher the peak temperature of each
SDO/AIA filter, the lower the filling factor of coronal loops. This tendency
is in agreement with the study of Sakamoto et al. (2009), which compared the
filling factors of the soft X-ray and EUV loops. This temperature dependence
implies that the hotter coronal loops have a sparser structure than the cooler
ones. Figure 16 presents the schematic picture of coronal loops in each
observation region. $\Delta S_{\mathrm{loop}}(T_{H})$ (red circle) and $\Delta
S_{\mathrm{loop}}(T_{L})$ (gray circle) represents the cross-sectional areas
of each loop composed of hotter and cooler plasmas, respectively. The red
square corresponds to that depicted in Figure 1.

The largest filling factor is greater than 100%. This can be caused by the
assumption that the observational region is regarded as a cube. Coronal loops
might be distributed along the line-of-sight direction deeper than
observational width $w$. On the other hand, in some cases, the filling factor
is very small, approximately 0.01%. Assuming $l$ of the elemental loops exist
in each observation region, the half radius of each elemental loop $d$ can be
described as:

$d=\sqrt{\frac{S_{\mathrm{obs}}\phi}{\pi l}}.$ (29)

Therefore, in case of $0.01\%<\phi<100\%$,

$10\lesssim d\lesssim 990~{}\mathrm{[km]}.$ (30)

This result roughly corresponds to the maximum half radius of the elemental
loop estimated by sounding rocket experiment, which is approximately
$15~{}\mathrm{km}$ (Peter et al., 2013). $10~{}\mathrm{km}$ is roughly
equivalent to $0.014~{}\mathrm{{}^{\prime\prime}}$ at the solar surface, which
is about 2.3 % of the resolution of SDO/AIA.

To make large amounts of simulation results for GA, we use the one-dimensional
coronal loop model which has a low computational cost. As mentioned in section
1, the evolution of the coronal loop can be described in 1D when ignoring some
processes. Strictly speaking, flares are magnetic reconnections between the
loops that cannot be described using the one-dimensional model. Instead of
this process, we directly substitute the flare heating term into the energy
equation. However, making thousands of simulation results by using
multidimensional model is unrealistic for current computational performance.

We use CANS solar flare package to calculate a coronal loop heated by small-
scale flares. Therefore we neglect some important processes in small-scale
flare analysis such as non-equilibrium of ionization (NEI), radiative
transfer, and non-thermal electron beams. Reale & Orlando (2008); Imada et al.
(2011, 2015) suggest that NEI should have effects on the emission from the
plasma which is impulsively heated by flares. Using the HYDRAD model (e.g.,
Bradshaw & Cargill, 2010), Bradshaw & Klimchuk (2011) studies a NEI effect on
the coronal loop emissions for the case that a loop is heated by nanoflares.
Because our simulation neglects NEI unlike HYDRAD, energy of small-scale
flares may be estimated smaller than the actual and some small-scale events
cannot be distinguished. Including this process to the simulation should make
the analysis more accurate. Polito et al. (2018) suggest that the Si IV
Doppler velocity at the chromosphere is a good indicator for detection of
nanoflares ($\simeq 10^{24}~{}\mathrm{erg}$) from the analysis based on the
RADYN model (e.g.,  Allred et al., 2015). It is necessary to validate the
result by analyzing both coronal loops and chromospheres. Testa et al. (2014)
also use the RADYN model and reveal that an observation of small and rapid
variability of intensity and velocity at the loop foot point is consistent
with heating by non-thermal electron beams generated by nanoflares. These
processes should also be included in our analysis, however.

We must derive the combination of simulations which reproduce the observations
the most. Therefore, we use GA, which is used to solve optimization and
searching problems such as the traveling salesman problem. GA is suitable for
statistical analysis, however, it has initial value dependence and may find
only a local minimum. Moreover, it is difficult to prove that the combined
synthetic light curve has only one combination because all combinations
(${}_{5000}\mathrm{C}_{100}>10^{200}$) must be calculated to prove that. In
addition, GA has some free parameters such as length of each gene and the
number of genes and iterations. The larger number of genes and iterations can
make the correlation better, and thus, they should be as large as possible. On
the other hand, the length of each gene defines the number of elemental loops
in each observation region. In this paper, we assume that the number of loops
in each macro pixel is 100 and try to estimate the best combination of
simulations under this condition. However, the number of elemental loops in
each identifiable loop is high and under discussion (e.g., Polito et al.,
2018; Tajfirouze et al., 2016). Therefore, it is necessary to analyze flares
in various cases of the length of each gene.

The authors thank K. Kusano for fruitful discussions. We would like to thank
T. Yokoyama for his efforts in developing CANS package. This work was
partially supported by the Grant-in-Aid for 17K14401 and 15H05816, and the
Program for Leading Graduate Schools, “PhD Professional: Gateway to Success in
Frontier Asia” by the Ministry of Education, Culture, Sports, Science and
Technology. The Solar Dynamics Observatory is a part of NASA’s Living with a
Star program. A part of this study was carried out by using the computational
resource of the Center for Integrated Data Science, Institute for Space-Earth
Environmental Research, Nagoya University.

## References

  * Allred et al. (2015) Allred, J. C., Kowalski, A. F., & Carlsson, M. 2015, ApJ, 809, 104, doi: 10.1088/0004-637X/809/1/104
  * Antolin & Shibata (2010) Antolin, P., & Shibata, K. 2010, ApJ, 712, 494, doi: 10.1088/0004-637X/712/1/494
  * Aschwanden et al. (2000) Aschwanden, M. J., Tarbell, T. D., Nightingale, R. W., et al. 2000, ApJ, 535, 1047, doi: 10.1086/308867
  * Benz & Krucker (2002) Benz, A. O., & Krucker, S. 2002, ApJ, 568, 413, doi: 10.1086/338807
  * Bradshaw & Cargill (2010) Bradshaw, S. J., & Cargill, P. J. 2010, ApJ, 717, 163, doi: 10.1088/0004-637X/717/1/163
  * Bradshaw & Klimchuk (2011) Bradshaw, S. J., & Klimchuk, J. A. 2011, ApJS, 194, 26, doi: 10.1088/0067-0049/194/2/26
  * Cargill (1994) Cargill, P. J. 1994, ApJ, 422, 381, doi: 10.1086/173733
  * Cargill et al. (2012) Cargill, P. J., Bradshaw, S. J., & Klimchuk, J. A. 2012, ApJ, 752, 161, doi: 10.1088/0004-637X/752/2/161
  * Christe et al. (2016) Christe, S., Glesener, L., Buitrago-Casas, C., et al. 2016, Journal of Astronomical Instrumentation, 5, 1640005, doi: 10.1142/S2251171716400055
  * Collura et al. (1988) Collura, A., Pasquini, L., & Schmitt, J. H. M. M. 1988, A&A, 205, 197
  * Datlowe et al. (1974) Datlowe, D. W., Elcan, M. J., & Hudson, H. S. 1974, Sol. Phys., 39, 155, doi: 10.1007/BF00154978
  * Dennis (1985) Dennis, B. R. 1985, Sol. Phys., 100, 465, doi: 10.1007/BF00158441
  * Drake (1971) Drake, J. F. 1971, Sol. Phys., 16, 152, doi: 10.1007/BF00154510
  * Freeland & Handy (1998) Freeland, S. L., & Handy, B. N. 1998, Sol. Phys., 182, 497, doi: 10.1023/A:1005038224881
  * Hori et al. (1997) Hori, K., Yokoyama, T., Kosugi, T., & Shibata, K. 1997, ApJ, 489, 426, doi: 10.1086/304754
  * Hudson (1991) Hudson, H. S. 1991, Sol. Phys., 133, 357, doi: 10.1007/BF00149894
  * Imada et al. (2015) Imada, S., Murakami, I., & Watanabe, T. 2015, Physics of Plasmas, 22, 101206, doi: 10.1063/1.4932335
  * Imada et al. (2011) Imada, S., Murakami, I., Watanabe, T., Hara, H., & Shimizu, T. 2011, ApJ, 742, 70, doi: 10.1088/0004-637X/742/2/70
  * Imada & Zweibel (2012) Imada, S., & Zweibel, E. G. 2012, ApJ, 755, 93, doi: 10.1088/0004-637X/755/2/93
  * Ishikawa et al. (2017) Ishikawa, S., Glesener, L., Krucker, S., et al. 2017, Nature Astronomy, 1, 771, doi: 10.1038/s41550-017-0269-z
  * Jess et al. (2019) Jess, D. B., Dillon, C. J., Kirk, M. S., et al. 2019, ApJ, 871, 133, doi: 10.3847/1538-4357/aaf8ae
  * Kawai et al. (2020) Kawai, T., Imada, S., Nishimoto, S., Watanabe, K., & Kawate, T. 2020, Journal of Atmospheric and Solar-Terrestrial Physics, 205, 105302, doi: 10.1016/j.jastp.2020.105302
  * Klimchuk (2006) Klimchuk, J. A. 2006, Sol. Phys., 234, 41, doi: 10.1007/s11207-006-0055-z
  * Klimchuk (2015) —. 2015, Philosophical Transactions of the Royal Society of London Series A, 373, 20140256, doi: 10.1098/rsta.2014.0256
  * Klimchuk et al. (2008) Klimchuk, J. A., Patsourakos, S., & Cargill, P. J. 2008, ApJ, 682, 1351, doi: 10.1086/589426
  * Kobayashi et al. (2014) Kobayashi, K., Cirtain, J., Winebarger, A. R., et al. 2014, Sol. Phys., 289, 4393, doi: 10.1007/s11207-014-0544-4
  * Lemen et al. (2012) Lemen, J. R., Title, A. M., Akin, D. J., et al. 2012, Sol. Phys., 275, 17, doi: 10.1007/s11207-011-9776-8
  * Lin et al. (1984) Lin, R. P., Schwartz, R. A., Kane, S. R., Pelling, R. M., & Hurley, K. C. 1984, ApJ, 283, 421, doi: 10.1086/162321
  * Ogawara et al. (1991) Ogawara, Y., Takano, T., Kato, T., et al. 1991, Sol. Phys., 136, 1, doi: 10.1007/BF00151692
  * Parker (1983) Parker, E. N. 1983, ApJ, 264, 635, doi: 10.1086/160636
  * Parker (1988) —. 1988, ApJ, 330, 474, doi: 10.1086/166485
  * Parnell & Jupp (2000) Parnell, C. E., & Jupp, P. E. 2000, ApJ, 529, 554, doi: 10.1086/308271
  * Pesnell et al. (2012) Pesnell, W. D., Thompson, B. J., & Chamberlin, P. C. 2012, Sol. Phys., 275, 3, doi: 10.1007/s11207-011-9841-3
  * Peter et al. (2013) Peter, H., Bingert, S., Klimchuk, J. A., et al. 2013, A&A, 556, A104, doi: 10.1051/0004-6361/201321826
  * Polito et al. (2018) Polito, V., Testa, P., Allred, J., et al. 2018, ApJ, 856, 178, doi: 10.3847/1538-4357/aab49e
  * Priest (1978) Priest, E. R. 1978, Sol. Phys., 58, 57, doi: 10.1007/BF00152555
  * Reale (2014) Reale, F. 2014, Living Reviews in Solar Physics, 11, 4, doi: 10.12942/lrsp-2014-4
  * Reale & Orlando (2008) Reale, F., & Orlando, S. 2008, ApJ, 684, 715, doi: 10.1086/590338
  * Rosner et al. (1978) Rosner, R., Tucker, W. H., & Vaiana, G. S. 1978, ApJ, 220, 643, doi: 10.1086/155949
  * Sakamoto et al. (2009) Sakamoto, Y., Tsuneta, S., & Vekstein, G. 2009, ApJ, 703, 2118, doi: 10.1088/0004-637X/703/2/2118
  * Schmelz et al. (2001) Schmelz, J. T., Scopes, R. T., Cirtain, J. W., Winter, H. D., & Allen, J. D. 2001, ApJ, 556, 896, doi: 10.1086/321588
  * Shakhovskaya (1989) Shakhovskaya, N. I. 1989, Sol. Phys., 121, 375, doi: 10.1007/BF00161707
  * Shimizu (1995) Shimizu, T. 1995, PASJ, 47, 251
  * Shimizu et al. (1994) Shimizu, T., Tsuneta, S., Acton, L. W., et al. 1994, ApJ, 422, 906, doi: 10.1086/173782
  * Shimizu et al. (1992) Shimizu, T., Tsuneta, S., Acton, L. W., Lemen, J. R., & Uchida, Y. 1992, PASJ, 44, L147
  * Snodgrass & Ulrich (1990) Snodgrass, H. B., & Ulrich, R. K. 1990, ApJ, 351, 309, doi: 10.1086/168467
  * Spitzer (1956) Spitzer, L. 1956, Physics of Fully Ionized Gases
  * Tajfirouze et al. (2016) Tajfirouze, E., Reale, F., Petralia, A., & Testa, P. 2016, ApJ, 816, 12, doi: 10.3847/0004-637X/816/1/12
  * Testa et al. (2013) Testa, P., De Pontieu, B., Martínez-Sykora, J., et al. 2013, ApJ, 770, L1, doi: 10.1088/2041-8205/770/1/L1
  * Testa et al. (2014) Testa, P., De Pontieu, B., Allred, J., et al. 2014, Science, 346, 1255724, doi: 10.1126/science.1255724
  * Tiwari et al. (2019) Tiwari, S. K., Panesar, N. K., Moore, R. L., et al. 2019, ApJ, 887, 56, doi: 10.3847/1538-4357/ab54c1
  * Tsuneta et al. (1991) Tsuneta, S., Acton, L., Bruner, M., et al. 1991, Sol. Phys., 136, 37, doi: 10.1007/BF00151694
  * Vekstein & Katsukawa (2000) Vekstein, G., & Katsukawa, Y. 2000, ApJ, 541, 1096, doi: 10.1086/309480
  * Vesecky et al. (1979) Vesecky, J. F., Antiochos, S. K., & Underwood, J. H. 1979, ApJ, 233, 987, doi: 10.1086/157462
  * Warren et al. (2008) Warren, H. P., Ugarte-Urra, I., Doschek, G. A., Brooks, D. H., & Williams, D. R. 2008, ApJ, 686, L131, doi: 10.1086/592960
  * Withbroe & Noyes (1977) Withbroe, G. L., & Noyes, R. W. 1977, ARA&A, 15, 363, doi: 10.1146/annurev.aa.15.090177.002051


# Life is not Always Depressing: Exploring the Happy Moments of People
Diagnosed with Depression

###### Abstract

In this work, we explore the relationship between depression and
manifestations of happiness in social media. While the majority of works
surrounding depression focus on symptoms, psychological research shows that
there is a strong link between seeking happiness and being diagnosed with
depression. We make use of Positive-Unlabeled learning paradigm to
automatically extract happy moments from social media posts of both controls
and users diagnosed with depression, and qualitatively analyze them with
linguistic tools such as LIWC and keyness information. We show that the life
of depressed individuals is not always bleak, with positive events related to
friends and family being more noteworthy to their lives compared to the more
mundane happy events reported by control users.

Keywords: happiness, depression, social media, positive-unlabeled learning,
LIWC, BERT

languageresourceLanguage Resources

Life is not Always Depressing: Exploring the Happy Moments of People Diagnosed
with Depression

Ana-Maria Bucur1,2, Adrian Cosma3, Liviu P. Dinu1
---
1 University of Bucharest, Romania 2 Universitat Politècnica de València,
Spain
2 University Politehnica of Bucharest, Romania
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>

Abstract content

## 1\. Introduction

”I have the choice of being constantly active and happy or introspectively
passive and sad. Or I can go mad by ricocheting in between”, wrote the famous
poet Sylvia Plath in her personal diaries [Plath and Hughes, 2000]. It is now
known that Sylvia Plath was clinically depressed most of her life, before
tragically committing suicide. Depression is the modern ailment of our times.
It is reported that more than 264 million people of all ages suffer from this
mental illness 111https://www.who.int/news-room/fact-sheets/detail/depression.
More than 700 thousand people take their own life every year, and even more
people attempt suicide 222https://www.who.int/news-room/fact-
sheets/detail/suicide. Thankfully, modern psychology has a fair understanding
of depression, its onset, its symptoms and potential treatments, more recently
aided by computational methods [Orabi et al., 2018, Burdisso et al., 2019,
Uban and Rosso, 2020, Inkpen et al., 2021, Bucur and Dinu, 2020, Bucur et al.,
2021a]. However, research in automatic understanding of manifestations of
moments of happiness is lacking. How much do we understand about the happy
moments of depressed individuals? In this preliminary work, we aim to bridge
this gap and develop a computational method for extracting and analyzing happy
moments from a large corpus of social media text.

Indeed, promoting happiness is an important aspect of our lives, and is one of
the main ways in which we can prevent disease and live a longer, healthier
life. While our happiness is determined largely by genetic factors [Diener et
al., 1999], it can be influenced in a significant manner through our behavior:
long-term relationships with friends and family, meaningful work and self-
care. Paradoxically, ?) showed that individuals in cultures that highly value
happiness are less likely to attain it. Moreover, the authors show that highly
valuing happiness is a key indicator and a potential risk factor of
depression. In the same vein, ?) found that happiness is positively related to
optimism and hope and that these traits are negatively associated with
depression and anxiety. Recent computational tools [Cobeli et al., 2022]
currently allow for accurate modelling of the optimism-pessimism continuum
from text, but very little research is performed on the specific instances of
happy moments, and on what aspect makes an individual happy. We posit that
while both depressed and control subjects experience happy moments, the
reasons for their happiness differ.

In this preliminary work, we explore through linguistic tools the qualitative
difference of happy moments between individuals diagnosed with depression and
controls in social media text. While the majority of work in the area of
exploring manifestations of depression in text focuses on specific symptoms
[Inkpen et al., 2021, Bucur and Dinu, 2020, Bucur et al., 2021a], we focus our
attention on exploring happiness manifestations in both depressed individuals
and control. For this goal, we make use of HappyDB [Asai et al., 2018], a
database of 100.000 expressions of happy moments, and conduct our analysis on
eRisk [Losada and Crestani, 2016], a popular dataset of social media posts
from Reddit for early detection of mental disorders. Different from sentiment
polarity, which determines the sentimental aspect of an opinion, happiness is
a broader concept that incorporates well-being and psychological health,
referring to a positive affective state. Since there is no available dataset
for modeling happy moments in depressed individuals and control, we use the
Positive-Unlabeled learning paradigm to automatically extract happy moments
from users’ posts, by combining HappyDB with eRisk to obtain a Positive-
Unlabeled dataset. We train an SVM in this paradigm [Elkan and Noto, 2008] on
top of BERT [Devlin et al., 2019] embeddings to extract happy moments, and use
LIWC [Pennebaker et al., 2001] and keyness statistics to make a qualitative
exploration of the happy moments. This work paves the way to a more in depth
understanding of depression through computational methods. In this way,
automatic tools can be developed to aid in early diagnosis, therapeutic
chatbots, and can provide more insights for developing suitable positive
psychological interventions (PPIs) [Parks and Layous, 2016] (e.g., feeling and
expressing gratitude, pursuing meaning, building hope) for people diagnosed
with depression. PPIs have been shown to improve physical and mental health in
patients with symptoms of depression [D’raven et al., 2015].

Through our analysis, we show that happy moments differ in people diagnosed
with depression compared to control subjects, users from the depression group
are more likely to consider happy moments those events related to family and
friends. Our analysis shows that, for users from the control group, the main
reasons for happiness are related to entertainment situations and financial
plans.

This paper has the following contributions:

  1. 1.

We propose a method for automatically extracting happy moments from social
media text using the Positive-Unlabeled learning paradigm.

  2. 2.

We make a qualitative exploration of happy moments present in the Reddit posts
of both depressed individuals and controls and found noteworthy differences
between the two groups, which are confirmed by psychological studies.

## 2\. Related Work

In psychology, there are different scales for measuring the happiness of an
individual, such as The Oxford Happiness Questionnaire [Hills and Argyle,
2002], The Satisfaction with Life Scale [Diener et al., 1985], The PNAS Scales
[Watson et al., 1988] and others. The PNAS Scales measures different feelings
or emotions (e.g., enthusiastic, alert, nervous, excited) in specific periods
of time (e.g., in this moment, today, past few days), while The Oxford
Happiness Questionnaire and The Satisfaction with Life Scale are comprised of
different statements on happiness, unhappiness and life satisfaction.

In the field of computational linguistics, ?) proposed a corpus-based approach
to finding happiness. The paper analyzes a corpus of blog posts from the
LiveJournal community annotated with sad and happy emotions. The authors’
findings reveal some patterns of happiness and sadness: 3 pm and 9-10 pm are
the happiest hours, and Wednesday is the saddest day of the week. After
analyzing the words from different semantic classes (e.g., actions, things,
socialness), the authors conclude the work with a recipe of happiness with the
following ingredients: ”something new”, ”lots of food that you enjoy”, ”your
favorite drink” and ”an interesting social place”.

In the realm of understanding manifestations of happiness, most notably, the
Workshop on Affective Content Analysis proposed to tackle the problem of
modeling happy moments in the CL-Aff Shared Task 2019 – in Pursuit of
Happiness [Jaidka et al., 2019]. It had two tasks: a) predicting the Agency
(if the author is in control, if the happy moment is influenced entirely by
their actions) and Sociality (if other people besides the author are involved
in the happy moment) of happy moments and b) suggesting insights to
characterize the happy moments automatically.

Even if many works exploring the social media discourse of people diagnosed
with depression [Orabi et al., 2018, Burdisso et al., 2019, Uban and Rosso,
2020, Bucur et al., 2021c] are paying attention to the emotions expressed in
their social media discourse [Aragon et al., 2021, Lara et al., 2021, Howes et
al., 2014], to the best of our knowledge, works focusing on happiness felt by
individuals diagnosed with depression are missing. Therefore, we aim to fill
this gap and explore the happy moments from the online discourse of users with
depression in comparison with control users.

## 3\. Data

Figure 1: Positive-Unlabeled training paradigm, applied to extracting happy
moments from eRisk depression dataset. We randomly sample 50,000 sentences
from eRisk (unlabeled positive and negative classes) and 50,000 sentences from
HappyDB (labeled positive class) to construct a Positive-Unlabeled dataset.
After encoding the sentences with BERT, we trained a linear SVM in the PU-
learning paradigm and use it to extract happy moments from eRisk validation
set.

In our experiments, we use two datasets comprised of English texts: HappyDB
and the eRisk 2018 depression dataset.

HappyDB [Asai et al., 2018] is a corpus of 100,000 happy moments, comprised of
sentences with crowdworkers’ responses to the questions ”what made you happy
in the past 24 hours?” or ”what made you happy in the past 3 months?”. It is a
diverse corpus, containing happy moments from a variety of topics such as
family (e.g., ”My daughter waved at me and said ’mama’ for the first time when
I came home this morning.”), food (e.g., ”I ate a juicy piece of pizza.”),
work (e.g., ”It made me happy to get a promotion at work.”), entertainment
(e.g., ”I went to a concert with friends.”) and others.

In this paper, we use the cleaned version of the HappyDB, in which spelling
errors were corrected and empty or single-word sentences were removed. An
important aspect of HappyDB is that it contains only positive examples (i.e.,
statements of happiness) and does not contain any negative examples, which
makes it difficult to use in standard classification settings.

The eRisk 2018 depression dataset [Losada and Crestani, 2016] contains posts
from Reddit and it was used for the early risk detection of depression task
from the Early risk prediction on the Internet (eRisk)
Workshop333https://early.irlab.org/. Users were annotated as having depression
by their mention of diagnosis (e.g., ”I was diagnosed with depression”) in
their posts. Expression such as ”I have depression” or ”I am depressed” were
not taken into account in annotating the users, only users with explicit
mentions of depression diagnosis were labeled as having depression. Users from
the control group are random users who do not have any mention of diagnosis in
their posts, some of them may be active in the depression
subreddit444https://www.reddit.com/r/depression/. The posts disclosing the
depression diagnosis were removed by the authors from the dataset. The
datasets contains 214 users diagnosed with depression and 1493 control users
with approximately 90,000 submissions from the depression group and over
985,000 submissions from the control group.

## 4\. Methods

To conduct our analysis, we first trained an SVM model on top of BERT sentence
embeddings that is able to extract happy moments from social media text. We
then analyzed the happy moments using LIWC [Pennebaker et al., 2001] and
keyness information.

### 4.1. Extracting Happy Moments

In order to extract relevant happy moments from the eRisk dataset, we train a
classifier according to the Positive-Unlabeled Learning framework [Elkan and
Noto, 2008]. In this setting, the goal is to train a binary classifier to
distinguish between positive classes (i.e., the happy moments) and negative
classes (i.e., any other type of text). However, during classifier training,
only some of the positive examples are labeled, and none of the negative
examples are labeled. Formally, we have a dataset of triples $(x,y,s)$, where
$x$ is an example, $y$ is the class label, $y\in\\{0,1\\}$, and $s$ is a
variable representing whether the pair $(x,y)$ is labeled or not. As such, if
the example is labeled, it belongs to the positive class: $P(y=1|s=1)=1$. When
variable $s=0$, the pair can belong to either class, 0 or 1 in our case. We
chose to adopt the PU-learning methodology, as opposed to a strictly
unsupervised approach, to include additional information in the computation
from the eRisk dataset, which contains Reddit posts qualitatively different
from crowd-sourced texts from HappyDB.

In our case, we split the posts from the eRisk dataset into sentences, and
construct a training dataset comprised of approximately 50,000 sentences from
HappyDB [Asai et al., 2018] and 50,000 sentences from the posts of the first
10 subjects from eRisk (from both depression and control classes). This
dataset is a Positive-Unlabeled dataset (PU dataset), since the sentences from
HappyDB are always the positive class, but the sentences from eRisk can
contain both happy moments and neutral (non-happy) text. This methodology is
presented in Figure 1. Due to limited computational resources, we choose to
use only a subset of the datasets for training our model. We use the rest of
the subjects as a validation set to conduct our analysis.

We train a linear SVM [Cortes and Vapnik, 1995] with high regularization using
the Elkanoto ?) training methodology on the PU dataset. We extract semantic
features using a pretrained SentenceBERT [Devlin et al., 2019] model and train
the SVM on the sentence representations. We used the
pulearn555https://github.com/pulearn/pulearn python library to train our
classifier. To conduct our analysis, we used the previously trained SVM model
on our eRisk validation set. Examples of extracted happy moments from eRisk
validation set are presented in Table 1.

Depression |  Control
---|---
Just found my first Rush dream setlist from long ago. |  Yesterday morning, we woke up to a remarkably pleasant surprise.
I just ordered two pairs of these…they look super comfy! |  I visited Japan recently.
I lost 120lbs through diet and exercise. |  This new vegan supermarket just opened near me!
Getting under 200 is the best feeling! |  I’m on vacation and there’s this amazing yarn shop nearby!
Did get a novella published with a small press, working on bigger stuff. |  I recently quit smoking weed.
I decided that I had to leave the job and make my new job the task of trying to improve my situation. |  About a month ago, my dad finally went into alcohol rehab.
I recently adopted a cat and he’s settled into his new home. |  We went to some amazing gigs and hung out.
My sister bought me The Last of Us Remastered and I’ve just been grinding through it. |  I get all excited when I receive an actual letter from someone I personally know.

Table 1: Selected happy moments automatically extracted from our validation
set of eRisk dataset using Positive-Unlabeled learning. Different from
sentiment polarity, happy moments express significant moments of positive
affective states.

### 4.2. Exploring Happiness

For exploring the happiness moments expressed in the eRisk 2018 dataset, we
use the Linguistic inquiry and word count (LIWC) lexicon [Pennebaker et al.,
2001] due to its psychologically meaningful categories. The LIWC 2001 lexicon
contains words from several categories such as linguistic dimensions (e.g.,
pronouns, articles, negations), psychological processes (e.g., affective
processes, cognitive processes), and personal concerns (e.g. occupation,
leisure). In our analysis, we aim to find the origin of people’s happy
moments, thus we analyze only the suitable categories from LIWC; we are not
interested in the linguistics dimensions. As such, we chose the categories of
sensory and perceptual processes (e.g. HEARING, FEELING), social processes
(e.g., FAMILY, FRIENDS), and personal concerns (e.g., SCHOOL, SPORTS,
RELIGION).

For comparing the happy moments from the two corpora (depression and control),
we compute the dominance scores [Mihalcea and Pulman, 2009] for the LIWC
features analyzed, thus comparing the relative difference in the use of words
from different categories between users diagnosed with depression and control
users. The coverage is computed by comparing the words from the foreground
corpus (F) with the words from the background corpus (B). The coverage of a
specific class (C) in a corpus (X) is computed as follows:

$Coverage_{X}(C)=\frac{\sum_{w_{i}\epsilon C}Frequency_{X}(W_{i})}{Size_{X}}$
(1)

where $C=W_{1},W_{2},...,W_{N}$ represents a class of words.

The dominance score of the foreground corpus with respect to the background
corpus is computed as a ratio between the coverages of the two corpora:

$Dominance_{F}(C)=\frac{Coverage_{F}(C)}{Coverage_{B}(C)}$ (2)

A dominance score equal or close to 1 for a given class indicates that the
foreground and background corpora contain a similar distribution of words from
the class. A score greater than 1 indicates that the foreground corpus
contains more words from a given class than the background corpus. In
contrast, a value lower than 1 suggests that words from a class are
underrepresented in the foreground corpus instead of the background corpus.

Depression |  Control
---|---
Category | Score |  Example words | Category | Score |  Example words
FRIEND S | 1.81 |  friend, love, boyfriend, girlfriend, mate, bud, neighbor | TV | 1.74 |  ad, show, movie, tv, actor, video
MUSIC | 1.48 |  song, sing, band, rap, listen, music | BODY | 1.21 |  eye, heal, breath, pregnant, cheek, skin
SLEEP | 1.47 |  night, bed, sleep, dream, woke, nap | MONEY | 1.41 |  fee, owe, rent, bought, pay, invest
FAMILY | 1.45 |  dad, kin, ex, son, great, family, mom | DEATH | 1.02 |  die, dead, decay, dying, ashes, terminate
SEXUAL | 1.38 |  sex, bi, love, stud, fuck, hug | RELIGION | 0.99 |  sin, christ, hell, god, angel, bless
SELF | 1.36 |  I, me, my, we, us, our, mine | SPORTS | 0.97 |  play, game, ski, team, running, sport, exercise
Table 2: Top classes and dominance scores for depression and control. Example
words for each category are the ones found in our two corpora.

## 5\. Results and Discussion

We computed dominance scores with each class in the foreground to reveal the
dominant LIWC categories in the two corpora (depression and control). In Table
2 we have on the left side the analysis with the happy moments from the
depression class as the foreground corpus and the happy moments of control
subjects as the background corpus. On the right side of Table 2, the moments
from the control users are the foreground corpus. From these results, it is
apparent that depressed individuals more frequently manifest happy moments
relating to family and friends (dominance scores $>$ 1.4): ”I’ve been training
regularly with my entire family for 3 months now and I love every bit of it.”,
”I’m with a woman I love and I have really supportive friends and family.”.
Such events are more likely to be considered of greater value when compared to
the happy moments of the control group. The control group, however, more
frequently exhibit happy moments in everyday situations, in contexts related
to entertainment, sports and financial plans: ”Knowing I make more money than
you does make me happy.”, ”Eating pop corn and seeing movie.”, ”I bought a new
controller and worked perfect.”. In line with psychology research, ?) found
that the support of family and friends is negatively correlated with
depression symptoms, in a study conducted on adolescents. Events related to
positive experiences with loved ones are more noteworthy in a depressed
subject’s life, more so than in the case of control subjects.

Figure 2: Keyness scores for nouns in the extracted happy moments from eRisk.
Individuals with depression are more likely to discuss subjects related to
family and friends than control users. Figure 3: Keyness scores for verbs in
the extracted happy moments from eRisk. Individuals with depression are more
likely to use emotion-related verbs, while users from the control group use
more action verbs.

Moreover, introspective discourse is more prevalent in depressed individuals
(i.e., SELF category), which is in accordance with other works [Bucur et al.,
2021b, Rude et al., 2004]. Happy moments related to sleep are also more
prevalent than those of the control category - poor sleep quality being one of
the depression symptoms [O’Leary et al., 2017]. Subjects diagnosed with
depression also have more happy moments related to music (e.g., ”Together we
crafted a legacy of music that has given so many people happiness and great
memories.”, ”I read a bit from my book and listened to music.”), which is
indicative of the positive effects music therapy has on mental well-being
[Tang et al., 2020].

A surprising result was to find the death category from LIWC in the reasons
for the happy moments of the control group. However, a closer manual analysis
of the moments containing keywords from the death category revealed that, even
if this category is comprised of death-related words (e.g., kill, dead,
dying), they may not be used with their primary sense (e.g., ”I have been able
to increase dead lifts from 70 kg to 100kg.”) or the words may be used in the
context of video games (e.g., ”I killed 5 King Black Dragons.”).

To further understand the differences in the happy moments reported by users
with depression and control, we compute the keyness score [Kilgarriff, 2009,
Gabrielatos, 2018] for the two corpora. We perform a keyness analysis by
comparing the frequencies of words from the happy moments of users diagnosed
with depression (target corpus) to the frequencies of words from the texts of
control users (reference corpus). We perform this analysis on nouns and verbs.
We present the top 20 nouns and verbs from the two corpora (depression and
control) ordered by their log-likelihood ratio (G2) [Dunning, 1993] in Figure
2 and 3.

In Figure 2 we show the most frequent nouns in the happy moments of
individuals diagnosed with depression in comparison with the control group. A
similar pattern of words related to family and friends is also found in the
keyness analysis of words from the happy moments of users from the depression
group.

In Figure 3 we present the most frequent 20 verbs from each group. Users
diagnosed with depression use more verbs related to emotions (e.g., feel, cry,
cuddle), as opposed to users from the control group, which use more action
verbs (e.g., build, lead, run) in the texts of their happy moments.

## 6\. Conclusion

In this work, we explored the reasons for happiness in social media text.
Using the HappyDB dataset of happy moments and the Positive-Unlabeled Learning
(PU-learning) framework, we automatically constructed a model to extract happy
moments from the eRisk dataset. Our linguistic analysis through LIWC and
keyness information showed that people diagnosed with depression more commonly
experience joy from moments related to their family and friends and attribute
greater significance to such events, while control subjects have happiness
moments related to common everyday activities such as subjects related to
money and entertainment.

This preliminary work is an exploration into manifestations of happiness of
depressed individuals. The datasets used may be positively biased towards
depressed individuals, which themselves attribute great value to their
experiences of happiness, and choose to express it online; it is not clear the
extent to which control users are expressing their happy moments online,
rather than just enjoying them in the real world. However, our preliminary
results are worthy of further development in future work, being in line with
research performed in psychological studies [Tang et al., 2020, Rude et al.,
2004], and could prove helpful in developing therapeutic chatbots and
developing suitable positive psychological interventions. Furthermore, we
believe that our results pave the way to a more in-depth analysis of
expressions of happiness, by analyzing aspects of agency and sociality in
happy moments. Moreover, due to computational limitations, we trained our PU-
learning model on a sub-sample of the corpus; in future work, we aim to
include the entire HappyDB dataset into the training data and perform our
analyses on larger datasets (i.e. RSDD [Yates et al., 2017]).

## 7\. Bibliographical References

## References

  * Alarcon et al., 2013 Alarcon, G. M., Bowling, N. A., and Khazon, S.  (2013).  Great expectations: A meta-analytic examination of optimism and hope.  Personality and Individual Differences, 54(7):821–827.
  * Aragon et al., 2021 Aragon, M. E., Lopez-Monroy, A. P., Gonzalez-Gurrola, L.-C. G., and Montes, M.  (2021).  Detecting mental disorders in social media through emotional patterns-the case of anorexia and depression.  IEEE Transactions on Affective Computing.
  * Asai et al., 2018 Asai, A., Evensen, S., Golshan, B., Halevy, A., Li, V., Lopatenko, A., Stepanov, D., Suhara, Y., Tan, W.-C., and Xu, Y.  (2018).  HappyDB: A corpus of 100,000 crowdsourced happy moments.  In Proceedings of the Eleventh International Conference on Language Resources and Evaluation (LREC 2018), Miyazaki, Japan, May. European Language Resources Association (ELRA).
  * Bucur and Dinu, 2020 Bucur, A.-M. and Dinu, L. P.  (2020).  Detecting early onset of depression from social media text using learned confidence scores.
  * Bucur et al., 2021a Bucur, A., Cosma, A., and Dinu, L. P.  (2021a).  Early risk detection of pathological gambling, self-harm and depression using BERT.  CoRR, abs/2106.16175.
  * Bucur et al., 2021b Bucur, A.-M., Podina, I. R., and Dinu, L. P.  (2021b).  A psychologically informed part-of-speech analysis of depression in social media.  In Proceedings of the International Conference on Recent Advances in Natural Language Processing (RANLP 2021), pages 199–207, Held Online, September. INCOMA Ltd.
  * Bucur et al., 2021c Bucur, A.-M., Zampieri, M., and Dinu, L. P.  (2021c).  An exploratory analysis of the relation between offensive language and mental health.  In Findings of the Association for Computational Linguistics: ACL-IJCNLP 2021, pages 3600–3606, Online, August. Association for Computational Linguistics.
  * Burdisso et al., 2019 Burdisso, S. G., Errecalde, M., and Montes-y Gómez, M.  (2019).  A text classification framework for simple and effective early depression detection over social media streams.  Expert Systems with Applications, 133:182–197.
  * Cobeli et al., 2022 Cobeli, S., Iordache, B., Yadav, S., Caragea, C., Dinu, L., and Iliescu, D.  (2022).  Detecting optimism in tweets using knowledge distillation and linguistic analysis of optimism.  In Proceedings of the Thirteenth International Conference on Language Resources and Evaluation (LREC 2022), Marseille, France. European Language Resources Association (ELRA).
  * Cortes and Vapnik, 1995 Cortes, C. and Vapnik, V.  (1995).  Support-vector networks.  Machine learning, 20(3):273–297.
  * Devlin et al., 2019 Devlin, J., Chang, M.-W., Lee, K., and Toutanova, K.  (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, Minneapolis, Minnesota, June. Association for Computational Linguistics.
  * Diener et al., 1985 Diener, E., Emmons, R. A., Larsen, R. J., and Griffin, S.  (1985).  The satisfaction with life scale.  Journal of personality assessment, 49(1):71–75.
  * Diener et al., 1999 Diener, E., Suh, E. M., Lucas, R. E., and Smith, H. L.  (1999).  Subjective well-being: Three decades of progress.  Psychological bulletin, 125(2):276.
  * Dunning, 1993 Dunning, T. E.  (1993).  Accurate methods for the statistics of surprise and coincidence.  Computational linguistics, 19(1):61–74.
  * D’raven et al., 2015 D’raven, L. T. L., Moliver, N., and Thompson, D.  (2015).  Happiness intervention decreases pain and depression, boosts happiness among primary care patients.  Primary health care research & development, 16(2):114–126.
  * Elkan and Noto, 2008 Elkan, C. and Noto, K.  (2008).  Learning classifiers from only positive and unlabeled data.  In Proceedings of the 14th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 213–220.
  * Ford et al., 2014 Ford, B. Q., Shallcross, A. J., Mauss, I. B., Floerke, V. A., and Gruber, J.  (2014).  Desperately seeking happiness: Valuing happiness is associated with symptoms and diagnosis of depression.  Journal of Social and Clinical Psychology, 33(10):890–905, December.
  * Gabrielatos, 2018 Gabrielatos, C.  (2018).  Keyness analysis: Nature, metrics and techniques.  In Corpus Approaches to Discourse, pages 225–258. Routledge.
  * Hills and Argyle, 2002 Hills, P. and Argyle, M.  (2002).  The oxford happiness questionnaire: a compact scale for the measurement of psychological well-being.  Personality and individual differences, 33(7):1073–1082.
  * Howes et al., 2014 Howes, C., Purver, M., and McCabe, R.  (2014).  Linguistic indicators of severity and progress in online text-based therapy for depression.  In Proceedings of the Workshop on Computational Linguistics and Clinical Psychology: From Linguistic Signal to Clinical Reality, pages 7–16.
  * Inkpen et al., 2021 Inkpen, D., Skaik, R., Gamaarachchige, P., Angelov, D., and Fredenburgh, M. T.  (2021).  uottawa at erisk 2021: Automatic filling of the beck’s depression inventory questionnaire using deep learning.  Working Notes of CLEF, pages 21–24.
  * Jaidka et al., 2019 Jaidka, K., Mumick, S., Chhaya, N., and Ungar, L.  (2019).  The cl-aff happiness shared task: Results and key insights.  In AffCon@ AAAI.
  * Kilgarriff, 2009 Kilgarriff, A.  (2009).  Simple maths for keywords.  In Proc. Corpus Linguistics.
  * Lara et al., 2021 Lara, J. S., Aragon, M. E., Gonzalez, F. A., and Montes-y Gomez, M.  (2021).  Deep bag-of-sub-emotions for depression detection in social media.  arXiv preprint arXiv:2103.01334.
  * Losada and Crestani, 2016 Losada, D. E. and Crestani, F.  (2016).  A test collection for research on depression and language use.  In International Conference of the Cross-Language Evaluation Forum for European Languages, pages 28–39. Springer.
  * Mihalcea and Liu, 2006 Mihalcea, R. and Liu, H.  (2006).  A corpus-based approach to finding happiness.  In AAAI Spring Symposium: Computational Approaches to Analyzing Weblogs, pages 139–144.
  * Mihalcea and Pulman, 2009 Mihalcea, R. and Pulman, S.  (2009).  Linguistic ethnography: Identifying dominant word classes in text.  In International Conference on Intelligent Text Processing and Computational Linguistics, pages 594–602. Springer.
  * Orabi et al., 2018 Orabi, A. H., Buddhitha, P., Orabi, M. H., and Inkpen, D.  (2018).  Deep learning for depression detection of twitter users.  In Proceedings of the Fifth Workshop on Computational Linguistics and Clinical Psychology: From Keyboard to Clinic, pages 88–97.
  * O’Leary et al., 2017 O’Leary, K., Bylsma, L. M., and Rottenberg, J.  (2017).  Why might poor sleep quality lead to depression? a role for emotion regulation.  Cognition and Emotion, 31(8):1698–1706.
  * Parks and Layous, 2016 Parks, A. C. and Layous, K.  (2016).  Positive psychological interventions. American Psychological Association.
  * Pennebaker et al., 2001 Pennebaker, J. W., Francis, M. E., and Booth, R. J.  (2001).  Linguistic inquiry and word count: Liwc 2001.  Mahway: Lawrence Erlbaum Associates, 71(2001):2001.
  * Plath and Hughes, 2000 Plath, S. and Hughes, T.  (2000).  The unabridged journals of Sylvia Plath.  Anchor Books.
  * Rude et al., 2004 Rude, S., Gortner, E.-M., and Pennebaker, J.  (2004).  Language use of depressed and depression-vulnerable college students.  Cognition & Emotion, 18(8):1121–1133.
  * Tang et al., 2020 Tang, Q., Huang, Z., Zhou, H., and Ye, P.  (2020).  Effects of music therapy on depression: A meta-analysis of randomized controlled trials.  PloS one, 15(11):e0240862.
  * Uban and Rosso, 2020 Uban, A.-S. and Rosso, P.  (2020).  Deep learning architectures and strategies for early detection of self-harm and depression level prediction.  In CEUR Workshop Proceedings, volume 2696, pages 1–12. Sun SITE Central Europe.
  * van Harmelen et al., 2016 van Harmelen, A.-L., Gibson, J. L., St Clair, M. C., Owens, M., Brodbeck, J., Dunn, V., Lewis, G., Croudace, T., Jones, P. B., Kievit, R. A., and Goodyer, I. M.  (2016).  Friendships and family support reduce subsequent depressive symptoms in at-risk adolescents.  PLOS ONE, 11(5):1–20, 05.
  * Watson et al., 1988 Watson, D., Clark, L. A., and Tellegen, A.  (1988).  Development and validation of brief measures of positive and negative affect: the panas scales.  Journal of personality and social psychology, 54(6):1063.
  * Yates et al., 2017 Yates, A., Cohan, A., and Goharian, N.  (2017).  Depression and self-harm risk assessment in online forums.  In Proceedings of the 2017 Conference on Empirical Methods in Natural Language Processing, pages 2968–2978, Copenhagen, Denmark, September. Association for Computational Linguistics.