# Quiescent Galaxies in the 3D-HST Survey: Spectroscopic confirmation of a
large number of galaxies with relatively old stellar populations at $z\sim 2$

Katherine E. Whitaker11affiliation: Astrophysics Science Division, Goddard
Space Flight Center, Code 665, Greenbelt, MD 20771, USA 88affiliation: NASA
Postdoctoral Program Fellow , Pieter G. van Dokkum22affiliation: Department of
Astronomy, Yale University, New Haven, CT 06520, USA , Gabriel
Brammer33affiliation: European Southern Observatory, Alonso de Córdova 3107,
Casilla 19001, Vitacura, Santiago, Chile , Ivelina G. Momcheva22affiliation:
Department of Astronomy, Yale University, New Haven, CT 06520, USA , Rosalind
Skelton22affiliation: Department of Astronomy, Yale University, New Haven, CT
06520, USA , Marijn Franx44affiliation: Sterrewacht Leiden, Leiden University,
NL-2300 RA Leiden, The Netherlands , Mariska Kriek55affiliation: Department of
Astronomy, University of California, Berkeley, CA 94720, USA , Ivo
Labbé44affiliation: Sterrewacht Leiden, Leiden University, NL-2300 RA Leiden,
The Netherlands , Mattia Fumagalli44affiliation: Sterrewacht Leiden, Leiden
University, NL-2300 RA Leiden, The Netherlands , Britt F.
Lundgren66affiliation: Department of Astronomy, University of Wisconsin,
Madison, WI 53706, USA , Erica J. Nelson22affiliation: Department of
Astronomy, Yale University, New Haven, CT 06520, USA , Shannon G.
Patel44affiliation: Sterrewacht Leiden, Leiden University, NL-2300 RA Leiden,
The Netherlands , Hans-Walter Rix77affiliation: Max Planck Institut fur
Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany

###### Abstract

Quiescent galaxies at $z\sim 2$ have been identified in large numbers based on
rest-frame colors, but only a small number of these galaxies have been
spectroscopically confirmed to show that their rest- frame optical spectra
show either strong Balmer or metal absorption lines. Here, we median stack the
rest-frame optical spectra for 171 photometrically-quiescent galaxies at
$1.4<z<2.2$ from the 3D-HST grism survey. In addition to H$\beta$
($\lambda$4861$\mathrm{\AA}$), we unambiguously identify metal absorption
lines in the stacked spectrum, including the G-band
($\lambda$4304$\mathrm{\AA}$), Mg I ($\lambda$5175$\mathrm{\AA}$), and Na I
($\lambda$5894$\mathrm{\AA}$). This finding demonstrates that galaxies with
relatively old stellar populations already existed when the universe was $\sim
3$ Gyr old, and that rest-frame color selection techniques can efficiently
select them. We find an average age of $1.3^{+0.1}_{-0.3}$ Gyr when fitting a
simple stellar population to the entire stack. We confirm our previous result
from medium-band photometry that the stellar age varies with the colors of
quiescent galaxies: the reddest 80% of galaxies are dominated by metal lines
and have a relatively old mean age of $1.6^{+0.5}_{-0.4}$ Gyr, whereas the
bluest (and brightest) galaxies have strong Balmer lines and a spectroscopic
age of $0.9^{+0.2}_{-0.1}$ Gyr. Although the spectrum is dominated by an
evolved stellar population, we also find [OIII] and H$\beta$ emission.
Interestingly, this emission is more centrally concentrated than the continuum
with $L_{\mathrm{OIII}}=1.7\pm 0.3\times 10^{40}$ erg s-1, indicating residual
central star formation or nuclear activity.

###### Subject headings:

galaxies: evolution — galaxies: formation — galaxies: high-redshift

††slugcomment: Accepted for publication in the Astrophysical Journal Letters

##  1\. Introduction

In the nearby universe, the most massive galaxies almost always have little
ongoing star-formation and red colors that exhibit a remarkably small
intrinsic scatter (e.g., Bower et al., 1992). These red, quiescent galaxies
form a well-defined color-mass relation known as the “red sequence”.
Photometric studies of large, representative samples of galaxies based on
broad-band and medium-band photometry have pushed the detection of the red
sequence out to $z\sim 2$ (e.g., Williams et al., 2009; Whitaker et al., 2011;
Brammer et al., 2011; Nicol et al., 2011). However, only a small fraction of
these distant quiescent galaxies have been spectroscopically confirmed (e.g.,
Cimatti et al., 2004; Daddi et al., 2005; Kriek et al., 2006; Cimatti et al.,
2008; Kriek et al., 2008, 2009; van Dokkum & Brammer, 2010; van de Sande et
al., 2011, 2012; Onodera et al., 2012; Toft et al., 2012; Bezanson et al.,
2013).

Interestingly, most distant quiescent galaxies with high quality rest-frame
optical spectra appear to exhibit young ages, showing strong Balmer absorption
lines (e.g., van de Sande et al., 2012; Bezanson et al., 2013). This could
imply that the bulk of high redshift quiescent galaxies were recently
quenched. However, most spectroscopic studies in the literature are biased
towards the brightest (and consequently youngest) of such galaxies. Therefore,
spectroscopic ages for a representative sample of quiescent galaxies remains
elusive.

Owing to the near-infrared (NIR) slitless spectroscopic capabilities provided
by the Wide-Field Camera 3 (WFC3) on the Hubble Space Telescope (HST), it is
now possible to obtain low-resolution spectroscopy of a mass-limited sample of
quiescent galaxies at $z\sim 2$. In this letter, we stack the spectra of 171
quiescent galaxies at $z\sim 2$ and demonstrate that they have absorption
features indicative of evolved stellar populations. Furthermore, with ages
derived from the stacked grism spectra, we are in a unique position to test
what drives the color spread of quiescent galaxies for the first time.

We assume a $\Lambda$CDM cosmology with $\Omega_{M}$=0.3,
$\Omega_{\Lambda}$=0.7, and $H_{0}$=70 km s-1 Mpc-1 throughout the paper. All
magnitudes are given in the AB system.

##  2\. Data

The 3D-HST treasury program (Brammer et al., 2012a), a 248–orbit NIR
spectroscopic survey with the HST/WFC3 G141 grism, provides spatially resolved
low-resolution spectra of all objects in five well-studied extra-galactic
fields to a $5\sigma$ depth for continuum magnitudes of
$H_{\mathrm{F140W}}\sim 23$. 3D-HST has targeted the AEGIS, COSMOS, GOODS-S
and UDS fields, as well as incorporated publicly-available data in the GOODS-N
field (GO:11600; PI:Weiner). WFC3 imaging for all five fields is available
from the CANDELS survey (Grogin et al., 2011; Koekemoer et al., 2011). With a
wavelength range of 1.10$\mu$m$<\lambda<$1.65$\mu$m, the prominent age-
dependent H$\beta\lambda$4861 and Mg$\lambda$5175 absorption features in
quiescent galaxies at $z\sim 2$ are observable. The first-order dispersion of
the G141 grism is 46 $\mathrm{\AA}$/pixel ($\mathrm{R}\sim 130$) with a
spatial resolution of 0.′′12, sampled with 0.′′06 pixels; as the spectra have
high spatial resolution and low spectral resolution, the line width almost
exclusively reflects the size of the galaxy in the dispersion direction.

The sample used in this letter is selected from WFC3-selected photometric
catalogs generated from the HST ACS and WFC3 images of the CANDELS and 3D-HST
survey fields, along with ancillary ground-based optical and NIR images and
mid-IR images with Spitzer/IRAC using a methodology similar to that described
by Whitaker et al. (2011). A full description of these catalogs is beyond the
scope of this letter and will be described in a forthcoming paper (Skelton et
al. in prep). We extract a two-dimensional flux-calibrated WFC3/G141 grism
spectrum for every object in the photometric catalog with
$H_{\mathrm{F140W}}<24$. Details of the reduction and extraction of the grism
spectra, including accounting for the contamination of overlapping objects,
are provided by Brammer et al. (2012a, b).

To determine the galaxy redshifts, we first compute a purely photometric
redshift from the photometry, using the EAZY code (Brammer et al., 2008). We
then fit the full two-dimensional grism spectrum separately with a combination
of the continuum template taken from the EAZY fit and a single emission-line-
only template with fixed line ratios taken from the SDSS composite star-
forming galaxy spectrum of Dobos et al. (2012). The final grism redshift,
z_grism, is determined on a finely-sampled redshift grid with the photometry-
only redshift probability distribution function used as a prior. This method
is somewhat more flexible than that originally described by Brammer et al.
(2012a), but the redshift precision is similar with $\sigma\sim 0.0035(1+z)$.
Finally, rest-frame colors and stellar population parameters are computed from
the photometry with the EAZY and FAST (Kriek et al., 2009) codes,
respectively, with the redshift fixed to z_grism and assuming a Chabrier
(2003) initial mass function.

##  3\. Sample Selection

A standard method for discriminating high redshift quiescent galaxies from
star-forming galaxies is selecting on the rest-frame $U$–$V$ and $V$–$J$
colors (e.g., Labbé et al., 2005; Wuyts et al., 2007; Williams et al., 2009;
Bundy et al., 2010; Cardamone et al., 2010; Whitaker et al., 2011; Brammer et
al., 2011; Patel et al., 2012); quiescent galaxies have strong Balmer/4000Å
breaks, characterized by red $U$–$V$ colors and relatively blue $V$–$J$
colors.

