# The role of Hi in regulating size growth of local galaxies

Zhizheng Pan<EMAIL_ADDRESS>Purple Mountain Observatory, Chinese Academy of
Sciences, 8 Yuan Hua Road, Nanjing, Jiangsu 210033, China  Jing Wang  Kavli
Institute for Astronomy and Astrophysics, Peking University, Yi He Yuan Lu 5,
Hai Dian District, Beijing 100871, China  Xianzhong Zheng  Purple Mountain
Observatory, Chinese Academy of Sciences, 8 Yuan Hua Road, Nanjing, Jiangsu
210033, China  Xu Kong  CAS Key Laboratory for Research in Galaxies and
Cosmology, Department of Astronomy,
University of Science and Technology of China, Hefei, Anhui 230026, China

###### Abstract

We study the role of atomic hydrogen (Hi) in regulating size growth of local
galaxies. The size of a galaxy, $D_{\rm r,~{}25}$, is characterized by the
diameter at which the $r-$band surface brightness reaches $\mu_{\rm
r}=25.0~{}\rm mag~{}arcsec^{-2}$. We find that the positions of galaxies in
the size ($D_{\rm r,~{}25}$)$-$stellar mass ($M_{\ast}$) plane strongly depend
on their Hi-to-stellar mass ratio ($M_{\rm HI}/M_{\ast}$). In the Hi-rich
regime, galaxies that are more rich in Hi tend to have larger sizes. Such a
trend is not seen in the Hi-poor regime, suggesting that size growth is barely
affected by the Hi content when it has declined to a sufficiently low level.
An investigation of the relations between size, $M_{\rm HI}/M_{\ast}$ and star
formation rate (SFR) suggests that size is more intrinsically linked with
$M_{\rm HI}/M_{\ast}$, rather than SFR. We further examine the Hi-to-stellar
disk size ratio ($D_{\rm HI}/D_{\rm r,~{}25}$) of galaxies and find that at
log($M_{\rm HI}/M_{\ast})>-0.7$, $D_{\rm HI}/D_{\rm r,~{}25}$ is weakly
correlated with $M_{\ast}$. These findings support a picture in which the Hi-
rich galaxies live in an inside-out disk growing phase regulated by gas
accretion and star formation. The angular momentum of the accreted materials
is probably the key parameter in shaping the size of an Hi-rich galaxy.

galaxies: evolution – galaxies: structure – galaxies: statistics

##  1 Introduction

The sizes of galaxies play a critical role in our understanding of how they
form and evolve. Observationally, the absence of a clear border of a galaxy
makes the measuring of its extent a non-trivial task. In the literature, there
are two popular approaches for characterizing the size of a galaxy. The first,
which is also the most popular approach, is defining the size of a galaxy as
the radial distance that enclosing half of its light, i.e., the effective
radius $R_{\rm e}$ (de Vaucouleurs, 1948). A second commonly used approach is
defining galaxy size at the radial location of a given isophote (for example
$R_{\rm 25}$, the radius at which the surface brightness reaches
$\mu=25.0~{}\rm mag~{}arcsec^{-2}$). Each approach has its superiority and
shortcoming. As already known, $R_{\rm e}$ is barely affected by the depth of
image but is quite sensitive to the light profile of a galaxy, making it
incapable of describing the global extension of a galaxy (Graham, 2019). The
second approach performs better in describing the global size of a galaxy but
requires the image to reach a certain depth.

Observational studies have demonstrated that the sizes of galaxies evolve
across cosmic time (e.g., van der Wel et al., 2014). On the one hand, mergers
are believed to significantly contribute to this evolution, especially at the
high mass end. Simulations suggest that the size of the emerged galaxy is
strongly dependent on the physical conditions (stellar mass, gas mass, etc) of
progenitor galaxies before merging (Hopkins et al., 2009). On the other hand,
in a merger-free case, theoretical studies predict that star-forming galaxies
(SFGs) grow from the inside out via gas accretion and star formation (e.g.,
Pichon et al., 2011). This inside-out galaxy formation picture is supported by
a growing body of observational studies (Wang et al., 2011; Dale et al., 2016;
Frankel et al., 2019; Chen et al., 2020).

The second scenario points out a potential link between size growth and the
cold gas content of a galaxy. Neutral atomic hydrogen (Hi) gas is the raw
material from which molecular gas and then star formation form. Studying the
scaling relations between galaxy properties and their Hi content thus provides
important insights on understanding galaxy formation. Based on the data
provided by the GALEX Arecibo SDSS Survey (GASS), Catinella et al. (2010)
showed that the concentration of galaxies is only weakly correlated with their
Hi content. Subsequent studies found that Hi-rich galaxies have bluer outer
disk than the Hi-normal ones, implying continuous growth of the outer disks
supplied by rich Hi gas reservoirs (Wang et al., 2011; Huang et al., 2014;
Kauffmann, 2015; Yıldız et al., 2017). Recently, Chen et al. (2020) study the
relations between the properties of bulge/disk component and the Hi content of
galaxies. They found that the color of disk is bluer in Hi rich galaxies,
whereas a similar trend is not found for the bulge component. This finding
suggests that Hi gas is closely related to the formation of disks but not
necessarily fuel the star formation of bulges in an efficient way.

Figure 1: Left: the correlation between mass-to-light ratio $M/L_{i}$ and the
SDSS $g-i$ color index. We fit the relation with a broken linear relation, as
shown in the red line. Right: the $\Delta~{}M$ distribution of the
spectroscopic-z sample, where $\Delta~{}M$ is the stellar mass difference
between that given by the JHU/MPA catalog and that estimated using the
$M/L_{i}$ vs.($g-i$) relation. The red dashed line shows a Gaussian function
with $\sigma=0.03$.  Figure 2: Left: the size$-$mass relation of the ALFALFA
sample. The large symbols represent the median log$(D_{\rm r,~{}25})$ in that
$M_{\ast}$ bin, with a bin size of $\Delta$log$M_{\ast}=0.2$. Error bars
represent the 1$\sigma$ standard deviation. We use a second-order polynomial
to fit this relation, and the best-fit result is shown in the gray dashed
line. Small symbols are color coded by log$(M_{\rm HI}/M_{\ast})$. The small
panel shows the residuals of the fit against stellar mass. Right: the
size$-$mass relation of the xGASS sample, with its best fit result shown in
black dash-dotted line.

Given the fact that Hi content is closely related to the formation of disk
component of galaxies, it is worthy to further investigate at what condition
the Hi gas could efficiently impact on size growth. In this paper, we aim to
study how size growth is regulated by the Hi gas of galaxies based on the data
from the extended GASS (xGASS, Catinella et al. 2018 ) and Arecibo Legacy Fast
ALFA (ALFALFA) survey. Throughout this paper, we adopt a concordance
$\Lambda$CDM cosmology with $\Omega_{\rm m}=0.3$, $\Omega_{\rm\Lambda}=0.7$,
$H_{\rm 0}=70$ $\rm km~{}s^{-1}$ Mpc-1 and a Kroupa (2001) initial mass
function (IMF).

##  2 Data

###  2.1 The ALFALFA sample

One of the Hi sample used in this work is drawn from the Arecibo Legacy Fast
ALFA (ALFALFA) survey (Giovanelli et al., 2005). We use the $\alpha.70$
catalog, which includes HI sources exacted from 70% of the final survey
area111http://egg.astro.cornell.edu/alfalfa/data/index.php. This catalog
contains $\sim$ 25,000 HI sources at $z<0.06$, with more than 95% having
identified optical counterparts. The Hi mass ($M_{\rm HI}$) was calculated by

$\frac{M_{\rm HI}}{\rm M_{\sun}}=(\frac{2.325\times 10^{5}}{1+z})(\frac{D_{\rm
L}^{2}}{\rm Mpc})(\frac{S_{\rm 21}}{\rm Jy~{}kms^{-1}}),$ (1)

where $D_{\rm L}$ is the luminosity distance of the source, $z$ is redshift
and $S_{\rm 21}$ is the integrated 21-cm line flux density. $D_{\rm L}$ is
determined using the Hubble Law for the sources that with $cz>6000~{}\rm
km~{}s^{-1}$. For sources with $cz<6000~{}\rm km~{}s^{-1}$, $D_{\rm L}$ is
derived using the local universe peculiar velocity model of Marsters (2005).
The $\alpha.70$ catalog also provides the coordinates of the matched optical
counterparts and the signal-to-noise ratio (S/N) of the Hi sources.

The optical data are drawn from the Sloan Digital Sky Survey (SDSS, York et
al. 2000). We cross-matched the ALFALFA sources that with confident (or
marginal) Hi detections (with Hi line detection flag 1 (or 2), see the
$\alpha.70$ catalog) with the SDSS data release 7 (DR7) database via the SDSS
online SQL tool 222http://skyserver.sdss.org/dr7/en/tools/crossid/crossid.asp,
yielding a sample of 17,134 galaxies.

Of the 17,134 SDSS-ALFALFA matched galaxies, 13,067 ones have spectroscopic
redshifts and stellar mass ($M_{\ast}$) estimates 333The $M_{\ast}$ of the
spec-z galaxies are available at http://www.mpa-garching.mpg.de/SDSS/DR7/.. In
this paper we wish to estimate the $M_{\ast}$ of the 17,134 galaxies in a
uniform way. We first derive the $i$-band absolute magnitude of individual
galaxies using the SDSS photometry and the distance information given in the
ALFALFA catalog. As confirmed by previous works, the optical broad-band color
is in good correlation with the mass-to-light ratio ($M/L$) of galaxies (Bell
et al., 2003; Taylor et al., 2011; Fang et al., 2013). In the left panel of
Figure 1, we show $M/L_{\rm i}$ as a function of $g-i$ color for the 13,067
galaxies that with $M_{\ast}$ estimates. The $g-i$ color has been corrected
for Galactic extinction based on the extinction map of Schlegel et al. (1998).
Since the log$(M/L_{\rm i})$ vs. ($g-i$) relation appears to have a slightly
shallower slope at red colors, we fit this relation with a broken linear
function:

$\displaystyle\mathrm{log}(M/L_{i})=\begin{cases}1.0(g-i)+0.93,&g-i\leq 0.8\\\
\\\ 0.85(g-i)+1.05,&g-i\geq 0.8\\\ \end{cases}$ (2)

We compare the stellar masses derived using equation (2) with those of the MPA
catalog and show the result in the right panel of Figure 1. As can be seen,
stellar masses derived from equation (2) are well consistent with those from
the JHU/MPA catalog. In the following sections, $M_{\ast}$ is referred to
which we estimated using equation (2). We further restrict galaxies to have
$M_{\ast}>10^{9.0}M_{\sun}$ and $z>0.01$, yielding a final sample of 11,378
galaxies.

Figure 3: Offset from the size$-$mass relation, $\Delta$log$D_{r,~{}25}$, as a
function of Hi-to-stellar mass ratio ($M_{\rm HI}/M_{\ast}$). In each panel,
the large solid circles represent the median $\Delta$log$D_{r,~{}25}$ in that
$M_{\rm HI}/M_{\ast}$ bin, with a bin size of $\Delta$log$(M_{\rm
HI}/M_{\ast})=0.3$. Error bars represent the 1$\sigma$ standard deviation.

###  2.2 The xGASS representative sample

The ALFALFA survey is relatively shallow and biased to Hi-rich galaxies. To
extend our study to more Hi-deficient galaxies, we supplement the Hi sample
with galaxies drawn from xGASS (extended GALEX Arecibo SDSS Survey, see
Catinella et al. 2018). xGASS is a census of the Hi content of 1179 galaxies
using the Arecibo telescope, and the galaxies are selected from the
intersected area of SDSS, GALEX and ALFALFA. The galaxy sample of xGASS spans
a mass range of $M_{\ast}>10^{9.0}M_{\sun}$ at redshift $0.01<z<0.05$. To
optimize the survey efficiency, galaxies that with reliable Hi detection in
ALFALFA were not observed again. The rest galaxies are observed until the Hi
line is detected, or a limit in $M_{\rm HI}/M_{\ast}$ is reached. The HI
spectrum of each galaxy is flagged according to its quality. A catalog
containing the Hi and other properties (such as star formation rates (SFRs),
concentration index, etc.) of the 1179 galaxies is compiled by Catinella et
al. (2018) 444The catalog is available at https://xgass.icrar.org/data.html.
We draw data from the xGASS catalog and convert $M_{\ast}$ from Chabrier
(2003) to Kroupa (2001) IMF by multiplying a factor of 1.06. A (1+$z$) factor
is multiplied to the xGASS $M_{\rm HI}$ to ensure that both the $M_{\rm HI}$
of ALFALFA and xGASS are calculated with a same method. In the following
analysis, we only use the 804 galaxies that with reliable Hi detections.

