# Spin-tensor Meissner currents of ultracold bosonic gas in an optical lattice

Xiaofan Zhou  State Key Laboratory of Quantum Optics and Quantum Optics
Devices, Institute of Laser spectroscopy, Shanxi University, Taiyuan 030006,
China  Collaborative Innovation Center of Extreme Optics, Shanxi University,
Taiyuan, Shanxi 030006, China  Suotang Jia  State Key Laboratory of Quantum
Optics and Quantum Optics Devices, Institute of Laser spectroscopy, Shanxi
University, Taiyuan 030006, China  Collaborative Innovation Center of Extreme
Optics, Shanxi University, Taiyuan, Shanxi 030006, China  Xi-Wang Luo
<EMAIL_ADDRESS>CAS Key Laboratory of Quantum Information, University of
Science and Technology of China, Hefei, Anhui 230026, China  Synergetic
Innovation Center of Quantum Information and Quantum Physics, University of
Science and Technology of China, Hefei, Anhui 230026, China

###### Abstract

We investigate the Meissner currents of interacting bosons subjected to a
staggered artificial gauge field in a three-leg ribbon geometry, realized by
spin-tensor–momentum coupled spin-1 atoms in a 1D optical lattice. By
calculating the current distributions using the state-of-the-art density-
matrix renormalization-group method, we find a rich phase diagram containing
interesting Meissner and vortex phases, where the currents are mirror
symmetric with respect to the middle leg (i.e., they flow in the same
direction on the two boundary legs opposite to that on the middle leg),
leading to the spin-tensor type Meissner currents, which is very different
from previously observed chiral edge currents under uniform gauge field. The
currents are uniform along each leg in the Meissner phase and form vortex-
antivortex pairs in the vortex phase. Besides, the system also support a
polarized phase that spontaneously breaks the mirror symmetry, whose ground
states are degenerate with currents either uniform or forming vortex-
antivortex pairs. We also discuss the experimental schemes for probing these
phases. Our work provides useful guidance to ongoing experimental research on
synthetic flux ribbons and paves the way for exploring novel many-body
phenomena therein.

##  I Introduction

Charged particles in a magnetic field showcase a remarkable variety of
macroscopic quantum phenomena, including quantized Hall resistance in
topological insulators Hasan2010 ; Qi2011 , Meissner effect in superconductors
Braunisch1992 ; Geim1998 . Recent experimental advances in realizing synthetic
gauge field in ultracold atomic system provide a powerful tool for exploring
these novel phenomena in a fully controllable, clean environment Lin2011 ;
Wang2012 ; Cheuk2012 ; Hamner2015 ; Wu2016 ; Huang2016 ; Li2017 ; SOCclock .
Chiral Meissner (topological edge) currents have been observed experimentally
with an atomic Bose (Fermi) gas in both ladder Boseladder1 ; Boseladder2  and
three-leg ribbon geometries Mancini2015 ; RamanRb , where the two-dimensional
(2D) lattice consists of the sites of an 1D optical lattice in the long
direction and the internal atomic spin states forming a synthetic lattice in
the short direction. Such synthetic dimensions enable investigating higher-
dimensional physics beyond the physical dimensions of the systems and bring
new opportunities in exploring novel quantum phenomena Boada2012 ; SDprl14 ;
Price2015 . The experimentally observed chiral currents result from an uniform
gauge field that is equivalent to ordinary spin-orbit coupling (i.e., spin-
vector–momentum coupling) with the form $qF_{z}$, where $q$ is the quasi-
momentum along the optical lattice direction and
$\mathbf{F}=\\{F_{x},F_{y},F_{z}\\}$ are the spin-vectors Zhai2015 .
Consequently, the corresponding Meissner currents also exhibit spin-vector
properties, that is, atoms with opposite spins propagate along opposite
directions Boseladder2 ; Mancini2015 ; RamanRb ; Livi2016 .

For higher spin ($\geq 1$) systems (e.g., the synthetic three-leg ribbon), it
is well known that there exist not only spin-vectors but also spin-tensors
tensor1 ; tensor2 . Spin-tensor–momentum coupling (STMC) in the form
$qF_{z}^{2}$ has been proposed and experimentally realized recently Luo2017 ;
Li2020 , which can significantly modify the band structures (e.g., dark-state
band, triply degenerate points) and lead to interesting many-body physics
(e.g., magnetic stripe phase) in the presence of interactions Luo2017 ;
Mai2018 ; Hu2018 ; Lei2022 ; Zhang2022 . It was also shown that STMC in a 1D
Mott bosonic lattice can support interesting spin-tensor magnetism orders
zhou2020 . So far, the studies of Meissner effects have been focused on the
aforementioned spin-vector types with uniform synthetic gauge fields.
Therefore, a natural and important question is to explore the Meissner effects
of interacting atoms in optical lattices with STMC.

In this paper, we investigate the ground-state Meissner currents of spin-
tensor–momentum coupled spin-1 bosons in a 1D optical lattice, using state-of-
the-art density-matrix renormalization-group (DMRG) numerical methods dmrg1 ;
dmrg2 . The system corresponds to a synthetic three-leg ribbon with a
staggered gauge field, where the magnetic domain wall is given by the middle
leg. We are interested in the Meissner effects and distinguish different
phases by examining the current and momentum distributions. In the non-
interaction limit, there are two phases depending on the minima of the lowest
single-particle band. The Meissner (M) phase with a single band minimum occurs
when the inter-leg couplings are strong, where the Meissner current is uniform
along each leg with its amplitude determined by the gauge field strength. As
the inter-leg coupling decreases across some critical value, the system
undergoes a transition to the vortex (V) phase with double band minima, where
atoms can occupy both minima and their interference can lead to vortex
structures of the currents, forming vortex-antivortex pairs. For both the M
and V phases, the current distributions are always mirror symmetric with
respect to the middle leg, they have a rank-2 spin-tensor form which is very
different from the spin-vector Meissner current. In particular, here the
current flows in the same direction on the two boundary legs, opposite to the
current direction on the middle leg. Meanwhile, the rank-0 scalar charge
current and rank-1 spin-vector current are both zero.

In the presence of interaction, the phase diagram is altered significantly. In
the weak interaction region, ferromagnetic spin-spin interaction can dominate
over density-density interaction and stabilize the V phase with equal
populations on the two band minima. In the strong interaction region, the wave
packet becomes broadened in the momentum space and the two wave packets at the
two band minima would merge into one, and thus the system favors the M phase
even when the lowest band has two minima. Interestingly, for intermediate
interaction strength, the spin-spin interaction can induce a new polarized (P)
phase that spontaneously breaks the mirror symmetry with respect to the middle
leg, though the STMC (i.e., synthetic gauge field) preserves such symmetry. In
the P phase, atoms start to occupy the dark middle band induced by the STMC,
and the ground-states are degenerate with currents either uniform or forming
vortex-antivortex pairs. The P phase resulting from the presence of dark
middle band is unique for the STMC systems, which is absent for spin-vector-
coupled systems. The Meissner currents in all phases persist no matter if the
system belongs to a superfluid or Mott insulator.

##  II Model and Hamiltonian

Figure 1: (a) Schematics of the system setup. The ultracold atoms are trapped
in a 1D optical lattice, which is generated by a pair of counterpropagating
lasers with the wavelength $\lambda_{\mathrm{L}}$ (green arrows). Three Raman
lasers (yellow, blue and red arrows) with wavelength $\lambda_{\Omega}$ and an
angle $\eta$ with respect to $y$-direction induce two Raman transitions. (b)
Raman transitions between the spin states $|0\rangle$ and
$|\\!\\!\uparrow\\!\\!(\downarrow)\rangle$ with detuning $\delta$. (c) The
effective three-leg ribbon with staggered gauge field along synthetic
dimension. Figure 2: (a) Single-particle band structure of the Hamiltonian
Eq.(3) for different values of $\Omega/t$ with flux $\phi/\pi=0.5$ and
$\Omega_{c}/t=2$. The system exhibits three bands, including upper band (blue
line), lowest band (red line) and the middle band (black dashed line). Blue
(red) lines with different markers correspond to the upper (lowest) band with
different $\Omega$, the square, dot, triangle, and star mark the bands with
$\Omega/t=0,1,2,3$, respectively. The middle band is the dark band with spin
state $|$-$\rangle$, and is unaffected by $\Omega$. The spin components
$|0\rangle$ and including $|$+$\rangle$ are indicated around the corresponding
band minima. The lowest band has a single minimum for $\Omega>\Omega_{c}$ and
two symmetric band minima for $\Omega<\Omega_{c}$. (b) The non-interacting
phase diagram in the plane $\Omega\\!\\!-\\!\\!\phi$, with quantum phases
vortex (V) and Meissner (M). The color bar denotes the current
$J_{\uparrow}^{\parallel}$, and the solid line is the critical line given by
$\Omega_{c}$. (c) (d) The profiles of the current $J_{\sigma}^{\parallel}$ as
functions of $\Omega$ and $\phi$, respectively, with $\phi/\pi=1/2$ in (c) and
$\Omega/t=2.0$ in (d). The filling factor is $\rho=1.0$ in (b)-(d). We set $t$
as energy unit.

We consider an experimental setup based on Bose-Einstein condensate (BEC) in a
1D optical lattice with STMC, as shown in Fig. 1(a). A pair of counter
propagating laser with wavelength $\lambda_{\mathrm{L}}$ (green arrows) is
used to realize the 1D optical lattices
$V_{\mathrm{lat}}\left(y\right)=-V_{0}\cos^{2}\left(k_{\mathrm{L}}y\right)$
along the $y$-direction, with wave number
$k_{\mathrm{L}}=2\pi/\lambda_{\mathrm{L}}$. Three Raman lasers $R_{1}$
(yellow), $R_{2}$ (blue) and $R_{3}$ (red) with wavelength $\lambda_{\Omega}$
and an incident angle $\eta$ with respect to lattice direction induce two
Raman transitions between the spin states $|0\rangle$ and
$|\\!\\!\uparrow\\!\\!(\downarrow)\rangle$ accompanied with the momentum
transfer $2k_{\mathrm{R}}$ Li2020 , in which
$k_{\mathrm{R}}=2\pi\cos(\eta)/\lambda_{\Omega}$, as shown in Fig. 1(b). The
tight-binding Hamiltonian reads

$\displaystyle H$ $\displaystyle=$ $\displaystyle-t\sum_{\left\langle
j,j^{\prime}\right\rangle,\sigma}\hat{b}_{j\sigma}^{{\dagger}}\hat{b}_{j^{\prime}\sigma}+\frac{\Omega}{\sqrt{2}}\sum_{j}\left(e^{i\phi
j}\hat{b}_{j0}^{{\dagger}}\hat{b}_{j+}+\mathrm{H.c.}\right)$ (1)
$\displaystyle+$
$\displaystyle\sum_{j}\frac{U_{0}}{2}\hat{n}_{j}\left(\hat{n}_{j}-1\right)+\frac{U_{2}}{2}\left(\mathbf{S}_{j}^{2}-2\hat{n}_{j}\right)+\delta
S_{j,z}^{2},$

where $\hat{b}_{j\sigma}^{{\dagger}}$ ($\hat{b}_{j\sigma}$) is the Bose
creating (annihilating) operator with spin basis
$\sigma=\\{\uparrow,0,\downarrow\\}$ and
$\hat{b}_{j\pm}=(\hat{b}_{j\uparrow}\pm\hat{b}_{j\downarrow})/\sqrt{2}$. The
particle number operator is $\hat{n}_{j}=\sum_{\sigma}\hat{n}_{j\sigma}$ with
$\hat{n}_{j\sigma}=\hat{b}_{j\sigma}^{{\dagger}}\hat{b}_{j\sigma}$. $t$ is the
tunneling amplitude between neighboring lattice sites, $\Omega$ is the Raman
coupling strength. The flux
$\phi=2\pi\cos(\eta)\lambda_{\mathrm{L}}/\lambda_{\Omega}=2\cos(\eta)k_{\mathrm{R}}a$
can be tuned by tuning the angle $\eta$. $\delta$ is the detuning for both
$\left|\uparrow\right\rangle$ and $\left|\downarrow\right\rangle$ states. For
simplicity, we set $\delta=0$ in the following. $U_{0}$ and $U_{2}$ are the
density-density and spin-spin on-site interactions, which are related to the
scattering lengths $a_{0,2}$ (corresponding to the channels with total spin 0
and 2 respectively) as $U_{0}=4\pi\hbar^{2}(a_{0}+2a_{2})/3M$ and
$U_{2}=4\pi\hbar^{2}(a_{2}-a_{0})/3M$, where $M$ is the mass of the atom
Ho1988 ; Ohmi1998 . The spin dependent interaction for 23Na is
antiferromagnetic with $U_{2}/U_{0}>0$, and ferromagnetic for 87Rb with
$U_{2}/U_{0}<0$.
$\mathbf{S}_{j}=\sum_{\sigma\acute{\sigma}}\hat{b}_{j\sigma}^{{\dagger}}\mathbf{F}_{\sigma\acute{\sigma}}\hat{b}_{j\acute{\sigma}}$
is the total spin at site $j$, where $\mathbf{F}_{\sigma\acute{\sigma}}$
represent the matrix elements of the spin-vector operator. We also set $t$ as
energy unit.

