# Tailoring Dynamical Fermionization: Delta kick cooling of a Tonks-Girardeau
gas

Léonce Dupays<EMAIL_ADDRESS>Department of Physics and Materials
Science, University of Luxembourg, L-1511 Luxembourg, G. D. Luxembourg  Jing
Yang<EMAIL_ADDRESS>Department of Physics and Materials Science, University
of Luxembourg, L-1511 Luxembourg, G. D. Luxembourg  Adolfo del Campo
<EMAIL_ADDRESS>Department of Physics and Materials Science, University
of Luxembourg, L-1511 Luxembourg, G. D. Luxembourg  Donostia International
Physics Center, E-20018 San Sebastián, Spain

###### Abstract

In one spatial dimension, quantum exchange statistics and interactions are
inextricably intertwined. As a manifestation, the expansion dynamics of a
Tonks-Girardeau gas is characterized by dynamical fermionization (DF), whereby
the momentum distribution approaches that of a spin-polarized Fermi gas. Using
a phase-space analysis and the unitary evolution of the one-body reduced
density matrix, we show that DF can be tailored and reversed, using a
generalization of delta kick cooling (DKC) to interacting systems,
establishing a simple protocol to rescale the initial momentum distribution.
The protocol applies to both expansions and compressions and can be used for
microscopy of quantum correlations.

In one spatial dimension, the interchange of particle coordinates inevitably
leads to their interactions. In the process, the quantum state accumulates a
phase shift that is the sum of two contributions, one stemming from quantum
statistics, and the other from scattering. This fact makes it possible to
relate the physics of some strongly correlated quantum gases to that of
noninteracting systems with different quantum statistics. This is the basis of
the Bose-Fermi duality introduced by Girardeau in 1960 to describe a gas of
one-dimensional hard-core bosons, which is now known as the Tonks-Girardeau
(TG) gas Girardeau (1960); Cazalilla  _et al._ (2011). The latter can be
described in terms of a spin-polarized one-dimensional Fermi gas with no
interactions. Pauli exclusion principle in the Fermi gas makes the
wavefunction vanish at contact, a feature shared by the TG gas due to the
hard-core interactions. The wavefunctions of the two systems are identical for
a given particle ordering and differ only in their symmetrization. The bosonic
TG wavefunction $\Psi_{\rm TG}$ can be obtained from that of the Fermi gas
$\Psi_{\rm F}$ by explicit symmetrization according to the Bose-Fermi mapping
$\Psi_{\rm TG}=\prod_{j<k}{\rm sgn}(x_{k}-x_{j})\Psi_{\rm F}$. Similar
relations exist in systems governed by strongly-attractive $p$-wave
interactions Granger and Blume (2004); Girardeau  _et al._ (2004); Girardeau
and Minguzzi (2006). One-dimensional dualities can be extended to general
exchange statistics, in and out of equilibrium Girardeau (2006); del Campo
(2008), paving the way to the experimental realization of hard-core anyons
Keilmann  _et al._ (2011). They can be further used to describe systems with
finite-interaction strength Buljan  _et al._ (2008); Batchelor  _et al._
(2008); Guan  _et al._ (2013a), mixtures Girardeau and Minguzzi (2007), and
spinor systems Deuretzbacher  _et al._ (2008), among other examples Cazalilla
_et al._ (2011); Guan  _et al._ (2013b).

The TG gas can be considered as the strong-coupling limit of the Lieb-Liniger
(LL) gas, which describes one-dimensional bosons subject to contact
interactions of finite-strength $c$ Lieb and Liniger (1963); Lieb (1963);
Olshanii (1998). This model is integrable and solvable by Bethe ansatz. The
relevance of the LL gas to ultracold atom physics was established by Olshanii,
who showed that ultracold atoms in tight waveguides are described by the LL
model with tunable coupling constant $c$ Olshanii (1998). The strongly
interacting limit $c\rightarrow+\infty$ leads to the TG regime, realized
experimentally by making use of an optical lattice Kinoshita  _et al._ (2004,
2005); Cazalilla (2004). The connection between the continuum and lattice
version of hard-core bosons is well-understood Cazalilla (2004), and dynamical
correlations are in one-to-one correspondence at low densities.

Local correlations such as the density profile are indistinguishable between
dual systems. By contrast, correlations depending on the off-diagonal elements
of the density matrix exhibit clear signatures of quantum statistics Penrose
and Onsager (1956); Lenard (1964). A prominent example is the momentum
distribution. While that of a 1D Fermi gas in the ground state exhibits a
characteristic flat profile, that of the TG gas is sharply peaked at $k=0$ and
has power-law tails decaying as $1/k^{4}$. Out of equilibrium, it was
predicted that a TG gas under free expansion in 1D exhibits dynamical
fermionization (DF), with the asymptotic momentum distribution of the TG gas
matching that of free fermions Rigol and Muramatsu (2005); Minguzzi and
Gangardt (2005). This phenomenon has been recently observed in the laboratory
for the first time Wilson  _et al._ (2020). DF also governs the asymptotic
behavior of an expanding LL gas, which enters the TG regime Jukić  _et al._
(2008, 2009). While it is conveniently described using scale-invariance, which
makes the density profile at different times self-similar, it does not rely on
it, and occurs whether the initial confinement is harmonic or not del Campo
and Muga (2006); del Campo (2008); Campbell  _et al._ (2015). Generalizations
of this phenomenon have been reported for a fermionic analog of the TG gas
Girardeau and Minguzzi (2006), hardcore anyons del Campo (2008), and spinor
quantum gases Alam  _et al._ (2021). DF is generally justified as a result of
free expansion along the axial direction: as the particle density decreases,
the asymptotic momentum distribution is that of the rapidities, which are the
conserved quantities in a many-body integrable quantum system Sutherland
(1998); Rigol and Muramatsu (2005); Minguzzi and Gangardt (2005); del Campo
and Muga (2006); del Campo (2008); Bolech  _et al._ (2012); Campbell  _et al._
(2015); Mei  _et al._ (2016).

In this Letter, we analyze DF in phase space in arbitrary scale-invariant
processes, showing that its appearance is not restricted to expansions but can
occur as well in an implosion protocol leading to a density increase. While
the momentum distribution and density profile of dual systems under DF become
equal, the one-body reduced density matrix is shown to evolve unitarily,
making the distinguishability of the corresponding quantum states independent
of time. As a result, DF can be reversed, making use of a generalization of
delta kick cooling (DKC) to interacting systems, pulsing an external
potential. This allows to engineer protocols that rescale the momentum
distribution for microscopy of quantum correlations.

TG gas in a time-dependent trap.— Consider a TG gas in a harmonic trap, dual
to an ideal Fermi gas in the same confinement Girardeau  _et al._ (2001). In
the ground state, the TG wavefunction is the absolute value of the fermionic
one, which is given by a Slater determinant, e.g., $\Psi_{\rm
F}(x_{1},\dots,x_{\rm N})=\frac{1}{\sqrt{{\rm N}!}}\det_{n=0,k=1}^{{\rm
N}-1,{\rm N}}[\phi_{n}(x_{k})]$ in terms of the single-particle harmonic
oscillator eigenstates. Both systems are scale invariant with dimensionality
$D=1$, and their time-dependent coherent states take the form Sutherland
(1998); Minguzzi and Gangardt (2005); Gritsev  _et al._ (2010); del Campo
(2011)

$\displaystyle\Psi\left(t\right)$ $\displaystyle=$
$\displaystyle\frac{1}{b^{\frac{{\rm N}}{2}}}\exp\left[i\frac{m\dot{b}}{2\hbar
b}\sum_{i=1}^{\rm N}x_{i}\,^{2}-i\int_{0}^{t}\frac{E(0)}{\hbar
b(t^{\prime})^{2}}dt^{\prime}\right]$ (1)
$\displaystyle\times\Psi\left(\frac{x_{1}}{b},\dots,\frac{x_{\rm
N}}{b},t=0\\!\right)\,,$

where the scaling factor $b(t)>0$ is the solution of the Ermakov equation
$\ddot{b}+\omega(t)^{2}b=\omega_{0}^{2}/b^{3}$ with the initial conditions
$b(0)=1$, $\dot{b}(0)=0$. Note that this scaling law is not restricted to the
ground-state but it is shared by any many-body eigenstate $\Psi(0)$ with
energy eigenvalue $E(0)$. Quantities derived from $|\Psi(t)|^{2}$ are shared
by dual systems related by the Bose-Fermi mapping, given that
$\prod_{j<k}[{\rm sgn}(x_{k}-x_{j})]^{2}=1$. By contrast, those dependent on
the coherence in real space generally differ. We focus on the one-body reduced
density matrix (OBRDM), that contains all the information required to analyze
one-body observables. It is defined as $\rho_{1}(x,x^{\prime},t)={\rm N}\int
dx_{2}\cdots dx_{\rm N}\Psi\left(x,x_{2}\cdots x_{\rm
N},t\right)\Psi^{*}\left(x^{\prime},x_{2}\cdots x_{\rm N},t\right)$. From it,
one can determine the density profile $\rho(x,t)=\rho_{1}(x,x,t)$, as well as
the momentum distribution, making use of the Fourier transform
$n(p,t)=\frac{1}{2\pi\hbar}\int
dxdx^{\prime}e^{-ip(x-x^{\prime})/\hbar}\rho_{1}(x,x^{\prime},t)$. Using (1),
the OBRDM evolves according to

$\displaystyle\rho_{1}(x,x^{\prime},t)=\frac{1}{b}\exp\left[i\frac{m\dot{b}}{2\hbar
b}(x^{2}-x^{\prime}\,{}^{2})\right]\rho_{1}\left(\frac{x}{b},\frac{x^{\prime}}{b},t=0\right).$
(2)

Similar relations hold for the time-evolution of higher-order reduced density
matrices. In the limit of adiabatic driving $\dot{b}/b\to 0$, the OBRDM is
rescaled as $\rho_{1}(x,x^{\prime},t)=\rho_{1}(x/b,x^{\prime}/b,0)/b$. The use
of controlled expansions involving time-dependent traps and engineered by
shortcuts to adiabaticity has been proposed for implementing such scaling
without the requirement of slow driving, but generally involve time-dependent
traps del Campo (2011); del Campo and Boshier (2012). Such protocols realize
in essence a dynamical microscope zooming in on correlations in the OBRDM.

Shared unitary evolution of the OBRDMs and its consequences.—Interestingly,
Eq. (2) indicates that the evolution of OBRDMs for both the TG gas and the
spin-polarized ideal Fermi gas is unitary. More precisely, we introduce a
generic label $A=\\{{\rm TG,F}\\}$ for any of the dual systems and define the
corresponding quantum state $\sigma^{A}=\frac{1}{{\rm N}}\int
dxdx^{\prime}\rho_{1}^{A}(x,x^{\prime})|x\rangle\langle x^{\prime}|$ such that
${\rm Tr}(\sigma^{A})=1$. Then Eq. (2) implies
$\sigma^{A}(t)=U(t)\sigma^{A}(0)U^{\dagger}(t)$, where

$\displaystyle U(t)=\exp\left[i\frac{m\dot{b}}{2\hbar
b}x^{2}\right]\exp\left[-i\frac{{\rm ln}b}{2\hbar}(xp+px)\right],$ (3)

and the rightmost term is the dilatation operator implementing a scaling
transformation in real space by a factor $b$. Note that $U(t)$ is the same for
both $A=\\{{\rm TG,F}\\}$. In fact, under scale invariant dynamics, the
evolution of the quantum state associated with any $s$-body reduced density
matrix is also unitary SM .

The identical unitary evolution of the OBRDMs has several consequences: (i)
The spectral decomposition of the OBRDM is of the form
$\rho_{1}(x,x^{\prime},t)=\sum_{\mu}\lambda_{\mu}(0)\phi_{\mu}(x,t)\phi_{\mu}(x^{\prime},t)$,
where the eigenvalues $\lambda_{\mu}(0)$ of the OBRDM, which correspond to the
occupation numbers of the natural orbitals $\phi_{\mu}(x,t)$, are constant in
time. It follows that the dynamics is isentropic, i.e., it preserves the von
Neumann entropy $S(\sigma^{A})=-{\rm Tr}[\sigma^{A}\log\sigma^{A}]$. The
situation in the continuum is thus in contrast with that reported for hard-
core bosons in an optical lattice Rigol and Muramatsu (2005). In addition, the
natural orbitals fulfill the relation
$\phi_{\mu}(x,t)=\exp\left[i\frac{m\dot{b}}{2\hbar
b}x^{2}\right]\phi_{\mu}(x/b,0)/\sqrt{b}$. (ii) Consider the Ulhmann fidelity
defined as $\mathcal{F}(\sigma,\sigma^{\prime})={\rm
Tr}\left(\sqrt{\sqrt{\sigma}\sigma^{\prime}\sqrt{\sigma}}\right)$ Uhlmann
(1992); Nielsen and Chuang (2000) as a distinguishablity measure between the
quantum states $\sigma$ and $\sigma^{\prime}$. Given that the Uhlmann fidelity
is invariant under conjugation of its arguments by a common unitary, it
follows that