Following the definition of Whitaker et al. (2012), our quiescent selection
box is shown in Figure 1, with the larger 3D-HST parent sample at $1.4<z<2.2$
shown in grey-scale. Using the dotted line in Figure 1
($(U-V)=-1.25\times(V-J)+2.85$), we further divide our quiescent sample into
“younger” (blue) and “older” (red) galaxies (Section 5). The nature of the
stellar populations of galaxies with rest-frame colors close to the
quiescent/star-forming division is not clear and we therefore restrict our
analysis to be more conservative by 0.08 mag (excluding those galaxies with
high transparency in Figure 1). We select galaxies at
$1.4<z_{\mathrm{grism}}<2.2$ with stellar masses log(M${}_{\star})>10.5$ M⊙,
further restricting our sample to require $>75\%$ wavelength coverage of the
grism spectrum with $<50\%$ of the pixels flagged as bad. The grism spectra of
quiescent galaxies do not have strong emission line features and therefore the
redshifts may not be reliable beyond certain magnitude limits; we additionally
require that we detect the continuum at $\gtrsim 5\sigma$ per resolution
element in the 1D spectrum ($H_{\mathrm{F140W}}<22.8$), removing the faintest
20% of galaxies.

Figure 1.— Rest-frame $U$–$V$ and $V$–$J$ color selection of quiescent
galaxies, further separating the sample into two bins of younger (blue) and
older (red) quiescent galaxies in Section 5. All galaxies with
$H_{\mathrm{F140W}}<24$ are shown in grayscale.

##  4\. Quiescent Stacks

Although the individual galaxies are too faint to discern spectral features,
by stacking the spectra of our final sample of 171 galaxies we can achieve the
necessary signal-to-noise ratio to robustly identify absorption features in a
large well-defined high redshift sample for the first time. To normalize the
spectra, we shift to the rest-frame and interpolate the flux values to a new
wavelength grid with 10$\mathrm{\AA}$ width bins. Next, we fit the continuum
in each spectrum with a third order polynomial, masking regions around
prominent absorption features. After dividing by the continuum, we determine
the median flux value at each wavelength bin and smooth with a boxcar of
20$\mathrm{\AA}$. The grey-scale error bars are derived from 100 bootstrap
iterations of the stacking analysis.

We show the stacked spectrum in Figure 2. In the top panel, we see many clear
absorption features: G-band ($\lambda 4304\mathrm{\AA}$) blended with
H$\gamma$($\lambda 4341\mathrm{\AA}$), Mg I ($\lambda 5175\mathrm{\AA}$), Fe I
($\lambda 5269\mathrm{\AA}$), and Na I ($\lambda 5894\mathrm{\AA}$). These
metal lines are indicative of an evolved stellar population. We caution that
systematic residuals in the blue part of the spectrum are of similar strength
as the emission lines. The residuals do not seem to be caused by photon noise,
and show up in all stacked spectra. We have not identified the cause of these
residuals. We cannot exclude errors in the stellar population synthesis
models, although that is unlikely in this well-studied wavelength range.

Figure 2.— Median rest-frame stack of 171 quiescent galaxies (top) and the
residuals after subtracting the best-fit model (red) with an age of 1.25 Gyr
(bottom), using the grism redshifts. The orange model shows an additional 5%
enhancement of centrally-concentrated H$\beta$ and [OIII] emission. The
accuracy of the grism redshifts enables us to resolve absorption features at
$z\sim 2$ with high signal-to-noise for the first time.

We perform a least-squares minimization using the Vazdekis et al. (2010)
models to find the best-fit age of the stack. The Vazdekis et al. (2010)
models provide moderately high resolution SEDs computed assuming an
instananeous burst, Solar metallicity and a Kroupa (2001) universal initial
mass function, based on the MILES (Sánchez-Blázquez et al., 2006) stellar
library. We simulate the expected morphological broadening by convolving the
high-resolution models with the grism response and the object morphology in
the $H_{\mathrm{F140W}}$ direct image (see Brammer et al., 2012b, for more
information on how the 2-D model spectra are generated). One-dimensional model
spectra are then extracted in the same way as was done for the observed G141
spectra. The red line in Figure 2 is the best-fit convolved model with an age
of 1.25 Gyr, fit to all wavelengths.

Figure 3.— The morphology of the residual [OIII] emission after subtracting a
1.25 Gyr model (histogram). The morphology is compared to the expectation for
a centrally-concentrated gas distribution (red line) and to a distribution
that follows the galaxy light (black line). The red shaded region indicates
the $1\sigma$ confidence interval resulting from the modeling analysis. The
resolved double-peak of the [OIII] feature is only expected for centrally-
concentrated emission with the low grism resolution.

The bottom panel of Figure 2 shows the residuals when subtracting the best-fit
model (red). Here, we see the [OIII] doublet ($\lambda 4959,5007\mathrm{\AA}$)
and H$\beta$ in emission. It is somewhat surprising to clearly resolve the
[OIII] doublet given the morphological broadening in the grism spectra. In
Figure 3, we demonstrate the expected resolution for the [OIII] emission
feature at a fixed line ratio of 3:1 when convolving with the interlaced PSF
(red) and the average galaxy $H_{\mathrm{F140W}}$ profile (black). The stacked
profile is constructed by re-centering, masking neighbors and normalizing by
the total flux of individual galaxies before taking the mean/median. The
[OIII] doublet is blended together when the emission is a global feature of
the galaxy, whereas centrally-concentrated emission (within roughly the
central pixel) shows two clear peaks. The point source nature of the [OIII]
emission lines also indicates that the high quality of the grism redshifts, as
redshift errors would broaden the observed [OIII] lines in the stack. We note
that although the relative ratio of the [OIII] doublet is mildly sensitive to
the assumed best-fit model, this effect is significantly less than the error
bars from the stacking analysis themselves. The residuals are consistent with
both centrally-concentrated [OIII] and H$\beta$ emission.

To model the central H$\beta$ and [OIII] enhancements, we simultaneously fit a
grid of Vazdekis et al. (2010) models convolved with the average galaxy
profile with enhancements ranging from 0–15% above the continuum flux for
H$\beta$ and [OIII] emission lines convolved with the PSF. We assume a fixed
line ratio of $\mathrm{[OIII]}\lambda 5007/\mathrm{H}\beta=1$, consistent with
the residuals in Figure 2 and line ratio diagnostics for post-starburst
galaxies in the Sloan Digital Sky Survey (SDSS; Mendel et al., 2013). The
best-fit model of 1.25 Gyr with a 5% [OIII] and H$\beta$ central enhancement
is shown in the top panel of Figure 2 (orange).

We measure the [OIII]$\lambda 5007$ line flux from the residuals in Figure 2
to be $1.4\pm 0.2\times 10^{-18}$ erg s-1 cm-2, which corresponds to
$L_{\mathrm{OIII}}$ of $1.7\pm 0.3\times 10^{40}$ erg s-1 for the median
redshift of 1.64. Similarly, we measure the H$\beta$ line flux to be $1.5\pm
0.2\times 10^{-18}$ erg s-1 cm-2. To determine the line flux errors, we add
the simulation errors in quadrature to the error in the average continuum
level. The resulting [OIII]$\lambda 5007$/H$\beta$ ratio is $0.9\pm 0.2$. If
we further correct $L_{\mathrm{OIII}}$ for the average dust extinction value
of $A_{V}=0.5$ determined from the photometry, we estimate the dust-corrected
[OIII] luminosity to be $L^{c}_{\mathrm{OIII}}=2.4\times 10^{40}$ erg s-1,
following Bassani et al. (1999). Adopting the conversion from
$L^{c}_{\mathrm{OIII}}$ to $L_{\mathrm{X}}$ from Lamastra et al. (2009) and a
factor of ten conversion from $L_{\mathrm{X}}$ to $L_{bol}$ (Lusso et al.,
2010), we estimate a typical active galactic nucleus (AGN) luminosity to be
$L_{bol}=2\times 10^{42}$ erg s-1. However, we note that it is not possible to
differentiate between centrally-concentrated residual star-formation or an AGN
without further information.

##  5\. What drives the spread in rest-frame colors?

Figure 4.— Median stacks of quiescent galaxies in two rest-frame color bins
(as defined in Figure 1). The top two panels show the number of galaxies and
median redshift of each wavelength bin. The observed stacks in the third panel
reveal that the bluest quiescent galaxies have younger ages, most notable by
the stronger H$\beta$ and weaker Mg I absorption. The fourth panel shows the
best-fit models, including a 5–6% centrally-concentrated enhancement of
H$\beta$ and [OIII] emission. The residuals are shown in the bottom panel.

With these high quality stacked spectra, we can now begin to assess the
driving factor behind the spread in rest-frame colors of quiescent galaxies.
The three possibilities include age, dust or metallicity variations amongst
galaxy populations. The difference in rest-frame $U$–$V$ and $V$–$J$ colors
corresponding to 1 Gyr of passive aging within the Bruzual & Charlot (2003)
models results in a 0.4 mag vector parallel to the quiescent sequence.
Similarly, increasing the metallicity from log(Z)=0.02 (Solar) to log(Z)=0.05
or adding 0.5 mag of dust reddening following the Calzetti et al. (2000) dust
law results in an almost identical vector in both direction and magnitude. Due
to degeneracies between age, dust and metallicity, it is notoriously difficult
to differentiate which of these effects dominates.

Metallicity variations will cause the color spread in galaxies to increase
over time as the galaxies passively evolve, reaching a roughly constant
scatter for the oldest stellar populations. As noted by Whitaker et al.
(2010), the trend for the intrinsic scatter in the rest-frame colors of
massive quiescent galaxies to steadily decrease from $z=2$ to the present
epoch (see also, e.g., Ruhland et al., 2009; Papovich et al., 2010) is
opposite to that predicted from metallicity variations alone. We therefore
expect that metallicity will not have a large effect at high redshifts.