##  III Phase diagram

###  III.1 Band structures

We firstly discuss the single-particle band structure of the system.
Effectively, the system corresponds to a three-leg ribbon with staggered gauge
field as shown in Fig. 1(c). The Hamiltonian without interaction can be
written as after a unitary transformation

$H_{0}=-te^{i\phi\gamma_{\sigma}}\\!\\!\\!\\!\\!\sum_{<j,j^{\prime}>,\sigma}\\!\\!\\!\hat{c}_{j\sigma}^{{\dagger}}\hat{c}_{j^{\prime}\sigma}+\frac{\Omega}{\sqrt{2}}\sum_{j}(\hat{c}_{j0}^{{\dagger}}\hat{c}_{j+}+\mathrm{H.c.})$
(2)

where
$(\hat{c}_{j\uparrow},\hat{c}_{j0},\hat{c}_{j\downarrow})^{T}=U(\hat{b}_{j\uparrow},\hat{b}_{j0},\hat{b}_{j\downarrow})^{T}$
with $U=\exp(i\phi jF_{z}^{2})$, and $\gamma_{\uparrow}=1$, $\gamma_{0}=0$,
$\gamma_{\downarrow}=1$. Then, we Fourier transform the Hamiltonian and obtain

$H_{0}=\sum_{q}\hat{C}_{q}^{{\dagger}}\left[-2t\cos(qa+\phi
F_{z}^{2})+\frac{\Omega}{\sqrt{2}}F_{x}\right]\hat{C}_{q},$ (3)

where
$\hat{C}_{q}=(\hat{c}_{q\uparrow},\hat{c}_{q0},\hat{c}_{q\downarrow})^{T}$ is
the corresponding operator in quasi-momentum $q$ space, and
$a=\pi/k_{\mathrm{L}}$ is the lattice constant, $F_{z}^{2}$ is the rank-2 spin
tensor. We see that the flux $\phi$ now represent the strength of STMC. The
system exhibits three bands after diagonalizing the Hamiltonian Eq.(3), as
shown in Fig. 2(a). The top and bottom bright-state bands exhibit the same
behavior as a spin-orbit-coupled spin-1/2 system with spin states $|0\rangle$
and
$|$+$\rangle=(\left|\uparrow\right\rangle+\left|\downarrow\right\rangle)/\sqrt{2}$
[see the color lines in Fig. 2(a)]. The middle band $E_{q}=-2t\cos(q+\phi)$ is
dark-state band which always has the spin state
$|$-$\rangle=(\left|\uparrow\right\rangle-\left|\downarrow\right\rangle)/\sqrt{2}$
and is independent from $\Omega$ [see the black dash line in Fig. 2(a)] since
it is decoupled from the Raman lasers. However, the dark-state band plays an
important role on both ground-state and Meissner current distributions in the
presence of interactions.

###  III.2 Order parameters

Before we calculate the phase diagram of the system, in the following, we
first introduce some order parameters that we will use to distinguish
different phases.

_Currents.—_ Since we are interested in the Meissner effects, the most
important property here would be the current along the legs and the rungs of
the ribbon. Local averaged currents along the leg-direction can be define as
Boseladder1 ; Boseladder2 ; helicalliquidnc ; zhou2017current ,

$\displaystyle
J_{j,\sigma}^{\parallel}=it\langle\hat{b}_{j+1\sigma}^{{\dagger}}\hat{b}_{j\sigma}-h.c.\rangle,J_{\sigma}^{\parallel}=\frac{1}{L}\sum_{j}J_{j,\sigma}^{\parallel}$
(4)

Similarly, we can define the rungs currents along the synthetic direction
Boseladder2 ; zhou2017current ,

$\displaystyle J_{j,0s}^{\perp}=i\Omega\langle e^{i\phi
j}\hat{b}_{j0}^{{\dagger}}\hat{b}_{js}-h.c.\rangle$ (5)

with $s=\\{\uparrow,\downarrow\\}$ . We then introduce the current orders

$\displaystyle
I_{0s}^{\perp}=\frac{1}{L}\sum_{j}\left|J_{j,0s}^{\perp}\right|\text{ and
}I_{z}^{\parallel}=\frac{1}{L}\sum_{j}|J_{j,\uparrow}^{\parallel}-J_{j,\downarrow}^{\parallel}|.$
(6)

The rungs current order $I_{0s}^{\perp}$ distinguishes the M phase (with
vanishing $I_{0s}^{\perp}$) from the V phase (with finite $I_{0s}^{\perp}$).
The leg current order $I_{z}^{\parallel}$ characterizes current polarization
associated to the mirror symmetry with respect to the middle leg (a non-
vanishing $I_{z}^{\parallel}$ represent the breaking of the mirror symmetry).
Different phases can be identified by different current distributions and
orders.

_Spin moments.—_ Experimentally, the spin moments of the atoms can be measured
by spin-selective imaging of the density distributions

$\displaystyle
n_{j,\sigma}=\langle\hat{b}_{j\sigma}^{{\dagger}}\hat{b}_{j\sigma}\rangle.$
(7)

The density distributions can be measured either in the basis
$\sigma=\\{\uparrow,0,\downarrow\\}$ or the basis $\sigma=\\{0,\pm\\}$. We can
define the spin moment orders as

$\displaystyle n_{\sigma}=\frac{1}{L}\sum_{j}n_{j,\sigma}$ (8)

which gives the population properties of the bands and can be used to
characterize the phase transition. We would like to mention that, $n_{-}$
corresponds to the occupation of the dark middle band.

_Entropy.—_ The many-body interacting phase diagram can be further confirmed
by entanglement entropy. Sharp features would emerge at the critical points in
the von Neumann entropy defined as Flammia2009 ; Hastings2010 ; Daley2012 ;
Abanin2012 ; jiang2012 ; Islam2015 ; Amico2008

$\displaystyle
S_{\mathrm{vN}}=-\mathrm{Tr_{A}[\hat{\rho}_{A}\log\hat{\rho}_{A}]}$ (9)

in which $\hat{\rho}_{A}=\mathrm{Tr_{B}|\psi\rangle\langle\psi|}$ is the
reduced density matrix, and $|\psi\rangle$ is the ground-state wavefunction
with $A,B$ corresponding to the left and right half of the 1D chain.

_Low energy level spacing.—_ The energy level spacing with respect to the
ground state energy, defined as

$\displaystyle\Delta_{i}=E_{i}-E_{0},$ (10)

may also be used to signal the phase transition. Where $E_{0}$ is the lowest
eigenenergy, and the $E_{i}$ is the $i$-th order eigenenergy with fixed
particle numbers.

###  III.3 Noninteracting phase diagram

In the noninteracting limit, the phase diagram is determined by the single-
particle band structures. The lowest band presents a single minimum in the M
phase $\Omega>\Omega_{c}$ and two symmetric minima in the V phase
$\Omega<\Omega_{c}$, where $\Omega_{c}$ is the critical coupling strength, as
shown in Fig. 2(a). According to the number of minima, we can plot the single-
particle phase diagram in the plane $\Omega\\!\\!-\\!\\!\phi$, as shown in
Fig. 2(b), which is similar as the phase diagram of spin-1/2 SOC system
Boseladder2 . However, the currents have a rank-2 spin-tensor form that are
very different from the spin-vector Meissner current Boseladder2 ; RamanRb .
In particular, here the current flows in the same direction on the boundary
legs, opposite to the current direction on the middle leg. We calculate the
ground-state currents of Hamiltonian Eq.(1), the averaged currents along each
leg are shown in Figs. 2(c) and 2(d). It is easy to check that the rank-0
scalar charge current $J^{\parallel}\equiv\sum_{\sigma}J_{\sigma}^{\parallel}$
and rank-1 spin-vector current $J_{z}^{\parallel}\equiv\sum_{\sigma}\sigma
J_{\sigma}^{\parallel}$ are both zero, while the rank-2 spin-tensor current
$J_{zz}^{\parallel}\equiv\sum_{\sigma}\sigma^{2}J_{\sigma}^{\parallel}=2J_{\uparrow}^{\parallel}$
is finite. As we increase $\Omega$ for a fixed $\phi$, $J_{zz}^{\parallel}$
increases and has a saturated value above the critical point, as shown in Fig.
2(c), the system undergoes a transition from V to M phase. When increasing the
flux $\phi$ with $\Omega/t=2.0$, $J_{zz}^{\parallel}$ increase firstly and
then decrease, reaching it maxima at the phase transition point from M to V
phase, as shown in Fig. 2(d).

Figure 3: The schematic of density distributions $n_{j,\sigma}$, currents
$J^{\parallel}_{j,\sigma}$ and $J^{\perp}_{j,0s}$ in different phases. The
points denote the density profile of the three-component with point size the
amount of atoms. The current direction and strength is denoted by the arrows.
(a) The M phase with $\Omega/t=2.5$ and $\phi/\pi=1/2$. (b1)-(b3) The V phase
with $\Omega/t=0.2$ and $\phi/\pi=1/4$. (b1) Left minimum ground-state [Blue
dot in Fig. 2(a)]. (b2) Right minimum ground-state [Red dot in Fig. 2(a)].
(b3) Equal superposition of the two minima.

In order to clearly show the local current properties of the V and M phases,
we plot the schematic of density distributions $n_{j,\sigma}$, currents
$J^{\parallel}_{j,\sigma}$ and $J^{\perp}_{j,0s}$, in Fig. 3. The current is
mirror symmetric with respect to the middle leg, analogous to the combination
of spin-1/2 (two-leg ladder) vector current and it’s mirror reflection. In the
M phase, both the currents and the densities are uniform, as shown in Fig.
3(a). In the V phase, the system exhibits two energy minima in the lowest
band. Atoms may populate on either minimum with the same tensor current as
shown in Figs. 3(b1) and 3(b2), or a linear combination of the two minima with
current vortex as shown in Fig. 3(b3).

Figure 4: (a) The current order parameters $I_{z}^{\parallel}$,
$I_{0s}^{\perp}$ versus interaction strength $U_{0}$. (b) Spin moment orders
$n_{\pm,0}$ versus $U_{0}$. (c) Entropy $S_{\text{vN}}$ and (d) Excited energy
gap $\Delta_{i}$ as functions of interaction strength $U_{0}$. In all
subfigures, we have $\Omega/t=0.2$, $\phi/\pi=1/4$, $U_{2}=-0.2U_{0}$,
$\rho=1.0$ and $L=64$.

###  III.4 Interacting phase diagram

The interplay between single-particle band structure and interaction is
fundamental to many areas of modern physics. The ability to engineer competing
interactions between atoms through lattice depth or Feshbach resonances
Chin2010  makes cold atomic system an ideal platform to generate a rich
variety of many-body phenomena Raghu2008 ; Bloch2008 ; Mueller2017 ;
Rachel2018 ; Junemann2017 ; zhou2017SPT ; Leseleuc2019 . Here, we perform the
state-of-the-art DMRG calculations to calculate the many-body ground-states of
the system under open boundary conditions. In our numerical simulations, we
set the cutoff of single-site atom number as $n_{\text{cutoff}}=4$. The effect
of such a cutoff can be neglected for strong interaction, while in the weak
interaction region, the phases and phase boundaries may be slightly affected
by the cutoff. We find that $n_{\text{cutoff}}=4$ is enough to determine the
phase boundaries (see Appendix A). We set lattice size up to $L=64$, for which
we retain 300 truncated states per DMRG block and perform 20 sweeps with a
maximum truncation error $\sim 10^{-7}$.