$\mathcal{F}[U(t)\sigma^{\rm TG}(0)U^{\dagger}(t),U(t)\sigma{{}^{\rm
F}}(0)U^{\dagger}(t)]=\mathcal{F}[\sigma^{\rm TG}(0),\sigma^{\rm F}(0)].$ (4)

As a result, even if the density profile, and the momentum distribution under
DF, are shared by both dual systems, their quantum states remain equally
distinguishable at all times. (iii) The unitary evolution in quantum mechanics
can be represented by a kernel in the phase space, convoluted with the initial
Wigner function of the system Garcia-Calderon and Moshinsky (1980). Under
scale-invariant dynamics, the convolution simplifies to a linear canonical
transformation on the initial Wigner function SM , which motivates the phase-
space analysis of DF. (iv) Finally, since the unitary evolution is invertible,
this paves the way to reverse DF, which we shall discuss in what follows.

Phase-space analysis of DF.— The Wigner function associated with the OBRDM can
be represented as a function of the coordinate $x$ and the canonically-
conjugated momentum $p$ Wigner (1932); Hillery  _et al._ (1984),

$\displaystyle
W(x,p,t)=\frac{1}{\pi\hbar}\int_{-\infty}^{\infty}\rho_{1}(x-y,x+y,t)e^{2ipy/\hbar}dy\,.$
(5)

The marginals of $W(x,p,t)$ correspond to the density profile
$\rho(x,t)=\rho_{1}(x,x,t)=\int W(x,p,t)dp$ and the momentum distribution
$n(p,t)=\rho_{1}(p,p,t)=\int W(x,p,t)dx$. From the dynamics (2) of the OBRDM,
following an arbitrary modulation of the trapping frequency $\omega(t)$, the
exact time-evolution of the Wigner function reads Shanahan  _et al._ (2018);
SM

$\displaystyle W(x,p,t)$ $\displaystyle=$ $\displaystyle
W\left(\frac{x}{b},bp-m\dot{b}x,t=0\right),$ (6)

where we note that $W$ does not need to be positive, i.e., it can describe a
non-classical state. Note that this evolution is common to all scale-invariant
systems, e.g., such as the single-particle time-dependent harmonic oscillator.
The Wigner function is rescaled, stretching (compressing) the density profile
along the $x$-axis, and compressing (stretching) the momentum distribution
along the $p$-axis if $b(t)>1$ ($b(t)<1$). The supplementary term $-m\dot{b}x$
involves a shift in phase-space, which induces DF.

Both for the TG and Fermi gas, the density profile exhibits explicitly the
scale invariance, while the asymptotic momentum distribution can be related to
the initial density profile SM :

$\displaystyle\rho(x,t)=\frac{1}{b}\rho\left(\frac{x}{b},0\right),n(p,t)\approx\frac{1}{m\dot{b}}\rho\left(\frac{p}{m\dot{b}},0\right).$
(7)

Note that the relation between $n(p,t)$ and $\rho(x,0)$ is different from what
is known in time-of-flight imaging, connecting $\rho(x,t)$ to $n(p,0)$. The
results in (7) are consistent with previous studies limited to sudden
expansions Jukić  _et al._ (2009); Campbell  _et al._ (2015). Note the first
equation of Eq. (7) is exact while the second is approximate, which requires
that the initial Wigner function decays over the characteristic spread $\Delta
p$ so that during the dynamics $\Delta p\ll x_{0}m\dot{b}b$ with
$x_{0}=\sqrt{\hbar/(m\omega_{0})}$, which is equivalent to the condition
$\dot{b}b\gg 2\omega_{0}$ taking $\Delta p\approx 2p_{0}$, where
$p_{0}=\hbar/x_{0}$. For a rigorous asymptotic analysis, see SM . In the
special case of a harmonic trap, the initial density profile of the TG gas
(and its dual system, the spin-polarized Fermi gas) can be expressed in terms
of the rescaled momentum distribution of the Fermi gas SM , making the term
“DF” natural in this setting Rigol and Muramatsu (2005); Minguzzi and Gangardt
(2005). For the sake of demonstration, we consider the expansion of a TG gas
initially confined in a harmonic trap with frequency $\omega_{0}$ that is
suddenly released in a wider trap with frequency $\omega_{1}$, illustrated in
Fig. 1. This process leads to periodic time-dependence of the scaling factor
$b(t)$ Minguzzi and Gangardt (2005); Dupays  _et al._ (2021); see SM  for
other protocols. The width of the cloud is controlled by $b(t)$ and oscillates
after the release of the TG gas into the wider trap. This behavior induces DF
periodically, with the nonequilibrium momentum distribution evolving between
that of an equilibrium TG gas and a Fermi gas, as predicted in Minguzzi and
Gangardt (2005) and recently observed experimentally Wilson  _et al._ (2020).
As emphasized, the phase-space dynamics in Eq. (6) holds for all scale-
invariant systems and is responsible for the asymptotic form of the momentum
distribution can be related to the initial density profile, whenever the width
of the atomic cloud varies swiftly. The fact that the asymptotic momentum
distribution can be also exactly related to that of the spin-polarized Fermi
gas is a specific feature of the harmonic confinement, as single-particle
eigenstates $\phi_{n}(x)$ are in this case expressed in terms of Hermite
polynomials, that are eigenstates of the Fourier transform. Note that despite
the coincidence of the marginals of the Wigner functions of the TG gas and the
Fermi gas at late times, the OBRDMs remain distinguishable as discussed above.

Figure 1: DF of a TG gas in a sudden expansion in phase space. Under scale-
invariant dynamics, DF results from a canonical transformation that describes
the evolution of the Wigner function. The density profile is scale-invariant
at all times. In addition, when the rate of change of the scaling factor is
large, the asymptotic momentum distribution can be related to the initial
density profile. The initial Wigner function is peaked along the axis $p=0$,
with riddles on both sides that take negative values. We chose
$\omega_{0}=5\omega_{1}$. At time $t_{1}=\frac{3\pi}{4\omega_{1}}$ the Wigner
function is rotated and dilated in the phase-space, so that the momentum
distribution corresponds to the rescaled density profile. We chose ${\rm
N}=10$.

Tailoring and reversing DF with kicks.— The momentum shift in (6) is
responsible for DF. Classically, one may expect to cancel it by applying a
conservative force for a short period of time $\tau_{k}$ inducing a momentum
change $\delta p=-\tau_{k}\partial_{x}V(x)$, i.e., pulsing an external
potential $V(x)$. This argument, limited to classical noninteracting systems,
is the basis of delta-kick cooling (DKC) Chu  _et al._ (1986); Ammann and
Christensen (1997); Morinaga  _et al._ (1999). In what follows, we make use of
the extension of DKC for scale-invariant interacting systems. Under Eq. (2),
excitations encoded in the phase factor proportional to $\dot{b}/b$, that are
responsible for DF, can be explicitly canceled in an interacting system by
applying a kick potential of appropriate strength. Canceling the phase allows
to tailor the momentum distribution and reverse DF. Given that the phase
oscillation in Eq. (1) is quadratic in the coordinates it can be canceled by
pulsing an external harmonic trap with a given frequency $\omega_{k}$. To this
end, consider the Hamiltonian with a $\delta$-kick applied at $t_{k}$

$\displaystyle
H_{k}(t)=H(t)+\delta(t-t_{k})\frac{1}{2}m\omega_{k}^{2}\sum_{i=1}^{{\rm
N}}x_{i}\,^{2}.$ (8)

The use of a delta function is justified when the duration of the pulse
$\tau_{k}$ is short with respect to other time scales Dupays  _et al._ (2021).
The corresponding time-evolution operator admits the factorization
$U_{\delta}(t_{F},0)=e^{-i\tau_{k}\frac{m\omega_{k}^{2}}{2\hbar}\sum_{i=1}^{{\rm
N}}x_{i}\,^{2}}U(t_{k},0)$, where $U(t,t^{\prime})$ is the propagator
associated with $H(t)$, and $\tau_{k}$ is a small time scale during which the
kick is applied. Considering the evolution from $t=0$ to time
$t_{k}+\tau_{k}$, one can choose the pulse parameters $\tau_{k}$ and
$\omega_{k}$ such that

$\displaystyle\tau_{k}\omega^{2}_{k}=\frac{\dot{b}(t_{k})}{b(t_{k})}\,.$ (9)

This requires pulsing a harmonic trap with $\omega_{k}>0$ in an expansion with
$\dot{b}(t_{k})>0$ and an inverted harmonic trap with purely imaginary
frequency $i\omega_{k}$ in a compression with $\dot{b}(t_{k})<0$. In either
case, the application of the kick cancels DF, and brings back the OBRDM to the
initial one up to a scaling of the coordinates with respect to $b$,

$\displaystyle\rho_{1}(x,x^{\prime},t_{k}+\tau_{k})$ $\displaystyle=$
$\displaystyle\frac{1}{b(t_{k})}\rho_{1}\left(\frac{x}{b(t_{k})},\frac{x\,^{\prime}}{b(t_{k})},t=0\right).$
(10)

Consequently, the momentum distribution after the kick is $n(p,t)=b\,n(pb,0)$,
and similarly the Wigner function reduces to $W(x,p,t)=W_{0}(x/b,bp,0)$ with
$b=b(t_{k})$. In turn, DKC prepares the same state that would have been
obtained under adiabatic dynamics, without the requirement of slow driving.

Imploding TG gas.— Intuitively the DF occurs for an expansion as the particle
density decreases. This is the case considered so far in theoretical and
experimental studies, in which it is possible to suppress DF by DKC as we show
in SM . Yet, the phase-space dynamics (6) also yields DF in a compression
process with $b(t_{F})<b(t=0)=1$ if the rate of change of the scaling factor
is fast enough so that $\dot{b}b$ is large. For the sake of illustration, we
consider a sudden compression protocol, where the trap of initial frequency
$\omega_{0}$ is compressed to a frequency $\omega_{1}>\omega_{0}$, leading to
the periodic scaling factor $b(t)$ displayed in Fig. 2. Large values of
$\dot{b}b$ induce a high-frequency phase modulation in the coordinate
representation.

Figure 2:  (a) The frequency modulation leading to the implosion protocol
relies on a sudden compression from the initial frequency $\omega_{0}$ to the
final frequency $\omega_{1}/5=\omega_{0}=1$. DF is reversed by making use of
pulsed attractive harmonic potential. For illustration, we chose the frequency
of the kick $\omega^{2}_{k}=1000\omega^{2}_{0}$, $t_{k}=3\pi/4\omega_{1}$ and
the final frequency determined by
$b_{F}=\sqrt{\omega_{0}/\omega_{F}}=\sqrt{1/5}$. (b) Evolution of the scaling
factor, its derivative, and the product $b\dot{b}$.

For the prescribed protocol, applying at the time $t_{k}$ a kick of duration
$\tau_{k}$, the required pulse parameters to reverse DF are set by

$\displaystyle\tau_{k}\omega^{2}_{k}=\frac{\omega_{1}\left(\omega^{2}_{0}-\omega^{2}_{1}\right)\sin(\omega_{1}t_{k})\cos(\omega_{1}t_{k})}{\left(\omega^{2}_{0}-\omega^{2}_{1}\right)\sin^{2}(\omega_{1}t_{k})+\omega^{2}_{1}}.$
(11)

The evolution of the the momentum distribution at different stages of the
protocol is shown in Fig. 3, for an initially confined TG gas undergoing an
implosion engineered by a sudden frequency increase. The oscillatory time-
dependence of $b(t)$ shown in Fig. 2 yields associated oscillations of the
momentum distribution, which exhibits DF exactly at the characteristic times
$t_{m}=\frac{(2m+1)\pi}{4\omega_{1}}$ with integer $m$. The required pulse
strength takes then the maximum value
$\left|\dot{b}(t_{m})b(t_{m})\right|=|\omega^{2}_{0}-\omega^{2}_{1}|$.
Application of a pulse satisfying the generalized DKC condition (11) gives
rise to the reversal of DF, rescaling the initial momentum distribution by a
factor $1/b(t_{k})$; see Fig. 3.

Figure 3: DKC of an imploding TG gas. (a) Momentum distribution of a trapped
TG gas at $t=0$, after the implosion at $t_{k}=\frac{3\pi}{4\omega_{1}}$,
after the implosion followed by a pulse of shorten duration $0.7\tau_{k}$
(green dashed-dotted line), and with the DKC duration $\tau_{k}$ satisfying
(11) (turquoise blue dashed line); ${\rm N}=7$. In dashed green the kick
strength is chosen as $\tau_{k}\omega^{2}_{k}$ in order to cancel DF, with
$\omega^{2}_{k}=1000\omega^{2}_{0}$. In choosing the strength of the kick
differently one can modulate DF. The implosion protocol is corresponds to the
sudden compression of the trap from $\omega_{0}$ to $\omega_{1}=5\omega_{0}$.