As the shape of the Balmer/4000Å break depends sensitively on the treatment of
dust reddening, one would ideally need to correlate direct, independent
measurements of the ages of massive galaxies with their colors. Surprisingly,
Bezanson et al. (2013) did not find a simple correlation between the rest-
frame color and the strength of the Balmer absorption lines, as some of the
reddest galaxies in the spectroscopic sample had strong Balmer absorption
lines. However, with a small sample size of 13 galaxies, it is difficult to
gauge if these results are representative of the overall quiescent population.

In Figure 4, we plot the median stacked spectra in two rest-frame color bins
for the 171 quiescent galaxies, as color-coded in Figure 1. The total number
of galaxies in each bin is labeled in Figure 1 and the number of galaxies per
wavelength bin is indicated in the top histogram of Figure 4. A minimum number
of 20 galaxies was required per wavelength bin to ensure a robust measurement.

Here, we assess the effects of age as a function of rest-frame color for a
mass-limited sample of quiescent galaxies by measuring the strength of
absorption features such as Mg I and H$\beta$. Both of these spectral features
are sensitive to the age of the stellar population, but not very sensitive to
metallicity or dust reddening. Given the narrow baseline of the spectral
features, dust extinction is not expected to have a significant effect. At the
grism resolution, the H$\beta$ absorption feature is expected to be about 4–5%
stronger in a 1 Gyr stellar population as compared to a 2 Gyr stellar
population, whereas the Mg I absorption feature is about 3–4% weaker.
Increasing the metallicity results in a slightly deeper Mg I feature, but only
on the order of $\lesssim 1\%$. To first order, any change measured in
H$\beta$ and Mg I is due to age differences between the galaxies.

The average age of galaxies with the bluest rest-frame colors is
$0.9^{+0.2}_{-0.1}$ Gyr, consistent with the expectation that they are in the
“post-starburst” phase (Kriek et al., 2010; Whitaker et al., 2010, 2012).
These galaxies clearly show significantly stronger H$\beta$ and H$\gamma$
absorption than the reddest galaxies. Note the difference in shape of the
blended G-band and H$\gamma$ between these two stacks. The redder galaxies are
0.7 Gyr older on average, with a best-fit age of $1.6^{+0.5}_{-0.4}$ Gyr. For
a 0.7 Gyr age difference, the Bruzual & Charlot (2003) models predict
$\Delta(U-V)=0.2$ and $\Delta(V-J)=0.3$, corresponding to exactly the observed
difference in the median rest-frame colors of the two bins.

The older quiescent galaxies show clear signs of either central residual star-
formation activity or an AGN, while the enhancement in the younger galaxies is
not well-constrained despite a well-determined spectroscopic age. We measure
the residual [OIII]$\lambda 5007$ and H$\beta$ line fluxes for the younger and
older populations, finding ratios of $0.3\pm 0.1$ and $1.1\pm 0.2$
respectively. The bluer galaxies lack a clear [OIII] feature, although the
difference with the redder galaxies is only marginally significant.

##  6\. Discussion

We present stacked grism spectroscopy of a mass-limited sample of 171
quiescent galaxies at $1.4<z<2.2$ from the 3D-HST treasury program (Brammer et
al., 2012a). We demonstrate that we can resolve absorption features such as
the G-band ($\lambda$4304$\mathrm{\AA}$), H$\beta(\lambda$4861$\mathrm{\AA}$),
Mg I ($\lambda$5175$\mathrm{\AA}$), and Na I ($\lambda$5894$\mathrm{\AA}$) in
the median stacks due to the high quality grism redshifts. We further detect
centrally-concentrated H$\beta$ and [OIII] in emission, indicating either
residual central star-formation (e.g., Bezanson et al., 2013) or an AGN.

The average L[OIII] derived from the stacks of $0.4\times 10^{7}$ L⊙ is
consistent with the results of Mendel et al. (2013), who find that 70-80% of
recently-quenched galaxies in the SDSS have LINER-like active nuclei
(L${}_{\mathrm{[OIII]}}\lesssim 10^{7}$ L⊙). Furthermore, the [OIII]$\lambda
5007$/H$\beta$ ratio of 0.9 for a median mass of $10^{10.83}$ M⊙ places the
average quiescent galaxy within the LINER region of the Mass-Excitation
diagnostic diagram (Juneau et al., 2011), albeit marginally consistent with
the “composite” crossover region within the error bars. From the average dust-
corrected [OIII] ($\lambda$5007$\mathrm{\AA}$) luminosity, we infer that
quiescent galaxies may have a typical AGN luminosity of $L_{bol}=2\times
10^{42}$ erg s-1. These results are consistent with the recent work of Olsen
et al. (2013), who find that 70–100% of massive quiescent galaxies at
$1.5<z<2.5$ contain a low- or high-luminosity AGN. The potential ubiquitous
nature of AGNs in massive, quiescent $z\sim 2$ galaxies provides observational
support that black hole accretion may be more effective at the high-mass end
(e.g., Kriek et al., 2007).

To understand the relationship between rest-frame color and our spectroscopic
age measurements, we further stack the data in two rest-frame color bins. We
find that the bluest quiescent galaxies have a more prominent H$\beta$
absorption line, with best-fit ages of $0.9^{+0.2}_{-0.1}$ Gyr, confirming the
idea that these galaxies are the most recently quenched. Furthermore, we find
that redder quiescent galaxies (80% of the population) have an older age of
$1.6^{+0.5}_{-0.4}$ Gyr. Although previous spectroscopic studies have measured
similarly old ages (e.g., Kriek et al., 2009; Onodera et al., 2012; van de
Sande et al., 2012), this is the first time we are able to probe the full
population. These results suggest that age varies significantly with the
colors of quiescent galaxies, consistent with Whitaker et al. (2010).

Galaxies with the reddest rest-frame colors likely have some degree of dust
extinction, as dust-free models cannot produce such red colors at these
redshifts. Dust extinction does not greatly affect the line-index measurements
for single stellar populations given the narrow baseline of the spectral
features, but its effect can be significant for the 4000Å break (MacArthur,
2005). Consequently, although the absorption features of the 3D-HST stacks
will be effectively independent of dust extinction, the rest-frame color can
change significantly in the presence of dust due to its sensitivity to the
continuum shape. We cannot rule out that the reddest galaxies have a
potentially significant contribution from younger, dusty quiescent galaxies.
Nonetheless, despite this caveat, we measure a significant age difference
between the two populations of 0.7 Gyr that is consistent with the observed
rest-frame color differences.

The key result from this paper is that we have conclusively demonstrated that
old, quiescent galaxies exist in large numbers at $z\sim 2$. This is not a
surprise given extensive previous results that were based on broadband and
medium-band photometry (e.g., Brammer et al., 2011), and ground-based
spectroscopy for small samples (e.g., Onodera et al., 2012), but our detection
of metal lines in the stacked spectrum puts any lingering concerns to rest.
Our result also implies that the “UVJ technique”, first described in Labbé et
al. (2005), is very effective in identifying quiescent galaxies at $z\sim 2$.

We sincerely thank the referee for the thorough review of the manuscript. This
research was supported by an appointment to the NASA Postdoctoral Program at
the Goddard Space Flight Center, administered by Oak Ridge Associated
Universities through a contract with NASA. Support from STScI grant GO-12177
and NASA ADAP grant NNX11AB08G is gratefully acknowledged.