##  3 Result

###  3.1 Impacts of Hi on the size-mass relation

In this paper, the size of a galaxy is defined as $D_{\rm
r,~{}25}=2R_{r,~{}25}$, where $R_{r,~{}25}$ is the radius at which the SDSS
$r-$band surface brightness reaches $\mu_{\rm r}=25.0~{}\rm
mag~{}arcsec^{-2}$. We use this size indicator because it gives a better
description to the global size of a galaxy than the effective radius $R_{\rm
e}$, as mentioned above. The quality of SDSS imaging data allows a robust
measurement of surface brightness profile of galaxies down to $\mu_{\rm
r}=26-27$ $\rm mag~{}arcsec^{-2}$ (Pohlen & Trujillo, 2006; Wang et al.,
2018). In this work, $R_{r,~{}25}$ is drawn from the SDSS DR7 photometric
database. We have checked the quality of the SDSS-measured $R_{r,~{}25}$ and
confirmed that it is suitable for our statistical study (see Appendix A).

In Figure 2, we show the $D_{\rm r,~{}25}-M_{\ast}$ relation for the ALFALFA
and xGASS galaxy samples. Symbols are color coded by Hi gas fraction, defined
as $f_{\rm HI}=M_{\rm HI}/M_{\ast}$. Since the $D_{\rm r,~{}25}-M_{\ast}$
relations of the two samples both show a shallower slope at the low-mass end,
we use a second-order polynomial to fit this relation over
log$(M_{\ast}/M_{\sun})=[9.0,11.5]$. As shown in the small panels of Figure 2,
the residual of the fitting is independent of stellar mass. The dispersion of
the $D_{\rm r,~{}25}-M_{\ast}$ relation is $\sigma\sim 0.12$ dex, which is
significantly smaller than that of the classical $R_{\rm e}-M_{\ast}$ relation
(Shen et al., 2003; van der Wel et al., 2014), as also reported by some
previous studies (Cortese et al., 2012; Trujillo et al., 2020). A most
important new feature revealed in Figure 2 is that, at fixed $M_{\ast}$,
galaxies that are more rich in Hi tend to have larger sizes. This feature also
holds when a new size indicator defined by the radial position of a given
isomass contour is utilized (see Appendix C).

As shown in previous studies, at fixed $M_{\ast}$, galaxies that with larger
$R_{\rm e}$ tend to have bluer colors and younger stellar ages (Shen et al.,
2003; Lange et al., 2015; Scott et al., 2017). Since blue galaxies also tend
to be rich in gas, the trend between the sizes and gas fraction of galaxies
shown in Figure 2 is expected.

We further divide our samples into 4 stellar mass bins and examine the
relation between $\Delta$log$D_{r,~{}25}$ and $f_{\rm HI}$.
$\Delta$log$D_{r,~{}25}$ is defined as the offset of a galaxy from the
size$-$mass relation, i.e., $\Delta$log$D_{\rm r,~{}25}$=log$D_{\rm
r,~{}25}-$log$\overline{D}_{\rm r,~{}25}$, where log$\overline{D}_{\rm
r,~{}25}$ is the median size of galaxies as fitted from the xGASS sample. The
results are shown in Figure 3. When replacing log$\overline{D}_{\rm r,~{}25}$
with that fitted from the ALFALFA sample, the conclusion is unchanged.
Overall, the trends between $\Delta$log$D_{r,~{}25}$ and $f_{\rm HI}$ are
similar in all panels. In the high $f_{\rm HI}$ regime, galaxies that are more
rich in Hi tend to have larger $\Delta$log$D_{r,~{}25}$. However, when $f_{\rm
HI}$ gets sufficiently low, this correlation breaks down, i.e.,
$\Delta$log$D_{r,~{}25}$ no longer depends on $f_{\rm HI}$.

The turn over gas fraction ($f_{\rm t}$) at which the
$\Delta$log$D_{r,~{}25}-f_{\rm HI}$ begins to break down seems to depend on
$M_{\ast}$. A visual inspection of Figure 3 suggests that at
log$(M_{\ast}/M_{\sun})<10.0$, log$f_{\rm t}=[-1.0,-0.7]$, while at
log$(M_{\ast}/M_{\sun})>10.0$, log$f_{\rm t}=[-1.5,-1.0]$. We speculate that
$f_{\rm t}$ could be somewhat related to the depth of the Hi survey used.
Nevertheless, for the mass range we probed, it should be safe to conclude that
$\Delta$log$D_{r,~{}25}$ positively correlates with $f_{\rm HI}$ at log$f_{\rm
HI}>-0.7$.

Figure 3 clearly demonstrates that there exists a positive correlation between
size and Hi fraction at log$f_{\rm HI}>-0.7$. If an isomass-defined size
indicator is used, this trend slightly weakens but still exists (see Appendix
C). This suggests that the $M/L$ effect also contributes to dispersion of
$D_{\rm r,~{}25}$ at fixed $M_{\ast}$. We do not use an isomass-defined size
indicator in the main body of this paper because this will significantly
narrow down the dynamical range of size at fixed $M_{\ast}$, making it
difficult to examine the trend between size and Hi fraction. However, if an
isomass-defined size indicator is used, the main conclusions of this work
should also hold.

For convenience, we refer galaxies with log$f_{\rm HI}>-0.7$ and log$f_{\rm
HI}<-0.7$ as Hi-rich and Hi-poor galaxies, respectively. For Hi-rich galaxies,
we interpret the positive correlation between $\Delta$log$D_{r,~{}25}$ and
$f_{\rm HI}$ as a natural consequence of the “inside-out” disk formation
scenario (e.g., Chiappini et al., 1997; Wang et al., 2011; Pérez et al., 2013;
Pan et al., 2015). In theory, normal SFGs are assumed to live in a quasi-
equilibrium phase balanced by gas accretion, star formation and outflows. In
this case, galaxies with a high $f_{\rm HI}$ preferentially have a higher gas
accretion rate (Wang et al., 2013). In the cold mode accretion, the accreted
materials usually have high specific angular momentum that help forming an
extended disk component (Stewart et al., 2011, 2013), which naturally results
in a positive $\Delta$log$D_{r,~{}25}-f_{\rm HI}$ correlation.

The weak correlation between $\Delta$log$D_{\rm r,~{}25}$ and $f_{\rm HI}$
exhibited in the Hi-poor regime implies that Hi-poor galaxies have a different
disk growth mode compared to their gas-rich counterparts. As shown below, Hi-
poor galaxies mainly distribute below the star formation main sequence (SFMS).
This suggests that many Hi-poor galaxies are quenching their star formation,
or have been quenched. During the quenching phase, external cold gas
replenishment is expected to be blocked and galaxies have left the star-
forming quasi-equilibrium phase (Peng et al., 2015). In this case, $f_{\rm
HI}$ simply reflects the amount of the residual gas reservoirs of galaxies. We
thus suggest that the weak correlation between size and Hi content shown in
the log$f_{\rm HI}<-0.7$ regime is not at odds with the inside-out disk
formation scenario.

Figure 4: The correlation between Hi-to-stellar size ratio ($D_{\rm HI}/D_{\rm
r,~{}25}$) and Hi-to-stellar mass ratio ($M_{\rm HI}/M_{\ast}$). Symbols are
color coded by $M_{\ast}$. The dashed line marks log$(M_{\rm
HI}/M_{\ast})=-0.7$, which separates the parameter plane into an Hi-poor and
an Hi-rich region. The large solid circles represent the median log$(D_{\rm
HI}/D_{\rm r,~{}25})$ in that $M_{\rm HI}/M_{\ast}$ bin, with a bin size of
$\Delta$log$(M_{\rm HI}/M_{\ast})=0.2$. Error bars represent the 1$\sigma$
standard deviation.

###  3.2 The Hi-to-stellar size ratio

Observationally, a SFG typically has an Hi gaseous disk that is more extended
than its stellar disk. As stars are formed from cold gas, a comparison between
the sizes of Hi disk and stellar disk may provide insights in understanding
the role of Hi in regulating size growth. In the literature, the Hi size of a
galaxy, $D_{\rm HI}$, is defined as the diameter at which the Hi surface
density reaches $\Sigma_{\rm HI}=1~{}\rm M_{\sun}\rm pc^{-2}$. Since both the
ALFALFA and xGASS do not provide resolved Hi images, in this work we use an
indirect method to estimate $D_{\rm HI}$ for individual galaxies. As first
reported by Broeils & Rhee (1997), there exists a tight correlation between
$D_{\rm HI}$ and $M_{\rm HI}$. With a sample of $\sim 500$ galaxies compiled
from literatures, Wang et al. (2016) revisited the $D_{\rm HI}-M_{\rm HI}$
relation and found

${\rm log}~{}D_{\rm HI}=(0.506\pm 0.003){\rm log}~{}M_{\rm HI}-(3.293\pm
0.009),$ (3)

which is very close to the one reported by Broeils & Rhee (1997).

With the large galaxy sample in hand, Wang et al. (2016) further quantified
the dispersion of the $D_{\rm HI}-M_{\rm HI}$ relation, finding $1\sigma\sim
0.06$ dex (or 14%). This scatter is independent of other galactic properties,
such as luminosity, $M_{\rm HI}$, Hi-to-luminosity ratio, etc. Remarkably,
galaxies across the Hubble types (from early-type, Sa-Sd, dIrr to even ultra-
diffuse galaxies) all follow a same $D_{\rm HI}-M_{\rm HI}$ relation (Wang et
al., 2016; Leisman et al., 2017). Since the galaxy samples of Wang et al.
(2016) and this work cover a similar range in redshift and $M_{\rm HI}$, in
this paper we use equation (3) to infer $D_{\rm HI}$ for the ALFALFA and xGASS
galaxies. We discuss in more detail in Appendix B to show that this
application should be valid.

Figure 4 shows $D_{\rm HI}/D_{\rm r,~{}25}$ as a function of $f_{\rm HI}$.
Inspired by Figure 3, we first divide the galaxies into 2 subsamples according
to $f_{\rm HI}$, with a demarcation line of log$f_{\rm HI}=-0.7$. For the Hi-
poor sample, $D_{\rm HI}/D_{\rm r,~{}25}$ is quite tightly correlated with
$f_{\rm HI}$, in the sense that more Hi-poor galaxies have a lower $D_{\rm
HI}/D_{\rm r,~{}25}$ ratio. For the Hi-rich sample, the correlation between
these two quantities significantly weakens. Note that there exists an up-
bending in the $(D_{\rm HI}/D_{\rm r,~{}25})-f_{\rm HI}$ relation at the high
$f_{\rm HI}$ end. This feature is due to selection effect. We have checked and
confirmed that a significant fraction of the galaxies with
$M_{\ast}<10^{9.0}M_{\sun}$ have similar $f_{\rm HI}$ but lower $D_{\rm
HI}/D_{\rm r,~{}25}$ compared to their massive counterparts. When galaxies
with $M_{\ast}<10^{9.0}M_{\sun}$ galaxies are included, the up-bending feature
shown at the high $f_{\rm HI}$ is absent.

At the low-$f_{\rm HI}$ end, Figure 4 should have missed some very low-$f_{\rm
HI}$ galaxies due to the Hi detection limits of xGASS and ALFALFA. In the
log$f_{\rm HI}<-1.0$ regime, the $D_{\rm HI}/D_{\rm r,~{}25}$ ratio is mainly
driven by $D_{\rm HI}$ (or $M_{\rm HI}$) since $D_{\rm r,~{}25}$ is no longer
sensitive to $f_{\rm HI}$ (see Figure 3). Therefore, the Hi-undetected
galaxies are expected to display in the extrapolation of the sequence defined
by the Hi-detected galaxies with log$f_{\rm HI}<-1.0$.

Figure 5: The correlation between $D_{\rm HI}/D_{\rm r,~{}25}$ and $M_{\ast}$
for the Hi-rich ALFALFA galaxies. Symbols are color coded by $M_{\rm
HI}/M_{\ast}$. The large solid circles represent the median log$(D_{\rm
HI}/D_{\rm r,~{}25})$ in that $M_{\ast}$ bin, with a bin size of
$\Delta$log$(M_{\ast}/M_{\sun})=0.2$. Error bars represent the 1$\sigma$
standard deviation. Figure 6: Top: the log($D_{\rm
r,~{}25}$)$-$log($M_{\ast}$) relation for the Hi-rich galaxies. Large symbols
represent the median log $D_{\rm r,~{}25}$ at each $M_{\ast}$ bin, with a bin
size of $\Delta$log$(M_{\ast}/M_{\sun})=0.2$. The dashed line represents the
best fit of the log($D_{\rm r,~{}25}$)$-$log($M_{\ast}$) relation. The small
panel show the size scatter distribution around the best fit relation, and the
red line shows the best fit gaussian function. Bottom: the log($D_{\rm
r,~{}25}$)$-$log($M_{\ast}$) relations for the Hi-poor galaxies, with a best
fit shown in solid line.