In summary, we have established the unitary character of the dynamics of the
OBRDM of a TG gas in a driven harmonic trap and point out its far-reaching
consequences. The time evolution is exactly isentropic and preserves the
occupation numbers of the natural orbitals. As a result, DF does not affect
the distinguishability between the OBRDMs of the dual systems, which is
independent of time, e.g., as quantified by the Uhlmann fidelity. For
arbitrary driving of the trap frequency, the unitary character of the dynamics
makes it possible to describe DF as a result of a canonical transformation in
phase space, that relates the asymptotic momentum distribution to the initial
density profile under rapid acceleration of the width of the atomic cloud.
This relation holds for expansions as well as compressions leading to an
increase of the interparticle density and is not restricted to the TG gas or
systems with hard-core interactions. Thanks to unitarity, DF can be further
tailored and completely reversed by applying a kick with a pulsed external
potential, generalizing DKC to interacting systems. This allows to rescale the
initial momentum distribution without the requirement of slow driving. Our
findings open the way to control nonequilibrium correlations in driven
ultracold gases and can be tested in laboratory settings used in recent
experiments Wilson  _et al._ (2020). They should be generalizable to atomic
mixtures, systems with $p$-wave interactions, finite coupling strength,
fractional exchange statistics, and spinor gases, among other examples.

Acknowledgements.— It is a pleasure to acknowledge discussions with Fernando
J. Gómez-Ruiz, Niklas Hörnedal, Federico Roccati, Naim Mackel, and Maxim
Olchanii.

## References

  * Girardeau (1960) M. Girardeau, Journal of Mathematical Physics 1, 516 (1960).
  * Cazalilla  _et al._ (2011) M. A. Cazalilla, R. Citro, T. Giamarchi, E. Orignac, and M. Rigol, Rev. Mod. Phys. 83, 1405 (2011).
  * Granger and Blume (2004) B. E. Granger and D. Blume, Phys. Rev. Lett. 92, 133202 (2004).
  * Girardeau  _et al._ (2004) M. Girardeau, H. Nguyen, and M. Olshanii, Optics Communications 243, 3 (2004).
  * Girardeau and Minguzzi (2006) M. D. Girardeau and A. Minguzzi, Phys. Rev. Lett. 96, 080404 (2006).
  * Girardeau (2006) M. D. Girardeau, Phys. Rev. Lett. 97, 100402 (2006).
  * del Campo (2008) A. del Campo, Phys. Rev. A 78, 045602 (2008).
  * Keilmann  _et al._ (2011) T. Keilmann, S. Lanzmich, I. McCulloch, and M. Roncaglia, Nature Communications 2, 361 (2011).
  * Buljan  _et al._ (2008) H. Buljan, R. Pezer, and T. Gasenzer, Phys. Rev. Lett. 100, 080406 (2008).
  * Batchelor  _et al._ (2008) M. T. Batchelor, X.-W. Guan, and A. Kundu, Journal of Physics A: Mathematical and Theoretical 41, 352002 (2008).
  * Guan  _et al._ (2013a) X.-W. Guan, M. T. Batchelor, and C. Lee, Rev. Mod. Phys. 85, 1633 (2013a).
  * Girardeau and Minguzzi (2007) M. D. Girardeau and A. Minguzzi, Phys. Rev. Lett. 99, 230402 (2007).
  * Deuretzbacher  _et al._ (2008) F. Deuretzbacher, K. Fredenhagen, D. Becker, K. Bongs, K. Sengstock, and D. Pfannkuche, Phys. Rev. Lett. 100, 160405 (2008).
  * Guan  _et al._ (2013b) X.-W. Guan, M. T. Batchelor, and C. Lee, Rev. Mod. Phys. 85, 1633 (2013b).
  * Lieb and Liniger (1963) E. H. Lieb and W. Liniger, Phys. Rev. 130, 1605 (1963).
  * Lieb (1963) E. H. Lieb, Phys. Rev. 130, 1616 (1963).
  * Olshanii (1998) M. Olshanii, Phys. Rev. Lett. 81, 938 (1998).
  * Kinoshita  _et al._ (2004) T. Kinoshita, T. Wenger, and D. S. Weiss, Science 305, 1125 (2004).
  * Kinoshita  _et al._ (2005) T. Kinoshita, T. Wenger, and D. S. Weiss, Phys. Rev. Lett. 95, 190406 (2005).
  * Cazalilla (2004) M. A. Cazalilla, Phys. Rev. A 70, 041604 (2004).
  * Penrose and Onsager (1956) O. Penrose and L. Onsager, Phys. Rev. 104, 576 (1956).
  * Lenard (1964) A. Lenard, Journal of Mathematical Physics 5, 930 (1964).
  * Rigol and Muramatsu (2005) M. Rigol and A. Muramatsu, Phys. Rev. Lett. 94, 240403 (2005).
  * Minguzzi and Gangardt (2005) A. Minguzzi and D. M. Gangardt, Phys. Rev. Lett. 94, 240404 (2005).
  * Wilson  _et al._ (2020) J. M. Wilson, N. Malvania, Y. Le, Y. Zhang, M. Rigol, and D. S. Weiss, Science 367, 1461 (2020).
  * Jukić  _et al._ (2008) D. Jukić, R. Pezer, T. Gasenzer, and H. Buljan, Phys. Rev. A 78, 053602 (2008).
  * Jukić  _et al._ (2009) D. Jukić, B. Klajn, and H. Buljan, Phys. Rev. A 79, 033612 (2009).
  * del Campo and Muga (2006) A. del Campo and J. G. Muga, Europhysics Letters (EPL) 74, 965 (2006).
  * Campbell  _et al._ (2015) A. S. Campbell, D. M. Gangardt, and K. V. Kheruntsyan, Phys. Rev. Lett. 114, 125302 (2015).
  * Alam  _et al._ (2021) S. S. Alam, T. Skaras, L. Yang, and H. Pu, Phys. Rev. Lett. 127, 023002 (2021).
  * Sutherland (1998) B. Sutherland, Phys. Rev. Lett. 80, 3678 (1998).
  * Bolech  _et al._ (2012) C. J. Bolech, F. Heidrich-Meisner, S. Langer, I. P. McCulloch, G. Orso, and M. Rigol, Phys. Rev. Lett. 109, 110602 (2012).
  * Mei  _et al._ (2016) Z. Mei, L. Vidmar, F. Heidrich-Meisner, and C. J. Bolech, Phys. Rev. A 93, 021607 (2016).
  * Girardeau  _et al._ (2001) M. D. Girardeau, E. M. Wright, and J. M. Triscari, Phys. Rev. A 63, 033601 (2001).
  * Gritsev  _et al._ (2010) V. Gritsev, P. Barmettler, and E. Demler, New Journal of Physics 12, 113005 (2010).
  * del Campo (2011) A. del Campo, Phys. Rev. A 84, 031606 (2011).
  * del Campo and Boshier (2012) A. del Campo and M. G. Boshier, Sci. Rep. 2, 648 (2012).
  * (38) See the Supplemental Material.
  * Uhlmann (1992) A. Uhlmann, Physics Letters A 161, 329 (1992).
  * Nielsen and Chuang (2000) M. Nielsen and I. Chuang,  _Quantum Computation and Quantum Information_ , Cambridge Series on Information and the Natural Sciences (Cambridge University Press, 2000).
  * Garcia-Calderon and Moshinsky (1980) G. Garcia-Calderon and M. Moshinsky, Journal of Physics A: Mathematical and General 13, L185 (1980).
  * Wigner (1932) E. Wigner, Phys. Rev. 40, 749 (1932).
  * Hillery  _et al._ (1984) M. Hillery, R. O’Connell, M. Scully, and E. Wigner, Physics Reports 106, 121 (1984).
  * Shanahan  _et al._ (2018) B. Shanahan, A. Chenu, N. Margolus, and A. del Campo, Phys. Rev. Lett. 120, 070401 (2018).
  * Dupays  _et al._ (2021) L. Dupays, D. C. Spierings, A. M. Steinberg, and A. del Campo, Phys. Rev. Research 3, 033261 (2021).
  * Chu  _et al._ (1986) S. Chu, J. E. Bjorkholm, A. Ashkin, J. P. Gordon, and L. W. Hollberg, Opt. Lett. 11, 73 (1986).
  * Ammann and Christensen (1997) H. Ammann and N. Christensen, Phys. Rev. Lett. 78, 2088 (1997).
  * Morinaga  _et al._ (1999) M. Morinaga, I. Bouchoule, J.-C. Karam, and C. Salomon, Phys. Rev. Lett. 83, 4037 (1999).
  * Schuch and Moshinsky (2008) D. Schuch and M. Moshinsky, SIGMA 4, 054 (2008).
  * Bleistein and Handelsman (1986) N. Bleistein and R. Handelsman,  _Asymptotic Expansions of Integrals_ , Dover Books on Mathematics Series (Dover Publications, 1986).
  * Pedri (2005) P. Pedri,  _Dynamical behavior of ultracold atomic gases_ , Ph.D. thesis, University of Trento (2005).
  * Pezer and Buljan (2007) R. Pezer and H. Buljan, Phys. Rev. Lett. 98, 240403 (2007).
  * Lapeyre  _et al._ (2002) G. J. Lapeyre, M. D. Girardeau, and E. M. Wright, Phys. Rev. A 66, 023606 (2002).
  * Forrester  _et al._ (2003) P. J. Forrester, N. E. Frankel, T. M. Garoni, and N. S. Witte, Phys. Rev. A 67, 043607 (2003).

Supplemental Material for
“Tailoring Dynamical Fermionization: Delta kick cooling of a Tonks-Girardeau
gas”

In this Supplemental Material, we present the time evolution of the OBRDM, the
Wigner function for the spin-polarized Fermi gas, the analysis of DF based on
saddle point approximation, the analysis of DF based on the Wigner function,
and protocols for tailoring DF and Delta kick cooling of a TG gas in a time-
dependent harmonic trap.

###### Contents

  1. I Time-evolution of the OBRDM
  2. II Representation of the scale-invariant unitary dynamics of the OBRDM by canonical transformation
  3. III Wigner function for the spin-polarized Fermi gas
  4. IV Analysis of Dynamical Fermionization in time-dependent traps
  5. V Dynamical fermionization in phase space: asymptotic marginal distributions of the Wigner function
  6. VI Tailoring Dynamical Fermionization and Delta kick cooling of a Tonks-Girardeau gas in a time-dependent harmonic trap
    1. VI.1 Ground-state OBRDM of a TG gas in a harmonic trap
    2. VI.2 Delta kick cooling of an expanding TG gas
    3. VI.3 Delta kick cooling of an imploding TG gas

##  I Time-evolution of the OBRDM

Consider a system in a harmonic trap with frequency $\omega_{0}$ which is
modulated as $\omega(t)$ for $t>0$. For scale-invariant dynamics, the
wavefunction of an energy eigenstate at $t=0$ evolves for $t>0$ according to
Gritsev  _et al._ (2010); del Campo (2011)

$\displaystyle\Psi(\vec{r}_{1},\dots,\vec{r}_{{\rm
N}},t)=\frac{1}{b^{\frac{DN}{2}}}\exp\left[i\frac{m\dot{b}}{2\hbar
b}\sum_{i=1}^{{\rm N}}\vec{r}_{i}\,^{2}-i\int_{0}^{t}\frac{E(0)}{\hbar
b^{2}(t^{\prime})}dt^{\prime}\right]\Psi\left(\frac{\vec{r}_{1}}{b},\cdots,\frac{\vec{r}_{{\rm
N}}}{b};0\right),$ (S1)

where ${\rm N}$ is the number of particles, $D$ is the dimension of the trap,
and $b(t)$ is the scaling factor determined by the Ermakov equation
$\ddot{b}+\omega(t)^{2}b=\omega_{0}^{2}/b^{3}$. Note that the coordinate-
dependent phase factor is common to all quantum states, while $E(0)$ is the
eigenvalue of the initial energy eigenstate $\Psi(0)$. In one dimension $D=1$,
the time-evolution of a many-particle eigenstate simply reads

$\displaystyle\Psi(t)=\frac{1}{b^{\frac{{\rm
N}}{2}}}\exp\left[i\frac{m\dot{b}}{2\hbar b}\sum_{i=1}^{{\rm
N}}x^{2}_{i}-i\int_{0}^{t}\frac{E(0)}{\hbar
b^{2}(t^{\prime})}dt^{\prime}\right]\Psi\left(\frac{x_{1}}{b},\cdots,\frac{x_{{\rm
N}}}{b},0\right).$ (S2)

The corresponding one-body reduced density matrix evolves according to

$\displaystyle\rho_{1}(x,x^{\prime};t)$ $\displaystyle=$ $\displaystyle{\rm
N}\int dx_{2}\cdots dx_{{\rm N}}\Psi(x,x_{2},\cdots,x_{{\rm
N}},t)\Psi^{*}(x^{\prime},x_{2},\cdots,x_{{\rm N}},t)$ (S3) $\displaystyle=$
$\displaystyle\exp\left[i\frac{m\dot{b}}{2\hbar b}(x^{2}-x^{\prime}\
{}^{2})\right]\frac{{\rm N}}{b}\int d\frac{x_{2}}{b}\cdots d\frac{x_{{\rm
N}}}{b}\Psi\left(\frac{x}{b},\frac{x_{2}}{b},\cdots,\frac{x_{{\rm
N}}}{b},t=0\right)\Psi^{*}\left(\frac{x^{\prime}}{b},\frac{x_{2}}{b},\cdots,\frac{x_{{\rm
N}}}{b},t=0\right)$ $\displaystyle=$
$\displaystyle\exp\left[i\frac{m\dot{b}}{2\hbar b}(x^{2}-x^{\prime}\
{}^{2})\right]\frac{1}{b}\rho_{1}\left(\frac{x}{b},\frac{x^{\prime}}{b},t=0\right).$