## References

  * Bassani et al. (1999) Bassani, L., Dadina, M., Maiolino, R., et al. 1999, ApJS, 121, 473
  * Bezanson et al. (2013) Bezanson, R., van Dokkum, P., van de Sande, J., et al. 2013, ApJ, 764, L8
  * Bower et al. (1992) Bower, R. G., Lucey, J. R., & Ellis, R. S. 1992, MNRAS, 254, 601
  * Brammer et al. (2012a) Brammer, G. B., van Dokkum, P. G., Franx, M., et al. 2012b, ApJS, 200, 13
  * Brammer et al. (2012b) Brammer, G. B., Sánchez-Janssen, R., Labbé, I., et al. 2012a, ApJ, 758, L17
  * Brammer et al. (2008) Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503
  * Brammer et al. (2011) Brammer, G. B., Whitaker, K. E., van Dokkum, P. G., et al. 2011, ApJ, 739, 24
  * Bruzual & Charlot (2003) Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000
  * Bundy et al. (2010) Bundy, K., Scarlata, C., Carollo, C. M., et al. 2010, ApJ, 719, 1969
  * Calzetti et al. (2000) Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682
  * Cardamone et al. (2010) Cardamone, C. N., Urry, C. M., Schawinski, K., et al. 2010, ApJ, 721, L38
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
  * Cimatti et al. (2008) Cimatti, A., Cassata, P., Pozzetti, L., et al. 2008, A&A, 482, 21
  * Cimatti et al. (2004) Cimatti, A., Daddi, E., Renzini, A., et al. 2004, Nature, 430, 184
  * Daddi et al. (2005) Daddi, E., Renzini, A., Pirzkal, N., et al. 2005, ApJ, 626, 680
  * Dobos et al. (2012) Dobos, L., Csabai, I., Yip, C.-W., et al. 2012, MNRAS, 420, 1217
  * Grogin et al. (2011) Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35
  * Juneau et al. (2011) Juneau, S., Dickinson, M., Alexander, D. M., et al. 2011, ApJ, 736, 104
  * Koekemoer et al. (2011) Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36
  * Kriek et al. (2010) Kriek, M., Labbé, I., Conroy, C., et al. 2010, ApJ, 722, L64
  * Kriek et al. (2008) Kriek, M., van der Wel, A., van Dokkum, P. G., et al. 2008, ApJ, 682, 896
  * Kriek et al. (2007) Kriek, M., van Dokkum, P. G., Franx, M., et al. 2007, ApJ, 669, 776
  * Kriek et al. (2006) —. 2006, ApJ, 649, L71
  * Kriek et al. (2009) Kriek, M., van Dokkum, P. G., Labbé, I., et al. 2009, ApJ, 700, 221
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231
  * Labbé et al. (2005) Labbé, I., Huang, J., Franx, M., et al. 2005, ApJ, 624, L81
  * Lamastra et al. (2009) Lamastra, A., Bianchi, S., Matt, G., et al. 2009, A&A, 504, 73
  * Lusso et al. (2010) Lusso, E., Comastri, A., Vignali, C., et al. 2010, A&A, 512, A34
  * MacArthur (2005) MacArthur, L. A. 2005, ApJ, 623, 795
  * Mendel et al. (2013) Mendel, J. T., Simard, L., Ellison, S. L., et al. 2013, MNRAS, 429, 2212
  * Nicol et al. (2011) Nicol, M.-H., Meisenheimer, K., Wolf, C., et al. 2011, ApJ, 727, 51
  * Olsen et al. (2013) Olsen, K. P., Rasmussen, J., Toft, S., et al. 2013, ApJ, 764, 4
  * Onodera et al. (2012) Onodera, M., Renzini, A., Carollo, M., et al. 2012, ApJ, 755, 26
  * Papovich et al. (2010) Papovich, C., Momcheva, I., Willmer, C. N. A., et al. 2010, ArXiv e-prints
  * Patel et al. (2012) Patel, S. G., Holden, B. P., Kelson, D. D., et al. 2012, ApJ, 748, L27
  * Ruhland et al. (2009) Ruhland, C., Bell, E. F., Häußler, B., et al. 2009, ApJ, 695, 1058
  * Sánchez-Blázquez et al. (2006) Sánchez-Blázquez, P., Peletier, R. F., Jiménez-Vicente, J., et al. 2006, MNRAS, 371, 703
  * Toft et al. (2012) Toft, S., Gallazzi, A., Zirm, A., et al. 2012, ApJ, 754, 3
  * van de Sande et al. (2012) van de Sande, J., Kriek, M., Franx, M., et al. 2012, ArXiv e-prints
  * van de Sande et al. (2011) —. 2011, ApJ, 736, L9+
  * van Dokkum & Brammer (2010) van Dokkum, P. G., & Brammer, G. 2010, ApJ, 718, L73
  * Vazdekis et al. (2010) Vazdekis, A., Sánchez-Blázquez, P., Falcón-Barroso, J., et al. 2010, MNRAS, 404, 1639
  * Whitaker et al. (2012) Whitaker, K. E., Kriek, M., van Dokkum, P. G., et al. 2012, ApJ, 745, 179
  * Whitaker et al. (2011) Whitaker, K. E., Labbé, I., van Dokkum, P. G., et al. 2011, ApJ, 735, 86
  * Whitaker et al. (2010) Whitaker, K. E., van Dokkum, P. G., Brammer, G., et al. 2010, ApJ, 719, 1715
  * Williams et al. (2009) Williams, R. J., Quadri, R. F., Franx, M., et al. 2009, ApJ, 691, 1879
  * Wuyts et al. (2007) Wuyts, S., Labbé, I., Franx, M., et al. 2007, ApJ, 655, 51


11institutetext: School of Computing, Informatics, and Decision Science
Engineering, Arizona State University, Tempe, USA
11email<EMAIL_ADDRESS>
22institutetext: Information Sciences Institute, University of Southern
California, Marina del Rey, USA
22email<EMAIL_ADDRESS>
# Discovering patterns of online popularity from time series

Mert Ozer 11 1234-5678-9012 Anna Sapienza 22 Andrés Abeliuk 22 Goran Muric 22
Emilio Ferrara 22

###### Abstract

How is popularity gained online? Is being successful strictly related to
rapidly becoming viral in an online platform or is it possible to acquire
popularity in a steady and disciplined fashion? What are other temporal
characteristics that can unveil popularity of online content? To answer these
questions, we leverage a multi-faceted temporal analysis of the evolution of
popular online contents. Here, we present dipm-SC: a multi-dimensional shape-
based time-series clustering algorithm with a heuristic to find the optimal
number of clusters. First, we validate the accuracy of our algorithm on
synthetic datasets generated from benchmark time series models. Second, we
show that dipm-SC can uncover meaningful clusters of popularity behaviors in a
real-world Twitter dataset. By clustering the multidimensional time-series of
popularity of contents coupled with other domain-specific dimensions, we
uncover two main patterns of popularity: bursty and steady temporal behaviors.
Moreover, we find that the way popularity is gained over time has no
significant impact on the final cumulative popularity.

###### Keywords:

multidimensional time series, shape-based clustering, online popularity,
social media

##  1 Introduction

How do contents get popular on the Web? Is being successful, strictly
constrained to a rapid and bursty virality [1] or is it possible to acquire
popularity in a steady and disciplined fashion? What are other temporal
characteristics that can unveil popularity of online contents?

Temporal dynamics of popularity have been studied extensively in Twitter [41],
news [24], and videos [9]. Thus far, previous work studying popularity either
explored topological patterns in the networks that lead to cascades [4, 25],
or focused on the short time period in which significant popularity is gained
[43, 23]. However, popularity is subject to endogenous and exogenous factors
[23] that affects its temporal dynamics. For instance, popularity of an item
may depend on the amount of external promotions it receives [35]. Analyzing
different temporal patterns that arise when considering both endogenous and
exogenous factors will further extend our understanding on temporal collective
behavior.

To take into account all these factors at a time and tackle the problem of
uncovering their relation with the dynamics of popularity, we focus on
multidimensional clustering methods. Multidimensional clustering can indeed
allow us to uncover the hidden relations between time series of events having
different evolutions in more than one dimension. For instance, the dynamics of
a user content can be not only described by its popularity but also by the
online activity and pace of the user, number of social connections in social
networks.

A wide range of studies have been focused on solving the problem of
efficiently clustering time series [2]. Some approaches have adopted methods
such as sub-sequence techniques [38, 46, 47], and dimensionality reduction
models [45, 7]. However, we aim at uncovering different patterns of popularity
on the basis of both their shapes and the time at which the main popularity
gain is achieved. Therefore, sub-sequence or subspace solutions are not
suitable for this analysis as they lose temporal granularity and order. In
this paper, we take inspiration from unidimensional clustering approaches [43,
30, 33, 12] that use a shape-based distance and a time series similarity
measure that is scale and shift invariant. We generalize these methods to
multidimensional time series and parametrized the temporal shifting, thus
allowing us to capture temporal patterns occurring at different times.

In the present work, we design a multi-faceted temporal analysis to study the
timeline of popular online content and provide insights on how popularity is
gained over time. To this aim, we extend K-Spectral Centroid (kSC) [43], a
unidimensional shape-based time series clustering algorithm into our
multidimensional version, called dipm-SC. It partitions multidimensional time
series data into clusters and uses a heuristic to find the optimal number of
clusters. Each cluster is characterized by a centroid, having a unique shape
per dimension identifying the average temporal behavior of the cluster
members.

We validate the accuracy of our algorithm on the datasets generated from
benchmark time series models with varying parameters of time series length,
dimensions and number of underlying clusters. To the best of our knowledge,
this is the first study which explores multidimensional use of k-SC.

Then, with the help of the dipm-SC, we study temporal multidimensional
behavior in two real world large datasets containing data from GitHub and
Twitter. In Github, we study the popularity of repositories considering a
series of other temporal dimensions that keep track of the repository
development, such as pushes and pull requests a repository receives, and
repository issues, such as issue entries and comments under them. Our
algorithm identifies 9 main clusters among most popular repositories on GitHub
in different shapes. In Twitter, we study how popular hashtags unfold in a
window of time where a major event happens. We keep track of number of times a
hashtag used hourly (as a measure of popularity) alongside with hourly
positive and negative sentiment rates and tweets containing the hashtag posted
by bot accounts. Our algorithm identifies 6 main clusters among most popular
hashtags on Twitter. We show that both bursty and steady clusters are evident
in both datasets. On Twitter, when a bursty popularity gain is present, we
observe increased sentiment and bot activity rates around the popularity
spike. We also note that sentiment and bot activity decay after the spike can
be characterized by the fat-tail of the popularity spike. Moreover, we find
that how repositories unfold over time shows no significant impact on the
cumulative popularity at the end of their lifetime. However, in the Twitter
data, our algorithm is able identify different cluster shapes leading to
different cumulative popularity gains.

##### Our contributions

  1. 1.

We propose a novel algorithm to cluster multidimensional shape-based time
series. Our approach detect the interplay of different dimensions that allow
to perform a non-supervised classification of different temporal patterns.

  2. 2.

We validate the accuracy of our algorithm on the datasets generated from
benchmark time series models, and compare its performance with other scalable
multidimensional time series methods.

  3. 3.

We present two real-world case studies, where we showcase the effectiveness of
our algorithm in both finding meaningful popularity dynamics and providing
clusters that are coherent and easy to interpret.

##  2 Methodology

In this section, we formally define the problem and discuss the need for
extensions to time series clustering algorithms in the literature. Then, we
develop our multidimensional extension based on a several unidimensional
shape-based time series clustering methods.

###  2.1 Problem Definition

In this work, we aim at finding temporal patterns of popular online content in
a multi-faceted manner. Let $N$ be the number of entities to be studied
described by $D$ temporal features. Each temporal feature is represented by a
discrete time series in $M$ time units (e.g., days or hours). This
representation of our data can be encoded in a tensor $\mathcal{X}$ of size
$N\times D\times M$.