We consider the repulsive density-density interaction with $U_{0}>0$ and
ferromagnetic spin-spin interaction with $U_{2}<0$ as for 87Rb or 7Li atoms.
In Fig. 4, we plot the order parameters as functions of the interacting
strength $U_{0}/t$, with fixed $\Omega/t=0.2$ and $\phi/\pi=1/4$, from which,
three phases (i.e., V, P and M) as well as their critical boundaries can be
identified. For weak interactions, the system stays in V phase with a large
rung current $I_{0s}^{\perp}$; for strong interactions, the system favors the
M phase with nearly vanishing $I_{0s}^{\perp}$. We want to mention that,
finite-size effects would induce a tiny averaged rungs current
$I_{0s}^{\perp}$ in the M phase, since the leg currents must form a closed
loop at the boundary through rungs. In both the V and M phase, we have
$I_{z}^{\parallel}=0$ due to the mirror symmetry, though
$J_{\sigma}^{\parallel}\neq 0$. For intermediate interaction strength, the P
phase emerges with $I_{z}^{\parallel}\neq 0$, as shown in Fig. 4(a), and the
mirror symmetry is spontaneously broken. The P phase may also support vortex
currents, leading to a large rungs current order $I_{0s}^{\perp}$. Different
phases have different spin occupations, therefore the spin moment orders also
exhibit sharp transitions at the phase boundary, as shown in Fig. 4(b).

Besides the currents and spin moments, the von Neumann entropy
$S_{\mathrm{vN}}$ as well as the low energy level spacing $\Delta_{i}$ can
also signal the transitions. When the entanglement entropy $S_{\mathrm{vN}}$,
an analytic function of correlations, is not analytic at some point, it must
correspond to a quantum phase transition (as long as the definition of the
entanglement entropy is analytic at that point). As demonstrated in Fig. 4(c),
sharp features in the entropy emerge at the critical phase transition points
($S_{\mathrm{vN}}$ or its derivatives show discontinuity) Amico2008 . The
degeneracy of the ground state can be seen from low energy level spacing which
signals the spontaneous symmetry breaking across the phase transition. As
demonstrated in Fig. 4(d), $\Delta_{i}$ vanish for $i<4$ in the P phase due to
ground-state degeneracy.

Figure 5: (a)-(d) The schematic of density distributions $n_{j,\sigma}$,
currents $J^{\parallel}_{j,\sigma}$ and $J^{\perp}_{j,0s}$ of the four-fold
degenerate ground-state in P phase, respectively. In all subfigure, we have
$\Omega/t=0.2$, $\phi/\pi=1/4$, $U_{0}/t=4.0$, $U_{2}=-0.2U_{0}$, $\rho=1.0$
and $L=64$ (only the central 17 sites are shown in the plot). Figure 6: The
phase diagrams in the planes of (a) $\Omega\\!\\!-\\!\\!U_{0}$ with
$\phi/\pi=1/4$, and (b) $\Omega\\!\\!-\\!\\!\phi$ with $U_{0}/t=2$. The red
dashed line in (a) is the phase transition points between superfluid (left
sides) and Mott insulator (right sides) phases, the black dashed line in (b)
is the single-particle phase boundary. (c) The finite-size scaling of the
critical points $U^{c}_{0}/t$ of phase transitions, with $\Omega/t=0.2$ and
$\phi/\pi=1/4$. In (a)-(c), $U_{2}=-0.2U_{0}$. (d) The current order
parameters $I_{z}^{\parallel}$, $I_{0s}^{\perp}$ versus interaction strength
$U_{2}$, with $\Omega/t=0.2$, $\phi/\pi=1/4$, $U_{0}/t=2.0$ and $L=64$. In all
subfigure, we have $\rho=1.0$.

The interaction can drive the system from V (at weak interaction regime) to M
(at strong interaction regime) phase though the lowest single-particle band
has two minima. Different from the noninteracting V phase where the ground-
state is an arbitrary linear combination of the two band minima, here the
interacting V phase is stabilized by ferromagnetic spin-spin interaction,
leading to an equal superposition of the two band minima. On one hand, the
interference between two band minima would induce density modulation and
increase density-density interaction energy. On the other hand, the
superposition also generate spin orders which lower the ferromagnetic spin-
spin interaction energy, and the system favors the V phase since the spin-spin
interaction dominates over the density-density interaction. We want to mention
that, for the V phase in the weak interacting region, the mean-field analysis
predicts that the relative phase between atoms at two band minima can either
be $0$ or $\pi$ Luo2017 , the ground-states should be two-fold degenerate,
which are related with each other through shifting the vortex by half period
along the leg direction. However, such degeneracy is lifted by the finite-size
and finite-$n_{\text{cutoff}}$ effects. In the strong interacting regime, the
atom distribution in momentum space is broadened and the two wave packets at
the two band minima would merge into one centered at $q=0$, therefore the
system favors the M phase even when the lowest band has two minima. The
current distributions of the interacting M and V phases are very similar as
that for non-interacting case shown in Figs. 3(a) and 3(b3).

Interestingly, a novel P phase emerges at intermediate interaction strength.
The typical current and density distributions are shown in Fig. 5 for flux
$\phi/\pi=1/4$. In this P phase, the mirror symmetry with respect to the
middle leg is spontaneously broken and atoms start to occupy the dark middle
band with $n_{-}\neq 0$. Besides the mirror symmetry, the $Z_{2}$ exchange
symmetry $|0\rangle\\!\\!\leftrightarrow\\!\\!|+\rangle$ is also broken,
leading to a four-fold degeneracy of the ground-states. The symmetry breaking
can be seen from the spin moments in Fig. 4(b) with $n_{0}\neq n_{+}$ and also
the density distributions in Fig. 5 which indicate that $n_{j\uparrow}\neq
n_{j\downarrow}$. Notice that the M and V phase preserve the mirror symmetry
and $Z_{2}$ exchange symmetry with $n_{0}=n_{+}$ and
$n_{\uparrow}=n_{\downarrow}$. The P phase can be a combined population of the
middle band with either left minimum [Figs. 5(a) and 5(b)] or right minimum
[Figs. 5(c) and 5(d)] of the lowest band, and for either combination, the
polarization $S_{j,z}=n_{j,\uparrow}-n_{j,\downarrow}$ is nonvanishing and can
be either positive or negative (related by the mirror symmetry), which is
determined by the relative phase between atoms on the two bands. For the two
ground-states with dominate population on the left minimum of the lowest band,
vortex-antivortex pairs are formed along the leg direction, leading to
modulations of spin densities and currents with period $2\pi/\phi$. We want to
mention that the order parmeters of the P phase in Figs. 4(a)-4(c) are
obtained from the vortex ground-state. Moreover, due to the breaking mirror
symmetry, the local spin-vector current $\sum_{\sigma}\sigma
J_{j,\sigma}^{\parallel}$ and spin-tensor current
$\sum_{\sigma}\sigma^{2}J_{j,\sigma}^{\parallel}$ are both finite in the P
phase. This dark middle band has spin state orthogonal to the lowest band, so
the superposition of the two bands will only induce spin orders without
generating density modulation. Notice that both the P and V phases have a
lower spin-spin interaction energy, the P phase has a higher single-particle
energy but a lower density-density interaction energy compared to V phase,
therefore a transition from V to P phase occurs as we increase the density-
density interaction strength.

In Figs. 6(a) and 6(b), we map out the phase diagrams in
$\Omega\\!\\!-\\!\\!U_{0}$ plane with $\phi/\pi=1/4$ and in
$\Omega\\!\\!-\\!\\!\phi$ plane with $U_{0}/t=2$, respectively. We see that
the M phase area is enlarged by interactions. Fig. 6(a) clearly shows that the
V and P phases are replaced eventually by the M phase as the interaction
strength increases. The system favors the P phase for weak inter-leg coupling
$\Omega$ and small flux $\phi$, where dark middle band can be occupied more
easily since the energy gap with the lowest band is small. For very small
$\Omega$, the mixing between states $|0\rangle$ and $|+\rangle$ is also weak,
so the density modulation and thereby the density-density interaction energy
in the V phase is weak, and the system favors the V phase instead of the P
phase that has higher single-particle energy. Similarly, when the flux is
large, density-density interaction energy for the V phase is small compared to
the high single-particle energy in the P phase, and the system favors V over P
phase, as shown in Fig. 6(b). As $\Omega$ increases across some critical
value, the density-density interaction energy for the V phase becomes stronger
compared to the ferromagnetic spin-spin interaction energy, atoms will
spontaneously populate around one of the two band minima with a single wave
packet centered around $q\neq 0$ and the system enters the M phase, further
increasing $\Omega$ will drive the wave packet center to $q=0$. For a small
flux, the system may directly enter the M phase with one wave packet centered
at $q=0$. In this paper, we will not distinguish the two different types of M
phase (centered around $q=0$ and $q\neq 0$) since they have similar current
distributions.

In obtaining the phase diagram, we have employed the finite-size scaling,
leading to the critical points of phase transitions in the thermodynamic
limit, which are almost the same as those of finite-size systems, as shown in
Fig. 6(c). The finite-size effects on the order parameters as well as the
energy level spacing and gap are discussed in Appendix B. As we discussed
previously, the V and P phases result from the ferromagnetic spin-spin
interaction. For antiferromagnetic spin-spin interactions $U_{2}/U_{0}>0$, one
only has the M phase. To see this, we plot the current order parameters as
function of spin-dependent interaction strength $U_{2}$ in Fig. 6(d), showing
that the P phase only exists in the ferromagnetic region. We want to emphasize
that the Meissner currents induced by the gauge field persist even when the
system enters the Mott insulator region, the red dashed line in Fig. 6(a)
shows the boundary between the superfluid and Mott insulator phases, which is
obtained by calculating the chemical potential gap gap1 ; gap2 ; gap3 ; gap4 .

##  IV Discussion and Conclusion

Experimentally, the spin moments can be measured by population imaging after
Stern-Gerlach separation SDprl14 ; RamanRb . The currents can be observed by
site-resolved detection of the quench dynamics Boseladder2 . The current
(motion of atoms) can also be proved by spin-selective time-of-flight imaging
of the lattice momentum distribution Boseladder2 ; Mancini2015

$\displaystyle
n_{\sigma}(k)=\sum_{i,j}e^{ik(i-j)}\langle\hat{b}_{i\sigma}^{{\dagger}}\hat{b}_{j\sigma}\rangle.$
(11)

The current is related to the lattice momentum unbalance, and we can define
the chirality of the atomic motion as

$\displaystyle\chi_{\sigma}=\frac{1}{L}\int_{0}^{\pi}h_{\sigma}(k)dk,$ (12)

which shows similar behaviors as the current $J^{\parallel}_{\sigma}$
Mancini2015 ; Livi2016 , where $h_{\sigma}(k)=n_{\sigma}(k)-n_{\sigma}(-k)$ is
the asymmetry function. Moreover, the entanglement entropy can be measured
using quantum interference of many-body twins of ultracold atoms in optical
lattices Islam2015 .

In summary, We studied the Meissner effects of interacting bosons in a one-
dimensional optical lattice with spin-tensor–momentum coupling. Using state-
of-the-art density-matrix renormalization-group numerical method, we obtained
the phase diagram with a rich variety of interesting phases, including the
Meissner, vortex and polarized phases. The current distributions show spin-
tensor properties and interesting vortex structures that are unique for the
spin-tensor–momentum coupled system of high spin system, which may find
possible applications in atomic spintronic devices. Our work reveals
nontrivial phases and transport properties resulting from the interplay
between spin-tensors, lattice physics and interactions, and thus paves the way
for exploring novel many-body phenomena of interacting particles in nonuniform
gauge fields.

##  Appendix A Effects of $n_{\text{cutoff}}$

As we discussed in the main text, the cutoff of single-site atom number
$n_{\text{cutoff}}$ would modify the phase and phase boundary at the weak
interacting region. A larger $n_{\text{cutoff}}$ leads to more accurate
results, but the calculation becomes more computational expensive. In Fig. 7,
we show the order parameter $n_{-}$ as a function of $U_{0}$ with different
atom number cutoff. We see that, the results are nearly the same in the region
$U_{0}/t>1.5$ for $n_{\text{cutoff}}=3,4,5$, while in the weak interacting
region $U_{0}/t<1$, the results with $n_{\text{cutoff}}=3$ deviate
significantly from that with $n_{\text{cutoff}}=4,5$. This can be easily
understood by noticing that, for weaker interaction, the probability of high
occupation is larger in the ground state. On the other hand, we find that the
results for $n_{\text{cutoff}}=4$ and $n_{\text{cutoff}}=5$ are nearly the
same (both give the phase boundary around $U_{0}/t=1$), indicating that
$n_{\text{cutoff}}=4$ is enough to determine the phase boundaries (notice that
most of the phase boundaries are located in the region $U_{0}/t\geq 1$ in the
phase diagram).