Using it, the momentum distribution is explicitly given by

$\displaystyle n(p,t)$ $\displaystyle=$ $\displaystyle\langle
p|\rho_{1}(t)|p\rangle$ (S4) $\displaystyle=$
$\displaystyle\frac{b}{2\pi\hbar}\int
dxdx\,^{\prime}\,\exp\left[i\frac{m\dot{b}b}{2\hbar}(x\,^{2}-x\,^{\prime}\,{}^{2})\right]e^{-i\frac{pb}{\hbar}(x-x\,^{\prime})}\rho_{1}\left(x,x\,^{\prime},t=0\right).$

Similar considerations apply to higher-order reduced density matrices that can
also be efficiently computed using scale-invariance. For compactness, we
introduce the notation $\vec{x}_{s}=(x_{1},\dots,x_{s})$, with
$s\in\mathbb{N}$. The time-evolution of the $s$-th reduced density matrix
reads

$\displaystyle\rho_{s}(\vec{x}_{s},\vec{y}_{s},t)$ $\displaystyle=$
$\displaystyle\frac{{\rm N}!}{({\rm N}-s)!}\int dx_{s+1}\dots dx_{{\rm N}}\
\Psi(\vec{x}_{s},x_{s+1},\dots,x_{{\rm
N}},t)\Psi(\vec{y}_{s},x_{s+1},\dots,x_{{\rm N}},t)^{*}.$ (S5)
$\displaystyle=$
$\displaystyle\frac{1}{b^{s}}\exp\left[i\frac{m\dot{b}}{2\hbar
b}(\vec{x}_{s}\,^{2}-\vec{y}_{s}\,^{2})\right]\rho_{s}\left(\frac{\vec{x}_{s}}{b},\frac{\vec{y}_{s}}{b},t=0\right).$
(S6)

Making use of DKC, upon application of the pulse at $t_{k}$ with duration
$\tau_{k}$, the dynamical phase is suppressed, and one recovers the same
density matrix that would be obtained in an adiabatic protocol

$\displaystyle\rho_{s}(\vec{x}_{s},\vec{y}_{s},t_{k}+\tau_{k})=\frac{1}{b^{s}}\rho_{s}\left(\frac{\vec{x}_{s}}{b},\frac{\vec{y}_{s}}{b},t=0\right),$
(S7)

with $b=b(t_{k})$. The joint characteristic function is obtained by the
$s$-dimensional Fourier transform,

$\displaystyle n_{s}(\vec{p}_{s},t)=\frac{1}{(2\pi\hbar)^{s}}\int
d\vec{x}_{s}d\vec{y}_{s}e^{-i\sum_{j=1}^{s}\frac{p_{j}(x_{j}-y_{j})}{\hbar}}\rho_{s}(\vec{x}_{s},\vec{y}_{s},t).$

After a DKC pulse, the characteristic function reads

$\displaystyle n_{s}(\vec{p}_{s},t_{k}+\tau_{k})$ $\displaystyle=$
$\displaystyle\frac{1}{(2\pi\hbar)^{s}}\int
d\vec{x}_{s}d\vec{y}_{s}\frac{1}{b^{s}}e^{-i\sum_{j=1}^{s}\frac{p_{j}(x_{j}-y_{j})}{\hbar}}\rho_{s}\left(\frac{\vec{x}_{s}}{b},\frac{\vec{y}_{s}}{b},t=0\right),$
(S9) $\displaystyle=$ $\displaystyle b^{s}n_{s}(b\vec{p}_{s},t=0),$ (S10)

with $b=b(t_{k})$.

From the equation of motions for the OBRDM and the $s$-order reduced density
matrix, it is possible to conclude that their evolution is described by the
unitary

$\displaystyle U^{(s)}(t)=\exp\left[i\frac{m\dot{b}}{2\hbar
b}\vec{x}_{s}^{2}\right]\exp\left[-i\frac{{\rm
ln}b}{2\hbar}(\vec{x}_{s}\cdot\vec{p}_{s}+\vec{p}_{s}\cdot\vec{x}_{s})\right].$
(S11)

##  II Representation of the scale-invariant unitary dynamics of the OBRDM by
canonical transformation

As shown in Ref. Garcia-Calderon and Moshinsky (1980), any classical canonical
transformation on the Wigner function corresponds to a unitary evolution in
quantum mechanics, which is represented by a kernel function in the phase
space. In particular, it was shown that for a linear canonical transformation,

$\displaystyle\begin{pmatrix}\vec{r}\\\
\vec{p}\end{pmatrix}=\begin{pmatrix}\alpha&\beta\\\
\gamma&\delta\end{pmatrix}\begin{pmatrix}\vec{r}\,^{\prime}\\\
\vec{p}\,^{\prime}\end{pmatrix},$ (S12)

with $\alpha\delta-\beta\gamma=1$, the kernel function in the quantum
mechanical evolution takes a very simple form,

$\displaystyle
K(\vec{r},\vec{p}|\vec{r}\,^{\prime},\vec{p}\,^{\prime})=\delta^{D}[\vec{r}\,^{\prime}-(\alpha\vec{r}+\beta\vec{p})]\delta^{D}[\vec{p}\,^{\prime}-(\gamma\vec{r}+\delta\vec{p})]\,,$

where $\delta^{D}$ the $D$-dimensional delta function. Thus, the evolution of
the Wigner function is dictated by Garcia-Calderon and Moshinsky (1980);
Schuch and Moshinsky (2008)

$\displaystyle W(\vec{r},\vec{p},t)=\iint
d^{D}\vec{r}\,^{\prime}d^{D}\vec{p}\,^{\prime}K(\vec{r},\vec{p}|\vec{r}\,^{\prime},\vec{p}\,^{\prime})W(\vec{r}\,^{\prime},\vec{p}\,^{\prime},0)=W(\alpha\vec{r}+\beta\vec{p},\gamma\vec{r}+\delta\vec{p},0).$
(S14)

Under scale-invariance, the time-dependent Wigner function can be solely
written in terms of the initial Wigner function:

$\displaystyle W(\vec{r},\vec{p},t)$ $\displaystyle=$
$\displaystyle\frac{1}{(\pi\hbar)^{D}}\int\frac{d^{D}\vec{y}}{b^{D}}\exp\left\\{i\frac{m\dot{b}}{2\hbar
b}[(\vec{r}-\vec{y})\,^{2}-(\vec{r}+\vec{y})\,^{2}]\right\\}\rho_{1}\left(\frac{\vec{r}-\vec{y}}{b},\frac{\vec{r}+\vec{y}}{b},0\right)e^{2i\vec{p}\cdot\vec{y}/\hbar},$
(S15) $\displaystyle=$
$\displaystyle\frac{1}{(\pi\hbar)^{D}}\int\frac{d^{D}\vec{y}}{b^{D}}\rho_{1}\left(\frac{\vec{r}-\vec{y}}{b},\frac{\vec{r}+\vec{y}}{b},0\right)e^{2i\left(b\vec{p}-m\dot{b}\vec{r}\right)\cdot\vec{y}/(b\hbar)},$
$\displaystyle=$ $\displaystyle
W\left(\frac{\vec{r}}{b},b\vec{p}-m\dot{b}\vec{r},0\right).$

Eqquation (S15) amounts to the following canonical transformation

$\displaystyle\begin{pmatrix}\alpha&\beta\\\
\gamma&\delta\end{pmatrix}=\begin{pmatrix}1/b&0\\\ -m\dot{b}&b\end{pmatrix}.$
(S16)

##  III Wigner function for the spin-polarized Fermi gas

It is of interest to compare the Wigner function of the TG gas in Fig. 1 to
the Wigner function of the dual, ideal Fermi gas plotted in Fig. 4. Given the
form of the OBRDM of the ideal Fermi gas $\rho^{\rm
F}_{1}(x,x^{\prime})=\sum_{n=0}^{{\rm
N}-1}\phi_{n}(x)\phi^{*}_{n}(x^{\prime})$, it follows that the Wigner function
is given as

$\displaystyle W^{\rm F}(x,p)=\sum_{n=0}^{{\rm N}-1}W_{n}(x,p),$ (S17)

in terms of the sum of the Wigner function of single-particle eigenstates of
the harmonic oscillator. The latter is well know and takes the form (see e.g.
Shanahan  _et al._ (2018))

$\displaystyle
W_{n}(x,p)=\frac{(-1)^{n}}{\pi\hbar}\exp\left[-\frac{2}{\hbar\omega_{0}}\left(\frac{p^{2}}{2m}+\frac{1}{2}m\omega^{2}x^{2}\right)\right]L_{n}\left[\frac{4}{\hbar\omega_{0}}\left(\frac{p^{2}}{2m}+\frac{1}{2}m\omega^{2}x^{2}\right)\right],$
(S18)

where $L_{n}(x)$ is the Laguerre polynomial. The Wigner function of an ideal
Fermi gas of ${\rm N}$ particles in a harmonic trap is shown in Fig. 4
together with the corresponding marginals: the density profile $\rho^{\rm
F}(x)$ and the momentum distribution $n^{\rm F}(p)$.

Figure 4: Wigner function of ground state of the Fermi gas, in a harmonic
trapping of frequency $\omega_{0}$. We note that the ellipsoïdal shape of the
Wigner function leads to riddles both in the density profile and in the
momentum distribution.

##  IV Analysis of Dynamical Fermionization in time-dependent traps

As we have seen in the preceding section, the momentum distribution can be
found in terms of the Fourier transform of the OBRDM
$n(p,t)=\frac{1}{2\pi\hbar}\int
dxdx\,^{\prime}\,e^{-ip(x-x\,^{\prime})/\hbar}\rho_{1}(x,x\,^{\prime};t)$. It
is useful to make the change of variables $x/b\to x$ in the distribution

$\displaystyle n(p,t)$ $\displaystyle=$ $\displaystyle\frac{b}{2\pi\hbar}\int
dxdx\,^{\prime}\,e^{i\phi(x,x^{\prime})}\rho_{1}\left(x,x\,^{\prime},t=0\right),$
(S19) $\displaystyle\phi(x,x^{\prime})$ $\displaystyle=$
$\displaystyle\frac{m\dot{b}b}{2\hbar}(x\,^{2}-x\,^{\prime}\,{}^{2})-\frac{pb}{\hbar}(x-x\,^{\prime}).$
(S20)

Previous literature analyzed DF Minguzzi and Gangardt (2005) making use of the
Stationary Phase Approximation (SPA) Bleistein and Handelsman (1986), even in
scenarios where the dynamics is not scale-invariant Pedri (2005); del Campo
(2008). Here, we review this approach. The phase $\phi(x,x^{\prime})$ (S20)
has stationary points $x^{*}=x^{\prime*}=\frac{p}{m\dot{b}}$ so that after
Taylor expansion

$\displaystyle\phi(x,x^{\prime})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\frac{m\dot{b}b}{\hbar}[(x-x^{*})^{2}-(x^{\prime}-x^{\prime*})^{2}].$
(S21)

We introduce the reference length $x_{0}=\sqrt{\frac{\hbar}{m\omega_{0}}}$ to
perform the integration in dimensionless units. Note that the OBRDM
$\rho_{1}(x,x^{\prime};t)$ scales as $1/x_{0}$, while the momentum
distribution scales as $x_{0}/\hbar$. One can insert the previous expression
(S21) in the momentum distribution (S19). Provided that
$\frac{1}{2}\frac{m\dot{b}b}{\hbar}\gg 1/x^{2}_{0}$ (i.e., $\dot{b}b\gg
2\omega_{0}$), one can make use of the SPA Bleistein and Handelsman (1986) as
proposed in Minguzzi and Gangardt (2005) to analyze DF. For instance, the free
expansion gives $b(t)=\sqrt{1+(\omega_{0}t)^{2}}$,
$\dot{b}=\frac{\omega^{2}_{0}t}{\sqrt{1+(\omega^{2}_{0}t)^{2}}}$, implying
$\dot{b}b=\omega^{2}_{0}t$ and $t_{DF}\gg\frac{2}{\omega_{0}}$, and leading to
the asymptotic momentum distribution

$\displaystyle n(p,t)$ $\displaystyle\approx$
$\displaystyle\frac{b}{2\pi\hbar}\int
dxdx\,^{\prime}\,e^{i\frac{1}{2}\frac{m\dot{b}b}{\hbar}[(x-x^{*})^{2}-(x^{\prime}-x^{\prime*})^{2}]}\rho_{1}\left(x^{*},x^{\prime*},t=0\right),$
(S22) $\displaystyle\approx$
$\displaystyle\frac{b}{2\pi\hbar}x^{2}_{0}\int\frac{dx}{x_{0}}\frac{dx\,^{\prime}}{x_{0}}\,e^{i\frac{1}{2}\frac{\dot{b}b}{\omega_{0}}\frac{[(x-x^{*})^{2}-(x^{\prime}-x^{\prime*})^{2}]}{x^{2}_{0}}}\rho_{1}\left(x^{*},x^{\prime*},t=0\right),$
$\displaystyle=$
$\displaystyle\frac{\omega_{0}}{\dot{b}}\frac{x^{2}_{0}}{\hbar}\rho\left(x^{*},t=0\right)=\frac{1}{m\dot{b}}\rho\left(\frac{p}{m\dot{b}},t=0\right).$