Given $\mathcal{X}$, we want to partition the $N$ multivariate time series
into an optimal number of $K$ clusters, by taking the temporal dynamics of
different $D$ dimensions into account at the same time. To this aim, we
perform a clustering process based on the cumulative similarity of shapes
among corresponding dimensions of the time series. Therefore, each cluster $k$
is represented by a multi-dimensional shape
$\mathbf{c}_{k}\in\mathbb{R}^{D\times M}$, which describes the overall
behavior of the cluster across multiple dimensions. In this work, solving the
aforementioned problem gives us the ability to uncover similar patterns of
popularity gain while differing in other dimensions, and vice versa.

###  2.2 Proposed Algorithm

The kSC algorithm proposed by Yang and Leskovec [43] stands out as a prominent
solution for analyzing time series data constructed from online human
behaviors. Along that line, there have been significant contributions to the
shape-based time series clustering literature [30, 31]. Here, we extend and
compare the following shape-based univariate time series clustering
algorithms: kShape[30], kDBA[33], kAVG+ED[12], and kSC[43]. We only report the
details about the methodology used to extend kSC[43] as it is the only
algorithm specifically developed for online social platforms and, as shown in
the supplementary material 111www.public.asu.edu/~mozer/dipm-SC, is the one
performing better in our synthetic scenario. However, we applied similar
extensions to the kShape, kDBA, and kAVG+ED algorithms.

Our extension focuses on two parts. First, we modify the existing kSC
algorithm to be applicable to multidimensional time series data. This process
consists of the two steps of updating the distance and averaging functions of
kSC. Second, we develop an iterative framework to find the optimal number of
multidimensional shape clusters inspired by a similar work developed for
kMeans [20].

####  2.2.1 kSC framework

K-Spectral Centroid (kSC) is an iterative kMeans based time series clustering
algorithm. The algorithm iteratively assigns time series data to clusters and
update cluster centroids until no updates to clusters is made. It adopts a
scale and shift invariant distance function [8] to measure the proximity of
two time series. Therefore, each cluster’s centroid represents the average
shape of its members. At each iteration, the method uses a spectral clustering
technique as its averaging function to update the centroids. However, as it is
originally proposed and designed for unidimensional time series data, we need
to extend its distance and averaging functions in order to apply it to a
multidimensional scenario.

####  2.2.2 Extending the Distance Function

There are several ways to extend the univariate distance function to its
multivariate version. First, we can compute distances separately on each
dimension and sum them. However, this option is not optimal, as usual time
series distances involve the manipulation of the signal through shifting and
warping process. Thus, the use of separate distances would lead to different
manipulations (shifting and warping) in each dimension. As our task requires
to observe their interplay instead, manipulating each dimension differently is
not desirable. Second, we can compute distances by using an optimal shared
manipulation across all dimensions. Here, we propose extending the kSC
algorithm’s distance function [8] as follows

$dist(\mathbf{c}_{k},\mathbf{x})=\underset{\alpha_{d},q}{\min}\frac{1}{D}\sum_{d=1}^{D}\dfrac{\left\|\mathbf{c}_{k}(d,:)-\alpha_{d}\mathbf{x}_{q}(d,:)\right\|}{\left\|\mathbf{c}_{k}(d,:)\right\|}\;\;,$
(1)

where $\mathbf{x}(d,:)$ is the $d$’th dimension of time series
$\mathbf{x}\in\mathbb{R}^{D\times M}$, $\mathbf{x}_{q}$ is the shifting of
$\mathbf{x}$ by $q$ steps in time, and $\mathbf{c}_{k}$ is the centroid of
cluster $k$. The above equation constitutes the essential part of calculating
the distance of multivariate time series $\mathbf{x}$ to a centroid. It indeed
enforces the same shifting parameter $q$ regardless of dimensions. Yet, it
keeps the scaling parameter $\alpha_{d}$ unique for each dimension. Here, we
assume that even if each dimension may experience trends in different scales,
they should have the same temporal order.

Notice that, for any fix $q$, the optimal $\alpha_{d}(q)$ that solves Eq. (1)
is

$\alpha_{d}(q)=\frac{\mathbf{x}_{q}(d,:)^{T}\mathbf{c}_{k}(d,:)}{\left\|\mathbf{x}_{q}(d,:)\right\|^{2}}\;\;,$
(2)

by setting the gradient to zero while q is fixed. Thus, minimizing the
distance function is reduced to a brute-force searching of the optimal $q$,
where $\alpha_{d}$ is calculated for each shifting option as in Eq. (2). The
algorithm can become really slow when full-length shifting is allowed. As we
are interested in temporal behavior of online content, we focus on moments in
time when the events happen. Therefore, we aim to keep the shifting parameter
to the minimum.

####  2.2.3 Extending the Averaging Function

As a second step in our method, we need to extend the averaging function of
kSC. The role of the averaging function is to find a suitable centroid for
each cluster, where the distance of its members is minimum to the centroid

$\displaystyle\mathbf{c}_{k}=\underset{\mathbf{c}}{\operatorname*{argmin}}\sum_{\mathbf{x}\in
S_{k}}dist(\mathbf{c},\mathbf{x})^{2}\;\;,$

where $S_{k}$ is the set of time series belonging to cluster $k$.

A possible way of extending the averaging function consists of finding a
univariate centroid which minimizes the distance to every dimension. We avoid
this option since there may be different behaviors on different dimensions and
our task involves observing them in parallel. Instead, we extract a unique
shape for each dimension separately. First, we replace the distance function
with its multivariate version

$\displaystyle\mathbf{c}_{k}=\underset{\mathbf{c}}{\operatorname*{argmin}}\sum_{\mathbf{x}\in
S_{k}}\underset{\alpha_{d},q}{\min}\frac{1}{D^{2}}\sum_{d}^{D}\dfrac{\left\|\mathbf{c}(d,:)-\alpha_{d}\mathbf{x}_{q}(d,:)\right\|^{2}}{\left\|\mathbf{c}(d,:)\right\|^{2}}\;\;.$
(3)

To be able to use the spectral properties of Rayleigh quotient in accordance
with the original kSC paper[43], we first need to find the optimal
$\alpha_{d}$ and $q$. $\alpha_{d}$ has a direct solution as proposed in Eq.
(2). However $q$ is dimension invariant. Here we set $q$ as the optimal $q$
found during the distance’s computation, and jointly shift every dimension of
the time series $\mathbf{x}$ accordingly during the new centroid calculation.
Finally, we invert the order of the two sums in Eq. (3) and get

$\displaystyle\mathbf{c}_{k}=\underset{\mathbf{c}}{\operatorname*{argmin}}\frac{1}{D^{2}}\sum_{d}^{D}\sum_{\mathbf{x}\in
S_{k}}\dfrac{\left\|\mathbf{c}(d,:)-\alpha_{d}\mathbf{x}(d,:)\right\|^{2}}{\left\|\mathbf{c}(d,:)\right\|^{2}}\;\;.$
(4)

Our minimization problem does not have dimension independent variables anymore
and we can solve each dimension’s centroid separately by following similar
steps to the univariate case

$\displaystyle\mathbf{c}_{k}(d,:)$
$\displaystyle=\underset{\mathbf{c}^{\prime}}{\operatorname*{argmin}}\sum_{\mathbf{x}\in
S_{k}}\dfrac{\left\|\alpha_{d}\mathbf{x}(d,:)-\mathbf{c}^{\prime}\right\|^{2}}{\left\|\mathbf{c}^{\prime}\right\|^{2}}$
$\displaystyle=\underset{\mathbf{c}^{\prime}}{\operatorname*{argmin}}\frac{1}{\left\|\mathbf{c}^{\prime}\right\|^{2}}\sum_{\mathbf{x}\in
S_{k}}\left\|\dfrac{\mathbf{x}(d,:)^{T}\mathbf{c}^{\prime}\mathbf{x}(d,:)}{\left\|\mathbf{x}(d,:)\right\|^{2}}-\mathbf{c}^{\prime}\right\|^{2}$
$\displaystyle=\underset{\mathbf{c}^{\prime}}{\operatorname*{argmin}}\frac{1}{\left\|\mathbf{c}^{\prime}\right\|^{2}}\sum_{\mathbf{x}\in
S_{k}}\left\|\left(\dfrac{\mathbf{x}(d,:)\mathbf{x}(d,:)^{T}}{\left\|\mathbf{x}(d,:)\right\|^{2}}-I\right)\mathbf{c}^{\prime}\right\|^{2}$
$\displaystyle=\underset{\mathbf{c}^{\prime}}{\operatorname*{argmin}}\frac{1}{\left\|\mathbf{c}^{\prime}\right\|^{2}}\mathbf{c}^{\prime
T}\sum_{\mathbf{x}\in
S_{k}}\left(I-\dfrac{\mathbf{x}(d,:)\mathbf{x}(d,:)^{T}}{\left\|\mathbf{x}(d,:)\right\|^{2}}\right)\mathbf{c}^{\prime}\;\;,$

and by replacing $\sum_{\mathbf{x}\in
S_{k}}(I-\dfrac{\mathbf{x}(d,:)\mathbf{x}(d,:)^{T}}{\left\|\mathbf{x}(d,:)\right\|^{2}})$
with $M$, we obtain the following minimization problem

$\displaystyle\mathbf{c}_{k}=\underset{\mathbf{c^{\prime}}}{\operatorname*{argmin}}\frac{\mathbf{c}^{\prime
T}M\mathbf{c}^{\prime}}{\left\|\mathbf{c}^{\prime}\right\|^{2}}\;\;,$

which achieves its minimum value with the eigenvector of $M$ corresponding to
the smallest eigenvalue. We refer further interested readers to the properties
of the Rayleigh quotient with parameters $\mathbf{c}^{\prime}$ and $M$ [16].

After successfully constructing the two building blocks of kSC for our
multivariate case, we present the pseudo code of our iterative algorithm
multivariate kSC (m-kSC), cf. Alg. 1.