Figure 7: Spin moment order $n_{-}$ as a function of interaction strength
$U_{0}/t$, for several $n_{\text{cutoff}}$, with $\Omega/t=0.2$,
$\phi/\pi=1/4$, $U_{2}=-0.2U_{0}$, $\rho=1.0$ and $L=32$. Figure 8: (a) The
current order parameter $I_{z}^{\parallel}$, (b) the current order parameters
$I_{0\downarrow}^{\perp}$, and (c) the spin moment orders $n_{0}$ versus
$U_{0}$ for $L=32$ and $64$. (d) The finite-size scaling of the
$I_{0\uparrow}^{\perp}$ for $U_{0}/t=4.0$ and $7.0$. In all subfigures, we
have $\Omega/t=0.2$, $\phi/\pi=1/4$, $U_{2}=-0.2U_{0}$, and $\rho=1.0$. Figure
9: (a) Entropy $S_{\text{vN}}$ versus $U_{0}/t$, for $L=32$ and $64$. (b)
$\Delta_{i}$ versus $U_{0}/t$ for $L=32$. The finite-size scaling of
$\Delta_{i}$ for (c) V phase with $U_{0}/t=0.6$ and M phase with $7.0$, and
(d) P phase with $U_{0}/t=3.0$ and $5.2$. (d). In all subfigure, we have
$\Omega/t=0.2$, $\phi/\pi=1/4$, $U_{2}=-0.2U_{0}$, $\rho=1.0$. The unit of
$\Delta_{i}$ is $t$.

##  Appendix B Finite-size effects

Now we briefly discuss the finite-size effects on the order parameters. As the
system size increases, the order parameters may slightly modified (especially
near the critical phase boundaries), and thus the phase boundaries (determined
by examining the spontaneous symmetry breaking) may shift slightly as we
increase the size. The phase boundaries at the thermodynamic limit are
obtained by finite-size scaling [see Fig. 6(c)]. In Figs. 8(a)-(c), we plot
the order parameters as functions of $U_{0}$ for different system size $L=32$
and $L=64$, away from the phase boundaries, we see that the order parameters
are almost the same for different sizes. In Fig. 8(d), we plot the typical
finite-size scaling of the order parameters and find that the thermodynamic
value is very close to the finite-size value with $L=64$. We would like to
mention that, the rung current in the M phase [see the red line in Fig. 8(d)]
decreases as the system size increase, this is because, the leg currents must
form a closed loop at the boundary through rungs, such boundary-rung current
lead to an averaged rung current scales as $1/L$. Numerically, we find that
the averaged rung current $I_{0s}^{\perp}$ does not vanish completely in the M
phase (though it is extremely small). For such interaction driven M phase,
though atoms occupy a single broadened wave packet in momentum space, the
double-well structure of the lowest band may result in residual interference
that leads to the tiny rung current.

As we discussed in the main text, the von Neumann entropy $S_{\mathrm{vN}}$ as
well as the low energy level spacing $\Delta_{i}$ can also signal the phase
transitions. For different system sizes, the entropy $S_{\text{vN}}$ exhibits
similar curves with peak positions always matching the phase boundaries given
by the order parameters (currents and spin moments), as shown in Fig. 9(a) for
$L=32$ and $64$. For the low energy level spacing, we have shown the results
for $L=64$ in the main text, while the results for $L=32$ are similar [see
Fig. 9(b)], the four fold degeneracy is clearly seen in the P phase. To see
how the energy gap scales with system size, we plot the finite-size scaling
results in Figs. 9(c) and 9(d), notice that the energy gap corresponds to
$\Delta_{1}$ ($\Delta_{4}$) for the V and M phase (P phase). The small gaps in
the P and M phases change slightly with system size and remain finite (of the
order of $\sim 10^{-2}t$) in the thermodynamic limit, which vanish at their
phase boundary. While for the V phase, the gap decreases significantly with
increasing system size, reaching to a small value (of the order of $10^{-3}t$)
in the thermodynamic limit. We also note that, in the V phase, higher excited-
state energy gaps $\Delta_{i}$ ($i>1$) also decrease with increasing system
size, the low energy levels are nearly equally spaced in the thermodynamic
limit. Since the V phase appears in the weak interaction region where the
effects of cutoff $n_{\text{cutoff}}$ is more significant, we expect that the
tiny gap in the V phase may be induced by the cutoff, and the V phase is
probably a gapless phase.

_Acknowledgments.—_ This work is supported by the National Key R&D Program of
China under Grant No. 2022YFA1404003, the National Natural Science Foundation
of China (NSFC) under Grant No. 12004230, 12174233 and 12034012, the Research
Project Supported by Shanxi Scholarship Council of China and Shanxi ’1331KSC’.
X.-W. Luo acknowledges support from the National Natural Science Foundation of
China (Grant No. 12275203) and the USTC start-up funding.

## References

  * (1) M. Hasan and C. Kane, Colloquium: Topological insulators, Rev. Mod. Phys. 82, 3045 (2010).
  * (2) X.-L. Qi and S.-C. Zhang, Topological insulators and superconductors, Rev. Mod. Phys. 83, 1057 (2011).
  * (3) W. Braunisch, N. Knauf, V. Kataev, S. Neuhausen, A. Grütz, A. Kock, B. Roden, D. Khomskii, and D. Wohlleben, Paramagnetic Meissner effect in Bi high-temperature superconductors, Phys. Rev. Lett. 68, 1908 (1992).
  * (4) A. K. Geim, S. V. Dubonos, J. G. S. Lok, M. Henini, and J. C. Maan, Paramagnetic Meissner effect in small superconductors, Nature (London) 369, 144 (1998).
  * (5) Y.-J. Lin, K. Jiménez-García, and I. B. Spielman, Spin-orbit coupled Bose-Einstein condensates, Nature (London) 471, 83 (2011).
  * (6) P. Wang, Z. Q. Yu, Z. Fu, J. Miao, L. Huang, S. Chai, H. Zhai, and J. Zhang, Spin-orbit coupled degenerate Fermi gases, Phys. Rev. Lett. 109, 095301 (2012).
  * (7) L. Cheuk, A. Sommer, Z. Hadzibabic, T. Yefsah, W. Bakr, and M. Zwierlein, Spin-injection spectroscopy of a spin-orbit coupled Fermi gas, Phys. Rev. Lett. 109, 095302 (2012).
  * (8) C. Hamner, Y. Zhang, M. A. Khamehchi, Matthew J. Davis, and P. Engels, Spin-orbit-coupled Bose-Einstein condensates in a one-dimensional optical lattice, Phys. Rev. Lett. 114, 070401 (2015).
  * (9) L. Huang, Z. Meng, P.Wang, P. Peng, S.-L. Zhang, L. Chen, D. Li, Q. Zhou, and J. Zhang, Experimental realization of two-dimensional synthetic spin-orbit coupling in ultracold Fermi gases, Nat. Phys. 12, 540 (2016).
  * (10) Z. Wu, L. Zhang, W. Sun, X.-T. Xu, B.-Z. Wang, S.-C. Ji, Y. Deng, S. Chen, X.-J. Liu, and J.-W. Pan, Realization of two dimensional spin-orbit coupling for Bose-Einstein condensates, Science 354, 83 (2016).
  * (11) J.-R. Li, J. Lee, W. Huang, S. Burchesky, B. Shteynas, F. C. Top, A. O. Jamison, and W. Ketterle, A stripe phase with supersolid properties in spin-orbit-coupled Bose-Einstein condensates, Nature (London) 543, 91 (2017).
  * (12) S. Kolkowitz, S. L. Bromley, T. Bothwell, M. L. Wall, G. E. Marti, A. P. Koller, X. Zhang, A. M. Rey, and J. Ye, Spin-orbit coupled fermions in an optical lattice clock, Nature (London) 542, 66 (2017).
  * (13) A. Petrescu and K. L. Hur, Bosonic Mott insulator with Meissner currents, Phys. Rev. Lett. 111, 150601 (2013).
  * (14) M. Atala, M. Aidelsburger, M. Lohse, J. T. Barreiro, B. Paredes, and I. Bloch, Observation of chiral currents with ultracold atoms in bosonic ladders, Nat. Phys. 10, 588 (2014).
  * (15) M. Mancini, G. Pagano, G. Cappellini, L. Livi, M. Rider, J. Catani, C. Sias, P. Zoller, M. Inguscio, M. Dalmonte, and L. Fallani, Observation of chiral edge states with neutral fermions in synthetic Hall ribbons, Science 349, 1510 (2015).
  * (16) B. K. Stuhl, H.-I. Lu, L. M. Aycock, D. Genkina, and I. B. Spielman, Visualizing edge states with an atomic Bose gas in the quantum Hall regime, Science 349, 1514 (2015).
  * (17) O. Boada, A. Celi, J. I. Latorre, and M. Lewenstein, Quantum simulation of an extra dimension, Phys. Rev. Lett. 108, 133001 (2012).
  * (18) A. Celi, P. Massignan, J. Ruseckas, N. Goldman, I. B. Spielman, G. Juzeliūnas, and M. Lewenstein, Synthetic gauge fields in synthetic dimensions, Phys. Rev. Lett. 112, 043001 (2014).
  * (19) H. M. Price, O. Zilberberg, T. Ozawa, I. Carusotto, and N. Goldman, Four-dimensional quantum Hall effect with ultracold atoms, Phys. Rev. Lett. 115, 195303 (2015).
  * (20) H. Zhai, Degenerate quantum gases with spin-orbit coupling: A review, Rep. Prog. Phys. 78, 026001 (2015).
  * (21) L. F. Livi, G. Cappellini, M. Diem, L. Franchi, C. Clivati, M. Frittelli, F. Levi, D. Calonico, J. Catani, M. Inguscio, and L. Fallani, Synthetic dimensions and spin-orbit coupling with an optical clock transition, Phys. Rev. Lett. 117, 220401 (2016).
  * (22) H. M. Bharath, Non-Abelian geometric phases carried by the spin fluctuation tensor, J. Math. Phys. (N.Y.) 59, 062105 (2018).
  * (23) H. M. Bharath, M. Boguslawski, M. Barrios, L. Xin, and M. S. Chapman, Exploring non-Abelian geometric phases in spin-1 ultracold atoms, Phys. Rev. Lett. 123, 173202 (2019).
  * (24) D. Li, L. Huang, P. Peng, G. Bian, P. Wang, Z. Meng, L. Chen, and J. Zhang, Experimental realization of a spin-tensor momentum coupling in ultracold Fermi gases, Phys. Rev. A 102, 013309 (2020).
  * (25) X.-W. Luo, K. Sun, and C. Zhang, Spin-tensor–momentum-coupled Bose-Einstein condensates, Phys. Rev. Lett. 119, 193001 (2017).
  * (26) X.-Y. Mai, Y.-Q. Zhu, Z. Li, D.-W. Zhang, and S.-L. Zhu, Topological metal bands with double-triple-point fermions in optical lattices, Phys. Rev. A 98, 053619 (2018).
  * (27) H. Hu, J. Hou, F. Zhang, and C. Zhang, Topological triply-degenerate points induced by spin-tensor-momentum couplings, Phys. Rev. Lett. 120, 240401 (2018).
  * (28) Z. Lei, Y. Deng, and C. Lee, Unpaired topological triply degenerate point for spin-tensor-momentum-coupled ultracold atoms, Phys. Rev. Res. 4, 033008 (2022).
  * (29) M. Zhang, X. Yuan, Y. Li, X.-W. Luo, C. Liu, M. Zhu, X. Qin, C. Zhang, Y. Lin, and J. Du, Observation of spin-tensor induced topological phase transitions of triply degenerate points with a trapped ion, Phys. Rev. Lett. 129, 250501 (2022).
  * (30) X. Zhou, X.-W. Luo, G. Chen, S. Jia, and C. Zhang, Quantum spiral spin-tensor magnetism, Phys. Rev. B 101, 140412(R) (2020).
  * (31) S. R. White, Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett. 69, 2863 (1992).
  * (32) U. Schollwök, The density-matrix renormalization group, Rev. Mod. Phys. 77, 259 (2005).
  * (33) T. L. Ho, Spinor Bose condensates in optical traps, Phys. Rev. Lett. 81, 742 (1998).
  * (34) T. Ohmi and K. Machida, Bose-Einstein condensation with internal degrees of freedom in alkali atom gases, J. Phys. Soc. Jpn., 67(6), 1822 (1998).
  * (35) S. Barbarino, L. Taddia, D. Rossini, L. Mazza, and R. Fazio, Magnetic crystals and helical liquids in alkaline-earth fermionic gases, Nat. Commun. 6, 8134 (2015).
  * (36) X. Zhou, J.-S. Pan, W. Yi, G. Chen, and S. Jia, Interaction-induced exotic vortex states in an optical lattice clock with spin-orbit coupling, Phys. Rev. A 96, 023627 (2017).
  * (37) S. T. Flammia, A. Hamma, T. L. Hughes, and X.-G. Wen, Topological entanglement Renyi entropy and reduced density matrix structure, Phys. Rev. Lett. 103, 261601 (2009).
  * (38) M. B. Hastings, I. González, A. B. Kallin, and R. G. Melko, Measuring Rényi entanglement entropy in quantum monte carlo simulations, Phys. Rev. Lett. 104, 157201 (2010).
  * (39) A. J. Daley, H. Pichler, J. Schachenmayer, and P. Zoller, Measuring entanglement growth in quench dynamics of bosons in an optical lattice, Phys. Rev. Lett. 109, 020505 (2012).
  * (40) D. A. Abanin and E. Demler, Measuring entanglement entropy of a generic many-body system with a quantum switch, Phys. Rev. Lett. 109, 020504 (2012).
  * (41) H.-C. Jiang, Z.-H. Wang, and L. Balents, Identifying topological order by entanglement entropy, Nat. Phys. 8, 902 (2012).
  * (42) R. Islam, R. Ma, P. M. Preiss, M. E. Tai, A. Lukin, M. Rispoli, and M, Greiner, Measuring entanglement entropy in a quantum many-body system, Nature 528, 77 (2015).
  * (43) L. Amico, R. Fazio, A. Osterloh, and V. Vedral, Entanglement in many-body systems, Rev. Mod. Phys. 80, 517 (2008).
  * (44) C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Feshbach resonances in ultracold gases, Rev. Mod. Phys. 82, 1225 (2010).
  * (45) S. Raghu, X.-L. Qi, C. Honerkamp, and S.-C. Zhang, Topological Mott insulators, Phys. Rev. Lett. 100, 156401 (2008).
  * (46) I. Bloch, J.n Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Rev. Mod. Phys. 80, 885 (2008).
  * (47) E. J. Mueller, Review of pseudogaps in strongly interacting Fermi gases, Rep. Prog. Phys. 80, 104401 (2017).
  * (48) S. Rachel, Interacting topological insulators: a review, Rep. Prog. Phys. 81, 116501 (2018).
  * (49) J. Jünemann, A. Piga, S.-J. Ran, M. Lewenstein, M. Rizzi, and A. Bermudez, Exploring interacting topological insulators with ultracold Atoms: The synthetic Creutz-Hubbard model, Phys. Rev. X 7, 031057 (2017).
  * (50) X. Zhou, J.-S. Pan, Z.-X. Liu, W. Zhang, W. Yi, G. Chen, and S. Jia, Symmetry-protected topological states for interacting fermions in alkaline-earth-like atoms, Phys. Rev. Lett. 119, 185701 (2017).
  * (51) S. de Léséleuc, V. Lienhard, P. Scholl, D. Barredo, S. Weber, N. Lang, H. P. Büchler, T. Lahaye, and A. Browaeys, Observation of a symmetry-protected topological phase of interacting bosons with Rydberg atoms, Science 365, 775 (2019).
  * (52) J. K. Freericks and H. Monien, Strong-coupling expansions for the pure and disordered Bose-Hubbard model, Phys. Rev. B 53, 2691 (1996).
  * (53) T. D. Kühner and H. Monien, Phases of the one-dimensional Bose-Hubbard model, Phys. Rev. B 58, R14741 (1998).
  * (54) T. D. Kühner, S. R. White, and H. Monien, One-dimensional Bose-Hubbard model with nearest-neighbor interaction, Phys. Rev. B 61, 12474 (2000).
  * (55) M. Rizzi, D. Rossini, G. De Chiara, S. Montangero, and R. Fazio, Phase diagram of spin-1 bosons on one-dimensional lattices, Phys. Rev. Lett. 95, 240404 (2005).