This equation can also be demonstrated without relying on scale-invariance
when the dynamics is described by the single-particle free-propagator (no
interactions) Pedri (2005).

For the harmonic trap, it is possible to go further and relate the density
profile of the TG gas to the momentum distribution of the Fermi gas. This is
possible as the OBRDM of the ideal Fermi gas can be directly expressed in
terms of the single-particle harmonic-oscillator eigenstates

$\displaystyle\rho^{\rm F}_{1}(x,x^{\prime})=\sum_{n=0}^{{\rm
N}-1}\phi_{n}(x)\phi^{*}_{n}(x^{\prime}),$ (S23)

whence it follows that the momentum distribution reads

$\displaystyle n^{\rm F}(p,0)=\sum_{n=0}^{{\rm
N}-1}\left|\tilde{\phi}_{n}(p)\right|^{2},$ (S24)

in terms of the Fourier transform $\tilde{\phi}_{n}(p)$ of the eigenstates
$\phi_{n}(x)$ of the Harmonic Oscillator
$H=\frac{p^{2}}{2m}+\frac{1}{2}m\omega^{2}_{0}x^{2}$. Given that
$\tilde{\phi}_{n}(p)=\frac{(-i)^{n}}{\sqrt{m\omega_{0}}}\phi_{n}\left(\frac{1}{m\omega_{0}}p\right)$,
we note that the density profile evaluated at the stationary points reads

$\displaystyle\rho^{\rm F}(x^{*},t=0)=\sum_{n=0}^{{\rm
N}-1}\left|\phi_{n}\left(\frac{p}{m\dot{b}},0\right)\right|^{2}=m\omega_{0}\
n^{\rm F}\left(\frac{\omega_{0}}{\dot{b}}p,t=0\right).$ (S25)

As the wavefunctions of the TG and Fermi gases only differ by the Bose-Fermi
mapping, the density profile of the TG and ideal Fermi gas are equal at all
times. In particular, they are equal at $t=0$, i.e., $\rho^{\rm
TG}(x,t=0)=\rho^{\rm F}(x,t=0)$. Thus, for the harmonic potential, the
asymptotic momentum distribution of the TG gas is given in terms of that the
Fermi gas, which is the central feature of DF Rigol and Muramatsu (2005);
Minguzzi and Gangardt (2005)

$\displaystyle n^{\rm TG}(p,t)$ $\displaystyle\approx$
$\displaystyle\frac{\omega_{0}}{\dot{b}}n^{\rm
F}\left(\frac{\omega_{0}}{\dot{b}}p,t=0\right).$ (S26)

##  V Dynamical fermionization in phase space: asymptotic marginal
distributions of the Wigner function

The marginals of the Wigner function read

$\displaystyle\rho(x,t)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{+\infty}dp\
W\left(\frac{x}{b},bp-m\dot{b}x,t=0\right)=\int_{-\infty}^{+\infty}\frac{dp}{b}\
W\left(\frac{x}{b},p,t=0\right),$ (S27) $\displaystyle n(p,t)$
$\displaystyle=$ $\displaystyle\int_{-\infty}^{+\infty}dx\
W\left(\frac{x}{b},bp-m\dot{b}x,t=0\right)=\int_{-\infty}^{+\infty}\frac{dk}{m\dot{b}}\
W\left(\frac{p}{m\dot{b}}-\frac{k}{m\dot{b}b},k,t=0\right),$ (S28)

where the last identity makes use of the change of variable $k=bp-m\dot{b}x$.
The first equation leads to

$\rho(x,t)=\frac{1}{b}\rho\left(\frac{x}{b},t=0\right).$ (S29)

We further assume $\dot{b}(t)b(t)$ is large and perform a Taylor expansion on
the integrand in r.h.s. of Eq. (S28), i.e.,

$W\left(\frac{p}{m\dot{b}}-\frac{k}{m\dot{b}b},k,0\right)=W\left(\frac{p}{m\dot{b}},k,0\right)+\sum_{n=1}^{\infty}\frac{1}{n!}\left(\frac{1}{mb\dot{b}}\right)^{n}W_{1}^{(n)}(x,k,\,0)\big{|}_{x=p/(m\dot{b})}k^{n},$
(S30)

we find

$n(p,\,t)=\frac{1}{m\dot{b}}\rho\left(\frac{p}{m\omega_{0}},t=0\right)+\frac{1}{m\dot{b}}\sum_{n=1}^{\infty}\frac{1}{n!}\left(\frac{1}{mb\dot{b}}\right)^{n}\int_{-\infty}^{\infty}dkW_{1}^{(n)}(x,k,\,0)\big{|}_{x=p/(m\dot{b})}k^{n},$
(S31)

where $W_{1}^{(n)}(x,\,y,t)$ denotes the $n$-th partial derivative with
respect to the first argument. Next, we show that when

$\displaystyle\frac{\partial^{p+q}\rho_{1}(x,\,y)}{\partial^{p}x\partial
y^{q}}$ $\displaystyle\to 0,\,\text{as}\,x\to\pm\infty,y\to\pm\infty,\,\forall
p,\,q,$ (S32)

one can find

$\left|\int_{-\infty}^{\infty}dkW_{1}^{(n)}(x,k,\,0)k^{n}\right|<\infty,\,\forall
n\geq 1.$ (S33)

Thus at late times, to the first-order approximation,

$n(p,t)\approx\frac{1}{m\dot{b}}\rho\left(\frac{p}{m\omega_{0}},t=0\right).$
(S34)

To show Eq. (S32) implies Eq. (S33), we first resort to the definition of the
Wigner transformation and find

$\displaystyle\int_{-\infty}^{\infty}dkW_{1}^{(n)}(x,k,\,0)k^{n}$
$\displaystyle=\frac{1}{\pi\hbar}\int_{-\infty}^{\infty}dy\int_{-\infty}^{\infty}dk\frac{\partial^{n}\rho_{1}(x-y,\,x+y)}{\partial
x^{n}}e^{2\text{i}ky/\hbar}k^{n}$
$\displaystyle=\frac{1}{\pi\hbar}\left(\frac{\hbar}{2i}\right)^{n}\int_{-\infty}^{\infty}dy\int_{-\infty}^{\infty}dk\frac{\partial^{n}\rho_{1}(x-y,\,x+y)}{\partial
x^{n}}\frac{\partial^{n}e^{2\text{i}ky/\hbar}}{\partial y^{n}}.$ (S35)

Integration by parts yields

$\int_{-\infty}^{\infty}dkW_{1}^{(n)}(x,k,\,0)k^{n}=\frac{1}{\pi\hbar}\left(\frac{\hbar}{2i}\right)^{n}\int_{-\infty}^{\infty}dy\int_{-\infty}^{\infty}dk\frac{\partial}{\partial
y}\frac{\partial^{n}\rho_{1}(x-y,\,x+y)}{\partial
x^{n}}\frac{\partial^{n-1}e^{2\text{i}ky/\hbar}}{\partial y^{n-1}},$ (S36)

where the boundary term vanishes since

$\lim_{y\to\pm\infty}\frac{\partial^{n}\rho_{1}(x-y,\,x+y)}{\partial
x^{n}}=0,$ (S37)

according to Eq. (S32). Iterating this process for $n-1$ times, we find

$\displaystyle\int_{-\infty}^{\infty}dkW_{1}^{(n)}(x,k,\,0)k^{n}$
$\displaystyle=\frac{1}{\pi\hbar}\left(\frac{\hbar}{2i}\right)^{n}\int_{-\infty}^{\infty}dy\int_{-\infty}^{\infty}dk\frac{\partial^{n}\rho_{1}(x-y,\,x+y)}{\partial
x^{n}\partial y^{n}}e^{2\text{i}ky/\hbar}$
$\displaystyle=\left(\frac{\hbar}{2i}\right)^{n}\int_{-\infty}^{\infty}dy\frac{\partial^{n}\rho_{1}(x-y,\,x+y)}{\partial
x^{n}\partial y^{n}}\delta(y)$
$\displaystyle=\left(\frac{\hbar}{2i}\right)^{n}\frac{\partial^{n}\rho_{1}(x-y,\,x+y)}{\partial
x^{n}\partial y^{n}}\bigg{|}_{y=0},$ (S38)

which immediately implies Eq. (S33).

Finally, we note for an initial eigenstate of the TG gas and the spin-
polarized fermi gas, according to the definition of the OBRDM,

$\rho_{1}(x,y,0)={\rm N}\int dx_{2}\dots dx_{{\rm N}}\
\Psi(x,x_{2},\dots,x_{{\rm N}},0)\Psi^{*}(y,x_{2},\dots,x_{{\rm N}},0).$ (S39)

In the case of harmonic confinement, Eq. (S32) is guaranteed by the property
of the many-body wave function,

$\partial^{n}\Psi(x,x_{2},\dots,x_{{\rm N}},0)/\partial x^{n}\to
0,\,\text{as},\,x\to\pm\infty,\,\forall n.$ (S40)

##  VI Tailoring Dynamical Fermionization and Delta kick cooling of a Tonks-
Girardeau gas in a time-dependent harmonic trap

###  VI.1 Ground-state OBRDM of a TG gas in a harmonic trap

To compute the evolution of the momentum distribution (S19) it is necessary to
determine the initial density matrix explicitly. A general expression for a
pure state was derived by Pezer and Buljan making use of the Bose-Fermi
mapping Pezer and Buljan (2007). Under harmonic confinement, further progress
is possible, recognizing that the Slater determinant is that of a Vandermonde
matrix. Closed formulas of the integrals can be found in terms of gamma and
hypergeometric functions Lapeyre  _et al._ (2002); Forrester  _et al._ (2003).
This expression is useful for the numerical evaluation of the momentum
distribution and we use it in our work. The ground state wavefunction of a
harmonically trapped TG gas can be simply written as

$\displaystyle\Psi_{\rm TG}(x_{1},\dots,x_{N})=\frac{1}{\sqrt{{\rm
N}!}}\left|\det_{n=0,k=1}^{{\rm N}-1,{\rm N}}[\phi_{n}(x_{k})]\right|,$ (S41)

in terms of the well-known eigenfunctions of the single-particle harmonic
oscillator $\phi_{n}(x_{k})$. Recognizing that Slater determinant involved is
that of a Vandermonde matrix, one finds that the ground-state wavefunction
admits the Jastrow form Girardeau  _et al._ (2001) taking $m=\omega_{0}=1$

$\displaystyle\Psi_{\rm TG}(x_{1},\dots,x_{N})=C_{{\rm
N}}\exp\left(-\sum_{j}\frac{x_{j}^{2}}{2x_{0}^{2}}\right)\prod_{j<k}|x_{j}-x_{k}|,$
(S42)

with normalization constant

$\displaystyle C_{{\rm N}}=\frac{2^{\frac{{\rm N}({\rm N}-1)}{4}}}{x_{0}^{{\rm
N}/2}\left(N!\pi^{\frac{{\rm N}}{2}}\prod_{n=0}^{{\rm N}-1}n!\right)^{1/2}}.$
(S43)

As a result, in units of $x_{0}$, the one-body density matrix reads Forrester
_et al._ (2003)