Since $D_{\rm HI}$ is inferred from $M_{\rm HI}$, Figure 4 thus also shows the
interplay between $M_{\ast}$, $D_{\rm r,~{}25}$ and $M_{\rm HI}$ similar to
Figure 2 and Figure 3. After transforming $M_{\rm HI}$ into $D_{\rm HI}$, one
can find that the gas disk shrinks more significantly compared to the stellar
disk once galaxies entering the gas-poor phase. This is likely due to the
different disk growth behaviors at different galaxy evolutional phases, as we
will argue below.

To investigate the Hi-to-stellar size ratio of Hi-rich galaxies more
specifically, we show $D_{\rm HI}/D_{\rm r,~{}25}$ as a function of $M_{\ast}$
in Figure 5. Interestingly, the Hi-to-stellar size ratio of these galaxies is
weakly correlated with $M_{\ast}$. With a sample of 108 galaxies observed by
the Westerbork Synthesis Radio Telescope, Broeils & Rhee (1997) found that
$D_{\rm HI}/D_{\rm r,~{}25}$ is only weakly correlated with luminosity. We
confirm the finding of Broeils & Rhee (1997) with a much larger sample and
emphasize that this trend only holds for Hi-rich galaxies.

At this point, it is worthy to revisit the size$-$mass relation separately for
the Hi-rich and Hi-poor galaxies. The results are shown in Figure 6. To be
consistent with the previous sections, we also use a second-order polynomial
to fit the size$-$mass relation. Overall, the size$-$mass relation of Hi-rich
galaxies has a shallower slope and a larger dispersion than the one of Hi-poor
galaxies. These features resemble those of the size$-$mass relations of star-
forming and quiescent galaxies (see Figure 3 of van der Wel et al. 2014). On
the one hand, this similarity is straightforward to interpret because the SFRs
and Hi richness of galaxies are correlated, in the sense that SFGs tend to be
more rich in Hi compared to quiescent galaxies. On the other hand, an Hi view
of the size$-$mass relation provides new insights in understanding the origin
of size scatter at fixed $M_{\ast}$ (also see Figure 3). In the next section,
we will show that the size scatter of SFGs is correlated with Hi richness,
rather than SFRs.

Figure 7: Top: the log(SFR)$-$log$M_{\ast}$ relation for xGASS galaxies. Hi-
rich and Hi-poor galaxies are indicated by solid circles and crosses,
respectively. The solid line indicate ridge line of the star formation main
sequence defined by Saintonge et al. (2016), while the dashed lines indicate
$\pm 0.3$ dex scatter around the main sequence. Bottom: the log($D_{\rm
r,~{}25}$)$-$log($M_{\ast}$) relation for the Hi-rich galaxies drawn form the
xGASS sample. Galaxies are color coded by their SFR offset from the star
formation main sequence ($\Delta SFMS$).

###  3.3 connections between Star formation rate, Hi fraction and Size

In a merger-free case, the size of a galaxy grows through in-situ star
formation. This naturally predicts a link between size growth and star
formation activity. To our knowledge, star formation is more physically
correlated with molecular gas, rather than Hi gas. Why the size of Hi-rich
galaxies strongly correlate with the Hi content? Is the size growth of Hi-rich
galaxies due to enhanced star formation activities?

To answer these questions, we first examine the distribution of xGASS galaxies
in the SFR$-M_{\ast}$ plane to explore the correlation between Hi-richness and
SFRs. The SFMS of xGASS galaxies is fitted by Saintonge et al. (2016), as
shown in the blue solid line of Figure 7. It can be seen that Hi-rich galaxies
distribute along the SFMS, i.e, they are mostly SFGs. In contrast, Hi-poor
galaxies distribute more broadly, and the majority of them distribute below
the SFMS. In low-mass regime of log$(M_{\ast}/M_{\sun})<10.2$, the Hi-poor
galaxies mainly distribute below the SFMS. At the high-mass end, a significant
fraction of the Hi-poor galaxies still distribute on the SFMS. This implies
that many massive SFGs may be short in Hi gas supply. Saintonge et al. (2016)
also argued that the flattening of the SFMS at the high mass end is primarily
due to the decrease of cold gas reservoirs.

Figure 8: Top: the log$(M_{\rm HI}/M_{\ast})-$log$M_{\ast}$ relation for xGASS
galaxies. Galaxies with $\Delta SFMS>-0.3$ and $\Delta SFMS<-0.3$ are
indicated by pink and black symbols, respectively. The solid line indicate
ridge line of the gas main sequence fitted from the galaxies with $\Delta
SFMS$¿-0.3, while the dashed lines indicate $\pm 0.4$ dex scatter around it.
Bottom: the log($D_{\rm r,~{}25}$)$-$log($M_{\ast}$) relation for the $\Delta
SFMS>-0.3$ galaxies drawn form the xGASS sample. Galaxies are color coded by
their offset from the Hi gas main sequence ($\Delta HIMS$).

To explore the correlation between size and SFRs, we divide the Hi-rich xGASS
galaxies into 3 subsamples according to their offset from the SFMS ridge line
($\Delta SFMS$, defined as $\Delta SFMS=\rm log(SFR/SFR_{\rm MS})$) and
examine their locations in the size$-$mass plane. We only include Hi-rich
galaxies since they are likely to live in an actively disk-growing phase, as
discussed above. The result is shown in the bottom panel of Figure 7. Clearly,
the location of galaxies in the size$-$mass plane is weakly dependent on
$\Delta SFMS$, i.e, there is little correlation between SFR and size at fixed
$M_{\ast}$. This result is broadly consistent with that reported in Lin et al.
(2020).

In the top panel of Figure 8, we show Hi fraction as a function of $M_{\ast}$.
Similar to the SFMS, we fit an Hi gas main sequence (HIMS) using xGASS
galaxies near the SFMS (with $\Delta SFMS$¿-0.3), as shown in the blue solid
line. It can be seen that galaxies with $\Delta SFMS$¿-0.3 generally have
higher Hi fraction than those with $\Delta SFMS$¡-0.3, but these two samples
have significantly overlap in the $M_{\rm HI}/M_{\ast}-M_{\ast}$ diagram.
Janowiecki et al. (2020) have investigated the HIMS based on the xGASS sample
and reported a similar feature. When investigating the molecular gas fraction
vs. stellar mass diagram, Janowiecki et al. (2020) found that the overlap is
significantly reduced (see their Fig.2 and Fig.3).

In the bottom panel of Figure 8, we examine the size$-$mass relation of SFGs
that with $\Delta SFMS>-0.3$. Symbols are color coded by the offset of
galaxies from the HIMS. As can be seen, at fixed $M_{\ast}$, Hi-rich SFGs tend
be more large in size. To conclude, at fixed $M_{\ast}$, the size scatter of
SFGs is more physically linked with their Hi content, rather than SFRs.

##  4 Discussion

###  4.1 size growth in the gas-rich phase

The central finding of this work is that the size of an Hi-rich galaxy is
closely linked to its Hi fraction, rather than SFR. Since the majority of star
formation occurs inside the stellar disk, a direct interpretation is that Hi-
rich galaxies have a more extended star-forming disk than the Hi-normal ones.
Previous works have studied the extension of star-forming disk through
ultraviolet (UV) observations and investigated its relation to the Hi
properties. Cortese et al. (2012) showed that the UV/optical size ratio of
galaxies is tightly correlated with their Hi fraction, in the sense that Hi-
rich galaxies tend to have larger UV/optical size ratios (see their Figure 6).
This is qualitatively consistent with our findings.

In an inside-out disk growth scenario, the materials that build a stellar disk
is primarily accreted from the surrounding environment. In this sense, the
properties of a galaxy would be to a large extent shaped by the properties of
the accreted material (such as specific angular momentum, accretion rate,
etc). Based on a suite of 30 cosmological magneto-hydrodynamical zoom
simulations, Grand et al. (2017) find that the sizes of galaxies are closely
linked to the angular momentum of halo material. Galaxies with the largest
disk size are produced by quiescent mergers that deposit high-angular momentum
material into the pre-existing disc, simultaneously increasing the spin of
dark matter and gas in the halo. In the observational work of Huang et al.
(2012), the authors studied the spin parameter $\lambda$ for the ALFALFA
galaxies, finding that Hi-rich galaxies tend to reside in high$-\lambda$
halos. A similar conclusion is reached by the recent work of Mancera Piña et
al. (2021), who show that gas-rich disk galaxies tend to have high specific
angular momentum. This can be interpreted since gas with high specific angular
momentum is more difficult to collapse and convert into stars, which naturally
results in a high $M_{\rm HI}/M_{\ast}$ ratio. We thus suggest that angular
momentum is the key parameter behind the relation between size and Hi
fraction.

Our study also suggests that cold gas accretion should exist in the majority
of Hi-rich galaxies. As shown in Figure 7, these galaxies are mostly normal
SFGs. If gas replenishment does not exist, the Hi disks of galaxies will
shrink due to continuous gas consumption via star formation. In addition,
without the acquisition of angular momentum from accreted materials (Stewart
et al., 2011, 2013), kinematic evolution will drive the remaining gas settling
towards galactic central regions, leading to a shrink of the gaseous disk. As
a result, $D_{\rm HI}/D_{\rm r,~{}25}$ will decrease towards lower $f_{\rm
HI}$. Figure 4 demonstrates that this is true at log$f_{\rm HI}<-0.7$,
suggesting that Hi-poor galaxies likely have their Hi supply shut down.

In the “bathtub” model (e.g., Bouché et al., 2010; Davé et al., 2011; Lilly et
al., 2013; Forbes et al., 2014; Peng & Maiolino, 2014), galaxies will approach
a quasi-equilibrium phase when cold gas replenishment, star formation and gas
outflows reach a balance. The weak correlation between $D_{\rm HI}/D_{\rm
r,~{}25}$ and $M_{\ast}$ at log$f_{\rm HI}>-0.7$ is likely a manifestation of
this phase. There is other observational evidence supporting that galaxies
will experience such a phase during their lifetime, for example, the existence
of a tight SFMS among SFGs from $z=0$ to $z=6$ (e.g., Brinchmann et al., 2004;
Noeske et al., 2007; Elbaz et al., 2007; Karim et al., 2011; Speagle et al.,
2014; Tasca et al., 2015). Galaxies on the ridge of the SFMS typically have
disk-like structures (Wuyts et al., 2011), suggesting that disk growth is
actively occurring in this phase.

###  4.2 size growth in the gas-poor phase

In the Hi-poor regime, we find that size is weakly correlated with Hi
fraction. As shown in Figure 7, galaxies that with log$f_{\rm HI}<-0.7$
typically distribute below the SFMS. Note that a significant fraction of
massive SFGs also have log$f_{\rm HI}<-0.7$, and we speculate many of them are
undergoing star formation quenching. By modeling the galaxy number density in
the NUV$-u$ color space, Lian et al. (2016) concluded that at
log($M_{\ast}/M_{\sun}$)$=$10.5 about 45% of the SFGs are undergoing
quenching. For Milky-Way mass SFGs, Pan et al. (2017) suggested that their
mass budget has been dominated by a quenched component, also implying on-going
quenching processes are actively taking place at the massive end. Although
many massive galaxies still appear “star-forming” at present, the quenching
progress will push them migrating onto the “dead sequence” with a time scale
of $2-4$ gigayears (Schawinski et al., 2014; Peng et al., 2015; Lian et al.,
2016; Hahn et al., 2017).

In the scenario of star formation quenching, Hi-poor galaxies grow their sizes
primarily via mergers. As shown in Catinella et al. (2010), Hi-poor galaxies
preferentially have bulge-dominated (early-type) morphologies. Shen et al.
(2003) show that the size$-$mass relation of SDSS early-type galaxies (ETGs)
is consistent with the assumption that they are the remnants of major mergers
of present-day disks. In this sense, the size scatter of ETGs should be more
physically related with their merger histories, rather than the current gas
content.