FE unequal | 0.479 | 0.406 | 0.340 | 0.303
RE unequal | NA | 0.406 | 0.336 | 0.311
random | 0.566 | 0.452 | 0.391 | 0.351
RE model under the equal prior
Design | Sample Size ($n$)
| 6 | 12 | 24 | 48
FE equal | NA | 0.506 | 0.412 | 0.367
RE equal | NA | 0.514 | 0.425 | 0.371
RE unequal | NA | 0.548 | 0.434 | 0.400
random | NA | 0.651 | 0.544 | 0.456
RE model under the unequal prior
Design | Sample Size ($n$)
| 6 | 12 | 24 | 48
FE unequal | NA | 0.415 | 0.346 | 0.312
RE equal | NA | 0.433 | 0.348 | 0.302
RE unequal | NA | 0.411 | 0.346 | 0.314
random | NA | 0.508 | 0.407 | 0.357

It is also of interest to see how well the optimal designs found under the
tree classification approach perform in terms of posterior model
probabilities. We conduct a simulation study under the FE structure using
either the equal or the unequal prior. For each design we want to assess, we
simulate 800 datasets from the joint distribution $p(\bm{y},m|\bm{d})$. For
each of the 800 datasets, we approximate the posterior model probability of
the model $m$ that generates the dataset $\bm{y}$ using IS with 100K prior
simulations. As for the classification results in Table 10, we are also
interested in the performance of optimal designs found under some wrongly
assumed scenarios. We also consider a ‘random’ setup where we select designs
randomly for each of the 800 datasets. Figure 10 shows the boxplots of the
estimated posterior model probabilities of the correct model for some of the
designs of interest when the true scenario is the FE structure with the equal
prior. The resulting boxplots when the true scenario is the FE structure with
the unequal prior are shown in Figure 11. It is again evident that the optimal
designs found are robust under the choice of the structure (FE or RE) and the
choice of the prior model probabilities (equal or unequal). We do not perform
a simulation study under the RE structure given the increasing difficulty of
estimating the posterior model probabilities under this structure.

It is important to note that the random forest-based validation results in
Table 10 were obtained in a small fraction of the time that it took to conduct
the simulation study used to produce the results in Figures 10 and 11.

(a) $n=6$

(b) $n=12$

(c) $n=24$

(d) $n=48$

Figure 10: Estimated posterior model probabilities for the correct model by
the validation study under the equal prior. Results based on sample sizes of
(a) $n=6$, (b) $n=12$, (c) $n=24$ and (d) $n=48$. Several designs are
considered: optimal design found under the correct (equal) prior, optimal
design found under the wrong (unequal) prior, optimal design found under the
wrong (RE) structure (no results for $n=6$) and randomly selected designs.

(a) $n=6$

(b) $n=12$

(c) $n=24$

(d) $n=48$

Figure 11: Estimated posterior model probabilities for the correct model by
the validation study under the unequal prior. Results based on sample sizes of
(a) $n=6$, (b) $n=12$, (c) $n=24$ and (d) $n=48$. Several designs are
considered: optimal design found under the correct (unequal) prior, optimal
design found under the wrong (equal) prior, optimal design found under the
wrong (RE) structure (no results for $n=6$) and randomly selected designs.

Figures 12 and 13 show misclassification matrices for the logistic regression
models under the FE structure for the equal and unequal priors, respectively.
To produce the results, 10K simulations from each model are used to train a
random forest with 100 trees. The misclassification matrices are then computed
based on a fresh test dataset of size $16\times 10K=160K$ with model
proportions in the dataset corresponding to the prior model probabilities
(under the equal prior, 10K simulations are taken from each model). The
improvement in classification accuracy is clear as the sample size is
increased, less so when increasing from $n=24$ to $n=48$. When the unequal
prior is selected, it is evident for small sample sizes that it is easier to
classify the models with higher prior probability. The misclassification
matrices for the RE structure are omitted because they are very similar.

(a) $n=6$

(b) $n=12$

(c) $n=24$

(d) $n=48$

Figure 12: Misclassification matrices obtained for the _FE structures_ of the
logistic regression example with the _equal prior_.

(a) $n=6$

(b) $n=12$

(c) $n=24$

(d) $n=48$

Figure 13: Misclassification matrices obtained for the _FE structures_ of the
logistic regression example with the _unequal prior_.

##  Appendix F Spatial Extremes Example

In this example, the goal is to place a fixed number of measuring sites in
space in order to maximise the ability to discriminate between different
spatial models for extreme outcomes (e.g., maximum annual temperatures). There
are many spatial models for extreme events, see Davison et al., (2012) for an
overview. For this example, we consider to discriminate between three
isotropic models: two _max-stable_ models and one _copula_ model.

###  F.1 Models

Max-stable processes are popular for modelling spatial extremes because they
are the only possible limits of renormalised pointwise maxima of infinitely
many independent copies of a stochastic process (de Haan and Ferreira,, 2006).
The advantage of working with the limiting process is that no knowledge about
the underlying true process is necessary. Inference for extreme outcomes based
on the true underlying process is fraught with high uncertainty and most often
not feasible because only the tails of the distribution are observed. If the
limiting assumption is (approximately) appropriate, it is much easier to model
the extreme data according to a max-stable process.

All the univariate marginal distributions of a max-stable process are members
of the family of _generalised extreme value (GEV)_ distributions. We assume
that all the univariate marginal distributions have a _unit Fréchet_
distribution ($\Pr\\{Y(\bm{x})\leq y\\}=\exp\\{-1/y\\},y>0$), so the focus is
on modelling the dependence structure of the process. The assumption of unit
Fréchet margins is not too restrictive from a practical perspective since a
simple transformation can be applied to the univariate margins to make them
unit Fréchet distributed, see Davison et al., (2012). The marginal parameters
needed for that transformation can be estimated in a separate step.
Alternatively, one may estimate the dependence and marginal parameters
together.

The _spectral representation_ of a max-stable process
$\\{Y(\bm{x}),\bm{x}\in\mathcal{X}\subseteq\mathbb{R}^{d}\\}$ with unit
Fréchet margins is given by

$\displaystyle Y(\bm{x})$ $\displaystyle=\max_{i\geq
1}\varphi_{i}(\bm{x}),\hskip 10.0pt\bm{x}\in\mathcal{X},$ (F.1)

where the _spectral functions_ $\varphi_{i}(\bm{x})=\zeta_{i}Z_{i}(\bm{x})$
are the products of the realisations $\\{\zeta_{i}\\}_{i=1}^{\infty}$ of a
Poisson point process on the positive real line with intensity
$d\Lambda(\zeta)=\zeta^{-2}d\zeta$ and of the independent realisations
$\\{Z_{i}(\bm{x}),\>\bm{x}\in\mathcal{X}\\}_{i=1}^{\infty}$ of a non-negative
stochastic process with continuous sample paths and
$\mathrm{E}[Z(\bm{x})]=1\;\forall\>\bm{x}\in\mathcal{X}$ (see, e.g., Ribatet,
(2013)).

Different max-stable processes are obtained by choosing different stochastic
processes $Z$. We consider two very popular stationary models, the _extremal-
$t$_ model (Opitz,, 2013) and the _Brown-Resnick_ model with power variogram
(Brown and Resnick,, 1977; Kabluchko et al.,, 2009). The specifications for
$Z_{i}(\bm{x})$ for each of the models are

$\displaystyle\text{Extremal-}t\text{:}\qquad$ $\displaystyle Z_{i}(\bm{x})$
$\displaystyle=\sqrt{\pi}\>2^{-(\nu-2)/2}\>\Gamma\left\\{(\nu+1)/2\right\\}^{-1}\>\max\\{0,\epsilon_{i}(\bm{x})\\}^{\nu},\;\nu>0,$
Brown-Resnick: $\displaystyle Z_{i}(\bm{x})$
$\displaystyle=\exp\left\\{\varepsilon_{i}(\bm{x})-\mathrm{Var}[\varepsilon_{i}(\bm{x})]/2\right\\},$

where $\epsilon_{i}$ and $\varepsilon_{i}$ are independent copies of Gaussian
processes.

In the case of the extremal-$t$ model, $\epsilon$ is a stationary Gaussian
process defined by the correlation function $\rho(h)$, where $h$ is the
Euclidean distance between two points. For our example, we assume the _powered
exponential_ or _stable_ correlation function:

$\displaystyle\rho(h)$
$\displaystyle=\exp\left[-\left(h/\lambda\right)^{\kappa}\right],\quad\lambda>0,\>0<\kappa\leq
2.$ (F.2)

The Brown-Resnick process is defined by its semi-variogram. If the process
$\varepsilon$ is a fractional Brownian motion centred at the origin, the
Brown-Resnick process is stationary and the semi-variogram has the form

$\displaystyle\gamma(h)$
$\displaystyle=\left(h/\lambda\right)^{\kappa},\quad\lambda>0,\>0<\kappa\leq
2,$

where $h$ denotes the distance between two locations.

Both models depend on two parameters governing the dependence between two
locations separated by a distance $h$: the _range_ parameter $\lambda$ and the
_smoothness_ parameter $\kappa$. In addition, the extremal-$t$ model has a
_degrees of freedom_ parameter denoted by $\nu$. We assume there is no
discontinuity of the correlation function at $h=0$ (i.e., no nugget effect).