$\displaystyle\rho_{1}(x,y)$ $\displaystyle=$ $\displaystyle{\rm
N}\int_{-\infty}^{+\infty}dx_{1}\cdots\int_{-\infty}^{+\infty}dx_{{\rm
N}-1}\Psi_{{\rm N}}(x_{1},\cdots,x_{{\rm N}-1},x)\Psi_{{\rm
N}}(x_{1},\cdots,x_{{\rm N}-1},y)^{*}$ (S47) $\displaystyle=$
$\displaystyle\frac{{\rm N}\ 2^{{\rm N}-1}}{({\rm
N}-1)!\sqrt{\pi}}e^{-(x^{2}/2)-(y^{2}/2)}\int_{-\infty}^{+\infty}dx_{1}\cdots\int_{-\infty}^{+\infty}dx_{{\rm
N}-1}\prod_{l=1}^{{\rm N}-1}e^{-x^{2}_{l}}|x-x_{l}|$
$\displaystyle\times|y-x_{l}|({\rm det}[\frac{1}{\sqrt{\sqrt{\pi}2^{j-1}}\
(j-1)!}x^{j-1}_{k}]_{j,k=1,\dots,{\rm N}-1})^{2}$ $\displaystyle=$
$\displaystyle\frac{{\rm N}\ 2^{{\rm N}-1}}{({\rm
N}-1)!\sqrt{\pi}}e^{-(x^{2}/2)-(y^{2}/2)}\int_{-\infty}^{+\infty}dx_{1}\cdots\int_{-\infty}^{+\infty}dx_{{\rm
N}-1}\prod_{l=1}^{{\rm N}-1}e^{-x^{2}_{l}}|x-x_{l}|$
$\displaystyle\times|y-x_{l}|\frac{1}{{\rm N}!}\sum_{P\in S_{{\rm
N}-1}}\epsilon(P)\sum_{Q\in S_{{\rm
N}-1}}\epsilon(Q)x^{P(l)-1}_{l}x^{Q(l)-1}_{l}$ $\displaystyle=$
$\displaystyle\frac{{\rm N}\ 2^{{\rm N}-1}}{({\rm
N}-1)!\sqrt{\pi}}e^{-(x^{2}/2)-(y^{2}/2)}\frac{1}{{\rm N}!}\sum_{P\in S_{{\rm
N}-1}}\epsilon(P)\sum_{Q\in S_{{\rm N}-1}}\epsilon(Q)\prod_{l=1}^{{\rm
N}-1}\int_{-\infty}^{+\infty}dx_{l}e^{-x^{2}_{l}}|x-x_{l}|$
$\displaystyle\times|y-x_{l}|\ \frac{1}{\sqrt{\sqrt{\pi}2^{P(l)-1}}\
(P(l)-1)!}\frac{1}{\sqrt{\sqrt{\pi}2^{Q(l)-1}}\ (Q(l)-1)!}\ x^{P(l)-1}_{l}\
x^{Q(l)-1}_{l}$ $\displaystyle=$ $\displaystyle\frac{2^{{\rm
N}-1}}{\Gamma({\rm N})\sqrt{\pi}}e^{-(x^{2}/2)-(y^{2}/2)}{\rm
det}\left[\frac{2^{(j+k)/2}}{2\sqrt{\pi}\sqrt{\Gamma(j)\Gamma(k)}}b_{j,k}(x,y)\right]_{j,k=1,\cdots,{\rm
N}-1}\,,$ (S48)

with

$\displaystyle b_{j,k}(x,y)=\int_{-\infty}^{+\infty}dt\ e^{-t^{2}}|x-t|\
|y-t|\ t^{j+k-2}.$ (S49)

Note that this provides the one-body reduced density matrix at $t=0$. Its
subsequent time evolution for $t>0$ can be described by making use of this
expression and scale invariance. The corresponding momentum distribution can
be readily obtained via Fourier transformation.

###  VI.2 Delta kick cooling of an expanding TG gas

Figure 5: DKC of an expanding TG gas. (a) Momentum distribution $n(p)/[{\rm
N}b(t)]$ where ${\rm N}$ is the number of particles and $b(t)$ the scaling
factor solution to the Ermakov equation. The initially peaked momentum
distribution broadens with the expansion as a result of DF, approaching that
of an ideal Fermi gas, at $2/\omega_{0}$. A general kick, such as that with
kick strength $0.7\tau_{k}$, further distorts the momentum distribution. Only
when the kick duration is chosen according to the DKC relation, DF is
completely reversed, and the momentum distribution of the final state matches
that of an adiabatic expansion and is given by that at $t=0$, up to scaling
(${\rm N}=7$).

We consider a cloud of atoms that is initially trapped in a harmonic potential
of frequency $\omega_{0}$. At time $t=0$ the trap is released so that
$\omega(t>0)=0$ and the solution of the Ermakov equation is given by
$b(t)=\sqrt{1+(\omega_{0}t)^{2}}$. Fig. 5 depicts the momentum distribution
for the TG gas trapped in a harmonic potential at initial time $t=0$, its
evolution after a free expansion governed by DF, and after a kick of duration
$\tau_{k}$, reversing DF. One can notice that the kick allows recovering a
distribution of bosonic shape. However, it is interesting to notice that the
distribution is rescaled compared to the initial momentum distribution,
following equation (S19), the expansion followed by the kick acting as a
microscope. This rescaling matches with the expression found in Minguzzi and
Gangardt (2005) predicting the evolution of the distribution for large $p$ as
$n(p,t)\approx p^{-4}b^{-3}$ that matches with the prediction of a decaying
tail upon completion of the DKC protocol $n_{\rm kick}(p,t)=b\times
n(pb,0)\approx b\times(pb)^{-4}$.

###  VI.3 Delta kick cooling of an imploding TG gas

Intuitively the DF occurs for a trap expansion when the particles get more
distant from each other and do not interact. We next show that it is also
possible to observe the DF in the less trivial case of a compression process
with $b(t_{F})<b(t=0)=1$ if the process is short enough so that the product
$\dot{b}b$ is large. The implosion can be engineered by a modulation of the
trapping potential. We proceed by reverse engineering. We first set a desired
time evolution of the scaling factor. For the sake of illustration, we
consider

$\displaystyle b(t)=1+(b_{F}-1)\left(\frac{t}{t_{F}}\right)^{3}.$ (S50)

Making use of the Ermakov equation, the time-dependent frequency
$\omega^{2}(t)$ that leads to thus time-dependence of $b(t)$ can be
determined, as shown in Fig. 7. Large values of $\dot{b}b$ induce a high-
frequency phase modulation in the coordinate representation.

Figure 6:  (a) The frequency modulation leading to the implosion protocol
(S50), involving a smooth quench at $t=0$ and a subsequently smooth
modulation. DF is reversed in an implosion making use of a pulsed inverted
oscillator. We chose the frequency of the kick
$\omega^{2}_{k}=1000\omega^{2}_{0}$, the compression time protocol as
$t_{k}=0.3/\omega_{0}$ and the final frequency determined by
$b_{F}=\sqrt{\omega_{0}/\omega_{F}}=\sqrt{1/3}$. (b) Strength of the phase
modulation shown as a function of the frequency ratio.

One may wonder whether it is possible to cancel DF by DKC. However, now the
derivative of the scaling parameter $\dot{b}$ is negative during the process
as shown in Fig. 7. Also, to cancel the DF it is necessary to apply a kick
with an inverted potential. One can now consider the “kicked” Hamiltonian
where the kick-pulse is repulsive $\omega_{k}=i\omega_{I}$, i.e.,
$H_{k}(t)=H(t)-\delta(t-t_{k})\frac{1}{2}m\omega_{I}^{2}\sum_{i=1}^{{\rm
N}}x_{i}\,^{2}$. The corresponding time-evolution operator reads as
$U_{\delta}(t_{k},0)=\exp\left(+i\frac{\tau_{k}m\omega_{I}^{2}}{2\hbar}\sum_{i=1}^{{\rm
N}}x_{i}\,^{2}\right)U(t_{k},0)$, where $U(t_{k},0)$ is the propagator
associated with $H(t)$, and $\tau_{k}$ is a small time scale during which the
kick is applied. As a result, the relation for the pulse parameters to reverse
DF is still given by Eq. (9). For the prescribed trajectory (S50), applying
the kick at the time $t_{k}$ and duration $\tau_{k}$, the required pulse
parameters to reverse DF are set by

$\displaystyle\tau_{k}\omega^{2}_{k}=3\frac{(b_{F}-1)}{t_{k}b_{F}}.$ (S51)

Figure 7: DKC of an imploding TG gas. (a) Momentum distribution $n(p)$ of a
trapped gas, after the implosion at $t_{k}=0.5/\omega_{0}$ and after the
delta-kick for $0.7\tau_{k}$ (green dashed-dotted line), after the delta-kick
$\tau_{k}$ (turquoise blue dashed line); ${\rm N}=7$. In dashed green the kick
strength is chosen as $\tau_{k}\omega^{2}_{k}$ in order to cancel DF, with
$\omega^{2}_{k}=1000\omega^{2}_{0}$. In choosing the strength of the kick
differently one can modulate DF. The implosion protocol is described by Eq.
(S50) with $b_{F}=\sqrt{1/5}$.


# Blind signal decomposition of various word embeddings based on join and
individual variance explained (JIVE)

Yikai Wang
Department of Biostatistics&Bioinformatics
Emory University
<EMAIL_ADDRESS>
Weijian Li
Department of Computer Science
Emory University
<EMAIL_ADDRESS>

###### Abstract

In recent years, natural language processing (NLP) has become one of the most
important areas with various applications in human’s life. As the most
fundamental task, the field of word embedding still requires more attention
and research. Currently, existing works about word embedding are focusing on
proposing novel embedding algorithms and dimension reduction techniques on
well-trained word embeddings. In this paper, we propose to use a novel joint
signal separation method - JIVE to jointly decompose various trained word
embeddings into joint and individual components. Through this decomposition
framework, we can easily investigate the similarity and difference among
different word embeddings. We conducted extensive empirical study on word2vec,
FastText and GLoVE trained on different corpus and with different dimensions.
We compared the performance of different decomposed components based on
sentiment analysis on Twitter and Stanford sentiment treebank. We found that
by mapping different word embeddings into the joint component, sentiment
performance can be greatly improved for the original word embeddings with
lower performance. Moreover, we found that by concatenating different
components together, the same model can achieve better performance. These
findings provide great insights into the word embeddings and our work offer a
new of generating word embeddings by fusing.

##  1 Introduction

Nowadays, with the advance of technology and information system, natural
language processing (NLP) has become one of the most important areas in both
industry and scientific research (Turney and Pantel, 2010; Cambria and White,
2014). Applications of NLP have been involved in many different areas, such as
language translation, keyword searching, machine comprehension, sentiment
analysis, user language interface, cross-language information retrieval,
speech recognition and other language related text mining tasks (Kim, 2014;
Salton and Buckley, 1988; Yin et al., 2017; Xu et al., 2015; Mikolov et al.,
2013; Lai et al., 2015; Dos Santos and Gatti, 2014). The core task of NLP is
to make computer system understand and manipulate natural human’s language to
perform various tasks. However, this is a very challenging task.

First of all, the fundamental question of any kind of NLP tasks is how to deal
with a large number of words (i.e. the smallest element in one sentence or
passage) in the vocabulary list (Mikolov et al., 2013). For example, Oxford
English Dictionary, published in 1989, contains full entries for 171,476 words
in current use (Murray et al., 1933). Initially, research built language
processing system based on N-gram language model or bag-of-word concept using
word counts or frequency, i.e. tf-idf (Aizawa, 2003; Cambria and White, 2014).
Obviously, this type of approach cannot fully capture the rich information of
human natural language and has limited performance. So far, one the most
fundamental and successful achievement in NLP is the concept of word
embedding, which aims at learning a vector space representation of each single
word in low dimensional space. Some famous and well-performed word embedding
algorithms includes Word2Vec (Mikolov et al., 2013), FastText (Bojanowski et
al., 2017) and GLoVE (Pennington et al., 2014). Each of these methods try to
learn a relatively low-dimensional vector representation of each word based on
the co-occurace pattern of words in the corpus in an unsupervised manner. For
example, Word2Vec trains the embedding of the word by using its surrounding
window of context words as predictor or outcomes. Therefore, the nearby words
tend to have similar word embeddings. Moreover, as for FastText, it further
makes use of the information in the character level (sub-element consists of
word), and GLoVE learns the embedding based on decomposing the global co-
occurrence matrix of words. Nowadays, these embedding algorithms are serving
as the foundation for NLP, and built upon these algorithm, a growing number of
machine learning and deep learning models are being developed in recent years
and achieve great success in many areas.

However, although word embedding algorithm facilitates the model development
of NLP, limited work or research has been conducted to analyze or understand
the output of these embedding algorithms. Most work simply treats them as
’hyper-parameter’ and changes the embedding algorithm based on their
performance on available datasets. However, given the importance of the word
embedding in NLP, we argue that more efforts are needed for investigating the
difference and the meaning of the output of word embedding algorithm.
Unfortunately, given the dimensionality, evaluating the difference or
similarity of various word embeddings is very challenging. Existing works are
only focusing on conducting dimension reduction on one word embedding output
instead of understanding the common and individual part of various word
embeddings (Levy and Goldberg, 2014; Iacobacci et al., 2015; Yin and Shen,
2018). For example, Yin and Shen (2018) provide a theoretical understanding of
word embedding to guide the selection the dimension. Iacobacci et al. (2015)
aims to further transfer the estimated embeddings to capture more information.
Different from these existing works, we consider multiple word embeddings and
utilize joint signal decomposition approach to map them into joint and
individual components for following analysis.

In this paper, we propose to jointly study the word embeddings learned from
different algorithms by decomposing them into common and individual components
based on the joint signal separation method - JIVE (Lock et al., 2013;
O’Connell and Lock, 2016; Feng et al., 2018). Specifically, JIVE can decompose
different word embeddings into the summation of shared column space and the
orthogonal embedding-specific space. Based on this outcome, we can answer
several fundamental questions: 1) what is the similarity among different word
embedding algorithms; 2) what is the difference between two word embedding
algorithms. In summary, we study 3 important questions for word embedding
based on JIVE analysis on word embeddings: 1) whether the information captured
by word embedding increases linearly with the dimension of embedding space; 2)
what is the difference and similarity for word embeddings trained with the
same algorithm but different corpus; 3) what is the difference and whether
there is similarity between different word embedding algorithms. Moreover,
there are other relevant approaches in the literature for joint decomposition
including canonical correlation analysis (CCA) (Hotelling, 1992), partial
least squares (PLS) (Helland, 2014), multi-level functional PCA (Di et al.,
2009), parallel and joint independent component analysis (Eichele et al.,
2008; Moosmann et al., 2008; Wang and Guo, 2019; Lukemire et al., 2018; Mejia
et al., 2020; Wang and Guo, 2020). These methods cannot provide joint and
individual information like JIVE and hence is not what we need. Furthermore,
our work can provide great insights for better understanding the word
embeddings and can potentially give directions for generating better word
embeddings or guidances on how to compress/control the information in word
embeddings.