However, it should be noted that not every Hi-poor galaxy is necessarily to be
quenched or undergoing quenching. For example, numerical studies show that
SFGs in the lower envelope of the SFMS typically have lower gas fraction and
longer gas depletion time scale ($\tau_{\rm dep}$) than their counterparts on
the SFMS ridge (Tacchella et al., 2016). These galaxies will quench star
formation if their gas replenishment time scale ($\tau_{\rm rep}$) is longer
than $\tau_{\rm dep}$. In the case of $\tau_{\rm rep}<\tau_{\rm dep}$, they
will evolve back to the SFMS ridge and become gas rich again. It is difficult
to quantify the fraction of galaxies with $\tau_{\rm rep}<\tau_{\rm dep}$ in
the Hi-poor regime, but the size growth of such galaxies should also follow
the scenario discussed for the Hi-rich galaxies.

##  5 Summary

In this paper, we study the role of Hi gas in regulating size growth of local
galaxies using galaxy samples from the xGASS and ALFALFA survey. In the Hi-
rich regime, galaxies that are more rich in Hi tend to have larger sizes,
which we interpret as a natural consequence of the “inside-out” disk assembly.
This trend is absent in the Hi-poor regime, indicating that size growth is
barely affected by Hi gas when it has declined to a sufficiently low level. We
also study the relations between size, Hi-fraction and SFR, finding that size
is more intrinsically linked with Hi fraction, rather than SFR. The Hi-to-
stellar size ratio of Hi-rich galaxies is found to be weakly dependent on
$M_{\ast}$. We conclude that in the Hi-rich phase, size growth is primarily
achieved by star formation. The size dispersion at fixed $M_{\ast}$ is
probably driven by the angular momentum of the accreted materials.

We are grateful to the anonymous referee for useful suggestions that helped
improving the presentation of this paper. This work is supported by the
National Key Research and Development Program of China (2017YFA0402703), the
National Science Foundation of China (11773076, 12073078), and the Chinese
Academy of Sciences (CAS) through a China-Chile Joint Research Fund (CCJRF
1809) administered by the CAS South America Center for Astronomy (CASSACA). We
wish to thank the entire ALFALFA and xGASS team for providing reduced Hi
catalog for the community. Funding for SDSS-III has been provided by the
Alfred P. Sloan Foundation, the Participating Institutions, the National
Science Foundation, and the U.S. Department of Energy Office of Science. The
SDSS-III web site is http://www.sdss3.org/

##  Appendix A The reliability of SDSS-measured $R_{25}$

It is known that the automatic SDSS photometry pipeline is not reliable for
angularly large sources, which is mainly due to the problems in background
substraction and blending large galaxies into multiple sources (West et al.,
2010). Since in this work we are going to use $D_{\rm r,~{}25}=2R_{\rm
r,~{}25}$ as a size indicator, we first need to check whether $R_{\rm
r,~{}25}$ is reliably measured by the SDSS photometry pipeline.

To do this, we compare the $g-$band isophotal radii at the surface brightness
of 25.0 $\rm mag~{}arcsec^{-2}$ ($R_{\rm 25}$) provided by the SDSS pipeline
with those measured by Wang et al. (2013) for the Blue disk sample, which
contains 50 galaxies with log($M_{\ast}/M_{\sun})=[10.0,11.0]$. The result is
shown in Figure 9. As can be seen, the SDSS-measured $g-$band $R_{\rm 25}$ is
slightly smaller than those measured by the Wang’s pipeline
($R_{g,25,SDSS}\approx 0.9R_{g,25,Bluedisk})$. This is due to the systematics
of these two pipelines. Since the $g-$band $R_{\rm 25}$ is similar to that
measured in the $r-$band, we conclude that the SDSS $R_{\rm 25}$ measurement
is reliable at least for galaxies with $R_{\rm 25}<45.0$ arcsec.

Figure 10 shows the $r-$band $R_{\rm 25}$ as a function of redshift for the
ALFALFA and xGASS samples. As can be seen, the majority of the galaxies have
$R_{\rm 25}<45.0$ arcsec, i.e., they are not very extended sources. We thus
conclude that the $R_{\rm r,~{}25}$ provided by the SDSS database can be used
in our statistical study.

Figure 9: Top: The comparison between the $g-$band $R_{\rm 25}$ measured by
SDSS and those measured by Wang et al. (2013) for the Blue disk sample. The
solid line indicates 1:1. Bottom: The distribution of the ratio between SDSS
$R_{\rm 25}$ and Wang’s $R_{\rm 25}$.  Figure 10: The distribution of galaxies
in the apparent $R_{\rm r,~{}25}$ vs. redshift $z$ space.

##  Appendix B Do ALFALFA and xGASS galaxies follow the $D_{\rm HI}-M_{\rm
HI}$ relation?

The $D_{\rm HI}$ presented in this work is inferred from the $D_{\rm
HI}-M_{\rm HI}$ relation, rather than from direct measurements. To make our
analysis work, we first need to verify that the ALFALFA and xGASS Hi-detected
galaxies also follow the $D_{\rm HI}-M_{\rm HI}$ relation. Since spatially-
resolved Hi imaging data are not available for ALFALFA and xGASS, a directly
investigation of the $D_{\rm HI}-M_{\rm HI}$ relation for ALFALFA and xGASS is
not feasible. We thus use an indirect approach to access whether the ALFALFA
and xGASS galaxies follow the $D_{\rm HI}-M_{\rm HI}$ relation. We do this by
placing the ALFALFA and xGASS galaxies together with Wang et al. (2016) sample
in the $M_{\ast}-M_{\rm HI}$ and $M_{\ast}-D_{\rm r,~{}25}$ plane to examine
the parameter distributions of these samples. At fixed $M_{\ast}$, if ALFALFA
and xGASS Hi-detected galaxies were more Hi-deficient than Wang’s sample, then
it may be problematic to apply the $D_{\rm HI}-M_{\rm HI}$ relation to ALFALFA
and xGASS. If not, the ALFALFA and xGASS galaxies should also follow the the
same $D_{\rm HI}-M_{\rm HI}$ relation established by Wang et al. (2016).

To derive the stellar mass for Wang’s sample, we cross-identified the sample
galaxies with SDSS footprints, yielding a sample of 195 galaxies. We then
estimated the stellar masses of these galaxies using equation (2). In left and
right panel of Figure 11, we show the parameter distribution of these three
samples in the $M_{\ast}-D_{\rm r,~{}25}$ and $M_{\ast}-M_{\rm HI}$ space,
respectively. Overall, the xGASS galaxies cover a similar region as Wang’s
sample in both parameter spaces. It can be seen that at
log$(M_{\ast}/M_{\sun})=9.0-10.0$, the ALFALFA galaxies appear to be
systematically more gas rich than Wang’s sample. This is because the ALFALFA
survey is biased to Hi rich galaxies. At high-masses, these two samples cover
a similar parameter region, which is due to the fact that Wang’s sample
includes 23 galaxies that with unusually high Hi mass fraction drawn from the
Bluedisk project (Wang et al., 2013). Since Wang’s sample well follow the
$D_{\rm HI}-M_{\rm HI}$ relation, we are thus confident that the Hi-detected
galaxies from ALFALFA and xGASS should also follow this relation since they do
not appear to be more Hi-deficient.

Figure 11: Left: the distribution of ALFALFA, xGASS and Wang et al. (2016)’s
sample in the $M_{\ast}-D_{\rm r,~{}25}$ plane. Right: the $M_{\ast}-M_{\rm
HI}$ plane. At a given stellar mass, it is clear that ALFALFA detected
galaxies appear more gas rich than Wang’s sample. At
log$(M_{\ast}/M_{\sun})>9.0$, both xGASS and Wang’s sample cover a similar
parameter region. Given this, the application of the $D_{\rm HI}-M_{\ast}$
relation of Wang et al. (2016) to ALFALFA and xGASS should be valid. Figure
12: Left:the log($D_{\rm r,~{}25}$)$-$log($M_{\ast}$) relation for face-on Hi-
rich galaxies. The large symbols represent the median log$(D_{\rm r,~{}25})$
in that $M_{\ast}$ bin, with a bin size of
$\Delta$log$(M_{\ast}/M_{\sun})=0.2$. Error bars represent the 1$\sigma$
standard deviation. We use a linear function to fit this relation, and the
best-fit result is shown in the solid line. Small symbols are color coded by
log$(M_{\rm HI}/M_{\ast})$. The small panel shows the size scatter around the
best fit relation. Right: the log($D_{\rm 5}$)$-$log($M_{\ast}$) relation.

##  Appendix C Replacing $D_{\rm r,~{}25}$ with $D_{\rm 5}$

In this work, the size of a galaxy is defined by $D_{\rm r,~{}25}$. Taking the
variation of mass-to-light ratio ($M/L$) among galaxies into account, a size
definition with fixed stellar mass surface density could be more physically
meaningful (Trujillo et al., 2020). We have selected a sample of $\sim$3700
ALFALFA galaxies with log$f_{\rm HI}>-0.7$ to compute their stellar mass
density profiles based on the SDSS public photometric data. The sample
galaxies are selected to have minor-to-major axis ratio of $b/a>0.7$, $R_{\rm
25}>10\farcs 0$ and $0.01<z<0.06$. The SDSS pipeline measures the surface
brightness profile of galaxies in five bands ($ugriz$) in a series of circular
annuli of fixed angular size. For each galaxy, we first corrected the annular
photometry for Galactic extinction and $k-$correction for each annulus. Then a
spline was fit to the cumulative light profile (to conserve flux), which was
then differentiated to derive the surface brightness profile interpolated over
a grid of angular sizes in $0\farcs 1$. The surface brightness profile was
then converted into stellar mass surface density profile using the mass-to-
light ratio expressed in equation (2).

We select a new sample rather than use the original sample presented in the
main body of the paper to do this test because: (1) only at log$f_{\rm
HI}>-0.7$ galaxy size shows a clear dependence on Hi richness; (2) it is
difficult to measure the stellar mass density profile for edge-on galaxies.
For this sample, we use $D_{\rm 5}$ to define the size of a galaxy, where
$D_{\rm 5}$ is the diameter at which the stellar mass surface density reaches
$5~{}M_{\sun}\rm pc^{-2}$. In Trujillo et al. (2020), the authors use $D_{\rm
1}$ as a size indicator. We do not use $D_{\rm 1}$ because our photometric
data is shallower than those used in Trujillo et al. (2020), which may result
in large uncertainties in the $D_{\rm 1}$ measurement. In left and right panel
of Figure 12, we show the $M_{\ast}-D_{\rm r,~{}25}$ and $M_{\ast}-D_{\rm 5}$
relations for our sample galaxies, respectively. To be consistent with
Trujillo et al. (2020), we use a linear function to fit both relations. As can
be seen, the $M_{\ast}-D_{\rm 5}$ relation is tighter and steeper than the
$M_{\ast}-D_{\rm r,~{}25}$ relation. A similar feature is also reported by
Trujillo et al. (2020). At fixed stellar mass, Hi-rich galaxies also tend to
have larger $D_{\rm 5}$, although this feature is less evident than that
revealed in the $M_{\ast}-D_{\rm r,~{}25}$ relation. We thus conclude that the
positive correlation between size and Hi richness for Hi-rich galaxies can not
be fully due to the $M/L$ effect.