The third model we consider is a copula model. Similar to the max-stable
models, the univariate marginal distributions of the copula model are unit
Fréchet. However, the extremal dependence between the locations is simply
modelled by a standard (non-extremal) copula. For an introduction to copulas
see Nelsen, (2006). We assume the multivariate Student-$t$ copula in our
example. The multivariate cumulative distribution function (CDF) at locations
$(\bm{x}_{1},\ldots,\bm{x}_{H})$ implied by the _non-extremal Student- $t$
copula_ model (Demarta and McNeil,, 2005) is

$\displaystyle\Pr\\{Y(\bm{x}_{1})\leq y_{1},\ldots,Y(\bm{x}_{H})\leq y_{H}\\}$
$\displaystyle=T_{H;\nu}\\{T_{1;\nu}^{-1}[F(y_{1})],\ldots,T_{1;\nu}^{-1}[F(y_{H})];\bm{\Sigma}\\},$

where $F(y)=\exp\\{-1/y\\}$ is the CDF of the unit Fréchet distribution,
$T_{1;\nu}^{-1}[\cdot]$ is the quantile function of the univariate Student-$t$
distribution with $\nu$ degrees of freedom, and
$T_{H;\nu}\\{\cdots;\bm{\Sigma}\\}$ is the CDF of the $H$-variate Student-$t$
distribution with $\nu$ degrees of freedom and dispersion matrix
$\bm{\Sigma}$. The diagonal elements of $\bm{\Sigma}$ are $1$ and the off-
diagonal elements contain the correlations between the locations. Therefore,
the entries of $\bm{\Sigma}$ are given by $\bm{\Sigma}_{ij}=\rho(h_{ij})$ for
$i,j=1,\ldots,H$, where $h_{ij}$ is the distance between locations $i$ and
$j$. As for the extremal-$t$ model, we assume the correlation function to be
the powered exponential correlation function (F.2). This also implies that the
non-extremal Student-$t$ copula model has the same set of parameters as the
extremal-$t$ model: range ($\lambda$), smoothness ($\kappa$), and degrees of
freedom ($\nu$).

###  F.2 Summary Statistics

If a reasonable amount of observations are collected at each location, the
data collected quickly becomes very high-dimensional, while each observation
is only marginally informative. This diminishes the classification power of
the classifiers we use. We therefore aim to reduce the dimension of the data
by generating informative summary statistics. Unfortunately, none of the
statistics we consider guarantee consistent model choice. This can potentially
result in large biases when estimating the posterior model probabilities
(Robert et al.,, 2011), which can also affect the estimates of the
misclassification error rates. However, trees and random forests work
reasonably well with a sizeable amount of moderately informative feature
variables. Therefore, we can include a wide variety of summary statistics,
where each contains some information about the process. Considering the
combined information of all the summary statistics, we expect that only a
small loss in information is incurred compared to the full dataset.

First, we include all the _F-madogram_ estimates for all the pairs of
locations. The F-madogram (Cooley et al.,, 2006) is similar to the semi-
variogram, but unlike the semi-variogram it also exists if the variances or
means of the random variables are not finite. Given $n$ observations
$\\{y_{1}(\bm{x}_{1}),\ldots,y_{n}(\bm{x}_{1})\\}$ and
$\\{y_{1}(\bm{x}_{2}),\ldots,y_{n}(\bm{x}_{2})\\}$ collected at locations
$\bm{x}_{1}$ as well as $\bm{x}_{2}$, the pairwise F-madogram between
locations $\bm{x}_{1}$ and $\bm{x}_{2}$ is estimated as

$\displaystyle\hat{\nu}_{F}(\bm{x}_{1},\,\bm{x}_{2})$
$\displaystyle=\frac{1}{2n}\sum_{i=1}^{n}|F\\{y_{i}(\bm{x}_{1})\\}-F\\{y_{i}(\bm{x}_{2})\\}|,$

where $F\\{y\\}=\exp\\{-1/y\\}$ is the CDF of the unit Fréchet distribution.

As a second set of summary statistics, we include estimates for all the
pairwise _extremal coefficients_ (Schlather and Tawn,, 2003). For a max-stable
process, the pairwise extremal coefficient between locations $\bm{x}_{1}$ and
$\bm{x}_{2}$ is defined as the value $\theta(\bm{x}_{1},\bm{x}_{2})$ for which

$\displaystyle\Pr(Y(\bm{x}_{1})\leq y,\>Y(\bm{x}_{2})\leq y)$
$\displaystyle=\Pr(Y(\bm{x}_{1})\leq
y)^{\theta(\bm{x}_{1},\bm{x}_{2})}=\exp\left(-\frac{\theta(\bm{x}_{1},\bm{x}_{2})}{y}\right).$
(F.3)

The pairwise extremal coefficient can assume values between $1$ and $2$. A
value of $\theta(\bm{x}_{1},\bm{x}_{2})=1$ indicates complete dependence
between the two locations. If $\theta(\bm{x}_{1},\bm{x}_{2})=2$, the two
locations are completely independent. We estimate it using the fast estimator
of Coles et al., (1999),

$\displaystyle\hat{\theta}(\bm{x}_{1},\bm{x}_{2})=\frac{n}{\sum_{i=1}^{n}{1/\max\\{y_{i}(\bm{x}_{1}),y_{i}(\bm{x}_{2})\\}}}.$
(F.4)

The extremal coefficient as defined by (F.3) only exists for max-stable
processes. In general, the coefficient also depends on the level $y$. However,
the quantities computed by Equation (F.4) might still provide useful
information about the dependence structure. For the $t$ copula model, Lee et
al., (2018) demonstrate by simulation that the estimates given by (F.4) are
indeed informative about the dependence structure.

The last set of summary statistics we consider is the set of _Kendall’s
$\tau$_ estimates between all pairs of locations. Kendall’s $\tau$ between
locations $\bm{x}_{1}$ and $\bm{x}_{2}$ is estimated by

$\displaystyle\hat{\tau}(\bm{x}_{1},\,\bm{x}_{2})$
$\displaystyle=\frac{2}{n\,(n-1)}\sum_{1\leq i<j\leq
n}\mathrm{sign}[y_{i}(\bm{x}_{1})-y_{j}(\bm{x}_{1})]\>\mathrm{sign}[y_{i}(\bm{x}_{2})-y_{j}(\bm{x}_{2})].$

Dombry et al., 2017b  show that for max-stable processes Kendall’s $\tau$ is
equal to the probability that the maxima at two locations occur concurrently
and are therefore attained for the same extremal function, so

$\displaystyle\tau(\bm{x}_{1},\,\bm{x}_{2})$
$\displaystyle=\Pr\left(\underset{i\geq
1}{\arg\max}\>\varphi_{i}(\bm{x}_{1})=\underset{i\geq
1}{\arg\max}\>\varphi_{i}(\bm{x}_{2})\right).$

All of the summary statistics we incorporate are also considered by Lee et
al., (2018), who perform ABC model selection using the summary statistic
projection method of Prangle et al., (2014) for a very similar set of models
as in this example. Therefore, a more detailed discussion of the summary
statistics can be found in Lee et al., (2018).

###  F.3 Bayesian Inference for Spatial Extremes Models

The likelihood functions of max-stable models are practically intractable for
most models for dimensions greater than two or three. Composite likelihood
methods have been the most popular way to conduct classical inference for max-
stable models, so model discrimination is usually based on the composite
likelihood information criterion (CLIC) (Padoan et al.,, 2010).

The observed extrema at several locations might occur at the same time, which
means that the extrema at these locations arise from the same extremal
function $\varphi_{i}(\bm{x})$ in Equation (F.1). The locations can then be
partitioned according to which extremal functions $\varphi_{i}(\bm{x})$
produce the extreme observations at the different locations. Stephenson and
Tawn, (2005) show that the joint likelihood of the extreme observations and
the partitions is substantially simpler than the likelihood of the extreme
observations without knowledge of the partitions. Thibaud et al., (2016) and
Dombry et al., 2017a  use this property to devise a Gibbs sampler with the
partitions as auxiliary variables to conduct Bayesian inference for max-stable
models. However, even the augmented likelihoods are expensive to evaluate for
the Brown-Resnick and extremal-$t$ model because they include multivariate
Gaussian (Brown-Resnick) and Student-$t$ (extremal-$t$) CDFs.

Due to the intractability of the likelihoods, ABC has also been a popular
method for Bayesian inference of max-stable models, see, e.g., Erhardt and
Smith, (2012) or the overview in Erhardt and Sisson, (2015). Lee et al.,
(2018) present an ABC application with the joint goal of model selection and
parameter estimation for the same set of models we consider. Hainy et al.,
(2016) seek to find optimal designs for parameter estimation for the
extremal-$t$ model with $\nu=1$ (called the ‘Schlather model’). They use ABC
to estimate the posterior variances, which they use as design criterion. Their
design algorithm is confined to very low-dimensional design spaces in order to
be able to store the reference table for all possible designs. They
sequentially select the best single location among a small set of possible
locations. With our classification approach, we are able to overcome these
limitations.

###  F.4 Settings and Results

In our example, we want to select $H$ ($H=3,\ldots,8$) locations on a regular
grid such that the ability to discriminate between the three models as
measured by the misclassification error rate is optimised. We search the $H$
optimal design points over a regular $6\times 6$ grid laid over a square with
edge length $10$. The data consist of $n=10$ independent realisations of the
process collected at all the locations. Due to the isotropic nature of the
processes, there are potentially many equivalent optimal solutions. With our
modification of the coordinate exchange algorithm using 20 random starts, we
seek to find one of these designs or at least a nearly optimal design.

We assume the following prior distributions:

$\displaystyle\log\left(\lambda\right)$ $\displaystyle\sim\mathcal{N}(1,\,4),$
$\displaystyle\kappa$ $\displaystyle\sim\mathcal{U}(0,\,2),$
$\displaystyle\log(\nu)$ $\displaystyle\sim\mathcal{N}(0,\,1)\text{ truncated
on }[-2.5,\,2.5].$

Furthermore, we assume equal prior model probabilities (= 1/3) for all models.

Simulating from the $t$ copula model is straightforward. It only requires
simulating from a multivariate $t$ distribution and then transforming the
margins with respect to the univariate $t$ CDF followed by the inverse unit
Fréchet CDF. For simulating from the max-stable models, we use the exact
simulation algorithm via extremal functions of Dombry et al., (2016).

During the design phase, we use classification trees evaluated on the training
set as well as random forests with $500$ trees using out-of-bag class
predictions to estimate the misclassification error rates. The training sets
for both methods contain 5K simulations per model. The optimal designs
obtained for these two methods are shown in Figures 14 (trees) and 15 (random
forests).

(a) 3 design points (b) 4 design points (c) 5 design points (d) 6 design
points (e) 7 design points (f) 8 design points

Figure 14: Optimal classification designs found using _trees_ for design sizes
from three to eight for the spatial extremes example. Selected design points
are marked by red triangles.

(a) 3 design points (b) 4 design points (c) 5 design points (d) 6 design
points (e) 7 design points (f) 8 design points

Figure 15: Optimal classification designs found using _random forests_ for
design sizes from three to eight for the spatial extremes example. Selected
design points are marked by red triangles.

To evaluate the designs found by our classification approach, we repeat
estimating the misclassification error rate via random forests with $500$
trees using out-of-bag class predictions on $100$ different simulated datasets
of size 15K (5K simulations per model). The distributions of the estimated
misclassification error rates are plotted in Figure 16. We also include the
distributions of the estimated misclassification error rates for $100$
simulated datasets generated on $100$ randomly selected designs. The optimal
classification designs found using random forests clearly perform best for all
design sizes. Using classification trees evaluated on the training set instead
of random forests leads to designs which are clearly worse (except for $H=7$).
Therefore, not considering the underestimation of the misclassification error
rate due to overfitting results in suboptimal designs with respect to the
misclassification error rate of the random forest classifier. However, the
average misclassification error rate of the classification tree designs is
still smaller than the average error rate of the random designs up until $7$
design points.

Figure 16: Spatial extremes example: distributions of the random forest-
estimated misclassification error rates over 100 random datasets of size 15K
generated at the optimal classification designs found using random forests
(rf) or trees (tr) for design sizes from three to eight. The distributions of
the misclassification error rates over 100 random datasets of size 15K
generated at 100 random designs (rd) are also shown for the same design sizes.

In addition to the misclassification error rate, we also compute the
misclassification matrix yielded by the random forest classifier for each of
the $100$ simulated datasets for each evaluated design. The average
misclassification matrices over the $100$ datasets are depicted in Figure 17
for the optimal designs obtained by the random forest classification approach.
They show that discriminating between the two max-stable models is more
difficult than discriminating between the $t$ copula model and either of the
max-stable models.

(a) 3 design points (b) 4 design points (c) 5 design points (d) 6 design
points (e) 7 design points (f) 8 design points