This paper is organized as follows: section 2 introduces the word embedding
algorithm being considered and the model specification of JIVE. Section 3
discusses the design of our experiments and section 4 presents the results as
well as the discussion. Section 5 discusses the finding in this paper and make
conclusion, and section 6 talks about limitation and the potential future
direction.

##  2 Methodology

In this section, we will first introduce three commonly used word embedding
algorithms and then formally define the model specification, hyper-parameter
selection and algorithm for JIVE.

###  2.1 Word Embedding

We will briefly introduce three commonly used word embedding algorithms being
considered in this paper: Word2Vec (Mikolov et al., 2013), fastText
(Bojanowski et al., 2017) and GloVe (Pennington et al., 2014).

Word2Vec is one of the first word embeddings that is applied in NLP tasks. It
is a predictive model trained on the context of a word in a large corpus.
During the training, either CBOW or skip-gram can be used to model the
objective function. CBOW uses the context to predict he center word, whereas
skip-gram uses the center word to predict he context. The objective function
aims to minimize the loss of this prediction given the word vectors. Either
negative sampling or hierarchical softmax can be used to reduce the
computational complexity.

FastText is trained in a very similar way as Word2Vec, except that fastText is
trained with character level n-grams instead of word level. With this
formation, the model can better capture sub-word information like prefix and
suffix, which also contains semantic meanings. Because of this nature,
fastText is good at handling unseen words and misspelling.

GloVe is a count-based model. It first construct a large co-occurrence matrix
based on the frequency of each word in its context. Matrix factorization then
is used to reduce the dimensionality of this large matrix to gain a vector
representation for each word. The objective function aims to minimize the
reconstruction loss based on this dimension-reduced matrix.

As we discussed, these word embedding methods all achieved great success in
different applications but little is known for the difference and similarity
among them. In the next section, we will introduce a novel statistical signal
decomposition technique which can jointly decompose different word embeddings
into joint and individual components for further analysis.

###  2.2 Joint and Individual Variance Explained

Joint and Individual Variance Explained (JIVE) is a statistical approach for
decomposing multiple data modalities into joint component and modality-
specific components. JIVE is originally proposed by Lock et al. (2013) for
decomposing different features on the same subject (gene and RNA). In this
section, we will provide the formal model specification of JIVE and its
assumption.

First, denote $X_{i}\in\mathcal{R}^{p_{i}\times n}$ to be the $i$th dataset
from same $n$ subjects with different feature spaces, $i=1,...,K$. For
example, in Lock et al. (2013), $X_{1}$ denotes the gene expression data same
and $X_{2}$ denotes miRNA data. In our case, $X_{i}$’s will denote different
word embeddings with $n$ to be the number of distinct words and $p_{i}$
denotes the dimension for word embedding. The unified model of JIVE will be:

$X_{i}=J_{i}+A_{i}+E_{i},$ (1)

Figure 1: Illustration of JIVE analysis for jointly decomposing word
embeddings

where $J_{i}$ captures the joint information among all data modalities and
$A_{i}$ denotes the modality-specific information, and $E_{i}$ denotes
residuals with expectation of zero. To ensure the identifiability of JIVE, we
have the following constraints: $J_{i}A_{i}^{\prime}=0$ for $i=1,...,K$. The
overal summary of JIVE is summarized in Figure 1.

Specifically, all joint components $J_{i}$ shares the same row space and can
be represented as $J_{i}=B_{i}J$ with $J\in\mathcal{R}^{r\times n}$ to be the
basis of the same row space and $B_{i}$ denotes the linear transformation and
is assumed to satisfy that $B_{i}^{\prime}B_{i}=I_{p_{i}}$. For individiual
component, we can also re-write it as $A_{i}=D_{i}H_{i}$ with
$D_{i}^{\prime}D_{i}=I_{p_{i}}$ and $JH_{i}^{\prime}=0$. Therefore, JIVE is a
natural extension of PCA for multiple modalities where $J$ and $H_{i}$ can be
viewed as the principle components for $X_{i}$. In summary, we have

$X_{i}=B_{i}J+D_{i}H_{i}+E_{i}.$ (2)

###  2.3 JIVE on Word Embeddings

Section 2.2 introduce the definion of JIVE. In this section, we will present
more details about JIVE analysis on word embeddings.

####  2.3.1 Joint Rank Selection

First, the hyper-parameter in JIVE is the rank of joint component $J_{i}$ or
the dimension of $J$, which is $r$. It controls the amount of information
shared acorss different word embeddings and its value can significantly affect
the conclusion. In the literature, there are three commonly used approaches
for selection of $r$: permutation test (Lock et al., 2013), BIC based
selection (O’Connell and Lock, 2016; Wang and Guo, 2020) and a theoretical
selection approach based on Wedin bound (Feng et al., 2018). Although
permutation based approach tends to be the most robust method, its application
for word embedding analysis is not desired because the computation burden is
very expensive because of the sample size, number of word in dictionary.
Moreover, BIC method requires the specification of probability distribution on
word embeddings. However, different embedding method has different
distributional assumption and hence the BIC based method can be very
complicated. Therefore, we select to use Wedin bound based selection approach
for JIVE analysis of word embeddings.

####  2.3.2 Algorithm for JIVE

Second, in this section, we present the algorithm for estimating the unknown
parameters in JIVE. Given the limited resources of JIVE in python comunity, we
followed the algorithm proposed in Lock et al. (2013) and implemented JIVE in
Python 3.0. We proposed an efficient approximate estimation method used for
initial starting point. The method for 2 word embeddings for estimating the
initial value and the iterative algorithm for JIVE are summarized in Algorithm
1. It can be easily extended to the case with more than 2 word embeddings.

Algorithm 1 JIVE with two word embeddings

Inputs: Concatenated embedding matrix
$\bm{X}=[\bm{X}_{1}^{\prime},\bm{X}_{2}^{\prime}]^{\prime}$ of dimension
$p_{1}+p_{2}\times N$, rank of joint component $r$, rank of individual terms
$r_{1},r_{2}$.

Initial values:

1\. Conduct SVD on $X=U$ diag
$(\Lambda_{1},...,\Lambda_{p_{1}+p_{2}})V^{\prime}$, set $J^{(0)}=U$ diag
$(\Lambda_{1},...,\Lambda_{r},0,..,0)V^{\prime}$.

2\. Conduct 2 separate SVD on first $p_{1}$ and last $p_{2}$ rows of
$X-J^{(0)}$, and only takes first $r_{i}$ terms for $A_{i}$,
$A^{(0)}=[A_{1}^{\prime},A_{2}^{\prime}]^{\prime}$.

3\. Calculate $R^{(0)}=\|X-J^{(0)}-A^{(0)}\|_{2}^{2}$.

Estimation

repeat

1\. Condition on $A^{(t)}$, update $J$ of rank $r$ to minimize
$\|X-J-A^{(t)}\|_{2}^{2}$.

2\. Condition on $J^{(t+1)}$, update $A_{1}$ of rank $r_{1}$ and $A_{2}$ of
rank $r_{2}$ to minimize $\|X-J^{(t+1)}-A\|_{2}^{2}$

3\. Calculate $R^{(t+1)}=\|X-J^{(t+1)}-A^{(t+1)}\|_{2}^{2}$

until convergence, i.e. $\frac{R^{(t+1)}-R^{(t)}}{R^{(t)}}<\epsilon$

The word embeddings $X_{i}$ in Algorithm 1 is assumed to be de-meaned at each
row/feature and re-scaled to have a standard deviation of 1 for the whole
embedding, i.e. $X_{i}=X_{i}/\|X_{i}\|$ where where $\|\cdot\|$ defines the
Frobenius norm.

The initial estimation follows a 2-step procedure which first conduct a single
PCA on concatenated data matrix to extract the joint component. In the second
step, it removes the estimated joint component from the data matrix and
conductes 2 separated PCA for each modalities. The esimation algorithm for
JIVE is adopted from Lock et al. (2013). We will check the algorithm
convergence in section 4.

##  3 Experimental Design

In this section, we will introduce the details about our experimental design.
We first introduce the datasets used in this paper and then the 3 different
downstream tasks being considered (summarized in Table 1). Finally, we will
discuss about the motivation and assumptions for each experiment settings.

Table 1. A brief summary of experimental design

###  3.1 Downstream Task and Datasets

For the downstream task, we wanted to choose a task of which the results and
performance depends largely on the quality of the word embeddings itself,
rather than on the structure of the neural network model. Our reasoning behind
this was that if the neural network itself is too complicated, we will not be
able to tell the performance differences between the word embeddings, because
a complicated model might compensate the negative effects caused by bad
embeddings. To this end, we chose sentiment analysis as our downstream task,
because the prediction can be affected by the meaning of a few keywords, in
which case the quality of word embeddings is decisive.

We chose two datasets for the sentiment analysis task. One is Stanford
Sentiment Treebank (Socher et al., 2013), the other is Twitter sentiment
analysis dataset. SST dataset contains movie reviews crawled from Rotten
Tomato. The language is comparatively more formal than Twitter. The training
set for SST has 8539 sentences, and testing set has 1101 sentences. The
vocabulary size is 18193. The target has 5 categories from 0 to 4, ranging
from most negative to most positive. The Twitter dataset contains 50248 tweets
crawled from Twitter, of which half were used for training and half for
testing. The vocabulary size for Twitter is 24604. Different from SST, the
Twitter dataset is a binary classification problem with only 2 sentiments,
positive and negative.

###  3.2 Same Embedding with Different Dimensionalities

We explored the effects brought by the difference in dimensionalities of the
same embedding model. Intuitively, the model with higher dimensionality should
encode more information than the model with lower dimensionality. Besides, the
model with higher dimensionality should encode all the information that is
already encoded by the model with lower dimensionality. With the joint and
individual matrices generated by JIVE, we were able to test our intuition.

For this experiment, we first ran JIVE on 2 word embeddings X and Y which have
different dimensionalities. The generated joint matrix represented the common
information encoded by the two embeddings, and had the shape of vocabulary
size by joint rank. The 2 individual matrices represented the distinctive
information that is not shared between the two embeddings. And this part of
information should be the potential cause of the different performances showed
by the two embeddings. With the joint and individual matrices, we were able to
directly use them as word embeddings to conduct downstream tasks. We also
concatenated 1) Joint matrix with individual matrix of X; 2) Joint matrix with
individual matrix of Y; 3) Individual matrix of X and that of Y and 4) Joint
matrix and the 2 individual matrices. This gave us 4 more word embeddings. The
performance of each word embedding gave us insights on what information are
encoded by the original embeddings.

Specifically, we conducted experiments on the following groups of embeddings
on the SST dataset: fastText 50d with fastText 100d, fastText 50d with
fastText 200d, fastText 100d with fastText 200d and GloVe 50d with GloVe 200d.

###  3.3 Same Embedding Trained on Different Corpus

We also explored the effects brought by the different training corpus of the
same embedding algorithms. The intuition was that the word embedding trained
on a formal language corpus, e.g. Wikipedia, might result in better
performance on formal language downstream tasks. On the other hand, the word
embedding trained on a colloquial corpus,e.g. Twitter, could perform better on
colloquial language downstream tasks. Even though the word embedding model was
generated by the same algorithm, the information they encoded could differ
drastically. Meanwhile, we wanted to explore the JIVE’s effects on these
embeddings. Similar to the precious setting, we will use the joint and the 2
individual matrices to generate 4 more new embeddings. Finally, we will use
this embeddings to conduct 2 downstream tasks: SST and Twitter sentiment
analysis. The reason behind this was that SST contains sentences that were
more formal than Twitter. The differences between the embeddings’ performances
can provide insights into these embeddings.

For this setting, we used the following embeddings: Word2Vec 400d trained on
400 million tweets, Word2Vec 300d trained on 6 billion tokens from Google
News, GloVe 200d trained on 2 billion tweets and GloVe 200d trained on 6
billion words from Wikipedia and Gigaword5. Both of the Word2Vec embeddings
were trained using skip-gram and negative sampling. We used the generated
embeddings to conduct experiments on both SST and Twitter sentiment analysis.

###  3.4 Different Embeddings Trained on the Same Corpus

Last but not least, we explored the different embeddings trained on the same
corpus. Since different embedding algorithms would encode the information
quite differently, JIVE analysis on those embeddings can provide us some
insights on what are the shared and individual information. We used Word2Vec
400d and GloVe 200d both trained on Twitter to conduct experiment on Twitter
sentiment analysis.

###  3.5 Neural Network Model