## References

  * Bell et al. (2003) Bell, E. F., McIntosh, D. H., Katz, N., & Weinberg, M. D. 2003, ApJS, 149, 289
  * Bouché et al. (2010) Bouché, N., Dekel, A., Genzel, R., et al. 2010, ApJ, 718, 1001
  * Broeils & Rhee (1997) Broeils, A. H., & Rhee, M.-H. 1997, A&A, 324, 877
  * Brown et al. (2015) Brown, T., Catinella, B., Cortese, L., et al. 2015, MNRAS, 452, 2479
  * Brinchmann et al. (2004) Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151
  * Catinella et al. (2010) Catinella, B., Schiminovich, D., Kauffmann, G., et al. 2010, MNRAS, 403, 683
  * Catinella et al. (2018) Catinella, B., Saintonge, A., Janowiecki, S., et al. 2018, MNRAS, 476, 875
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
  * Chiappini et al. (1997) Chiappini, C., Matteucci, F., & Gratton, R. 1997, ApJ, 477, 765
  * Chen et al. (2020) Chen, X., Wang, J., Kong, X., et al. 2020, MNRAS, 492, 2393
  * Cortese et al. (2012) Cortese, L., Boissier, S., Boselli, A., et al. 2012, A&A, 544, A101
  * Dale et al. (2016) Dale, D. A., Beltz-Mohrmann, G. D., Egan, A. A., et al. 2016, AJ, 151, 4.
  * Davé et al. (2011) Davé, R., Finlator, K., & Oppenheimer, B. D. 2011, MNRAS, 416, 1354
  * de Vaucouleurs (1948) de Vaucouleurs, G. 1948, Annales d’Astrophysique, 11, 247
  * Elbaz et al. (2007) Elbaz, D., Daddi, E., Le Borgne, D., et al. 2007, A&A, 468, 33
  * Fang et al. (2013) Fang, J. J., Faber, S. M., Koo, D. C., & Dekel, A. 2013, ApJ, 776, 63
  * Frankel et al. (2019) Frankel, N., Sanders, J., Rix, H.-W., et al. 2019, ApJ, 884, 99
  * Forbes et al. (2014) Forbes, J. C., Krumholz, M. R., Burkert, A., & Dekel, A. 2014, MNRAS, 438, 1552
  * Giovanelli et al. (2005) Giovanelli, R., Haynes, M. P., Kent, B. R., et al. 2005, AJ, 130, 2598
  * Graham (2019) Graham, A. W. 2019, PASA, 36, e035
  * Grand et al. (2017) Grand, R. J. J., Gómez, F. A., Marinacci, F., et al. 2017, MNRAS, 467, 179
  * Hahn et al. (2017) Hahn, C., Tinker, J. L., & Wetzel, A. 2017, ApJ, 841, 6
  * Huang et al. (2012) Huang, S., Haynes, M. P., Giovanelli, R., & Brinchmann, J. 2012, ApJ, 756, 113
  * Huang et al. (2014) Huang, S., Haynes, M. P., Giovanelli, R., et al. 2014, ApJ, 793, 40
  * Hopkins et al. (2009) Hopkins, P. F., Hernquist, L., Cox, T. J., et al. 2009, ApJ, 691, 1424
  * Janowiecki et al. (2020) Janowiecki, S., Catinella, B., Cortese, L., et al. 2020, MNRAS, 493, 1982
  * Karim et al. (2011) Karim, A., Schinnerer, E., Martínez-Sansigre, A., et al. 2011, ApJ, 730, 61
  * Kauffmann (2015) Kauffmann, G. 2015, MNRAS, 450, 618
  * Yim & van der Hulst (2016) Yim, K. & van der Hulst, J. M. 2016, MNRAS, 463, 2092
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231
  * Lange et al. (2015) Lange, R., Driver, S. P., Robotham, A. S. G., et al. 2015, MNRAS, 447, 2603
  * Leisman et al. (2017) Leisman, L., Haynes, M. P., Janowiecki, S., et al. 2017, ApJ, 842, 133
  * Lian et al. (2016) Lian, J., Yan, R., Zhang, K., & Kong, X. 2016, ApJ, 832, 29
  * Lilly et al. (2013) Lilly, S. J., Carollo, C. M., Pipino, A., Renzini, A., & Peng, Y. 2013, ApJ, 772, 119
  * Lin et al. (2020) Lin, L., Faber, S. M., Koo, D. C., et al. 2020, ApJ, 899, 93
  * Mancera Piña et al. (2021) Mancera Piña, P. E., Posti, L., Pezzulli, G., et al. 2021, arXiv:2107.02809
  * Marsters (2005) Masters,K. L., 2005, PhD thesis, Cornell Univ.
  * Noeske et al. (2007) Noeske, K. G., Weiner, B. J., Faber, S. M., et al. 2007, ApJ, 660, L43
  * Pan et al. (2015) Pan, Z., Li, J., Lin, W., et al. 2015, ApJ, 804, L42
  * Pan et al. (2017) Pan, Z., Zheng, X., & Kong, X. 2017, ApJ, 834, 39
  * Peng & Maiolino (2014) Peng, Y.-j., & Maiolino, R. 2014, MNRAS, 443, 3643
  * Peng et al. (2015) Peng, Y., Maiolino, R., & Cochrane, R. 2015, Nature, 521, 192
  * Pérez et al. (2013) Pérez, E., Cid Fernandes, R., González Delgado, R. M., et al. 2013, ApJ, 764, L1
  * Pichon et al. (2011) Pichon, C., Pogosyan, D., Kimm, T., et al. 2011, MNRAS, 418, 2493
  * Pohlen & Trujillo (2006) Pohlen, M. & Trujillo, I. 2006, A&A, 454, 759
  * Saintonge et al. (2016) Saintonge, A., Catinella, B., Cortese, L., et al. 2016, MNRAS, 462, 1749.
  * Schawinski et al. (2014) Schawinski, K., Urry, C. M., Simmons, B. D., et al. 2014, MNRAS, 440, 889
  * Scott et al. (2017) Scott, N., Brough, S., Croom, S. M., et al. 2017, MNRAS, 472, 2833
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525
  * Shen et al. (2003) Shen, S., Mo, H. J., White, S. D. M., et al. 2003, MNRAS, 343, 978
  * Speagle et al. (2014) Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15
  * Stewart et al. (2011) Stewart, K. R., Kaufmann, T., Bullock, J. S., et al. 2011, ApJ, 738, 39
  * Stewart et al. (2013) Stewart, K. R., Brooks, A. M., Bullock, J. S., et al. 2013, ApJ, 769, 74
  * Taylor et al. (2011) Taylor, E. N., Hopkins, A. M., Baldry, I. K., et al. 2011, MNRAS, 418, 1587
  * Tacchella et al. (2016) Tacchella, S., Dekel, A., Carollo, C. M., et al. 2016, MNRAS, 457, 2790.
  * Tasca et al. (2015) Tasca, L. A. M., Le Fèvre, O., Hathi, N. P., et al. 2015, A&A, 581, A54
  * Trujillo et al. (2020) Trujillo, I., Chamba, N., & Knapen, J. H. 2020, MNRAS, 493, 87
  * van der Wel et al. (2014) van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2014, ApJ, 788, 28
  * Wang et al. (2011) Wang, J., Kauffmann, G., Overzier, R., et al. 2011, MNRAS, 412, 1081
  * Wang et al. (2013) Wang, J., Kauffmann, G., Józsa, G. I. G., et al. 2013, MNRAS, 433, 270
  * Wang et al. (2016) Wang, J., Koribalski, B. S., Serra, P., et al. 2016, MNRAS, 460, 2143
  * Wang et al. (2018) Wang, J., Zheng, Z., D’Souza, R., et al. 2018, MNRAS, 479, 4292
  * West et al. (2010) West, A. A., Garcia-Appadoo, D. A., Dalcanton, J. J., et al. 2010, AJ, 139, 315
  * Wuyts et al. (2011) Wuyts, S., Förster Schreiber, N. M., van der Wel, A., et al. 2011, ApJ, 742, 96
  * Yıldız et al. (2017) Yıldız, M. K., Serra, P., Peletier, R. F., Oosterloo, T. A., & Duc, P.-A. 2017, MNRAS, 464, 329
  * York et al. (2000) York, D. G., Adelman, J., Anderson, J. E., Jr., et al. 2000, AJ, 120, 1579


# Exploiting Emotional Dependencies with Graph Convolutional Networks for
Facial Expression Recognition

Panagiotis Antoniadis, Panagiotis Paraskevas Filntisis, Petros Maragos
School of E.C.E., National Technical University of Athens, Greece This work
was not supported by any organization.

###### Abstract

Over the past few years, deep learning methods have shown remarkable results
in many face-related tasks including automatic facial expression recognition
(FER) in-the-wild. Meanwhile, numerous models describing the human emotional
states have been proposed by the psychology community. However, we have no
clear evidence as to which representation is more appropriate and the majority
of FER systems use either the categorical or the dimensional model of affect.
Inspired by recent work in multi-label classification, this paper proposes a
novel multi-task learning (MTL) framework that exploits the dependencies
between these two models using a Graph Convolutional Network (GCN) to
recognize facial expressions in-the-wild. Specifically, a shared feature
representation is learned for both discrete and continuous recognition in a
MTL setting. Moreover, the facial expression classifiers and the valence-
arousal regressors are learned through a GCN that explicitly captures the
dependencies between them. To evaluate the performance of our method under
real-world conditions we perform extensive experiments on the AffectNet and
Aff-Wild2 datasets. The results of our experiments show that our method is
capable of improving the performance across different datasets and backbone
architectures. Finally, we also surpass the previous state-of-the-art methods
on the categorical model of AffectNet.

##  I INTRODUCTION

Facial expressions are one of the most powerful nonverbal ways for human
beings to convey their emotional state [11]. Facial expression recognition
(FER) has been a topic of study for decades due to its potential applications
in various fields including human-computer interaction, digital entertainment,
advertisement, health care and intelligent robot systems [10], [3], [43],
[46], [20]. However, recognizing facial expressions in the wild is still very
challenging due to variations, occlusions and the ambiguity of human emotion
[37], [55].

While the cultural and ethnic background of a person can affect his expressive
style, Ekman indicated that humans perceive certain basic emotions in the same
way regardless of their culture [18], [17]. These six universal facial
expressions (happiness, sadness, surprise, fear, disgust and anger) constitute
the categorical model. Contempt was subsequently added as one of the basic
emotions [42]. Due to its direct and intuitive definition of facial
expressions, the categorical model is used in the majority of FER algorithms
([4], [1], [28], [26], etc) and large-scale databases (AffectNet [44], RAF-DB
[36], SFEW [13], FER-2013 [22], EmotionNet [19], etc). However, the
subjectivity and ambiguity of restricting human emotion to discrete categories
result in large intra-class variations and small inter-class differences.

Recently, the dimensional model proposed by Russell [51] has gained a lot of
attention where emotion is described using a set of 2 latent dimensions that
are valence (how pleasant or unpleasant a feeling is) and arousal (how likely
is the person to take action under the emotional state). Another dimension
called dominance is used sometimes to know whether the person is controlling
the situation or not. Since a continuous representation can distinguish
between subtly different displays of affect and encode small changes in the
intensity, some recent algorithms [47], [7], [31] and databases (Aff-Wild [60]
and AFEW-VA [32]) have utilized the dimensional model for uncontrolled FER.
Even so, predicting a 2-dimensional continuous value instead of a category
increases a lot the task complexity and lacks intuitiveness.

Altogether, the categorical and the dimensional model have their respective
benefits and drawbacks [28]. Therefore, recent studies try to leverage both
representations, along with Action Units (AU) detection, through multi-task
learning (MTL) [7], [31], [58]. However, the strong dependence between the
categorical and the dimensional model is not fully exploited when they only
share a common feature representation. Fig. 1 illustrates this relation using
the validation set of AffectNet.

Figure 1: Distribution of the basic expressions in the VA space using the
validation set of AffectNet that illustrates the emotional dependencies
between the categorical and the dimensional model. The basic emotions are
located around the neutral emotion that appears when valence and arousal are
close to zero.

In multi-label classification, there has been a lot of research on how to
properly capture and explore the correlation between labels [38], [35]. In
[9], Chen proposed ML-GCN model for multi-label image recognition that
explicitly learns the labels correlation by generating the object classifiers
via a Graph Convolutional Network (GCN) [30]. Inspired by this architecture,
we propose Emotion-GCN a novel GCN based MTL framework for FER in the wild.
The main idea in this paper is to to generate dependent expression classifiers
and valence-arousal (VA) regressors though a GCN based mapping function
instead of learning them as separate parameter vectors. The generated vectors
are then applied to an extracted image representation to enable end-to-end
training. Hence, the dependence between the categorical and the dimensional
emotion models is explicitly captured through both a shared feature
representation and the dependent classifiers and regressors. Experiments on
the AffectNet and Aff-Wild2 datasets indicate that Emotion-GCN increases the
performance across different datasets and backbone networks managing to
achieve state-of-the-art results on the categorical model of AffectNet.

The rest of this paper is organized as follows. Section II presents the recent
work on FER systems and MTL. Section III describes the deep learning method
that this paper proposes. Section IV discusses the experimental results,
providing a clear view of the accuracy improvements introduced by our method.
This section also includes a description of the data used during the
experiments. Finally, Section V summarizes the key aspects of our work and
concludes the paper.

##  II RELATED WORK

Image-based FER has been extensively studied for many years. Typically, a FER
system consists of three stages: face detection, feature extraction and
classification. Traditional approaches tend to conduct FER by using
handcrafted features, such as Local Binary Patterns [53], Gabor wavelets [40],
[45] and Histogram of Oriented Gradients [5]. While there is a lot of
intuition behind these features and their performance on several lab-
controlled databases is impressive, they lack generalizability and sufficient
learning capacity [37].