Figure 17: Average misclassification matrices over 100 simulated datasets
obtained for the _random forest classification designs_ for the spatial
extremes example. Design sizes from 3 – 8 design points are considered.


#  A Thorough Review on Recent Deep Learning Methodologies for Image
Captioning

Ahmed Elhagry
<EMAIL_ADDRESS>Karima Kadaoui
<EMAIL_ADDRESS>
###### Abstract

Image Captioning is a task that combines computer vision and natural language
processing, where it aims to generate descriptive legends for images. It is a
two-fold process relying on accurate image understanding and correct language
understanding both syntactically and semantically. It is becoming increasingly
difficult to keep up with the latest research and findings in the field of
image captioning due to the growing amount of knowledge available on the
topic. There is not, however, enough coverage of those findings in the
available review papers. We perform in this paper a run-through of the current
techniques, datasets, benchmarks and evaluation metrics used in image
captioning. The current research on the field is mostly focused on deep
learning-based methods, where attention mechanisms along with deep
reinforcement and adversarial learning appear to be in the forefront of this
research topic. In this paper, we review recent methodologies such as UpDown,
OSCAR, VIVO, Meta Learning and a model that uses conditional generative
adversarial nets. Although the GAN-based model achieves the highest score,
UpDown represents an important basis for image captioning and OSCAR and VIVO
are more useful as they use novel object captioning. This review paper serves
as a roadmap for researchers to keep up to date with the latest contributions
made in the field of image caption generation.

Keywords: Image Captioning, Computer Vision, NLP, Reinforcement Learning,
Attention Mechanism, GANs

##  1 Introduction

Image captioning has a variety of applications. From automatic image indexing
to assistive technology, there is more and more incentive to build solutions
that don’t rely on human annotation since it quickly becomes cost inefficient.
A lot of work has been done in the field already and some image captioning
survey papers are available, but none of them cover some important models that
came out between 2018 and 2020. In a field that evolves at this rate, it is
imperative to keep track of the latest models and frameworks. That is why
we’re presenting this paper, in an attempt to cover some of the most relevant
papers that have recently been published and compare the performances of their
models and describing how they work.

One of the techniques that have a crucial role in image captioning today is
the use of attention mechanisms. Ever since transformers [26] were introduced,
a number of different tasks such as machine translation and language modeling
saw significant improvements thanks to them. Image captioning is no different:
there is an extensive use of top-down visual attention in different models as
will be presented in the paper. Another technique that sparked researchers’
interest is deep reinforcement learning. It has shown to be particularly good
with unusual images (e.g. bed in a forest). The way it’s used in image
captioning is through optimizing the reward function by maximizing its
expected value, which cannot be done using MLE since metrics are not
differentiable.

###  1.1 Related Work

As mentioned in the review paper [9], the authors presented a comprehensive
review of the state-of-the-art deep learning-based image captioning techniques
by late 2018. The paper gave a taxonomy of the existing techniques, compared
the pros and cons, and handled the research topic from different aspects
including learning type, architecture, number of captions, language models and
feature mapping. They also discussed both the strengths and weaknesses of
different datasets and evaluation metrics. According to another review paper
[25] published later mid 2019, the authors compared the image captioning
methodologies from 2016 to 2019 including newer ones on two datasets: MSCOCO
and Flickr30k. An investigation was done on different feature extractors
including AlexNet, VGG-16 Net, ResNet, GoogleNet with all the nine Inception
models, and DenseNet. In addition, language models were covered such as LSTM,
RNN, CNN, GRU and TPGN. This is while comparing several evaluation metrics
including BLEU (1 to 4), CIDEr and METEOR. In [18], the paper put the
spotlight on some of the advancement on the image captioning task until early
2020, where various approaches were discussed including N-cut, color-based
segmentation and hybrid engine. It also discussed how model engineering and
incorporating more hyper-parameters improve the overall pipeline and result in
the best accuracy for such models. In the same year, another study [6] covered
a review about the literature from 2017 to 2019, where they discussed
different datasets and architectures. They stated that the CNN-LSTM models
outperformed the CNN-RNN ones, and the most evaluation metric used was BLEU (1
to 4). They also found that the best methods for such a model implementation
are encode-decoder and attention mechanism. Furthermore, they mentioned that a
combination of both methods can help in improving the results in such a task.
Image captioning remains an active research area, and new methodologies keep
being published up until this moment. That was one of the main motivations to
write this review paper in order to tackle all the recent advances in the past
few years including 2020.

##  2 Methods

###  2.1 UpDown

Most common mechanisms that rely on visual attention today are of the top-down
kind. They are fed their partially finished caption at each time step to gain
context. The issue with these models however is that there is no deliberation
as to what regions of an image will receive attention. This has an effect on
the quality of the captions as focusing on salient object regions will provide
descriptions that are similar to ones given by humans [23].

[4] introduce Up-Down, a model that joins an entirely visual bottom-up
mechanism and a task-specific context top-down one. The former gives proposals
on image regions that it deems salient, while the latter uses context to
compute an attention distribution over them, thus allowing attention to be
directed to the important objects in the input image.

Implementation Details The bottom-up mechanism employs the Faster R-CNN [20]
object detection model, responsible for recognizing class objects and
surrounding them with bounding boxes. For pre-training, it is initialized with
Resnet-101 [8] and trained on the Visual Genome dataset [14]. The top-down
mechanism uses a visual attention LSTM and a language one. The attention LSTM
is fed the previous language LSTM outputs, the word generated at time t-1 and
mean-pooled image features to decide which regions should receive attention.
The generated caption up until that point is then used to compute the
conditional distribution over potential output words, and the product of all
of the conditional distributions gives the distribution over the complete
captions.

###  2.2 OSCAR

Vision-language pre-training (VLP) is widely used for learning cross-modal
representations. It suffers, however, from 2 issues [16]: a difficulty in
discerning features due to the overlap of their image regions and a lack of
alignment caption words and their corresponding image regions. [16] remedy to
this by using object tags as ”anchor points”. More specifically, they use
triples as inputs composed of image region features, object tags and word
sequence (caption). Doing this helps, because when one channel is incomplete
or noisy, the other might complete the information (you can describe an object
both through image and language). It is therefore simple to make the
alignments, because the most important elements in the image appear in the
matching caption and are also the ones that are expected to receive the most
attention.

Figure 1: The semantic space used by OSCAR [16]. In the example of a dog
sitting on a couch, ”couch” and ”dog” are close in region features since
they’re roughly in the same area of the image, but they farther apart in word
embeddings because of their different meanings.

Implementation Details OSCAR detects object tags using Faster R-CNN [20] and
presents a 2-view perspective:
(1) A dictionary-view with a linguistic semantic space encompassing the tags
and caption tokens, and a visual semantic space where image regions lay (Fig.
1).
(2) A modality-view that consists of an image modality containing image
features and tags, and a language modality with caption tokens. The total loss
of pretraining is defined by the addition of a masked token loss of predicting
masked tokens form the linguistic semantic space and the contrastive loss of
predicting if an image-tag sequence is polluted (from replaced tags).

At inference time, the input consists of image regions and tags At each time
step of the generation, a [MASK] token is appended to the sequence before
being replaced by a token from the vocabulary until the [STOP] token is
generated.

###  2.3 VIVO

In the nocaps [1] challenge, the only allowed image-caption dataset is the MS
COCO [17] one, making conventional VLP methods inapplicable [10]. For that
reason, [10] came up with VIVO (VIsual VOcabulary pre-training) What it does
differently is defining a ”visual vocabulary”, which is a joint embedding
space of tags and image region features, where vectors of semantically close
objects (e.g. accordion and instrument) are close to each other (Fig. 2).
After pre-training the vocabulary, the model is fine-tuned with image-caption
pairs using the MS COCO dataset [17]. The key difference between VIVO and
other VLP models is that VIVO is only pre-trained on image-tag pairs and no
captions are involved before fine-tuning. This can prove very useful since
tags are easier to be generated automatically which allows the use of a huge
number of them for no annotation costs.

Figure 2: Visual Vocabulary used by VIVO. [10] Objects that are similar
semantically are closer together. o represents regions and + represents tags.
Yellow objects and tags are novel.

Implementation Details VIVO uses a multi-layer Transformer responsible for
aligning tags with their corresponding image region features followed by a
linear layer and softmax. During pre-training, image region features are
extracted from the input image using Updown’s object detector [4] and fed to
the Transformer along with a set of pairs of images and tags. One or more tags
are randomly masked and the model makes predictions based on the remaining
tags and the image regions.

In fine-tuning, the model is fed a triplet of image regions, tags and a
caption where some of the caption’s tokens are randomly masked and the model
learns to predict them. A uni-directional attention mask is applied and the
parameters are optimized using a cross-entropy loss.

At inference time, image region features are extracted from the input image
and tags are detected. A caption is then generated one token at a time in an
auto-regressive manner (using the previous tokens as input) until the end
token is generated or the caption reached its maximum length.

###  2.4 Meta Learning

One of the drawbacks of reinforcement learning is the reward hacking problem,
which, in other words, is overfitting on the reward function which occurs when
the agent finds a way to maximize the score without generating captions of a
better quality. When using a CIDEr optimization [21] for example, common
phrases are given less weight and punishment is given to a caption that is too
short. As a result, when a short caption is generated, common phrases are
added to it to make it longer, ending up with unnatural sentence endings such
as ”a little girl holding a cat in a of a.” [15]

[15] introduce meta learning, which is learning a meta model that is able to
optimize and adapt to several different tasks [7]. In this case, the model
simultaneously optimizes the reward function (reinforcement task) and uses
supervision from the ground truth (supervision task) by taking gradient steps
in both directions. This guarantees the distinctiveness of the captions and
their propositional correctness and results in sound human-like sentences.

Additionally, they import the SPICE [3] metric and add it to the CIDEr [27]
reward term since it performs semantic propositional evaluation using a scene
graph. What this means is that an unusual caption ending will be an object-
relation pair in the scene graph without a match (Fig. 3). Unfortunately,
SPICE has a reward hacking issue of its own since it allows duplicate tuples.
It is therefore not easy to develop an ideal evaluation metric.

Figure 3: SPICE scene graph for the caption ”A young girl standing on top of a
tennis court”. The objects are marked red, the relations blue and the
attributes green. [3]

Implementation Details [15] use the UpDown architecture [4] as outlined above.
The two tasks that the model needs to optimize are the maximum likelihood
estimate task and the reinforcement learning task. In other words, it needs to
take 2 gradient steps in order to update the parameter $\theta$. In the first
step, the model adapts $\theta$ to the two tasks and calculates their
respective losses. Then, $\theta$ updates itself in what is named a ”meta
update”. Doing things this way, the model learns the parameter $\theta$ that
optimizes both tasks instead of simply taking a step in between the two
gradients when adding up the losses.

###  2.5 Conditional GAN-Based Model

To overcome reward hacking, [5] use discriminator networks to decide whether a
generated caption is from a human or a machine. Since they do not give their
model a name, we will call it IC-GAN (Image Captioning GAN) for the sake of
practicality.

Implementation Details [5] experimented with 2 different architectures for the
discriminator: one of them using a CNN with a fully connected layer and a
sigmoid transformation, and the other an RNN (LSTM) with a fully connected
layer and a softmax. They also experiment with an ensemble of 4 CNNs and 4
RNNs. For the generator, a number of different architectures were used but we
will focus in the results on the generator that uses the UpDown architecture
[4]. In all cases, the generator and discriminator need to be pre-trained
before being alternatively fine-tuned.

###  2.6 Evaluation Metrics

To compare the quality of the generated captions to the ground-truth, a number
of evaluation metrics are used. The most common used ones being CIDEr, SPICE,
BLEU and METEOR. The common metrics across all the covered literature are
CIDEr and SPICE, which is why we will be using them. CIDEr [27] is an image
classification evaluation metric that uses term frequency-inverse document
frequency (TF-IDF) [22] to achieve human consensus. SPICE [3] is a new
semantic concept-based caption assessment metric based on scene-graph, a
graph-based semantic representation (Fig. 3) [11][24].

###  2.7 Benchmarks

[1] develop a benchmark called nocaps, which, in addition to the image
captioning dataset Microsoft COCO Captions [17], makes use of the Open Images
object detection dataset [13] to introduce novel objects not seen in the
former. The nocaps benchmark is made up of 166,100 captions that describe
15,100 images from the OpenImages validation and test sets. OSCAR [16] and
VIVO [10] methods are evaluated on the nocaps validation set [1]. Karpathy
splits [12] are used in the evaluation of the meta learning model [15] and the
IC-GAN one [5]. Finally, UpDown [4] is evaluated on both benchmarks.

##  3 Results