We used the sipmle CNN model proposed by Kim (2014). This model used CNN
kernels of size $N$ by $L$, where $N$ was the kernel size and $L$ was the
sentence length to simulate n-grams and sequential information of a sentence.
In our case, we first padded all sentences to the max length and applied a
kernel of size 5 to the sentences.

##  4 Empirical Results

In this section, we will present the results corresponding to section 3. First
of all, all JIVE iterations are converged.

Table 2. Variance explained percentage for the same embedding algorithm with
different dimensions

###  4.1 Setting 1: Same embedding different dimension - large dimension
always captures

In this section, we present the results of using same embedding method
(FastText and GLoVE) with different number of dimension. First of all, among
all settings discussed in section 3.2, the embeddings with higher dimension
can always capture all the information in lower dimensional embedding based on
JIVE as shown in Table 2, which is as expected and also justifies that method
works properly. However, the among of information stored in embeddings is not
linearly increased with the number of dimension. For example, the 50
dimensional FastText takes 84.6% information in 100 dimensional FastText.
Table 5 summarizes the performance on SST-based sentiment analysis for
different components. The joint component of rank 50 for GLoVE has a better
performance than the original GLoVE with 50 features. Moreover, for FastText,
the concatenated data between FastText50 and FastText200 achieved a higher
performance than the original FastText200. These findings suggest that JIVE
could help improve the performance under some situations.

Table 3. Variance explained percentage for the same embedding algorithm
trained with different corpus Table 4. Variance explained percentage for the
GLoVE and Word2Vec trained with on the same Twitter corpus Table 5. Accuracy
on SST sentiment analysis for the same embedding and different dimensions
Table 6. Accuracy on Twitter sentiment analysis for the GLoVE and Word2Vec
trained on different corpus

###  4.2 Setting2: Same embedding trained on different corpus

Table 3 summarizes the percentage of information captured by different
components from word2vec and GLoVE trained on the different corpus. Although
trained on different corpus, the information in joint component is higher than
70% in all situations which indicates a strong similarity. But different from
setting 1, all trained embedding methods still have a considerable amount of
information in their individual components. For example, word2vec trained on
Twitter corpus still keeps a 22.0% information in the individual component.
This indicates that different corpus did contribute to the final estimated
embeddings significantly.

For sentiment analysis, Table 6 summarizes the results on Twitter sentiment
analysis for GLoVe trained on Twitter or Wiki corpus and for Word2Vec trained
on Twitter and Google News corpus. First of all, for the original trained
embedding (before JIVE), embeddings trained on Twitter corpus achieved a
higher accuracy for Twitter sentiment analysis which is as expected. The joint
component for word2vec trained on Twitter and Google News has a better
accuracy than original performance of word2vec on Google News corpus. Although
the joint part has slightly lower accuracy than original for word2vec on
Twitter, the concatenated data (Joint+Ind) for word2vec-Twitter achieved a
much higher accuracy than its original one (by 5%). Similar for GLoVE on
Twitter and Wiki, the (Joint+Ind) for GLoVE-Twitter has a better accuracy than
its original performance. In either cases, the final concatenated (Ind X+Y+
Joint) can always provide much higher sentiment classification accuracy, which
indicates the possibility of generating better word embedding by JIVE-based
concatenation. We also note that Wiki corpus has relatively lower quality
compared with Google News. This might be reason that the joint component for
GLoVE-Wiki and GLove-Twitter has worse accuracy than their orginal ones.

Moreover, we also included the results on SST sentiment analysis for word2vec
trained on Twitter and Google News corpus in Table 7 which serves as a further
sanity check since both embeddings are trained from the different corpus than
the corpus used for testing. We find similar results where the joint component
improves the original ones with lower accuracy.

Table 7. Accuracy on SST sentiment analysis for the Word2Vec trained on
Twitter and Google News corpus Table 8. Accuracy on Twitter sentiment analysis
for the GLoVE and Word2Vec trained on Twitter corpus

###  4.3 Setting3: Different embeddings trained on the same corpus

Table 4 summarizes the percentage of information captured by different
components from word2vec and GLoVE trained on the same Twitter corpus. The
joint rank of GLoVE and Word2Vec trained on Twitter data is 172, and although
trained by different algorithm, the joint component captures 95.3% information
for GLoVE. This finding indicates that corpus is very important in training
word embedding method. Table 8 summarizes the performance on Twitter sentiment
analysis for different components. Again, the joint component has higher
accuracy than the original lower one - GLoVE-Twitter, and the JIVE-based
concatenated embeddings always have better performance than their original
ones, which further justifies the usage of JIVE in analyzing different word
embeddings.

##  5 Discussion and Conclusion

In this paper, we analyzed the trained word embeddings from fastText, word2vec
and Glove with different dimensionalities and corpus based on JIVE. With the
help of JIVE method, we decompose different word embedding results into joint
and individual components. Based on sentiment analysis on 2 public available
datasets, we evaluate the information contained in different part of the word
embedding. We found that for the same embedding algorithm trained with the
same corpus, higher dimensional embeddings can always capture all the
information in the lower dimensional embeddings but the amount of information
is not increased linearly with dimension. This finding justifies the mechanism
of word embedding algorithm and provides hints for selection of dimensionality
by checking the amount of information captured in different cases.
Furthermore, we find that corpus and embedding algorithm all contribute to the
final trained embeddings significantly where the individual component always
capture considerable amount of information. However, in many situations, by
mapping 2 word embeddings into the same space, the joint component can improve
the original one with lower performance even with a smaller dimension in the
Joint component. This finding is very interesting and provides a possibility
for improve the performance of a word embedding. Finally, based on Twitter
sentiment analysis, JIVE-based concatenated word embeddings always have a
better performance than their original ones. This suggests 2 things: 1) JIVE
can find important features and remove noises in word embedding; 2) combining
different word embeddings could generate a better and more informative word
embeddings. Overall, our work provides great insights for understanding the
similarity and difference of different word embedding algorithms from a
numerical perspective, and provides many opportunities for researchers in NLP
field to further study and improve the word embeddings.

##  6 Limitation and Future Work

Currently, we only evaluated the performance based on sentiment analysis. In
the future, we plan to consider other types of NLP tasks such as sequence
tagging (NER) to better evaluate the estimated joint and individual components
from JIVE. Another limitation of the current work is on the embeddings.

## References

  * Aizawa (2003) Akiko Aizawa. 2003.  An information-theoretic perspective of tf–idf measures.  _Information Processing & Management_, 39(1):45–65.
  * Bojanowski et al. (2017) Piotr Bojanowski, Edouard Grave, Armand Joulin, and Tomas Mikolov. 2017.  Enriching word vectors with subword information.  _Transactions of the Association for Computational Linguistics_ , 5:135–146.
  * Cambria and White (2014) Erik Cambria and Bebo White. 2014.  Jumping nlp curves: A review of natural language processing research.  _IEEE Computational intelligence magazine_ , 9(2):48–57.
  * Di et al. (2009) Chong-Zhi Di, Ciprian M Crainiceanu, Brian S Caffo, and Naresh M Punjabi. 2009.  Multilevel functional principal component analysis.  _The annals of applied statistics_ , 3(1):458.
  * Dos Santos and Gatti (2014) Cicero Dos Santos and Maira Gatti. 2014.  Deep convolutional neural networks for sentiment analysis of short texts.  In _Proceedings of COLING 2014, the 25th International Conference on Computational Linguistics: Technical Papers_ , pages 69–78.
  * Eichele et al. (2008) Tom Eichele, Vince D Calhoun, Matthias Moosmann, Karsten Specht, Marijtje LA Jongsma, Rodrigo Quian Quiroga, Helge Nordby, and Kenneth Hugdahl. 2008.  Unmixing concurrent eeg-fmri with parallel independent component analysis.  _International Journal of Psychophysiology_ , 67(3):222–234.
  * Feng et al. (2018) Qing Feng, Meilei Jiang, Jan Hannig, and JS Marron. 2018.  Angle-based joint and individual variation explained.  _Journal of multivariate analysis_ , 166:241–265.
  * Helland (2014) Inge Helland. 2014.  Partial least squares regression.  _Wiley StatsRef: Statistics Reference Online_.
  * Hotelling (1992) Harold Hotelling. 1992.  Relations between two sets of variates.  In _Breakthroughs in statistics_ , pages 162–190. Springer.
  * Iacobacci et al. (2015) Ignacio Iacobacci, Mohammad Taher Pilehvar, and Roberto Navigli. 2015.  Sensembed: Learning sense embeddings for word and relational similarity.  In _Proceedings of the 53rd Annual Meeting of the Association for Computational Linguistics and the 7th International Joint Conference on Natural Language Processing (Volume 1: Long Papers)_ , volume 1, pages 95–105.
  * Kim (2014) Yoon Kim. 2014.  Convolutional neural networks for sentence classification.  _arXiv preprint arXiv:1408.5882_.
  * Lai et al. (2015) Siwei Lai, Liheng Xu, Kang Liu, and Jun Zhao. 2015.  Recurrent convolutional neural networks for text classification.  In _Twenty-ninth AAAI conference on artificial intelligence_.
  * Levy and Goldberg (2014) Omer Levy and Yoav Goldberg. 2014.  Neural word embedding as implicit matrix factorization.  In _Advances in neural information processing systems_ , pages 2177–2185.
  * Lock et al. (2013) Eric F Lock, Katherine A Hoadley, James Stephen Marron, and Andrew B Nobel. 2013\.  Joint and individual variation explained (jive) for integrated analysis of multiple data types.  _The annals of applied statistics_ , 7(1):523.
  * Lukemire et al. (2018) Joshua Lukemire, Yikai Wang, Amit Verma, and Ying Guo. 2018.  Hint: A toolbox for hierarchical modeling of neuroimaging data.  _arXiv preprint arXiv:1803.07587_.
  * Mejia et al. (2020) Amanda F Mejia, Mary Beth Nebel, Yikai Wang, Brian S Caffo, and Ying Guo. 2020.  Template independent component analysis: Targeted and reliable estimation of subject-level brain networks using big data population priors.  _Journal of the American Statistical Association_ , 115(531):1151–1177.
  * Mikolov et al. (2013) Tomas Mikolov, Ilya Sutskever, Kai Chen, Greg S Corrado, and Jeff Dean. 2013.  Distributed representations of words and phrases and their compositionality.  In _Advances in neural information processing systems_ , pages 3111–3119.
  * Moosmann et al. (2008) Matthias Moosmann, Tom Eichele, Helge Nordby, Kenneth Hugdahl, and Vince D Calhoun. 2008.  Joint independent component analysis for simultaneous eeg–fmri: principle and simulation.  _International Journal of Psychophysiology_ , 67(3):212–221.
  * Murray et al. (1933) James AH Murray, Henry Bradley, William A Craigie, and Charles T Onions. 1933.  _The Oxford english dictionary_ , volume 1.  Clarendon Press Oxford.
  * O’Connell and Lock (2016) Michael J O’Connell and Eric F Lock. 2016.  R. jive for exploration of multi-source molecular data.  _Bioinformatics_ , 32(18):2877–2879.
  * Pennington et al. (2014) Jeffrey Pennington, Richard Socher, and Christopher Manning. 2014.  Glove: Global vectors for word representation.  In _Proceedings of the 2014 conference on empirical methods in natural language processing (EMNLP)_ , pages 1532–1543.
  * Salton and Buckley (1988) Gerard Salton and Christopher Buckley. 1988.  Term-weighting approaches in automatic text retrieval.  _Information processing & management_, 24(5):513–523.
  * Socher et al. (2013) Richard Socher, Alex Perelygin, Jean Wu, Jason Chuang, Christopher D Manning, Andrew Ng, and Christopher Potts. 2013.  Recursive deep models for semantic compositionality over a sentiment treebank.  In _Proceedings of the 2013 conference on empirical methods in natural language processing_ , pages 1631–1642.
  * Turney and Pantel (2010) Peter D Turney and Patrick Pantel. 2010.  From frequency to meaning: Vector space models of semantics.  _Journal of artificial intelligence research_ , 37:141–188.
  * Wang and Guo (2019) Yikai Wang and Ying Guo. 2019.  A hierarchical independent component analysis model for longitudinal neuroimaging studies.  _NeuroImage_ , 189:380–400.
  * Wang and Guo (2020) Yikai Wang and Ying Guo. 2020.  Locus: A novel decomposition method for brain network connectivity matrices using low-rank structure with uniform sparsity.  _arXiv preprint arXiv:2008.08915_.
  * Xu et al. (2015) Kelvin Xu, Jimmy Ba, Ryan Kiros, Kyunghyun Cho, Aaron Courville, Ruslan Salakhudinov, Rich Zemel, and Yoshua Bengio. 2015.  Show, attend and tell: Neural image caption generation with visual attention.  In _International conference on machine learning_ , pages 2048–2057.
  * Yin et al. (2017) Zi Yin, Keng-hao Chang, and Ruofei Zhang. 2017.  Deepprobe: Information directed sequence understanding and chatbot design via recurrent neural networks.  In _Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pages 2131–2139. ACM.
  * Yin and Shen (2018) Zi Yin and Yuanyuan Shen. 2018.  On the dimensionality of word embedding.  In _Advances in Neural Information Processing Systems_ , pages 895–906.