Later, many in-the-wild facial expression databases were developed that
enabled the research of FER in more challenging environments. Deep
Convolutional Neural Networks (CNNs) have achieved promising recognition
performance by learning powerful high-level features [61], [15], [59], [29],
[24], [14]. In [39] and [57] region-based attention networks were designed for
pose and occlusion aware FER, where the regions are either cropped from
landmark points or fixed positions. Facial Motion Prior Networks were proposed
in [8] that generate a facial mask to focus on facial muscle moving regions.
In [25] deep Siamese neural networks equipped with a supervised loss function
were used to reduce the high intra-class variation of the task. In [21] deep
and handcrafted features were combined and a local learning framework was used
at test time based on SVMs. Recently, Self-Cure network [56] was proposed that
suppresses the uncertainties caused by ambiguous expressions and the
suggestiveness of the annotators.

MTL was first explored in [6] based on the idea that learning complementary
tasks in parallel while using a shared representation improves the
generalization performance. Since then, MTL has been used in many areas of
computer vision, such as classification and detection [52] or geometry and
regression tasks [16]. In the facial analysis domain, a multi-purpose
algorithm for seven face-related tasks was proposed in [50]. FATAUVA-Net [7]
performs sequential facial attribute recognition, AU detection, and VA
estimation on videos. In [31] a holistic framework was proposed that jointly
learns three facial behavior tasks (recognition of basic emotions, VA
estimation and AUs detection) and two simple strategies were used for coupling
the tasks during training. Closer to our work, in [58] a two-level attention
with a two-stage MTL (2Att-2Mt) framework was proposed for facial emotion
estimation on static images. However, the work was focused mainly on the
estimation of VA and the dependencies between the emotion representations was
not further explored.

Apart from MTL, research on the dependencies between the categorical and the
dimensional emotion representations is limited. Recently, in CAKE [28] the
link between the two representations was explored and a 3-dimensional
representation of emotion learned in a multi-domain fashion was proposed.

##  III OUR PROPOSED METHOD

In this section the architecture of the proposed network, the dependent
classifiers and regressors and the MTL setting for recognition are introduced.
The overall architecture of the proposed Emotion-GCN is shown in Fig. 2.

Figure 2: Overall architecture of our Emotion-GCN model for FER in the wild.
Our graph contains the seven expression labels and the two VA dimensions that
are connected to each other based on the adjacency matrix $\mathbf{A}$.
Stacked GCNs are learned over the graph to map the word embeddings of the
nodes $\mathbf{Z}$ into a set of dependent classifiers and regressors,
$\mathbf{W^{c}}$ and $\mathbf{W^{r}}$ respectively. These vectors are then
applied to an image representation $\mathbf{x}$ extracted from the input image
$\mathbf{I}$ via a DenseNet $\mathbf{f_{cnn}}$ followed by a global max-
pooling function $\mathbf{f_{gmp}}$. The whole network is trained end-to-end
for both basic expression classification ($\mathbf{L^{c}}$) and VA regression
($\mathbf{L^{r}}$).

###  III-A Feature extraction

Given an input image $\mathbf{I}$ of size $227\times 227$ pixels, we obtain
$1024\times 7\times 7$ feature maps using a Dense Convolutional Network
(DenseNet) [27]. Each layer in DenseNet obtains additional inputs from all
preceding layers and passes on its own feature-maps to all subsequent layers
to ensure maximum information flow between layers. Then, a global max-pooling
function is applied to get the feature vector $\mathbf{x}$ of the image:

$\mathbf{x}=f_{gmp}(f_{cnn}(\mathbf{I}))\in\mathbb{R}^{D}$ (1)

where $f_{cnn}$ corresponds to the convolution layers of the DenseNet,
$f_{gmp}$ to the the global max-pooling function and $D=1024$.

###  III-B Dependent classifiers and regressors

Given the feature vector $\mathbf{x}$ of the image, we want to learn one
classifier for each facial expression (classification task) and one regressor
for each dimension in the VA space (regression task). Each classifier acts as
a weight vector that is multiplied with the feature $\mathbf{x}$ and then
passes though a softmax activation function assigning a probability score to
its emotion category. In parallel, each regressor follows the same procedure
without passing through an activation function since its output is a predicted
continuous value for either valence or arousal. We denote the expression
classifiers by $\mathbf{W^{c}}\in\mathbb{R}^{7\times D}$ and the VA regressors
by $\mathbf{W^{r}}\in\mathbb{R}^{2\times D}$ where each row of the matrices
corresponds to a classifier and a regressor respectively. Traditionally, these
vectors are optimized as individual parameters of the deep learning network.
To capture the correlation between the categorical and the dimensional model,
we propose to learn these vectors using a GCN that maps their word embeddings
to dependent classifiers and regressors retaining the shared information
between the two tasks.

####  III-B1 GCN Overview

Generally, the goal of a GCN is to learn a function $f$ on a graph $G=(V,E)$
that takes as input a feature description for each node of the graph
$\mathbf{H^{l}}\in\mathbb{R}^{n\times d}$ ($n=|V|$) and a correlation matrix
$\mathbf{A}\in\mathbb{R}^{n\times n}$ and produces new node-level features
$\mathbf{H^{l+1}}\in\mathbb{R}^{n\times d^{\prime}}$. The update rule is
formulated as follows:

$\mathbf{H^{l+1}}=h(\mathbf{\hat{A}}\mathbf{H^{l}}\mathbf{W^{l}})$ (2)

where $\mathbf{W^{l}}\in\mathbb{R}^{d\times d^{\prime}}$ is the weight matrix
for the $l$-th GCN layer, $\mathbf{\hat{A}}$ is the normalized version of
matrix $\mathbf{A}$ such that all rows sum to one and $h(\cdot)$ is LeakyRELU
[41]. For more information on GCN we refer readers to [30].

In our case, the nodes of the graph correspond to the seven expression labels
and the two VA dimensions i.e. V = {Neutral, Happy, Sad, Surprise, Fear,
Disgust, Anger, Valence, Arousal}. Hence, the input is
$\mathbf{Z}\in\mathbb{R}^{9\times d}$ that contains the word embedding of each
node ($d$ is the dimensionality of the embeddings). Each GCN layer $l$ takes
the node representations from the previous layer $\mathbf{H^{l}}$ as inputs
and outputs new node representations $\mathbf{H^{l+1}}$ using (2). Finally,
the output representation of the last layer $\mathbf{W}\in\mathbb{R}^{9\times
D}$ contains the dependent classifiers and regressors. The first seven rows of
the output matrix $\mathbf{W}$ constitute the classification part
$\mathbf{W^{c}}$ and the rest two the regression part $\mathbf{W^{r}}$.

####  III-B2 Adjacency Matrix Design

According to the update rule of a GCN (2) the feature of a node in the graph
is the weighted sum of its own feature and the adjacent nodes’ features. Since
the purpose behind the use of a GCN is to exploit the dependencies between the
categorical and the dimensional model, we should design the adjacency matrix
$\mathbf{A}$ to this direction. As before, we assume that the first seven rows
of $\mathbf{A}$ correspond to the basic expressions and the last two are the
continuous dimensions i.e. valence and arousal. Since we deal with a multi-
task and not a multi-label recognition problem, we are interested only in the
correlation between the categorical and the dimensional model. Hence, we set
the other pairs of $\mathbf{A}$ to zero except for the diagonal to allow self-
loops. Also, we take the absolute value of the correlation to ignore its type
(positive or negative) and focus on its amplitude. The correlation matrix
$\mathbf{A}\in\mathbb{R}^{9\times 9}$ can be written as:

$A_{ij}=\begin{cases}1,&\ \text{if}\ \ i=j\\\ |c_{ij}|,&\ \text{if}\ \
i\in\textit{Cat}\ \land\ j\in\textit{Dim}\\\ |c_{ij}|,&\ \text{if}\ \
j\in\textit{Cat}\ \land\ i\in\textit{Dim}\\\ 0,&\ \text{else}\end{cases}$ (3)

where Cat and Dim are the set of indices of the categorical and the
dimensional labels respectively. As a correlation metric, we use the
Spearman’s rank correlation coefficient [54] that for two variables
$\mathbf{X}=\\{x_{1},...,x_{N}\\}$ and $\mathbf{Y}=\\{y_{1},...,y_{N}\\}$ is
defined as:

$c_{xy}=\frac{\sum_{k=1}^{N}x_{k,r}y_{k,r}}{\sqrt{\sum_{k=1}^{N}x_{k,r}^{2}\sum_{k=1}^{N}y_{k,r}^{2}}}$
(4)

where $x_{k,r}$ and $y_{k,r}$ are the rank transformation of the initial
values $x_{k}$ and $y_{k}$. Following the ideas in ML-GCN [9], we use a
threshold $\tau$ to filter the noisy edges as follows:

$A^{\prime}_{ij}=\begin{cases}1,&\ \text{if}\ A_{ij}\geq\tau\\\ 0,&\
\text{if}\ A_{ij}<\tau\end{cases}$ (5)

where $\tau=0.1$ to enable the propagation of information between weakly
correlated nodes. As shown in [34], the learned features of each node may be
over-smoothed and become indistinguishable when a binary correlation matrix is
used. To alleviate the over-smoothing problem, we apply the re-weighted scheme
of ML-GCN that is defined as:

$A^{{}^{\prime\prime}}_{ij}=\begin{cases}(p/\sum_{\begin{subarray}{c}j=1\\\
i\neq j\end{subarray}}^{9}A^{\prime}_{ij})\times A^{\prime}_{ij},&\ \text{if}\
i\neq j\\\ 1-p,&\ \text{if}\ i=j\end{cases}$ (6)

where the variable $p$ determines the weights assigned to a node itself and
its adjacent nodes. We set $p=0.7$ to increase the influence of the
neighborhood nodes. Fig. 3 shows the output adjacency matrix
$\mathbf{A^{{}^{\prime\prime}}}$ of Emotion-GCN using the training set of
AffectNet. As we expected, pairs like happy-valence and surprise-arousal are
strongly connected since the presence of the one usually denotes the presence
of the other (Fig. 1).

Figure 3: Adjacency matrix of our Emotion-GCN.

###  III-C Multi-task learning setting

After we generate the dependent matrices $\mathbf{W^{c}}$, $\mathbf{W^{r}}$
and the feature vector $\mathbf{x}$, we perform FER simultaneously on the
categorical and the dimensional model. For the classification task, the
predicted scores are computed as:

$\hat{y}_{i}=\frac{e^{\mathbf{w_{i}^{c}}\cdot\mathbf{x}}}{\sum_{k=1}^{7}e^{\mathbf{w_{k}^{c}}\cdot\mathbf{x}}}$
(7)

where $\hat{y}_{i}$ is the probability of emotion $i$ and $\mathbf{w_{i}^{c}}$
is the classifier of emotion $i$. The network is trained using a weighted
version of the traditional categorical cross entropy loss $\mathbf{L^{c}}$
since the dataset is highly imbalanced [44]. In other words, the network is
penalized more for misclassifying samples from under-represented classes than
from well-represented classes as follows:

$L^{c}=-\sum_{i=1}^{7}w_{i}\ y_{i}\log(\hat{y}_{i})$ (8)

$w_{i}=\frac{f_{i}}{f_{min}}$ (9)

where $y_{i}=1$ if class $i$ is the ground truth expression, $f_{i}$ is the
number of samples of the $i$-th class and $f_{min}$ is the number of samples
in the most under-represented class i.e. Disgust. For the regression task, the
predicted values are computed as:

$\hat{p}_{i}=\mathbf{w^{r}_{i}}\cdot\mathbf{x}$ (10)

where $\hat{p}_{1},\hat{p}_{2}$ are the predicted values of valence and
arousal respectively and $\mathbf{w^{r}_{1}},\mathbf{w^{r}_{2}}$ correspond to
the VA regressors. Most studies in the literature use the Mean Squared Error
(MSE) as a loss function to train regression models. However, more and more
works on emotion recognition use a Correlation-based loss function to measure
the agreement between the true emotion dimension and the predicted emotion
degree [31], [23]. The Concordance Correlation Coefficient (CCC) is often used
since it takes the bias into Pearson’s correlation coefficient [2]. It is
defined as:

$\rho_{c}=\frac{2s_{xy}}{s_{x}^{2}+s_{y}^{2}+(\bar{x}-\bar{y})^{2}}$ (11)

where $s_{x}$ and $s_{y}$ denote the variance of the predicted and ground
truth values respectively, $\bar{x}$ and $\bar{y}$ are the corresponding mean
values and $s_{xy}$ is the respective covariance value. The range of CCC is
from -1 (perfect disagreement) to 1 (perfect agreement). Hence, in our case we
define $\mathbf{L^{r}}$ as:

$L^{r}=1-\frac{\rho_{v}+\rho_{a}}{2}$ (12)

where $\rho_{v}$ and $\rho_{a}$ are the respective CCC of valence and arousal.
The overall loss function of our network training is defined as
$L=L^{c}+L^{r}$.

##  IV EXPERIMENTS

In this section, we present the experimental evaluation of the proposed
Emotion-GCN model (Fig. 2). First, we briefly describe the facial expression
datasets that we use and the details of our implementations. Then, an ablation
study of our models follows and finally a comparison with the current state-
of-the-art methods.

###  IV-A Datasets

####  IV-A1 AffectNet

AffectNet is by far the largest database of facial expressions that provides
both categorical and VA annotations. It contains more than one million facial
images collected from the Internet by querying three major search engines
using 1,250 emotion-related keywords in six different languages [44]. About
half of the retrieved images are manually annotated for the presence of seven
discrete facial expressions (categorical model) and the intensity of valence
and arousal (dimensional model). It is a very challenging database as it
contains images of people from different races and ethnic groups as well as
high variety in the background, lighting, pose, point of view, etc. To be
consistent with the majority of past studies [61], [39], [8], [28], [14],
[21], [25], we exclude the emotion of contempt and train our models on the
images with neutral and the 6 basic emotions (approximately 280,000 training
samples). Since the test set of AffectNet is not publicly available, we
evaluate our approaches on the validation set that contains 500 images for
each emotion. The mean class accuracy is used as the evaluation metric for the
classification task because the validation set is balanced and the CCC for the
regression task.

####  IV-A2 Aff-Wild2

Aff-Wild2 is the first ever database annotated for valence-arousal estimation,
action unit detection and basic expression classification [31]. It consists of
548 videos collected from YouTube and shows both subtle and extreme human
behaviours in real-world settings. We trained and evaluated our models on a
subset of the database that contains only the frames with both categorical and
VA annotations, as required by our method. As proposed by the authors of Aff-
Wild2 a weighted average between the F1 score and mean class accuracy is used
as the evaluation metric for the classification task
($0.67\times\textit{F1}+0.33\times\textit{Acc}$).

###  IV-B Implementation details

Before training the network, the face region of each image is cropped using
the provided bounding box by the database and scaled to $227\times 227$
pixels. Also, we perform face alignment based on the position of the eyes to
obtain a normalized representation of each face. Specifically, we compute the
location of the eyes by taking the mean value of the six detected landmarks in
each eye. Then, we rotate the image by $\theta$ that is defined as
$\theta=tan^{-1}\left(\frac{r_{y}-l_{y}}{r_{x}-l_{x}}\right)$, where $(r_{x},\
r_{y})$ and $(l_{x},\ l_{y})$ are the coordinates of the right and the left
eye respectively. To augment the data, six types of augmentation techniques
are used (flip, rotation and changes in brightness, contrast, hue and
saturation) [25]. The GCN module consists of two layers defined with output
dimensionality of 512 and 1024 respectively. For the word representations, we
use 300-dimensional GloVe [49] trained on the Wikipedia dataset. The networks
are trained for 10 epochs using a batch size of 35 and a learning rate of
0.001. Stochastic Gradient Descent is adopted as the optimization algorithm
with a momentum of 0.9 and PyTorch [48] is used as our deep learning
framework.

###  IV-C Ablation Study

First, we investigate the performance of four different networks on the
categorical and the dimensional model of affect to present the improvements
that Emotion-GCN introduces.

####  IV-C1 Categorical Model

TABLE I: Performance of our models on the categorical model of AffectNet and Aff-Wild2. Method | AffectNet | Aff-Wild2
---|---|---
Single-task | 64.37 | 45.06
Multi-task + MSE | 64.8 | 43.1
Multi-task + CCC | 65.69 | 43.33
Emotion-GCN | 66.46 | 48.92

We trained (i) a single-task network for discrete FER using the weighted CE
loss of (8). Then, two multi-task networks were trained for discrete and
continuous FER using a weighted CE loss for the classification task and a (ii)
MSE or (iii) CCC loss (12) for the regression task. Finally, we trained (iv)
the proposed Emotion-GCN model that generates dependent classifiers and
regressors using a 2-layer GCN as presented in Fig. 2. Table I shows the
results of our experiments. In AffectNet we can see that learning to predict
the VA values as an additional task in a MTL framework increases the accuracy
by 1.32% since the shared representation improves the generalization of the
network. Also, in agreement with recent studies, using the CCC loss for the
regression task increases the classification accuracy of the model verifying
that in MTL a correlation-based regression loss performs better. Finally,
Emotion-GCN achieves a total accuracy of 66.46%, since the dependent
classifiers manage to effectively capture the dependencies between the two
emotion representations. The confusion matrices for these models are shown in
Fig. 4. It can be seen that our proposed method increases the accuracy for
most classes while the single-task network performs better only in the neutral
class. Similarly, in Aff-Wild2 the total evaluation metric increases by 3.86%
indicating that the benefits of Emotion-GCN generalize in more datasets. To
investigate how different choices for the parameters of the GCN affect the
performance, we perform additional experiments in Table II. We observe that
the number of GCN layers ($L$) is the most crucial parameter since using only
one layer decreases the accuracy a lot.

(a) Single-task

(b) Multi-task + MSE

(c) Multi-task + CCC

(d) Emotion-GCN

Figure 4: Confusion matrices of our models on the validation set of AffectNet.
TABLE II: Classification accuracy of Emotion-GCN on the categorical model of
AffectNet using different values for $\tau$, $p$ and $L$ (number of GCN
layers). In the first table $p=0.7$ and in the second table $\tau=0.1$.

$L$ $\tau$ | 0 | 0.05 | 0.1 | 0.15
---|---|---|---|---
1 | 53.09 | 43.85 | 53.71 | 54.94
2 | 65.29 | 65.31 | 66.46 | 65.69
3 | 65.77 | 64.94 | 64.74 | 65.74

$L$ $p$ | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8
---|---|---|---|---|---|---|---
1 | 49.6 | 50.94 | 52.57 | 52.66 | 52.86 | 53.71 | 59.31
2 | 65.71 | 66.14 | 65.29 | 66.29 | 66.09 | 66.46 | 65.06
3 | 65.29 | 65.54 | 65.89 | 66.09 | 65.54 | 64.74 | 65.83

####  IV-C2 Dimensional Model

TABLE III: Performance of our models on the dimensional model of AffectNet and Aff-Wild2. Method | AffectNet | Aff-Wild2
---|---|---
CCC-V | CCC-A | CCC-V | CCC-A
Single-task | 0.761 | 0.628 | 0.416 | 0.501
Multi-task + MSE | 0.752 | 0.572 | 0.435 | 0.378
Multi-task + CCC | 0.768 | 0.651 | 0.408 | 0.481
Emotion-GCN | 0.767 | 0.649 | 0.457 | 0.514

For the evaluation of our models on the VA space, the networks (ii), (iii) and
(iv) are the same since they are trained for both discrete and continuous FER.
Additionally, a single-task network for VA regression was trained using the
CCC loss of (12). Table III shows the performance evaluation of our
experiments on the dimensional model. In AffectNet we can see that the multi-
task network with the CCC loss improves the regression performance along with
the classification one verifying that both tasks benefit from the shared
feature representation. Actually, we observe that the increase in the
performance of arousal prediction is higher that that of valence. This is due
to the fact that most emotions have positive value of arousal (Fig. 1) and a
simultaneous emotion classification provides more useful information in the
task of arousal regression. Finally, our GCN based approach achieves similar
performance with that of the multi-task network on the dimensional model. In
Aff-Wild2, our proposed model surpasses the performance of both the single-
task and multi-task networks. Overall, our method presents significant
improvements in both the categorical and the dimensional model. In Fig. 5 we
present some positive and negative results of our Emotion-GCN model on
AffectNet.

Figure 5: Predictions of our models on samples of AffectNet. The first column
indicates the ground truth values. The second and the third column present the
predictions of the multi-task network trained with CCC loss and our Emotion-
GCN respectively. To examine which network better exploits the emotional
dependencies, we selected samples where the predictions of the networks on the
dimensional model are close. In the first four samples our Emotion-GCN model
successfully recognizes the depicted emotion while the multi-task network
fails indicating that our proposed model effectively captured the dependencies
presented in the VA space (Fig. 1). However, there are cases where our network
fails (last two samples).

###  IV-D Visualization

To further analyze the effectiveness of our approach, we investigate the
similarity between the dependent classifiers and regressors. In Fig. 6 the
cosine similarity of the learned vectors on AffectNet by our single-task
networks, our best multi-task network and our Emotion-GCN are presented. As we
can see, learning a shared representation for both tasks through MTL slightly
increases the similarity between the expression classifiers and the VA
regressors. Our proposed method increases their similarity even more in
consistence with the dependencies presented in Fig. 1. Specifically, the
regressor of valence comes closer to the classifier of happy and the regressor
of arousal closer to the classifiers of anger, disgust, fear and surprise.
These emotions appear in regions of VA space where the values of valence or
arousal are high. Therefore, the proposed network has successfully captured
the dependence between the categorical and the dimensional emotion
representation.

We also observe an increase in the similarity between the pairs valence-
surprise and happy-surprise while their respective nodes are not adjacent in
the graph (Fig. 3). These dependencies are successfully captured by our
network since in a 2-layer GCN a node incorporates information from a 2-hop
neighborhood [12]. To further examine whether the dependence between the
emotions of happiness and surprise is reasonable, we compute their co-
occurrence in a multi-label emotion dataset. The EMOTIC dataset [33] is a
collection of images of people in unconstrained environments annotated
according to their apparent emotional states. Each person is annotated for 26
discrete categories, with multiple labels assigned to each image. About 25% of
the samples labeled as happy are labeled as surprise too that indicates that
these two emotions are strongly related indeed.

(a) Single-task

(b) Multi-task + CCC

(c) Emotion-GCN

Figure 6: Visualization of the cosine similarity between the learned
classifiers and regressors by our models on AffectNet.

###  IV-E Dependence between classifiers

Inspired by the fact that there are dependencies between the basic emotions,
we trained a similar GCN based network but we enabled this time the direct
propagation of information between the seven emotion classifiers. Instead of
setting their correlation to zero as in Emotion-GCN, we compute the
conditional probability matrix of the basic emotions like in ML-GCN. The model
achieves a total accuracy of 66.23% on the categorical model of AffectNet.
Despite the intuition behind this approach, the recognition performance is
slightly lower than that of our proposed method. We believe that by enabling
the propagation of information between the basic emotions the dependence
between the categorical and dimensional model is ignored. Also, we deal with a
single-label recognition task and the possible benefits of this approach are
more suitable in a multi-label recognition task.

###  IV-F Comparison with the State of the Art

In Table IV, we compare the performance of our Emotion-GCN model with several
state-of-the-art methods for FER on the categorical model of AffectNet.
Regarding Aff-Wild2, we only used the subset which contains both categorical
and continuous annotations, while methods in the bibliography typically report
results on the whole dataset, making direct comparisons not possible. In
IPA2LT [61], the authors designed LTNet to discover the latent truths from the
human annotations and the machine annotations trained from different FER
datasets. In gACNN [39] and OADN [14] attention networks were proposed for
occlusion aware FER. Facial Motion Prior Network in [8] generates a facial
mask so as to focus on facial muscle moving regions. In [21] deep (CNNs) and
handcrafted features (BOVW) were combined and in [25] a deep Siamese network
along with a supervised loss function was used to reduce the intra-class
variation of the task. Closer to our work, CAKE [28] proposed a 3-dimensional
representation of emotion learned in a multi-domain fashion. Our Emotion-GCN
model considerably outperforms these recent state-of-the-art methods,
achieving an accuracy of 66.46%.

In [24] BReG-NeXt achieved state-of-the-art accuracy in 8-way classification
on AffectNet (including contempt) by introducing a residual-based network
architecture. To investigate the ability of Emotion-GCN to generalize over
other model architectures as well, we replaced our DenseNet backbone network
with BReG-NeXt using the publicly available
code111https://github.com/behzadhsni/BReG-NeXt. Since the provided code does
not include the data preprocessing strategy, we followed our preprocessing and
training pipeline described before (same with DenseNet), which explains the
different results compared to [24]. As we can see in Table V, our proposed
model outperforms the single-task and multi-task models when using BReG-NeXt
as the backbone network indicating that Emotion-GCN generalizes across
different model architectures as well.