Method |  CIDEr |  SPICE
---|---|---
Resnet Baseline |  111.1 |  20.2
UpDown |  120.1 |  21.4
MLE Maximization |  110.2 |  20.3
RL Maximization |  120.4 |  21.3
MLE + RL Maximization |  119.3 |  21.2
Meta Learning |  121.0 |  21.7
IC-GAN (Updown/CNN-GAN) |  123.2 |  22.1
IC-GAN (Updown/RNN-GAN) |  122.2 |  22.0
IC-GAN (Updown/ensemble) |  125.9 |  22.3
Table 1: Results of the overall performance on MS COCO Karpathy test split
[15] [5]. Methods with a * symbol use reinforcement learning with CIDEr
optimization.

Method | in-domain | near-domain | out-of-domain | overall
---|---|---|---|---
CIDEr |  SPICE |  CIDEr |  SPICE |  CIDEr |  SPICE |  CIDEr |  SPICE
Validation Set
UpDown (2019) |  78.1 |  11.6 |  57.7 |  10.3 |  31.3 |  8.3 |  55.3 |  10.1
UpDown + CBS |  80.0 |  12.0 |  73.6 |  11.3 |  66.4 |  9.7 |  73.1 |  11.1
UpDown + ELMo + CBS |  79.3 |  12.4 |  73.8 |  11.4 |  71.7 |  9.9 |  74.3 |  11.2
OSCAR (2020) |  79.6 |  12.3 |  66.1 |  11.5 |  45.3 |  9.7 |  63.8 |  11.2
OSCAR + CBS |  80.0 |  12.1 |  80.4 |  12.2 |  75.3 |  10.6 |  79.3 |  11.9
OSCAR+SCST+CBS |  83.4 |  12.0 |  81.6 |  12.0 |  77.6 |  10.6 |  81.1 |  11.7
VIVO (2020) |  88.8 |  12.9 |  83.2 |  12.6 |  71.1 |  10.6 |  81.5 |  12.2
VIVO + CBS |  90.4 |  13.0 |  84.9 |  12.5 |  83.0 |  10.7 |  85.3 |  12.2
VIVO+SCST+CBS |  92.2 |  12.9 |  87.8 |  12.6 |  87.5 |  11.5 |  88.3 |  12.4
Human |  84.4 |  14.3 |  85.0 |  14.3 |  95.7 |  14.0 |  87.1 |  14.2

Table 2: Evaluation on nocaps validation set [10]

###  3.1 MS COCO Karpathy Splits Benchmark

[4] initially ran experiments on both an ablated Resnet [8] baseline model and
the UpDown model to measure the impact of the bottom-up attention. Since the
meta learning model [15] and IC-GAN [5] use the UpDown architecture
themselves, we have decided to include the Resnet baseline to the comparison
(Table 1) to get an idea about the individual effect of the bottom-up + top-
down approach.

We notice that UpDown shows an important gain in performance going from 111.1
to 120.1 and from 20.2 to 20.4 in the CIDEr and SPICE metrics respectively.
This represents a relative improvement of 8% in CIDEr and 3% in SPICE.
Therefore, adding bottom-up attention has an important positive impact on
image captioning.

Concerning the experiments in [15], we observe that the model that uses meta
learning receives a CIDEr score of 121.0 and a SPICE score of 21.7. It is the
most performant on both evaluation metrics compared to maximizing using the
maximum likelihood estimate, reinforcement learning and the MLE+RL
maximization (which relies on simply adding up the gradients from the
supervision and reinforcement tasks). It also has a slight improvement on the
UpDown model without the use of meta learning.

IC-GAN on the other hand, shows the highest performance with a significant
improvement on all 3 models. The relative improvements compared to the UpDown
model range between 1.7% (RNN-GAN) and 4.6% (ensemble). Additionally,
comparing with conventional reinforcement learning approaches, the proposed
adversarial learning method boosts the performance from 8.9% to 16.3% [5].

It is important to note that although using a CNN-GAN slighlty improves the
score compared to using an RNN-GAN, the latter can save up to 30% of training
time compared to the former [5].

Evaluation scores aside, IC-GAN is able to generate human-like captions and
avoids mistakes commonly made by tradition reinforcement learning methods,
such as duplicated words and logical errors like ”a group of people standing
on top of a clock” [5].

###  3.2 nocaps Benchmark

The OSCAR model [16] is characterized by being highly efficient when it comes
to parameters due to the anchor points making the semantic alignments learning
easier. When used on its own, it outperfoms the UpDown model on all in-domain,
near-domain and out-of-domain subsets. By adding Constrained Beam Search [2]
and Self Critical Sequential Training (SCST) [21] the performance improves
tremendously, particularly in the out-of-domain subset going from 45.3 to 77.6
in CIDEr and from 9.7 to 10.6 in SPICE.

OSCAR does not score however, as high as VIVO as shown in Table 2. In the in-
domain case, VIVO on its own outperforms the combinations OSCAR+SCST+CBS and
UpDown+ELMo+CBS [19] by a difference ranging from 5.4 to 9.5 in CIDEr and from
0.5 to 0.9 on SPICE.

The VIVO+SCST+CBS version shows the highest performance with CIDEr scores that
even surpass the human ones across in-domain and near-domain. The out-of-
domain subset results still show great scores that surpass all other models.

##  4 Discussion

It seems that the current research on image captioning is heavily focused on
deep-learning techniques, and for a good reason. Image captioning is a very
complex task that combines both computer vision and natural language
processing and it therefore needs powerful techniques that could handle that
level of complexity. Attention mechanisms along with deep reinforcement and
adversarial learning, appear to be actively-researched methods for this task,
as showcased in this paper. Faster R-CNN is a popular network choice, along
with the LSTM. In particular, the UpDown model seems to be used as a basis for
multiple papers that were published between 2018 and 2020 which gives it a key
role and notable impact on the advances in the field. Novel object captioning
also seems to be gathering a lot of interest after proving its usefulness.

Although the VIVO and OSCAR models don’t show scores as high as the ones that
use meta learning and adversarial learning, they are superior in terms of
their usability ”in the wild”, since the MS COCO dataset [17] on which all of
these models are trained contains only a small part of the objects that we run
into in real life. We should therefore not take evaluation scores at face
value, especially after demonstrating the reward hacking problem in
reinforcement learning. In the future, additional efforts should be put into
making more robust reward functions and into more research on novel object
captioning, especially replacing human-annotated object detection datasets
with fully machine-generated tags.

##  5 Conclusions

Image captioning is an active research topic that creates the space for
competition among researchers. It is noticed that the existing review papers
do not cover some of the important recent advances, although there are new
methodologies with outperforming performance. The outgoing research on image
captioning is focused on deep learning-based methods, where attention
mechanisms exist along with deep reinforcement and adversarial learning. Our
paper discusses the recent methods and their implementations. State-of-the-art
techniques include Updown, OSCAR, VIVO, Meta Learning and a GAN-based model.
The GAN-based model is the most performant, UpDown has the most impact and
OSCAR and VIVO are the more useful. We hope this review will provide the
community with a complementary guideline along with the existing review papers
for further pursued research in the image captioning research topic.

## References

  * [1] H. Agrawal, K. Desai, Y. Wang, X. Chen, R. Jain, M. Johnson, D. Batra, D. Parikh, S. Lee, and P. Anderson.  nocaps: novel object captioning at scale.  CoRR, abs/1812.08658, 2018.
  * [2] P. Anderson, B. Fernando, M. Johnson, and S. Gould.  Guided open vocabulary image captioning with constrained beam search.  CoRR, abs/1612.00576, 2016.
  * [3] P. Anderson, B. Fernando, M. Johnson, and S. Gould.  Spice: Semantic propositional image caption evaluation.  Computer Vision – ECCV 2016 Lecture Notes in Computer Science, page 382–398, 2016.
  * [4] P. Anderson, X. He, C. Buehler, D. Teney, M. Johnson, S. Gould, and L. Zhang.  Bottom-up and top-down attention for image captioning and VQA.  CoRR, abs/1707.07998, 2017.
  * [5] C. Chen, S. Mu, W. Xiao, Z. Ye, L. Wu, and Q. Ju.  Improving image captioning with conditional generative adversarial nets.  Proceedings of the AAAI Conference on Artificial Intelligence, 33:8142–8150, 2019.
  * [6] M. Chohan, A. Khan, M. Saleem, S. Hassan, A. Ghafoor, and M. Khan.  Image captioning using deep learning: A systematic literature review.  International Journal of Advanced Computer Science and Applications, 11(5), 2020.
  * [7] C. Finn, P. Abbeel, and S. Levine.  Model-agnostic meta-learning for fast adaptation of deep networks.  In International Conference on Machine Learning, pages 1126–1135. PMLR, 2017.
  * [8] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  CoRR, abs/1512.03385, 2015.
  * [9] M. Z. Hossain, F. Sohel, M. F. Shiratuddin, and H. Laga.  A comprehensive survey of deep learning for image captioning, 2018.
  * [10] X. Hu, X. Yin, K. Lin, L. Wang, L. Zhang, J. Gao, and Z. Liu.  Vivo: Visual vocabulary pre-training for novel object captioning, 2021\.
  * [11] J. Johnson, R. Krishna, M. Stark, L.-J. Li, D. A. Shamma, M. S. Bernstein, and L. Fei-Fei.  Image retrieval using scene graphs.  2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2015.
  * [12] A. Karpathy and F. Li.  Deep visual-semantic alignments for generating image descriptions.  CoRR, abs/1412.2306, 2014.
  * [13] I. Krasin, T. Duerig, N. Alldrin, V. Ferrari, S. Abu-El-Haija, A. Kuznetsova, H. Rom, J. Uijlings, S. Popov, A. Veit, S. Belongie, V. Gomes, A. Gupta, C. Sun, G. Chechik, D. Cai, Z. Feng, D. Narayanan, and K. Murphy.  Openimages: A public dataset for large-scale multi-label and multi-class image classification.  Dataset available from https://github.com/openimages, 2017.
  * [14] R. Krishna, Y. Zhu, O. Groth, J. Johnson, K. Hata, J. Kravitz, S. Chen, Y. Kalantidis, L. Li, D. A. Shamma, M. S. Bernstein, and F. Li.  Visual genome: Connecting language and vision using crowdsourced dense image annotations.  CoRR, abs/1602.07332, 2016.
  * [15] N. Li, Z. Chen, and S. Liu.  Meta learning for image captioning.  Proceedings of the AAAI Conference on Artificial Intelligence, 33:8626–8633, 2019.
  * [16] X. Li, X. Yin, C. Li, P. Zhang, X. Hu, L. Zhang, L. Wang, H. Hu, L. Dong, F. Wei, et al.  Oscar: Object-semantics aligned pre-training for vision-language tasks.  In European Conference on Computer Vision, pages 121–137. Springer, 2020.
  * [17] T. Lin, M. Maire, S. J. Belongie, L. D. Bourdev, R. B. Girshick, J. Hays, P. Perona, D. Ramanan, P. Dollár, and C. L. Zitnick.  Microsoft COCO: common objects in context.  CoRR, abs/1405.0312, 2014.
  * [18] K. C. Nithya and V. V. Kumar.  A review on automatic image captioning techniques.  In 2020 International Conference on Communication and Signal Processing (ICCSP), pages 0432–0437, 2020.
  * [19] M. E. Peters, M. Neumann, M. Iyyer, M. Gardner, C. Clark, K. Lee, and L. Zettlemoyer.  Deep contextualized word representations.  CoRR, abs/1802.05365, 2018.
  * [20] S. Ren, K. He, R. B. Girshick, and J. Sun.  Faster R-CNN: towards real-time object detection with region proposal networks.  CoRR, abs/1506.01497, 2015.
  * [21] S. J. Rennie, E. Marcheret, Y. Mroueh, J. Ross, and V. Goel.  Self-critical sequence training for image captioning.  CoRR, abs/1612.00563, 2016.
  * [22] S. Robertson.  Understanding inverse document frequency: on theoretical arguments for idf.  Journal of Documentation, 60(5):503–520, 2004.
  * [23] B. J. Scholl.  Objects and attention: The state of the art.  Cognition., 80(1).
  * [24] S. Schuster, R. Krishna, A. Chang, L. Fei-Fei, and C. D. Manning.  Generating semantically precise scene graphs from textual descriptions for improved image retrieval.  Proceedings of the Fourth Workshop on Vision and Language, 2015\.
  * [25] R. Staniūtė and D. Šešok.  A systematic literature review on image captioning.  Applied Sciences, 9(10):2024, 2019.
  * [26] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and I. Polosukhin.  Attention is all you need.  CoRR, abs/1706.03762, 2017.
  * [27] R. Vedantam, C. L. Zitnick, and D. Parikh.  Cider: Consensus-based image description evaluation.  2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2015.