Algorithm 1 m-kSC Algorithm

1:Input:$\\{\mathcal{X},K\\}$ where $\mathcal{X}\in\mathbb{R}^{N\times D\times
M}$ is the tensor containing $N$ multivariate time series and $K$ is number of
clusters.

2:Output: $\\{\mathcal{C},S\\}$ where $\mathcal{C}\in\mathbb{R}^{K\times
D\times M}$ is the tensor of cluster centroids and $S$ contains each cluster
assignments.

3:Initialize cluster assignments $S$ randomly

4:while $S$ changes on every iteration do

5: for $k=1:K$ do

6: for $d=1:D$ do

7: $M=\sum_{\mathbf{x}_{n}\in
S_{k}}(I-\dfrac{\mathbf{x}_{n}(d,:)\mathbf{x}_{n}(d,:)^{T}}{\left\|\mathbf{x}_{n}(d,:)\right\|^{2}})$

8: $\mathbf{C}(k,d,:)=$ Smallest eigenvector of $M$.

9: end for

10: end for

11: for $n=1:N$ do

12:
$k=\underset{k=1,...,K}{\operatorname*{argmin}}dist(\mathbf{c}_{k},\mathbf{x}_{n})$
using Eq. 1

13: $S(n)=k$

14: end for

15:end while

####  2.2.4 Finding the optimal K

As with many kMeans based approaches, m-kSC also needs the number of clusters
as a parameter. A common practice to find the optimal $K$ involves post-
processing, where we run the clustering algorithm several times for different
$k$ values. In this case, a trade off between the final loss function value
and model complexity decides the optimal $K$. However, to avoid re-running the
algorithm for every possible $k$, here we adopt dip-dist, the iterative
strategy proposed by Kalogeratos and Likas [20]. We choose the dip-dist
technique over others [32, 18] as the only assumption dip-dist makes is the
unimodality of pairwise distances of cluster members. It suggests that an
optimal cluster structure should not have more than a single mode among
pairwise distances of its members (multimodality). The density of pairwise
distance values should reach to maximum around the mode and decay while moving
away.

Multimodality of each cluster is checked by its members’ distance to each
other. First, we calculate the distance of each time series to the other
members in the cluster. Then, we sort the distance vector in decreasing order.
If the null hypothesis of unimodality in the distribution of the sorted
distances is rejected by Hartigan’s dip test [19] with a significance level
$p$-value$<\alpha$, it is considered to be a splitter in the cluster. Finally,
if the number of splitters is higher than a given threshold $v$, we split the
cluster into two disjoint sub-clusters. Our splitting strategy involves a
local search with m-kSC where k is set to 2. Heuristically, we choose the
least squared error objective function centroids over 10 runs. We iteratively
split clusters until there is no clusters left with ratio of splitters greater
than the given threshold. As suggested by Kalogeratos and Likas [20], at each
iteration, we only split the cluster with maximum number of splitters. We
report our complete dip-test based multidimensional spectral centroid
algorithm dipm-SC, cf. Alg. 2. We refer further interested readers to
Kalogeratos and Likas [20] for kMeans extension and to Hartigan and Hartigan
[19] for the Hartigan’s dip test. We make the source code and the experimental
settings available at www.public.asu.edu/~mozer/dipm-SC/dipm_source_code.zip.

The computational cost of m-kSC is dominated by the calculation of matrix M
for each $x_{n}$ and for d dimensions ($\mathcal{O}(m^{2}dn)$) and the
eigendecomposition of M for k clusters ($\mathcal{O}(m^{3}dk)$). Multimodality
test in dipm-SC for each cluster k has a complexity of $\mathcal{O}(k(bn\log
n+n^{2}))$. So, the total complexity of the dipm-SC becomes the maximum of the
$\\{\mathcal{O}(m^{2}dn),\mathcal{O}(m^{3}dk),\mathcal{O}(k(bn\log
n+n^{2}))\\}$.

Algorithm 2 dipm-SC Algorithm

1:Input:$\\{\mathcal{X},\alpha,s\\}$ where $\alpha$ is the significance level,
and $v$ is the split threshold.

2:Output: $\\{\mathcal{C}^{*},S^{*}\\}$.

3:Assign time series to one cluster $S=ones(N,1)$

4:for $i=1:N$ do

5: for $j=1:N$ do

6: $\mathbf{d}(i,j)=dist(\mathcal{X}(i,:,:),\mathcal{X}(j,:,:))$ Equation 1

7: end for

8:end for

9:while $max(score_{k})>v$ do

10: for $k=1:K$ do

11: $score_{k}=checkClusterModality(\mathbf{d}(S_{k},S_{k}),S_{k},\alpha)$

12: end for

13: if $max(score_{k})>v$ then

14: $k=\underset{k}{\operatorname*{argmax}}(score_{k})$

15: $[\mathbf{c}_{1},\mathbf{c}_{2}]=$splitCluster$(\mathcal{X}(S_{k},:,:))$

16: $\mathcal{C}(k,:,:)=\mathbf{c}_{1}$

17: $\mathcal{C}(K+1,:,:)=\mathbf{c}_{2}$

18: $K=K+1$

19: $[\mathcal{C},S]=$m-kSC$(\mathcal{X},K,\mathcal{C})$

20: end if

21:end while

##  3 Evaluation

We evaluate the results of our algorithm applied to one synthetically
generated and two real-world multivariate time series datasets: GitHub and
Twitter. We refer interested readers to
https://www.public.asu.edu/~mozer/dipm-SC/supp_material.pdf for the comparison
of the dipm-SC with extension of other shape-based time series clustering
algorithms on synthetically generated time series dataset and the GitHub
analysis. For the sake of brevity of this paper, here we present our
experimental results for Twitter hashtag dataset. First, we perform the
qualitative analysis by studying different behavioral patterns of the detected
multivariate centroids. Secondly, we characterize the clusters based on their
”periodicity”, locating them on a scale from periodic to viral. We expect that
periodic time series will behave steadily and that the viral time series will
exhibit significant spikes over short time periods. To this aim, we use two
previously introduced metrics; burstiness and memory [15]. These two metrics
are computed from the distribution of the inter-event time $P(\tau)$ between
consecutive events of time series data. In particular, given the inter-event
time distribution $P(\tau)$ having mean $m_{\tau}$ and standard deviation
$\sigma_{\tau}$, we calculate burstiness as

$B=\frac{\sigma_{\tau}-m_{\tau}}{\sigma_{\tau}+m_{\tau}}\;\,$ (5)

where $B\in[-1,1]$: $B=1$ corresponds to the most bursty sequence of events,
and $B=-1$ is a completely regular (periodic) sequence.

This measure however is not enough to characterize the correlations taking
place between consecutive events. Correlation of consecutive events can be
thought as a memory process. A measure that can be used to compute such memory
coefficient between consecutive events $(\tau_{i},\tau_{i+lag})$ is defined as

$M=\frac{1}{n_{\tau}-1}\sum^{n_{\tau}-1}_{i=1}\frac{\left(\tau_{i}-m_{1}\right)\left(\tau_{i+lag}-m_{2}\right)}{\sigma_{1}\sigma_{2}}\;\,$
(6)

where, $n_{\tau}$ is the number of intervals between events in the time
series, while $m_{1,2}$ and $\sigma_{1,2}$ are respectively the sample mean
and standard deviation of $\tau_{i,i+lag}$’s. We set the value of lag based on
the empirically observed value of periodicity for each case.

##  4 Datasets description

We collected tweets starting from February 14th until March 6th, 2018 related
to the Parkland school shooting event from GNIP. To this aim, we provided 140
terms (hashtags, words, phrases) related to the event and the broader gun
debate in the United States.222https://www.public.asu.edu/~mozer/dipm-
SC/GNIP_query_list.txt GNIP provides any tweet that contains at least one of
the queried terms, while not suffering from any bias implications of free API
streams [29]. Our data contains $~{}23.7$M tweets coming from $~{}3.7$M users.

We then identify the top $1,000$ most occurring (popular) hashtags in our
dataset and build hourly time series for each of them. Alongside with hourly
counts of each hashtag, we construct the hourly ratio of tweets sent by bot
accounts to the total hourly volume. To identify the most obvious bot
accounts, we used the free API provided by Botometer [10]. We tested about
$5\%$ of the most active accounts ($193$K approximately) and found that
$~{}19$K of them are likely bots—in line with other research on bot frequency
[39].

As the third and fourth dimensions, we consider the hourly rate of positive
and negative sentiment tweets. We use an off-the-shelf short text
classification tool SentiStrength [37] to detect positive and negative
sentiment of each tweet where hashtags occur. SentiStrength provides scores
for positive and negative sentiments separately. We aggregate each sentiment’s
score on an hourly basis by simply summing them. Finally, we build the hourly
sentiment ratio time series for each hashtag by normalizing the sentiment
volumes by total tweet volumes.

Therefore, the timeline of each hashtag over time is represented by 4 separate
dimensions: total count, bot tweet rate, positive sentiment rate, and negative
sentiment rate.

##  5 Popularity Gain Patterns on Twitter