TABLE IV: Comparison with state-of-the-art methods on AffectNet (7-way classification). Method | Accuracy
---|---
IPA2LT [61] | 57.31
gACNN [39] | 58.78
Facial Motion Prior Network [8] | 61.52
CAKE [28] | 61.7
OADN [14] | 61.89
CNNs and BOVW + global SVM [21] | 63.31
Siamese [25] | 64
Emotion-GCN (ours) | 66.46
TABLE V: Performance of Emotion-GCN using BReG-NeXt as the backbone network on AffectNet. Method | Network | Accuracy
---|---|---
Single-task | BReG-NeXt | 60.49
Multi-task+CCC | BReG-NeXt | 61.14
Emotion-GCN | BReG-NeXt | 61.94

##  V CONCLUSION

In this work, a novel GCN based MTL framework is proposed for in-the-wild FER.
Specifically, our Emotion-GCN model learns a shared feature representation for
both discrete and continuous expression recognition to exploit the
dependencies between the categorical and the dimensional model of affect. To
further capture these dependencies, the expression classifiers and the VA
regressors are learned though a GCN that maps their word representation to
dependent vectors inspired by recent work in multi-label image recognition.
Experimental results on AffectNet, the largest facial expression database,
have demonstrated that our Emotion-GCN outperforms the performance of the
recent state-of-the-art methods for discrete FER.

## References

  * [1] D. Acharya, Z. Huang, D. Pani Paudel, and L. Van Gool.  Covariance pooling for facial expression recognition.  In Proc. CVPRW, 2018.
  * [2] B. T. Atmaja and M. Akagi.  Evaluation of error and correlation-based loss functions for multitask learning dimensional speech emotion recognition.  arXiv preprint arXiv:2003.10724, 2020.
  * [3] P. M. Blom, S. Bakkes, C. Tan, S. Whiteson, D. Roijers, R. Valenti, and T. Gevers.  Towards personalised gaming via facial expression recognition.  In Proc. AIIDE, 2014.
  * [4] J. Cai, Z. Meng, A. S. Khan, Z. Li, J. O’Reilly, and Y. Tong.  Island loss for learning discriminative features in facial expression recognition.  In Proc. FG, 2018.
  * [5] P. Carcagnì, M. Del Coco, M. Leo, and C. Distante.  Facial expression recognition and histograms of oriented gradients: a comprehensive study.  SpringerPlus, 4, 2015.
  * [6] R. Caruana.  Multitask learning.  Machine learning, 28, 1997.
  * [7] W.-Y. Chang, S.-H. Hsu, and J.-H. Chien.  Fatauva-net: An integrated deep learning framework for facial attribute recognition, action unit detection, and valence-arousal estimation.  In Proc. CVPRW, 2017.
  * [8] Y. Chen, J. Wang, S. Chen, Z. Shi, and J. Cai.  Facial motion prior networks for facial expression recognition.  In Proc. VCIP, 2019.
  * [9] Z.-M. Chen, X.-S. Wei, P. Wang, and Y. Guo.  Multi-label image recognition with graph convolutional networks.  In Proc. CVPR, 2019.
  * [10] C. A. Corneanu, M. O. Simón, J. F. Cohn, and S. E. Guerrero.  Survey on rgb, 3d, thermal, and multimodal approaches for facial expression recognition: History, trends, and affect-related applications.  IEEE Trans. on PAMI, 38, 2016.
  * [11] C. Darwin.  The expression of the emotions in man and animals.  University of Chicago press, 2015.
  * [12] T. Derr, Y. Ma, and J. Tang.  Signed graph convolutional networks.  In Proc. ICDM, 2018.
  * [13] A. Dhall, R. Goecke, S. Lucey, and T. Gedeon.  Static facial expression analysis in tough conditions: Data, evaluation protocol and benchmark.  In Proc. ICCVW, 2011.
  * [14] H. Ding, P. Zhou, and R. Chellappa.  Occlusion-adaptive deep network for robust facial expression recognition.  In Proc. IJCB, 2020.
  * [15] H. Ding, S. K. Zhou, and R. Chellappa.  Facenet2expnet: Regularizing a deep face recognition net for expression recognition.  In Proc. FG, 2017.
  * [16] D. Eigen and R. Fergus.  Predicting depth, surface normals and semantic labels with a common multi-scale convolutional architecture.  In Proc. ICCV, 2015.
  * [17] P. Ekman.  Strong evidence for universals in facial expressions: a reply to russell’s mistaken critique.  Psychological bulletin, 1994.
  * [18] P. Ekman and W. V. Friesen.  Constants across cultures in the face and emotion.  Journal of personality and social psychology, 17, 1971.
  * [19] C. Fabian Benitez-Quiroz, R. Srinivasan, and A. M. Martinez.  Emotionet: An accurate, real-time algorithm for the automatic annotation of a million facial expressions in the wild.  In Proc. CVPR, 2016.
  * [20] P. P. Filntisis, N. Efthymiou, P. Koutras, G. Potamianos, and P. Maragos.  Fusing body posture with facial expressions for joint recognition of affect in child–robot interaction.  Robotics and Automation Letters, 4, 2019.
  * [21] M.-I. Georgescu, R. T. Ionescu, and M. Popescu.  Local learning with deep and handcrafted features for facial expression recognition.  IEEE Access, 7, 2019.
  * [22] I. J. Goodfellow, D. Erhan, P. L. Carrier, A. Courville, M. Mirza, B. Hamner, W. Cukierski, Y. Tang, D. Thaler, D.-H. Lee, et al.  Challenges in representation learning: A report on three machine learning contests.  In ICONIP, 2013.
  * [23] J. Han, Z. Zhang, M. Schmitt, M. Pantic, and B. Schuller.  From hard to soft: Towards more human-like emotion recognition by modelling the perception uncertainty.  In Proc. ACM international conference on Multimedia, 2017.
  * [24] B. Hasani, P. S. Negi, and M. Mahoor.  Breg-next: Facial affect computing using adaptive residual networks with bounded gradient.  IEEE Trans. on Affective Computing, 2020.
  * [25] W. Hayale, P. Negi, and M. Mahoor.  Facial expression recognition using deep siamese neural networks with a supervised loss function.  In Proc. FG, 2019.
  * [26] P. Hu, D. Cai, S. Wang, A. Yao, and Y. Chen.  Learning supervised scoring ensemble for emotion recognition in the wild.  In Proc. ICMI, 2017.
  * [27] G. Huang, Z. Liu, L. Van Der Maaten, and K. Q. Weinberger.  Densely connected convolutional networks.  In Proc. CVPR, 2017.
  * [28] C. Kervadec, V. Vielzeuf, S. Pateux, A. Lechervy, and F. Jurie.  Cake: Compact and accurate k-dimensional representation of emotion.  CoRR, 2018.
  * [29] B.-K. Kim, H. Lee, J. Roh, and S.-Y. Lee.  Hierarchical committee of deep cnns with exponentially-weighted decision fusion for static facial expression recognition.  In Proc. ICMI, 2015.
  * [30] T. N. Kipf and M. Welling.  Semi-supervised classification with graph convolutional networks.  In Proc. ICLR, 2017.
  * [31] D. Kollias and S. Zafeiriou.  Expression, affect, action unit recognition: Aff-wild2, multi-task learning and arcface.  BMVC, 2019.
  * [32] J. Kossaifi, G. Tzimiropoulos, S. Todorovic, and M. Pantic.  Afew-va database for valence and arousal estimation in-the-wild.  Image and Vision Computing, 65, 2017.
  * [33] R. Kosti, J. M. Alvarez, A. Recasens, and A. Lapedriza.  Emotic: Emotions in context dataset.  In Proc. CVPRW, 2017.
  * [34] Q. Li, Z. Han, and X.-M. Wu.  Deeper insights into graph convolutional networks for semi-supervised learning.  In Proc. AAAI Conference on Artificial Intelligence, 2018.
  * [35] Q. Li, M. Qiao, W. Bian, and D. Tao.  Conditional graphical lasso for multi-label image classification.  In Proc. CVPR, 2016.
  * [36] S. Li and W. Deng.  Reliable crowdsourcing and deep locality-preserving learning for unconstrained facial expression recognition.  IEEE Trans. on Image Processing, 28, 2018.
  * [37] S. Li and W. Deng.  Deep facial expression recognition: A survey.  IEEE Trans. on Affective Computing, 2020.
  * [38] X. Li, F. Zhao, and Y. Guo.  Multi-label image classification with a probabilistic label enhancement model.  In Proc. UAI, 2014.
  * [39] Y. Li, J. Zeng, S. Shan, and X. Chen.  Occlusion aware facial expression recognition using cnn with attention mechanism.  IEEE Trans. on Image Processing, 28, 2018.
  * [40] C. Liu and H. Wechsler.  Gabor feature based classification using the enhanced fisher linear discriminant model for face recognition.  IEEE Trans. on Image processing, 11, 2002.
  * [41] A. L. Maas, A. Y. Hannun, and A. Y. Ng.  Rectifier nonlinearities improve neural network acoustic models.  In Proc. ICML, 2013.
  * [42] D. Matsumoto.  More evidence for the universality of a contempt expression.  Motivation and Emotion, 16, 1992.
  * [43] D. McDuff, R. El Kaliouby, J. F. Cohn, and R. W. Picard.  Predicting ad liking and purchase intent: Large-scale analysis of facial responses to ads.  IEEE Trans. on Affective Computing, 6, 2014.
  * [44] A. Mollahosseini, B. Hasani, and M. H. Mahoor.  Affectnet: A database for facial expression, valence, and arousal computing in the wild.  IEEE Trans. on Affective Computing, 10, 2017.
  * [45] S. Moore and R. Bowden.  Local binary patterns for multi-view facial expression recognition.  Computer vision and image understanding, 115, 2011.
  * [46] G. Muhammad, M. Alsulaiman, S. U. Amin, A. Ghoneim, and M. F. Alhamid.  A facial-expression monitoring system for improved healthcare in smart cities.  IEEE Access, 5, 2017.
  * [47] M. A. Nicolaou, H. Gunes, and M. Pantic.  Continuous prediction of spontaneous affect from multiple cues and modalities in valence-arousal space.  IEEE Trans. on Affective Computing, 2, 2011.
  * [48] A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, and A. Lerer.  Automatic differentiation in pytorch.  In NIPS-W, 2017.
  * [49] J. Pennington, R. Socher, and C. D. Manning.  Glove: Global vectors for word representation.  In Proc. EMNLP, 2014.
  * [50] R. Ranjan, S. Sankaranarayanan, C. D. Castillo, and R. Chellappa.  An all-in-one convolutional neural network for face analysis.  In Proc. FG, 2017.
  * [51] J. A. Russell.  A circumplex model of affect.  Journal of personality and social psychology, 39, 1980.
  * [52] P. Sermanet, D. Eigen, X. Zhang, M. Mathieu, R. Fergus, and Y. LeCun.  Overfeat: Integrated recognition, localization and detection using convolutional networks.  CoRR, 2014.
  * [53] C. Shan, S. Gong, and P. W. McOwan.  Facial expression recognition based on local binary patterns: A comprehensive study.  Image and vision Computing, 27, 2009.
  * [54] C. Spearman.  The proof and measurement of association between two things.  International journal of epidemiology, 1961.
  * [55] M. F. Valstar, B. Jiang, M. Mehu, M. Pantic, and K. Scherer.  The first facial expression recognition and analysis challenge.  In Proc. FG, 2011.
  * [56] K. Wang, X. Peng, J. Yang, S. Lu, and Y. Qiao.  Suppressing uncertainties for large-scale facial expression recognition.  In Proc. CVPR, 2020.
  * [57] K. Wang, X. Peng, J. Yang, D. Meng, and Y. Qiao.  Region attention networks for pose and occlusion robust facial expression recognition.  IEEE Trans. on Image Processing, 29, 2020.
  * [58] W. Xiaohua, P. Muzi, P. Lijuan, H. Min, J. Chunhua, and R. Fuji.  Two-level attention with two-stage multi-task learning for facial emotion recognition.  Journal of Visual Communication and Image Representation, 62, 2019\.
  * [59] Z. Yu and C. Zhang.  Image based static facial expression recognition with multiple deep network learning.  In Proc. ICMI, 2015.
  * [60] S. Zafeiriou, D. Kollias, M. A. Nicolaou, A. Papaioannou, G. Zhao, and I. Kotsia.  Aff-wild: valence and arousal’in-the-wild’challenge.  In Proc. CVPRW, 2017.
  * [61] J. Zeng, S. Shan, and X. Chen.  Facial expression recognition with inconsistently annotated datasets.  In Proc. ECCV, 2018.