In this section we present Twitter hashtag timeline analysis with the help of
multidimensional shape clusters identified by our dipm-SC algorithm. We use 24
hours as our shifting parameter. We also apply gaussian smoothing of window
size 24 to eliminate noise. In total we are able to identify 13 clusters in
different timeline shapes. Here we present 6 clusters which spans nearly 78%
percent of popular hashtags under study. We leave other 7 smaller clusters
reachable externally for brevity of the paper
(http://www.public.asu.edu/~mozer/dipm-SC). We report these 6 most popular
timeline shapes in Fig. 1.

(a) Cluster-P (b) Cluster-B1 (c) Cluster-B2

(d) Cluster-B3 (e) Cluster-Boycott1 (f) Cluster-Boycott2

Figure 1: Shape of the uncovered clusters of Twitter

First, we observe that cluster-P in Figure 1a is significantly different in
shape when compared with others. In particular, it does not show any burst of
activity in each of the four dimensions considered (popularity, sentiment
ratios, and bot tweet ratio). On the contrary, this cluster shows a steady
pattern with daily periodic activities.

On the other hand, the remaining 5 clusters are characterized by a sudden
spike in the popularity signal at a certain point in time and different trends
among the other dimensions. The popularity burst of Cluster-B1 and B2 occurs
at the same time. Spike in popularity gain considerably differs in the tail
part. Cluster-B1 shows a longer tail while Cluster-B2 dies off sooner.
Cluster-B2 temporal pattern shows not only a decreasing trend in popularity
but also a faster decay in both sentiment and bot ratios then the one we
observe in cluster-B1. We notice an analogous discrepancy among the last two
clusters (i.e., Boycott1 and Boycott2).

Finally, we analyze the temporal coherence of hashtags clustered together. We
note that our algorithm does not consider semantics while clustering timeline
of hashtags. The detected clusters are indeed identified because they share
similar temporal patterns among different dimensions. We report the top 8
closest hashtags to the cluster centroids in the top part of each panel in
Figure 1.

Cluster-P contains broader issue related hashtags such as #NRA, #2A(second
amendment), #MAGA, #Trump, which show steady and periodic behavior among
dimensions. Cluster B1 and B2 exhibits spikes around the time when the tragic
shooting takes place. In cluster-B1 we observe event-related long lasting
hashtags, while in cluster-B2 we observe event-related expiring hashtags such
as #prayfor variants, #MentalHealth. Cluster B3 is unique in terms of its
spike time when student survivors start the #neveragain movement and appear in
a town hall organized by CNN. Lastly, we notice two boycott related hashtags
dominated clusters spiking in popularity gain around the same time. It
overlaps with the time of students’ call to boycott NRA and companies that
shares financial interests with it. It is easy to see these two clusters
differing in their tail parts. Cluster-Boycott2 (longer popularity gain tail)
shows more NRA & Amazon related hashtags while Cluster-Boycott1 involves other
companies names. We will later present statistically significant difference in
cumulative popularity gains between these two clusters.

###  5.1 Interplay of Dimensions

The key piece of information we acquire through this analysis is how
popularity spike of a hashtag correlates with increased sentiment ratios and
bot involvement. For every bursty popularity gain, we observe an increase in
all other dimensions. Previous studies on bursty attention mechanisms of
online content usually focus on the characteristic of the fat-tail [44, 23,
34]. Here, we observe that both sentiment and bot tweet ratios stay more
steady with a longer tail of popularity (e.g., B1, B3, and Boycott2) than
dying ones (e.g., B2, Boycott1). We evidence it by fitting a linear curve and
measuring the average slope of cluster members. When compared, hashtags
belonging to B1 have a higher slope than the ones in B2 in sentiment and bot
involvement dimensions after the spike. We observe the same pattern between
Boycott2 and Boycott1 clusters, where Boycott2 has a higher slope in all three
dimensions after the spike.

Another interesting behavior we direct attention to is the sudden increase in
bot tweet ratio after the popularity spike in cluster B2, B3 and Boycott1.
These clusters are characterized by a weaker tail in popularity gain after a
burst in popularity. This signals a steady bot involvement for a while
although abandonment of the overall activity takes place.

(a)

(b)

Figure 2: Burstiness and Memory of Clusters. Figure (a): Burstiness and memory
of Twitter clusters. Figure (b): Cumulative popularity of cluster members at
the end of the data timeline. Each blue dot represents a hashtag.

###  5.2 Burstiness and Memory of Clusters

First, we quantify popularity gain’s temporal behavior with the described
burstiness and memory metrics. We present average burstiness and memory of
each cluster in Figure 2a. As expected, cluster-P has the lowest burstiness
while enjoying higher values of memory. Lowest memory and highest burstiness
belongs to Boycott1 cluster.

Next, we investigate if different multi-dimensional timeline shapes lead to
significant difference in hashtags’ cumulative popularity gain at the end of
the time window we observe. Since this analysis is based on the timeline of a
hashtag rather than a lifespan, we do not align hashtags in their creation
date. That is why, for example, comparing Cluster-B1 against Cluster-B3 is not
a meaningful test. However, it is noteworthy that two-sample Kolmogorov-
Smirtnov test rejects the null hypothesis for Cluster-Boycott1 and Cluster-
Boycott2 of which we know the beginning of campaigns. Differing distributions
for these two shape clusters can be also observed from Figure 2b. Moreover, it
is also appealing to see no cumulative popularity distribution difference
between hashtags of steady cluster P and bursty cluster B1 (p-value 0.056).

##  6 Related Work

Our work overlaps with the current literature among two branches. First is the
task of clustering time series, and the second is the task of uncovering
temporal patterns of online content.

Clustering time series has been a challenging task and research area that
produced vast amount of literature [2]. Usually, when external feature
learning or modeling [42, 26, 3] is not imposed, clustering of time series
involves choosing a proper time series distance function and an algorithm.
Numerous works have introduced various distance measures to calculate
proximity of time series data [40, 11]. There have been partitional [30, 43,
33, 12], subsequence matching [6, 22], hierarchical [21] approaches using
these different time series distance measures. Moreover, a major body of work
[38, 47] exists in subsequence matching based time series clustering where
they identify shorter most identifying portions of time series data also known
as shapelets to group them. For the multivariate time series data, same
categorizations can be made as modeling based [17, 14], and variants of
generalizing univariate solutions to multivariate cases [36, 13]. Our approach
falls into the second category where we extend existing distance functions
available for univariate time series data and update centroid finding (i.e.
averaging function) accordingly.

As the second line of work from literature, we compare our effort with
analysis of temporal patterns of online content. There is significant amount
of literature focusing on characterizing and modeling bursts and long tail
dynamics [5, 27, 23]. First and foremost, our study does not focus on
modeling, rather it finds clusters of shapes inherent to online temporal
behavior. However, we show that our cluster shapes can be identified
distinctively by their burstiness and periodicity through postprocessing. Our
work also differentiates itself from others [28, 27, 43] by analyzing lifetime
or a longer timeline of online entities rather than focusing on a window where
activity peaks occur and omitting the rest.

##  7 Conclusion

In the present work, we set to study the temporal patterns that lead a user’s
content to gain popularity in online platforms. In particular, our study
tackled the following problems: (i) uncovering the different temporal patterns
of popularity in online platforms; (ii) studying how different dimensions’
interplay; and (iii) proposing dipm-SC, a novel algorithm that extends the
state-of-the-art models allowing to cluster multidimensional time series of
events at a time.

First, we compared our method with other extensions of univariate time series
clustering algorithms on synthetic data, where we successfully demonstrated
its higher accuracy. We then applied our framework on the multivariate time
series data deriving from two major online platforms: GitHub and Twitter.
Through these applications, we showcased the efficacy of our method in
uncovering fundamental pattern of popularity in online contexts. Our method
can indeed be used not only to find the overall difference between time series
shapes, such as bursty versus steady behaviors, but it can also uncover
differences in these two main trends depending on the time in which popularity
is acquired. Moreover, the results provided by our approach are easily
interpretable and sheds light on analyzing the interplay of multiple
dimensions. In the Twitter scenario, we found that the uncovered clusters are
temporally coherent in the hashtags used and related to certain types of
events, such as the Parkland school shooting and boycott campaigns.

In conclusion, we devised a methodology that extend the state-of-the-art
literature in the area of multivariate time series clustering and uses a
similarity metric based on the shape of the time series. Moreover, our method
can be used both to uncover fundamental patterns based solely on the shapes of
the time series and on the temporal occurrence of the events. Our approach
also provide a way of studying popularity of online content and understanding
their dynamics over time.

####  7.0.1 Acknowledgements

The authors are grateful to the Defense Advanced Research Projects Agency
(DARPA), contract W911NF-17-C-0094, for their support.

## References

  * Adar et al. [2004] Adar, E., Zhang, L., Adamic, L.A., Lukose, R.M.: Implicit structure and the dynamics of blogspace.  In: Workshop on the weblogging ecosystem. vol. 13, pp. 16989–16995 (2004)
  * Aghabozorgi et al. [2015] Aghabozorgi, S.R., Shirkhorshidi, A.S., Teh, Y.W.: Time-series clustering - a decade review.  Inf. Syst. 53, 16–38 (2015)
  * Alon et al. [2003] Alon, J., Sclaroff, S., Kollios, G., Pavlovic, V.: Discovering clusters in motion time-series data.  In: 2003 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2003. Proceedings. vol. 1 (June 2003)
  * Bakshy et al. [2011] Bakshy, E., Hofman, J.M., Mason, W.A., Watts, D.J.: Everyone’s an influencer: quantifying influence on twitter.  In: Fourth ACM international conference on Web search and data mining. pp. 65–74. ACM (2011)
  * Barabási [2005] Barabási, A.L.: The origin of bursts and heavy tails in human dynamics.  Nature 435, 207–211 (2005)
  * Bergroth et al. [2000] Bergroth, L., Hakonen, H., Raita, T.: A survey of longest common subsequence algorithms.  In: Proceedings of the Seventh International Symposium on String Processing Information Retrieval (SPIRE’00). pp. 39– (2000)
  * Chakrabarti et al. [2002] Chakrabarti, K., Keogh, E., Mehrotra, S., Pazzani, M.: Locally adaptive dimensionality reduction for indexing large time series databases.  ACM Trans. Database Syst. 27(2), 188–228 (Jun 2002)
  * Chu and Wong [1999] Chu, K.K.W., Wong, M.H.: Fast time-series searching with scaling and shifting.  In: Proc. of the 18th ACM SIGMOD-SIGACT-SIGART Symposium on Principles of Database Systems. pp. 237–248. PODS ’99, ACM (1999)
  * Crane and Sornette [2008] Crane, R., Sornette, D.: Robust dynamic classes revealed by measuring the response function of a social system.  Proceedings of the National Academy of Sciences 105(41), 15649–15653 (2008)
  * Davis et al. [2016] Davis, C.A., Varol, O., Ferrara, E., Flammini, A., Menczer, F.: Botornot: A system to evaluate social bots.  In: Proceedings of the 25th International Conference Companion on World Wide Web. pp. 273–274 (2016)
  * Ding et al. [2008] Ding, H., Trajcevski, G., Scheuermann, P., Wang, X., Keogh, E.: Querying and mining of time series data: Experimental comparison of representations and distance measures.  Proc. VLDB Endow. 1(2), 1542–1552 (Aug 2008)
  * Faloutsos et al. [1994] Faloutsos, C., Ranganathan, M., Manolopoulos, Y.: Fast subsequence matching in time-series databases.  In: Proc. of the 1994 ACM SIGMOD International Conference on Management of Data. pp. 419–429 (1994)
  * Fontes and Pereira [2016] Fontes, C.H., Pereira, O.: Pattern recognition in multivariate time series – a case study applied to fault detection in a gas turbine.  Engineering Applications of Artificial Intelligence 49, 10 – 18 (2016)
  * Ghassempour et al. [2014] Ghassempour, S., Girosi, F., Maeder, A.: Clustering multivariate time series using hidden markov models.  International Journal of Environmental Research and Public Health 11(3), 2741–2763 (2014)
  * Goh and Barabási [2008] Goh, K.I., Barabási, A.L.: Burstiness and memory in complex systems.  EPL (Europhysics Letters) 81(4), 48002 (2008)
  * Golub and van Loan [2013] Golub, G.H., van Loan, C.F.: Matrix Computations.  JHU Press, fourth edn. (2013)
  * Hallac et al. [2017] Hallac, D., Vare, S., Boyd, S., Leskovec, J.: Toeplitz inverse covariance-based clustering of multivariate time series data.  In: Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. pp. 215–223 (2017)
  * Hamerly and Elkan [2003] Hamerly, G., Elkan, C.: Learning the k in k-means.  In: Proceedings of the 16th International Conference on Neural Information Processing Systems. pp. 281–288. NIPS’03, MIT Press, Cambridge, MA, USA (2003)
  * Hartigan and Hartigan [1985] Hartigan, J.A., Hartigan, P.M.: The dip test of unimodality.  The Annals of Statistics 13(1), 70–84 (1985)
  * Kalogeratos and Likas [2012] Kalogeratos, A., Likas, A.: Dip-means: An incremental clustering method for estimating the number of clusters.  In: Proceedings of the 25th International Conference on Neural Information Processing Systems. pp. 2393–2401. NIPS’12, USA (2012)
  * Kaufman and Rousseeuw [2008] Kaufman, L., Rousseeuw, P.: Finding Groups in Data: an Introduction to Cluster Analysis.  Wiley (2008)
  * Keogh and Lin [2005] Keogh, E., Lin, J.: Clustering of time-series subsequences is meaningless: Implications for previous and future research.  Knowl. Inf. Syst. 8(2), 154–177 (Aug 2005)
  * Lehmann et al. [2012] Lehmann, J., Gonçalves, B., Ramasco, J.J., Cattuto, C.: Dynamical classes of collective attention in twitter.  In: Proceedings of the 21st international conference on World Wide Web. pp. 251–260. ACM (2012)
  * Lerman and Ghosh [2010] Lerman, K., Ghosh, R.: Information contagion: An empirical study of the spread of news on digg and twitter social networks.  Icwsm 10, 90–97 (2010)
  * Leskovec et al. [2007] Leskovec, J., McGlohon, M., Faloutsos, C., Glance, N., Hurst, M.: Patterns of cascading behavior in large blog graphs.  In: Proceedings of the 2007 SIAM international conference on data mining. pp. 551–556. SIAM (2007)
  * Längkvist et al. [2014] Längkvist, M., Karlsson, L., Loutfi, A.: A review of unsupervised feature learning and deep learning for time-series modeling.  Pattern Recognition Letters pp. 11–24 (2014)
  * Matsubara et al. [2012] Matsubara, Y., Sakurai, Y., Prakash, B.A., Li, L., Faloutsos, C.: Rise and fall patterns of information diffusion: Model and implications.  In: Proceedings of the 18th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. pp. 6–14. KDD ’12, ACM (2012)
  * McGlohon et al. [2007] McGlohon, M., Leskovec, J., Faloutsos, C., Hurst, M., Glance, N.: Finding patterns in blog shapes and blog evolution (2007)
  * Morstatter et al. [2013] Morstatter, F., Pfeffer, J., Liu, H., Carley, K.M.: Is the sample good enough? comparing data from Twitter’s streaming API with Twitter’s Firehose.  Proceedings of ICWSM (2013)
  * Paparrizos and Gravano [2016] Paparrizos, J., Gravano, L.: k-shape: Efficient and accurate clustering of time series.  SIGMOD Rec. 45(1), 69–76 (Jun 2016)
  * Paparrizos and Gravano [2017] Paparrizos, J., Gravano, L.: Fast and accurate time-series clustering.  ACM Trans. Database Syst. 42(2), 8:1–8:49 (Jun 2017)
  * Pelleg and Moore [2000] Pelleg, D., Moore, A.: X-means: Extending k-means with efficient estimation of the number of clusters.  In: In Proceedings of the 17th International Conf. on Machine Learning. pp. 727–734. Morgan Kaufmann (2000)
  * Petitjean et al. [2011] Petitjean, F., Ketterlin, A., Gançarski, P.: A global averaging method for dynamic time warping, with applications to clustering.  Pattern Recogn. 44(3), 678–693 (Mar 2011)
  * Ratkiewicz et al. [2010] Ratkiewicz, J., Fortunato, S., Flammini, A., Menczer, F., Vespignani, A.: Characterizing and modeling the dynamics of online popularity.  Physical review letters 105(15), 158701 (2010)
  * Rizoiu et al. [2017] Rizoiu, M.A., Xie, L., Sanner, S., Cebrian, M., Yu, H., Van Hentenryck, P.: Expecting to be hip: Hawkes intensity processes for social media popularity.  In: Proceedings of the 26th International Conference on World Wide Web. pp. 735–744 (2017)
  * Shokoohi-Yekta et al. [2017] Shokoohi-Yekta, M., Hu, B., Jin, H., Wang, J., Keogh, E.: Generalizing dtw to the multi-dimensional case requires an adaptive approach.  Data Mining and Knowledge Discovery 31(1), 1–31 (Jan 2017)
  * Thelwall et al. [2010] Thelwall, M., Buckley, K., Paltoglou, G., Cai, D., Kappas, A.: Sentiment strength detection in short informal text.  J. Am. Soc. Inf. Sci. Technol. 61(12), 2544–2558 (Dec 2010)
  * Ulanova et al. [????] Ulanova, L., Begum, N., Keogh, E.: Scalable clustering of time series with u-shapelets.  In: Proceedings of the 2015 SIAM International Conference on Data Mining. pp. 900–908 (????)
  * Varol et al. [2017] Varol, O., Ferrara, E., Davis, C., Menczer, F., Flammini, A.: Online human-bot interactions: Detection, estimation, and characterization (2017)
  * Wang et al. [2013] Wang, X., Mueen, A., Ding, H., Trajcevski, G., Scheuermann, P., Keogh, E.: Experimental comparison of representation methods and distance measures for time series data.  Data Mining and Knowledge Discovery 26(2), 275–309 (Mar 2013)
  * Wu and Huberman [2007] Wu, F., Huberman, B.A.: Novelty and collective attention.  Proceedings of the National Academy of Sciences 104(45), 17599–17601 (2007)
  * Xiong and Yeung [2002] Xiong, Y., Yeung, D.Y.: Mixtures of arma models for model-based time series clustering.  In: Proceedings of the 2002 IEEE International Conference on Data Mining. pp. 717–. ICDM ’02, IEEE Computer Society, Washington, DC, USA (2002)
  * Yang and Leskovec [2011a] Yang, J., Leskovec, J.: Patterns of temporal variation in online media.  In: Proceedings of the Fourth ACM International Conference on Web Search and Data Mining. pp. 177–186. WSDM ’11, ACM (2011a)
  * Yang and Leskovec [2011b] Yang, J., Leskovec, J.: Patterns of temporal variation in online media.  In: Proceedings of the fourth ACM international conference on Web search and data mining. pp. 177–186. ACM (2011b)
  * Yang and Shahabi [2004] Yang, K., Shahabi, C.: A pca-based similarity measure for multivariate time series.  In: Proceedings of the 2Nd ACM International Workshop on Multimedia Databases. pp. 65–74. MMDB ’04, ACM, New York, NY, USA (2004)
  * Ye and Keogh [2009] Ye, L., Keogh, E.: Time series shapelets: A new primitive for data mining.  In: Proceedings of the 15th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. pp. 947–956. KDD ’09, ACM (2009)
  * Zhang et al. [2016] Zhang, Q., Wu, J., Yang, H., Tian, Y., Zhang, C.: Unsupervised feature learning from time series.  In: Proceedings of the Twenty-Fifth International Joint Conference on Artificial Intelligence. pp. 2322–2328. IJCAI’16, AAAI Press (2016)