# Stability of supercurrents in a superfluid phase of spin-1 bosons in an
optical lattice

Shion Yamashika1, Ryosuke Yoshii2, and Shunji Tsuchiya1 1Department of
Physics, Chuo University, 1-13-27 Kasuga, Tokyo 112-8551, Japan
2Center for Liberal Arts and Sciences, Sanyo-Onoda City University, 1-1-1
Daigaku-Dori, Sanyo-Onoda, Yamaguchi 756-0884, Japan

###### Abstract

We study collective modes and superfluidity of spin-1 bosons with
antiferromagnetic interactions in an optical lattice based on the time-
dependent Ginzburg-Landau (TDGL) equation derived from the spin-1 Bose-Hubbard
model. Specifically, we examine the stability of supercurrents in the polar
phase in the vicinity of the Mott insulating phase with even filling factors.
Solving the linearized TDGL equation, we obtain gapless spin-nematic modes and
gapful spin-wave modes in the polar phase that arise due to the breaking of
$S^{2}$ symmetry in spin space. Supercurrents exhibit dynamical instabilities
induced by growing collective modes. In contrast to the second-order phase
transition, the critical momentum for mass currents is finite at the phase
boundary of the first-order superfluid-Mott insulator (SF-MI) phase
transition. Furthermore, the critical momentum remains finite throughout the
metastable SF phase and approaches zero towards the phase boundary, at which
the metastable SF state disappears. We also study the stability of spin
currents motivated by recent experiments for spinor gases. The critical
momentum for spin currents is found to be zero, where a spin-nematic mode
causes the dynamical instability. We investigate the origin of the zero
critical momentum for spin currents and find it attributed to the fact that
the polar state becomes energetically unstable even in the presence of an
infinitesimal spin current. We discuss implications of the zero critical
momentum for spin currents for the stability of the polar state.

##  I Introduction

Ultracold atoms in an optical lattice have opened a versatile research field
that lies at the interface of condensed matter physics, statistical physics,
and atomic, molecular, and optical physics [1, 2]. Spinor bosons in an optical
lattice provide with an ideal platform for studying quantum magnetism, quantum
phase transitions, and non equilibrium quantum dynamics [3, 4, 5, 6, 7, 8, 9].
In particular, the superfluid-Mott insulator (SF-MI) phase transition of
spin-1 bosons in an optical lattice has attracted much attention. The rich
structure of the phase diagram that arises from the interplay between strong
correlations and spin degrees of freedom has been a main focus of intense
theoretical studies. It exhibits many interesting features including the
parity effect of the MI phase [10, 11, 12, 13], the first-order SF-MI
transition [14, 15, 16, 17], and the spin-nematic order in the MI phase [18,
19, 20, 21, 12, 22, 23, 24, 25, 26]. On the other hand, recent progress of
experimental studies have made possible to observe the signature of the first-
order phase transition [27] and quantum critical dynamics [28] in this system.

One of the remarkable features of bosonic superfluids in an optical lattice is
dynamical instabilities of supercurrents [29, 30, 31, 32, 33, 34, 35, 36, 37,
38, 39]. In fact, this phenomenon is considered to be deeply related with the
SF-MI phase transition. It has been predicted for spineless bosons that the
critical superfluid velocity approaches zero at the second-order SF-MI phase
transition reflecting the diverging healing length [39, 38] and later it was
experimentally confirmed [33]. Given the rich physics associated with the SF-
MI phase transition and the recent experimental progress, therefore, it is of
particular interest to investigate dynamical instabilities of spin-1 bosons in
an optical lattice. Recently, the critical momentum for mass currents and the
stability phase diagram have been calculated numerically using the dynamical
Gutzwiller approximation [40]. The critical momentum in the vicinity of the
first-order SF-MI phase transition, however, has not been fully investigated.
Meanwhile, critical dynamics of spin currents has been observed in an
antiferromagnetic spin-1 Bose-Einstein condensate [41]. In this experiment,
spin currents were induced by counter flow of two spin components. This
experiment inspires us to investigate the stability of spin currents in an
optical lattice.

In this paper, we study spin-1 bosons with antiferromagnetic interactions in
an optical lattice. Specifically, we examine the stability of supercurrents in
the polar phase in the vicinity of the MI phase with even filling factors
based on the time-dependent Ginzburg-Landau (TDGL) equation, which we derive
from the spin-1 Bose-Hubbard model. Our focus is on the superfluidity in the
vicinity of the first-order phase transition. We calculate the critical
momentum for mass currents and find that it has a finite value not only at the
first-order SF-MI phase boundary, but throughout the metastable SF phase, in
contrast to the second-order case. We also study the stability of spin
currents motivated by the recent experiments for spinor gases [41]. The
critical momentum for spin currents is found to be zero. We clarify the origin
of the instability of spin currents and briefly discuss about its implications
for the stability of the polar phase.

The organization of this paper is as follows: In Sec. II, we introduce the
spin-1 Bose-Hubbard model and calculate the metastable phase diagram within
the perturbative mean-field theory. In Sec. III, we study collective modes in
the SF phase using the TDGL equation. In Sec. IV, we study the stability of
supercurrents and calculate the critical momenta for mass and spin currents.
In Sec. V, we investigate the origin of the instability of spin currents and
discuss the stability of the polar phase. Finally, we summarize our results in
Sec. VI. Details of the calculations in the perturbative mean-field theory are
given in Appendix A and B. The derivations of the TDGL equation and the
conservation laws are given in Appendix C and D, respectively. We set
$\hbar=k_{\rm B}=1$ throughout this paper.

##  II spin-1 Bose-Hubbard model and SF-MI phase transition

We consider spin-1 bosons trapped in a $d$-dimensional cubic optical lattice
at zero temperature. In this paper, we neglect effect of a harmonic trapping
potential for simplicity. If the lattice potential is sufficiently deep, the
system is well described by the spin-1 Bose-Hubbard model [11, 19]:

$\displaystyle\hat{H}$ $\displaystyle=$ $\displaystyle-t\sum_{\langle
i,j\rangle,\alpha}(\hat{b}^{\dagger}_{i\alpha}\hat{b}_{j\alpha}+\mathrm{h.c.})-\mu\sum_{i}\hat{n}_{i}$
(1) $\displaystyle+$
$\displaystyle\frac{U_{0}}{2}\sum_{i}\hat{n}_{i}(\hat{n}_{i}-1)+\frac{U_{2}}{2}\sum_{i}(\hat{\bm{S}}_{i}^{2}-2\hat{n}_{i}),$

where $\hat{b}_{i\alpha}\,(\hat{b}_{i\alpha}^{\dagger})$ is the annihilation
(creation) operator for an atom in the hyperfine state
$|F=1,m=\alpha\rangle\,(\alpha=1,0,-1)$ at site $i$. Here, $\langle
i,j\rangle$ denotes a summation over nearest-neighbor sites, $t$ the hopping
matrix element, $\mu$ the chemical potential, and $U_{0}$ and $U_{2}$ the on-
site spin-independent and spin-dependent interactions, respectively. In this
paper, we assume $U_{2}>0$, i.e., an antiferromagnetic interaction.
$U_{2}/U_{0}$ is fixed for each atomic species. For example,
$U_{2}/U_{0}=0.04$ for ${}^{23}{\rm Na}$ [11]. In this paper, we set
$U_{2}/U_{0}$ optimal values for validity of the theory. The operators
$\hat{n}_{i}=\sum_{\alpha}\hat{b}_{i\alpha}^{\dagger}\hat{b}_{i\alpha}$, and
$\hat{\bm{S}}_{i}=\sum_{\alpha,\beta}\hat{b}_{i\alpha}^{\dagger}\bm{F}_{\alpha\beta}\hat{b}_{i\beta}$
represent the number of particles and the spin at site $i$, respectively.
Here, $\bm{F}$ denotes the spin-1 matrices

$\displaystyle F_{x}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}\left[\begin{array}[]{ccc}0&1&0\\\ 1&0&1\\\
0&1&0\end{array}\right],\
F_{y}=\frac{i}{\sqrt{2}}\left[\begin{array}[]{ccc}0&-1&0\\\ 1&0&-1\\\
0&1&0\end{array}\right],$ (8) $\displaystyle F_{z}$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{ccc}1&0&0\\\ 0&0&0\\\
0&0&-1\end{array}\right].$ (12)

In the following subsections, we briefly review the perturbative mean-field
theory developed in Ref. [11] and then extend it to describe the first-order
SF-MI phase transition as well as the metastable SF and MI phases.

###  II.1 Mott insulating phase

We first study the MI state in the limit of $t=0$, where the Hamiltonian (1)
can be written as

$\displaystyle\hat{H}$ $\displaystyle=$
$\displaystyle\sum_{i}\hat{H}^{0}_{i},$ (13) $\displaystyle\hat{H}^{0}_{i}$
$\displaystyle=$
$\displaystyle-\mu\hat{n}_{i}+\frac{U_{0}}{2}\hat{n}_{i}(\hat{n}_{i}-1)+\frac{U_{2}}{2}(\hat{\bm{S}}_{i}^{2}-2\hat{n}_{i}).$
(14)

Note that $\hat{\bm{S}}^{2}_{i}$, $\hat{S}_{iz}$ and $\hat{n}_{i}$ commute
with each other. The simultaneous eigenstates for these operators
$|S_{i},m_{i},n_{i}\rangle\ (-S_{i}\leq m_{i}\leq S_{i})$ satisfy

$\displaystyle\hat{\bm{S}}_{i}^{2}|S_{i},m_{i},n_{i}\rangle$ $\displaystyle=$
$\displaystyle S_{i}(S_{i}+1)|S_{i},m_{i},n_{i}\rangle,$ (15)
$\displaystyle\hat{S}_{iz}|S_{i},m_{i},n_{i}\rangle$ $\displaystyle=$
$\displaystyle m_{i}|S_{i},m_{i},n_{i}\rangle,$ (16)
$\displaystyle\hat{n}_{i}|S_{i},m_{i},n_{i}\rangle$ $\displaystyle=$
$\displaystyle n_{i}|S_{i},m_{i},n_{i}\rangle.$ (17)

$|S_{i},m_{i},n_{i}\rangle$ is thus an energy eigenstate of $\hat{H}^{0}_{i}$:

$\displaystyle\hat{H}^{0}_{i}|S_{i},m_{i},n_{i}\rangle=E^{0}(S_{i},n_{i})|S_{i},m_{i},n_{i}\rangle,$
(18)

where the energy eigenvalue $E^{0}(S_{i},n_{i})$ is given by

$\displaystyle E^{0}(S_{i},n_{i})$ $\displaystyle=$ $\displaystyle-\mu
n_{i}+\frac{U_{0}}{2}n_{i}(n_{i}-1)$ (19)
$\displaystyle+\frac{U_{2}}{2}[S_{i}(S_{i}+1)-2n_{i}].$

$S_{i}$ takes the minimum value in the ground state due to $U_{2}>0$.
Since the orbital part of the wave function in a single lattice site is
symmetric under permutation of any two atoms, the spin part has to be also
symmetric under permutation of atoms due to Bose statistics. As a result,
$S_{i}$ is even $(S_{i}=0,2,4,...,n_{i})$ for even $n_{i}$, while $S_{i}$ is
odd $(S_{i}=1,3,5,...,n_{i})$ for odd $n_{i}$ [42]. The ground state of the
single-site Hamiltonian (14) is thus $|0,0,n_{i}\rangle$ for even filling
factors, while it is $|1,m_{i},n_{i}\rangle\,(m_{i}=1,0,-1)$ for odd filling
factors that is triply degenerate. For even filling factors, the ground state
$|0,0,n_{i}\rangle$ is a spin-singlet insulator, in which all atoms form spin-
singlet pairs [11] as

$\displaystyle|0,0,n_{i}\rangle=\frac{1}{\sqrt{f(\frac{n_{i}}{2};0)}}(\hat{\Theta}_{i}^{\dagger})^{\frac{n_{i}}{2}}|\mathrm{vac}\rangle,$
(20)

where
$\hat{\Theta}^{\dagger}_{i}=(\hat{b}^{\dagger}_{i0})^{2}-2\hat{b}^{\dagger}_{i1}\hat{b}^{\dagger}_{i-1}$
is the creation operator of a spin-singlet pair. The normalization factor is
given by [43]

$\displaystyle f(Q;S)=S!Q!2^{Q}\frac{(2Q+2S+1)!!}{(2S+1)!!}.$ (21)

For odd filling factors, there remains a single atom that cannot form a spin-
singlet pair in the ground state as

$\displaystyle|1,m_{i},n_{i}\rangle=\frac{1}{\sqrt{f(\frac{n_{i}-1}{2};1)}}\hat{b}^{\dagger}_{i\,m_{i}}(\hat{\Theta}^{\dagger}_{i})^{\frac{n_{i}-1}{2}}|\mathrm{vac}\rangle.$
(22)

###  II.2 Perturbative mean-field theory

We employ the perturbative mean-field theory to study the SF-MI phase
transition [11]. We introduce the superfluid order parameter
$\psi_{\alpha}\equiv\langle\hat{b}_{i\alpha}\rangle\equiv\sqrt{n_{\rm
c}}\zeta_{\alpha}$, where $n_{\rm c}=\sum_{\alpha=\pm 1,0}|\psi_{\alpha}|^{2}$
is the number of condensate atoms per site and $\zeta_{\alpha}$ is the
normalized spinor: $\sum_{\alpha=\pm 1,0}\zeta^{*}_{\alpha}\zeta_{\alpha}=1$.
Linearizing the hopping term with respect to fluctuation
$\delta\hat{b}_{i\alpha}=\hat{b}_{i\alpha}-\psi_{\alpha}$, we obtain

$\displaystyle-t\sum_{\langle
i,j\rangle,\alpha}(\hat{b}^{\dagger}_{i\alpha}\hat{b}_{j\alpha}+\mathrm{h.c.})$
$\displaystyle\sim-
zt\sum_{i,\alpha}[(\psi_{\alpha}^{*}\hat{b}_{i\alpha}+\psi_{\alpha}\hat{b}_{i\alpha}^{\dagger})]+ztN_{\rm
s}\sum_{\alpha}\psi_{\alpha}^{*}\psi_{\alpha},$ (23)

where $z=2d$ denotes the number of nearest-neighbor sites and $N_{\rm s}$ the
total number of lattice sites. The Hamiltonian (1) thus reduces to

$\displaystyle\hat{H}$ $\displaystyle=$
$\displaystyle\sum_{i}\hat{H}^{\mathrm{mf}}_{i},$ (24)
$\displaystyle\hat{H}^{\mathrm{mf}}_{i}$ $\displaystyle=$
$\displaystyle\hat{H}_{i}^{0}+\hat{V}_{i}+zt\sum_{\alpha}\psi_{\alpha}^{*}\psi_{\alpha},$
(25) $\displaystyle\hat{V}_{i}$ $\displaystyle=$ $\displaystyle-
zt\sum_{\alpha}(\psi_{\alpha}\hat{b}_{i\alpha}^{\dagger}+\psi_{\alpha}^{*}\hat{b}_{i\alpha}).$
(26)

Here, $\hat{V}_{i}$ describes transfer of atoms between the $i$-th site and
the condensate $\psi_{\alpha}$. We treat $\hat{V}_{i}$ as a perturbation
assuming small t. We omit the site index in the rest of this section.
We focus on the MI phase with even filling factors and the SF phase around it
in this paper. In this case, since the MI phase has a unique ground state
(20), the standard perturbative mean-field theory is applicable. For odd
filling factors, extensions are required for the degenerate ground state (22).
To describe the first-order SF-MI phase transition predicted in Ref. [14], it
is necessary to expand the ground-state energy up to sixth-order. Applying the
standard perturbation theory, the ground-state energy per site can be
calculated as

$\displaystyle E$ $\displaystyle=$ $\displaystyle E^{0}(0,n)+C_{2}n_{\rm
c}+C_{4}n_{\rm c}^{2}+C_{4}^{\prime}n_{\rm c}^{2}\langle\bm{F}\rangle^{2}$
(27) $\displaystyle+C_{6}n_{\rm c}^{3}+C_{6}^{\prime}n_{\rm
c}^{3}\langle\bm{F}\rangle^{2},$

where $\vec{\zeta}=(\zeta_{1},\zeta_{0},\zeta_{-1})^{T}$ and
$\langle\bm{F}\rangle=\sum_{\alpha,\beta}\zeta_{\alpha}^{*}\bm{F}_{\alpha\beta}\zeta_{\beta}$.
The details of the calculations and the explicit forms of $C_{2}$, $C_{4}$,
$C_{4}^{\prime}$, $C_{6}$, and $C_{6}^{\prime}$ are given in Appendix B. The
fourth-order term with $C_{4}^{\prime}>0$ lifts the degeneracy in spin-space
and the polar state is realized in the SF phase [11]:
$\langle\bm{F}\rangle^{2}$ takes its minimum in the ground state, i.e.,
$\langle\bm{F}\rangle=\bm{0}$.
Since all spinors are related to each other by gauge transformation
$e^{i\theta}$ and spin rotations $R(\alpha,\beta,\gamma)=e^{-i\alpha
F_{z}}e^{-i\beta F_{y}}e^{-i\gamma F_{z}}$, where $(\alpha,\beta,\gamma)$ are
the Euler angles, the general form of the order parameter for the polar state
is given as [44]

$\displaystyle\vec{\zeta}=e^{i\theta}R(\alpha,\beta,\gamma)\left[\begin{array}[]{c}0\\\
1\\\
0\end{array}\right]=e^{i\theta}\left[\begin{array}[]{c}-\frac{e^{-i\alpha}}{\sqrt{2}}\sin\beta\\\
\cos\beta\\\ \frac{e^{i\alpha}}{\sqrt{2}}\sin\beta\end{array}\right].$ (34)

Since $\vec{\zeta}$ in Eq. (34) is invariant under shift of $\gamma$ as well
as the discrete transformation
$(\beta,\theta)\rightarrow(\beta+\pi,\theta+\pi)$, the polar state has the
symmetry characterized by the isotropy group $H={\rm
U}(1)_{\gamma}\otimes({\bf Z}_{2})_{\beta,\theta}$ [8]. The symmetry of the
spin-singlet state in Eq. (20) is characterized by the isotropy group $G={\rm
U}(1)_{\theta}\otimes{\rm SO}(3)_{\alpha,\beta,\gamma}$. The order parameter
manifold for the polar state is thus given by [45]

$\displaystyle M=G/H$ $\displaystyle=$ $\displaystyle\frac{{\rm
U}(1)_{\theta}\otimes{\rm SO}(3)_{\alpha,\beta,\gamma}}{{\rm
U}(1)_{\gamma}\otimes({\mathbf{Z}}_{2})_{\beta,\theta}}$ (35) $\displaystyle=$
$\displaystyle\frac{{\rm U}(1)_{\theta}\otimes
S^{2}_{\alpha,\beta}}{({\mathbf{Z}}_{2})_{\beta,\theta}}.$

It is convenient to introduce the $d$ vector defined as [46]

$\displaystyle\bm{d}=\left[\begin{array}[]{c}d_{x}\\\ d_{y}\\\
d_{z}\end{array}\right]=\left[\begin{array}[]{c}\frac{1}{\sqrt{2}}(-\psi_{1}+\psi_{-1})\\\
\frac{1}{\sqrt{2}i}(\psi_{1}+\psi_{-1})\\\ \psi_{0}\end{array}\right].$ (42)

Note that the $d$ vector transforms as a vector under rotations in spin-space.
From Eq. (34), the polar state is characterized by the $d$ vector

$\displaystyle\bm{d}=\sqrt{n_{\rm
c}}e^{i\theta}\left[\begin{array}[]{c}\sin\beta\cos\alpha\\\
\sin\beta\sin\alpha\\\ \cos\beta\end{array}\right],$ (46)

where $\alpha$ and $\beta$, respectively, represent the polar and azimuthal
angles that parametrize the surface of the unit sphere $S^{2}$. Due to the
invariance of $\vec{\zeta}$ under
$(\beta,\theta)\rightarrow(\beta+\pi,\theta+\pi)$, $\bm{d}$ specifies a
preferred axis, not a preferred direction along that axis.

###  II.3 Second-order SF-MI phase transition

Figure 1: (Color online) (a) Phase diagram of the second-order SF-MI phase
transition. We set $U_{2}/U_{0}=0.33$. The dotted line represents the phase
boundary between the spin-singlet insulator with $n=2$ and the polar state
determined by the perturbative mean-field theory. The red crosses represent
the phase boundary obtained by numerically diagonalizing Eq. (25). The blue
dash-dotted line indicates the particle-hole symmetry line determined by $K=0$
(see Appendix D). (b) Ground state energy as a function of $\sqrt{n_{\rm c}}$
in the MI phase (squares), the SF phase (triangles), and at the phase boundary
(PB) (diamonds). The solid lines plot Eq. (47). The dotted lines are obtained
by diagonalizing Eq. (25).

In the case of the second-order SF-MI phase transition, where
$U_{2}/U_{0}>0.32$ [14], it is enough to expand the ground-state energy up to
fourth order as

$\displaystyle E=E^{0}(0,n)+C_{2}n_{\rm c}+C_{4}n_{\rm c}^{2}.$ (47)

Note that we have assumed the polar state and set
$\langle\bm{F}\rangle=\bm{0}$ in Eq. (47). The right-hand side of Eq. (47) is
plotted as a function of $\sqrt{n_{\rm c}}$ in Fig. 1 (b). From the standard
Landau theory for second-order phase transitions [47], the phase boundary
between the MI and SF phases is determined by the condition $C_{2}=0$. Figure
1 (a) shows the resulting phase diagram [11].
The ground-state energy can be calculated by numerically diagonalizing Eq.
(25), as shown in Fig. 1 (b). The phase boundary can be determined by $n_{\rm
c}$ that minimizes the numerically calculated ground-state energy. It
precisely agrees with that of the perturbative mean-field theory, as shown in
Fig. 1 (a).

###  II.4 First-order SF-MI phase transition

Figure 2: (Color online) (a) Phase diagram of the first-order SF-MI phase
transition. We set $U_{2}/U_{0}=0.31$. The dotted line represents the phase
boundary between the spin-singlet insulator with $n=2$ and the polar state
determined by the perturbative mean-field theory. The red crosses represent
the phase boundary obtained by numerically diagonalizing Eq. (25). The blue
dash-dotted line indicates the particle-hole symmetry line determined by $K=0$
(see Appendix D). The solid black line indicates the phase boundary of the
first-order phase transition. (b) Magnification of the tip of the Mott lobe in
(a). The solid line is the phase boundary. The metastable SF (MI) state exists
in the thick (thin) gray region, in which the condition
$3C_{2}C_{6}<C_{4}^{2}<4C_{2}C_{6}$ $(0<4C_{2}C_{6}<C_{4}^{2})$ is satisfied.
(c) Ground-state energy as a function of $\sqrt{n_{\rm c}}$. The solid lines
plot Eq. (48). The dotted lines are calculated by diagonalizing Eq. (25).
$zt/U_{0}$ and $\mu/U_{0}$ for each line are set to the values for the same
symbols in (b).

It has been found that the SF-MI phase transition is of first order in the
vicinity of the tip of the Mott lobe (see Fig. 2) when $U_{2}/U_{0}<0.32$
[14]. The metastable SF and MI phases appear inside the MI and SF regions,
respectively, in the phase diagram [17]. We need in this case the ground-state
energy expanded up to sixth order,

$\displaystyle E=E^{0}(0,n)+C_{2}n_{\rm c}+C_{4}n_{\rm c}^{2}+C_{6}n_{\rm
c}^{3}.$ (48)

Note that $C_{4}$ is negative on the first-order part of the phase boundary.
The right-hand side of Eq. (48) is plotted as a function of $\sqrt{n_{\rm c}}$
in Fig. 2 (c).
In the SF phase, the equation for $n_{\rm c}$,

$\displaystyle C_{2}n_{\rm c}+C_{4}n_{\rm c}^{2}+C_{6}n_{\rm c}^{3}=0,$ (49)

has a positive solution as shown in Fig. 2 (c) (see the curves with green
inverted triangles and red diamonds). The ground state is thus in the SF phase
if $C_{4}^{2}-4C_{2}C_{6}>0$. In the MI phase, i.e.,
$C_{4}^{2}-4C_{2}C_{6}<0$, Eq. (49) has a single solution $n_{\rm c}=0$ [see
the curves with orange circles and blue squares in Fig. 2(c)]. The phase
boundary is, therefore, determined by the condition $C_{4}^{2}-4C_{2}C_{6}=0$.
If $C_{2}>0$ in the SF phase, since the ground-state energy has a local
minimum at $\sqrt{n_{\rm c}}=0$, the metastable MI state appears. On the other
hand, the ground-state energy has a local minimum at $\sqrt{n_{\rm c}}>0$ in
the MI phase, if the metastable SF state exists. From the condition for a
local minimum $\partial E/\partial\sqrt{n_{\rm c}}=0$, we obtain

$\displaystyle n_{\rm
c}=0,\,\frac{-C_{4}\pm\sqrt{C_{4}^{2}-3C_{2}C_{6}}}{3C_{6}}.$ (50)

The positive solution ($n_{\rm c}>0$) in Eq. (50) corresponds to the
metastable SF state. The metastable SF state thus appears if
$C_{4}^{2}-3C_{2}C_{6}>0$. Consequently, the conditions for the existence of
the metastable SF and MI states are given as
$3C_{2}C_{6}<C_{4}^{2}<4C_{2}C_{6}$ and $0<4C_{2}C_{6}<C_{4}^{2}$,
respectively. Figures 2 (a) and 2 (b) show the regions of the metastable SF
and MI phases in the phase diagram when $U_{2}/U_{0}=0.31$. The phase boundary
obtained by diagonalizing Eq. (25) agrees well with that of the perturbative
mean-field theory, as shown in Fig. 2 (a), as far as $U_{2}/U_{0}\lesssim
0.32$. If $0<U_{2}/U_{0}\ll 0.32$, since the SF order parameter is large even
in the vicinity of the phase boundary, the perturbative expansion of the
ground-state energy by the SF order parameter breaks down. The correct phase
boundary is obtained by diagonalizing Eq. (25) [17].

##  III collective modes

Figure 3: (a) Dispersion relations for the collective modes in Eqs. (96),
(97), and (98). We set $U_{2}/U_{0}=0.33$, $zt/U_{0}=0.3834$, and
$\mu/U_{0}=0.9831$. Figure 4: (Color online) Schematic representation of the
spin-nematic modes $\omega_{\alpha}$ and $\omega_{\beta}$ in Eqs. (101) and
(102). $\omega_{\alpha}$ ($\omega_{\beta}$) indicates fluctuation of $\alpha$
($\beta$).

We study collective modes in the SF phase based on the TDGL equation,

$\displaystyle iK\partial_{t}\Psi_{\alpha}-J\partial_{t}^{2}\Psi_{\alpha}$
$\displaystyle=$
$\displaystyle-\frac{\nabla^{2}}{2m^{*}}\Psi_{\alpha}+C_{2}\Psi_{\alpha}+2c_{4}({\bf\Psi}^{\dagger}{\bf\Psi})\Psi_{\alpha}$
(51)
$\displaystyle+2c_{4}^{\prime}\langle\\!\langle\bm{F}\rangle\\!\rangle\cdot\sum_{\beta}\bm{F}_{\alpha\beta}\Psi_{\beta}$
$\displaystyle+3c_{6}({\bf\Psi}^{\dagger}{\bf\Psi})^{2}\Psi_{\alpha}+c_{6}^{\prime}\langle\\!\langle\bm{F}\rangle\\!\rangle^{2}\Psi_{\alpha}$
$\displaystyle+2c_{6}^{\prime}({\bf\Psi}^{\dagger}{\bf\Psi})\langle\\!\langle\bm{F}\rangle\\!\rangle\cdot\sum_{\beta}\bm{F}_{\alpha\beta}\Psi_{\beta},$

where ${\bf\Psi}=(\Psi_{1},\Psi_{0},\Psi_{-1})^{T}$ denotes the SF order
parameter. (The derivation of the TDGL equation is summarized in Appendix C).
Here
$\langle\\!\langle\bm{F}\rangle\\!\rangle=\sum_{\alpha,\beta}\Psi_{\alpha}^{*}\bm{F}_{\alpha\beta}\Psi_{\beta}$
is the spin average.

Neglecting the sixth-order terms, we first examine the basic characters of
collective modes within the fourth-order TDGL equation. We assume a
commensurate filling and set $K=0$ for simplicity [see Eq. (312) in Appendix
D].
We set $(\alpha,\beta,\theta)=(0,\pi/2,0)$ for the static order parameter,

$\displaystyle{\bf\Psi}^{\rm
0}=\sqrt{\frac{\rho}{2}}\left[\begin{array}[]{c}-1\\\ 0\\\
1\end{array}\right],$ (55)

where $\rho=({{\bf\Psi}^{\rm 0}})^{\dagger}{\bf\Psi}^{\rm 0}=-C_{2}/2c_{4}$ is
the SF density. The $d$ vector for Eq. (55) $\bm{d}=\sqrt{\rho}(1,0,0)^{T}$ is
shown in Fig. 4. We introduce fluctuation of the order parameter around the
static solution ${\bf\Psi}^{\rm 0}$ as

$\displaystyle{\bf\Psi}$ $\displaystyle=$
$\displaystyle\sqrt{\rho+\delta\rho}e^{i\delta\theta}R(\delta\alpha,\pi/2+\delta\beta,\delta\gamma)\left[\begin{array}[]{c}0\\\
1\\\ 0\end{array}\right]$ (59) $\displaystyle\sim$
$\displaystyle{\bf\Psi}^{\rm 0}+\delta{\bf\Psi},$ (60)

where

$\displaystyle\delta{\bf\Psi}=\left[\begin{array}[]{c}\delta\Psi_{1}\\\
\delta\Psi_{0}\\\
\delta\Psi_{-1}\end{array}\right]=\sqrt{\rho}\left[\begin{array}[]{c}-\frac{1}{\sqrt{2}}\left(\frac{\delta\rho}{2\rho}+i\delta\theta-i\delta\alpha\right)\\\
-\delta\beta\\\
\frac{1}{\sqrt{2}}\left(\frac{\delta\rho}{2\rho}+i\delta\theta+i\delta\alpha\right)\end{array}\right].$
(67) (68)

The $d$ vector that corresponds to Eq. (60) is given by

$\displaystyle\bm{d}=\sqrt{\rho}\left[\begin{array}[]{c}1\\\ 0\\\
0\end{array}\right]+\sqrt{\rho}\left[\begin{array}[]{c}\frac{\delta\rho}{2\rho}+i\delta\theta\\\
\delta\alpha\\\ -\delta\beta\end{array}\right].$ (75)

Linearizing the fourth-order TDGL equation by fluctuation $\delta{\bf\Psi}$
and Fourier transforming by
$\delta\Psi_{\mu}=\sum_{\bm{k},\omega}e^{i(\bm{k}\cdot\bm{r}-\omega
t)}\delta\Phi_{\mu}({\bm{k}},\omega)\ (\mu=0,\pm 1)$, we obtain

$\displaystyle\left[\begin{array}[]{cccc}\varepsilon_{k}+c_{+}\rho&c_{+}\rho&c_{-}\rho&c_{-}\rho\\\
c_{+}\rho&\varepsilon_{k}+c_{+}\rho&c_{-}\rho&c_{-}\rho\\\
c_{-}\rho&c_{-}\rho&\varepsilon_{k}+c_{+}\rho&c_{+}\rho\\\
c_{-}\rho&c_{-}\rho&c_{+}\rho&\varepsilon_{k}+c_{+}\rho\\\
\end{array}\right]\left[\begin{array}[]{c}\delta\Phi_{1}\\\
\delta\Phi_{1}^{*}\\\ \delta\Phi_{-1}\\\ \delta\Phi_{-1}^{*}\\\
\end{array}\right]$ $\displaystyle=$ $\displaystyle
J\omega^{2}\left[\begin{array}[]{c}\delta\Phi_{1}\\\ \delta\Phi_{1}^{*}\\\
\delta\Phi_{-1}\\\ \delta\Phi_{-1}^{*}\\\ \end{array}\right],$ (88)
$\displaystyle\left[\begin{array}[]{cc}\varepsilon_{k}+2c_{4}^{\prime}\rho&-2c_{4}^{\prime}\rho\\\
-2c_{4}^{\prime}\rho&\varepsilon_{k}+2c_{4}^{\prime}\rho\end{array}\right]\left[\begin{array}[]{c}\delta\Phi_{0}\\\
\delta\Phi_{0}^{*}\end{array}\right]$ $\displaystyle=$ $\displaystyle
J\omega^{2}\left[\begin{array}[]{c}\delta\Phi_{0}\\\
\delta\Phi_{0}^{*}\end{array}\right],$ (95)

where $\varepsilon_{k}\equiv\frac{\bm{k}^{2}}{2m^{*}}$ and $c_{\pm}\equiv
c_{4}\pm c_{4}^{\prime}$. Solving Eqs. (88) and (95), we obtain the three
degenerate gapless modes with the same dispersion,

$\displaystyle\omega_{\theta}(\bm{k})=\omega_{\alpha}(\bm{k})=\omega_{\beta}(\bm{k})=\sqrt{\varepsilon_{k}/J}.$
(96)

In addition, we obtain the three gapful modes that have the dispersions,

$\displaystyle\omega_{\rho}(\bm{k})$ $\displaystyle=$
$\displaystyle\sqrt{[\varepsilon_{k}+4c_{4}{\rho}]/J},$ (97)
$\displaystyle\omega_{y}(\bm{k})=\omega_{z}(\bm{k})$ $\displaystyle=$
$\displaystyle\sqrt{\left[\varepsilon_{k}+4c_{4}^{\prime}{\rho}\right]/J}.$
(98)

Figure 3 shows the dispersions in Eqs. (96), (97), and (98). The amplitudes of
the normal modes for $\omega_{\theta}$ and $\omega_{\rho}$ are given,
respectively, as

$\displaystyle\delta\Phi_{1}-\delta\Phi_{1}^{*}-\delta\Phi_{-1}+\delta\Phi_{-1}^{*}$
$\displaystyle\propto$ $\displaystyle\delta\theta,$ (99)
$\displaystyle\delta\Phi_{1}+\delta\Phi_{1}^{*}-\delta\Phi_{-1}-\delta\Phi_{-1}^{*}$
$\displaystyle\propto$ $\displaystyle\delta\rho.$ (100)

They represent the phase and amplitude modes that arise due to the spontaneous
breaking of ${{\rm U}(1)}_{\theta}$ in Eq. (35) [48].
The amplitudes of the normal modes for
$\omega_{\alpha}(\bm{k}),\,\omega_{\beta}(\bm{k})$, $\omega_{y}(\bm{k})$, and
$\omega_{z}(\bm{k})$ are given, respectively, as

$\displaystyle\delta\Phi_{1}-\delta\Phi_{1}^{*}+\delta\Phi_{-1}-\delta\Phi_{-1}^{*}$
$\displaystyle\propto$ $\displaystyle\delta\langle\\!\langle
Q_{xy}\rangle\\!\rangle\propto\delta\alpha,$ (101)
$\displaystyle\delta\Phi_{0}+\delta\Phi_{0}^{*}$ $\displaystyle\propto$
$\displaystyle\delta\langle\\!\langle
Q_{zx}\rangle\\!\rangle\propto\delta\beta,$ (102)
$\displaystyle\delta\Phi_{0}-\delta\Phi_{0}^{*}$ $\displaystyle\propto$
$\displaystyle\delta\langle\\!\langle F_{y}\rangle\\!\rangle,$ (103)
$\displaystyle\delta\Phi_{1}+\delta\Phi_{1}^{*}+\delta\Phi_{-1}+\delta\Phi_{-1}^{*}$
$\displaystyle\propto$ $\displaystyle\delta\langle\\!\langle
F_{z}\rangle\\!\rangle.$ (104)

Here, we introduce the nematic tensor [49],

$\displaystyle
Q_{\mu\nu}\equiv\frac{1}{2}(F_{\mu}F_{\nu}+F_{\nu}F_{\mu})-\delta_{\mu\nu}\frac{\bm{F}^{2}}{3}.$
(105)

$Q_{\mu\nu}$ characterizes the spin-nematic order [50] and $\langle\\!\langle
Q_{\mu\nu}\rangle\\!\rangle=\sum_{\alpha,\beta}\Psi_{\alpha}^{*}(Q_{\mu\nu})_{\alpha\beta}\Psi_{\beta}$.
$\omega_{y}$ and $\omega_{z}$ represent spin-wave modes associated with
fluctuations of magnetization, while $\omega_{\alpha}$ and $\omega_{\beta}$
represent spin-nematic modes associated with fluctuations of the nematic
tensor [51]. They induce fluctuations of the polar and azimuthal angles of the
$d$ vector $\alpha$ and $\beta$, as shown in Fig. 4. The spin-nematic modes
$\omega_{\alpha}$ and $\omega_{\beta}$ arise due to the spontaneous breaking
of ${S}^{2}_{\alpha,\beta}$ in Eq. (35).
In the case of the first-order phase transition, dynamics of the SF order
parameter is governed by the sixth-order TDGL equation (51). In the same
manner as the fourth-order TDGL equation, the dispersions of the collective
modes can be derived as

$\displaystyle\omega_{\rho}(\bm{k})$ $\displaystyle=$
$\displaystyle\sqrt{[\varepsilon_{k}+4(c_{4}+3c_{6}{\rho}^{\prime}){\rho}^{\prime}]/J},$
(106) $\displaystyle\omega_{y}(\bm{k})$ $\displaystyle=$
$\displaystyle\omega_{z}(\bm{k})=\sqrt{[\varepsilon_{k}+4(c_{4}^{\prime}+c_{6}^{\prime}{\rho}^{\prime}){\rho}^{\prime}]/J},$
(107) $\displaystyle\omega_{\theta}(\bm{k})$ $\displaystyle=$
$\displaystyle\omega_{\alpha}(\bm{k})=\omega_{\beta}(\bm{k})=\sqrt{\varepsilon_{k}/J},$
(108)

where the superfluid density is given by

$\displaystyle{\rho}^{\prime}=\frac{-c_{4}+\sqrt{c_{4}^{2}-3C_{2}c_{6}}}{3c_{6}}.$
(109)

Figures 5 (a) and 5 (b) plot the energy gap of the gapful modes
$\Delta_{\rho}=\omega_{\rho}(0)$ and $\Delta_{y}=\omega_{y}(0)=\omega_{z}(0)$
as functions of $U_{0}$ near the phase boundary. $\Delta_{\rho}$ and
$\Delta_{y}$ vanish at the second-order phase boundary, whereas they are
finite at the first-order phase boundary. The finite energy gap at the phase
boundary reflects the existence of the metastable SF phase. $\Delta_{\rho}$
vanishes at $U_{0}=U_{\rm 0c}^{\prime}$ simultaneously with the disappearance
of the metastable SF state.

The TDGL method is valid in the vicinity of the MI phase, where the order
parameter is so small that the perturbative mean-field expansion is allowed.
It describes low energy superfluid dynamics, in which the order parameter
varies slowly in space and time. Thus, the analytical expressions for the
collective modes derived by the TDGL equation should be correct in the low
energy regime for small momenta.

To check the validity of the TDGL equation, we also evaluate the excitation
spectrum of collective modes numerically using the time-dependent Gutzwiller
ansatz [52, 7]. The time-dependent Gutzwiller ansatz approximates the many-
body wave function as a product of single-site wave functions, each of which
is expanded by finite number of Fock states. The validity of the time-
dependent Gutzwiller ansatz depends crucially on the number of Fock states
involved in the expansion. This restriction to the number of Fock states as
well as the approximation for the many-body wave function as a product of
single-site wave functions limits the validity of the time-dependent
Gutzwiller ansatz within the low energy regime near the phase boundary,
because large numbers of Fock states are required to describe the SF phase in
highly excited states and/or away from the phase boundary. The collective
modes calculated by the TDGL equation and the time-dependent Gutzwiller
ansatz, therefore, should agree in the low energy regime for small momenta in
the vicinity of the MI phase.

Figure 5 shows a comparison of the excitation spectrum calculated by the TDGL
equation and the time-dependent Gutzwiller ansatz. They agree well in the low
energy regime for small $k$ as expected.

Figure 5: (a) Energy gap of the gapful modes in Eqs. (97) and (98) normalized
by the spin-independent interaction at the second-order SF-MI phase boundary
$U_{\rm 0c}$. We set $U_{2}/U_{0}=0.33$ and $\mu/U_{0}=0.9831$. (b) Energy gap
of the gapful modes in Eqs. (106) and (107) normalized by $U_{\rm 0c}$. The
metastable SF phase disappears at $U_{0}=U_{0{\rm c}}^{\prime}$. We set
$U_{2}/U_{0}=0.31$ and $\mu/U_{0}=1.008$. Figure 6: (a) Dispersion relations
for the collective modes in Eqs. (106), (107), and (108). We set
$U_{2}/U_{0}=0.50$, $zt/U_{0}=0.451$, and $\mu/U_{0}=0.775$. The red crosses
are the dispersion relations calculated using the time-dependent Gutzwiller
ansatz. We expand the single-site wave function by the Fock state
$|n_{1},n_{0},n_{-1}\rangle$ with $0\leq n_{1}+n_{0}+n_{-1}\leq 10$.

##  IV Dynamical instability of supercurrents

###  IV.1 Mass currents

To examine the stability of mass currents, we consider a current carrying
state, in which $m=\pm 1$ components are flowing in the same direction with
momentum $\bm{p}$, as

$\displaystyle{\bf\Psi}^{\rm 0}_{\rm
m}=\sqrt{\frac{\bar{\rho}}{2}}e^{i\bm{p}\cdot\bm{r}}\left[\begin{array}[]{c}-1\\\
0\\\ 1\end{array}\right],$ (113)

where $\bar{\rho}=-(\varepsilon_{p}+C_{2})/2c_{4}$ is the SF density. Note
that $p\equiv|\bm{p}|$ should be smaller than the maximum value $p_{\rm
max}=\sqrt{-2m^{*}C_{2}}$ for $\bar{\rho}>0$. The mass and spin currents for
Eq. (113) are given by

$\displaystyle\bm{j}_{\rm m}$ $\displaystyle=$
$\displaystyle\bar{\rho}\frac{\bm{p}}{m^{*}},$ (114)
$\displaystyle\bm{j}^{\bm{n}}_{\rm s}$ $\displaystyle=$ $\displaystyle{\bf
0},$ (115)

respectively. Here, mass and spin currents are defined, respectively, as

$\displaystyle\bm{j}_{\rm m}$ $\displaystyle=$
$\displaystyle\frac{1}{2im^{*}}[{\bf\Psi}^{\dagger}\nabla{\bf\Psi}-(\nabla{\bf\Psi}^{\dagger}){\bf\Psi}],$
(116) $\displaystyle\bm{j}^{\mu}_{\rm s}$ $\displaystyle=$
$\displaystyle\frac{1}{2im^{*}}[{\bf\Psi}^{\dagger}F_{\mu}\nabla{\bf\Psi}-(\nabla{\bf\Psi}^{\dagger})F_{\mu}{\bf\Psi}].$
(117)

$\bm{j}^{\mu}_{\rm s}$ denotes a spin current for the $\mu$-th component of
magnetization $\langle\\!\langle F_{\mu}\rangle\\!\rangle\ (\mu=x,y,z)$. [We
derive Eqs. (116) and (117) in Appendix D.] We introduce fluctuation of the
order parameter around the static solution ${\bf\Psi}^{\rm 0}_{\rm m}$ as

$\displaystyle{\bf\Psi}$ $\displaystyle=$ $\displaystyle
e^{i\bm{p}\cdot\bm{r}}\sqrt{\bar{\rho}+\delta\rho}e^{i\delta\theta}R(\delta\alpha,\pi/2+\delta\beta,\delta\gamma)\left[\begin{array}[]{c}0\\\
1\\\ 0\end{array}\right]$ (121) $\displaystyle\sim$
$\displaystyle{\bf\Psi}^{\rm 0}_{\rm m}+\delta{\bf\Psi},$ (122)

where

$\displaystyle\delta{\bf\Psi}$ $\displaystyle=$ $\displaystyle
e^{i\bm{p}\cdot\bm{r}}\left[\begin{array}[]{c}\delta\Psi_{1}\\\
\delta\Psi_{0}\\\ \delta\Psi_{-1}\end{array}\right]$ (126) $\displaystyle=$
$\displaystyle\sqrt{\bar{\rho}}e^{i\bm{p}\cdot\bm{r}}\left[\begin{array}[]{c}-\frac{1}{\sqrt{2}}\left(\frac{\delta\rho}{2\bar{\rho}}+i\delta\theta-i\delta\alpha\right)\\\
-\delta\beta\\\
\frac{1}{\sqrt{2}}\left(\frac{\delta\rho}{2\bar{\rho}}+i\delta\theta+i\delta\alpha\right)\end{array}\right].$
(130)

Substituting Eq. (122) into the fourth-order TDGL equation with $K=0$ and
linearizing with respect to $\delta{\bf\Psi}$, we obtain

$\displaystyle\left[\begin{array}[]{cccc}\varepsilon_{k}+\frac{\bm{p}\cdot\bm{k}}{m^{*}}+c_{+}\bar{\rho}&c_{+}\bar{\rho}&c_{-}\bar{\rho}&c_{-}\bar{\rho}\\\
c_{+}\bar{\rho}&\varepsilon_{k}-\frac{\bm{p}\cdot\bm{k}}{m^{*}}+c_{+}\bar{\rho}&c_{-}\bar{\rho}&c_{-}\bar{\rho}\\\
c_{-}\bar{\rho}&c_{-}\bar{\rho}&\varepsilon_{k}+\frac{\bm{p}\cdot\bm{k}}{m^{*}}+c_{+}\bar{\rho}&c_{+}\bar{\rho}\\\
c_{-}\bar{\rho}&c_{-}\bar{\rho}&c_{+}\bar{\rho}&\varepsilon_{k}-\frac{\bm{p}\cdot\bm{k}}{m^{*}}+c_{+}\bar{\rho}\\\
\end{array}\right]\left[\begin{array}[]{c}\delta{\Phi}_{1}\\\
\delta{\Phi}_{1}^{*}\\\ \delta{\Phi}_{-1}\\\ \delta{\Phi}_{-1}^{*}\\\
\end{array}\right]$ $\displaystyle=$ $\displaystyle
J\omega^{2}\left[\begin{array}[]{c}\delta{\Phi}_{1}\\\ \delta{\Phi}_{1}^{*}\\\
\delta{\Phi}_{-1}\\\ \delta{\Phi}_{-1}^{*}\\\ \end{array}\right],$ (143)

$\displaystyle\left[\begin{array}[]{cc}\varepsilon_{k}+\frac{\bm{p}\cdot\bm{k}}{m^{*}}+2c_{4}^{\prime}\bar{\rho}&-2c_{4}^{\prime}\bar{\rho}\\\
-2c_{4}^{\prime}\bar{\rho}&\varepsilon_{k}-\frac{\bm{p}\cdot\bm{k}}{m^{*}}+2c_{4}^{\prime}\bar{\rho}\end{array}\right]\left[\begin{array}[]{c}\delta{\Phi}_{0}\\\
\delta{\Phi}_{0}^{*}\end{array}\right]$ $\displaystyle=$ $\displaystyle
J\omega^{2}\left[\begin{array}[]{c}\delta{\Phi}_{0}\\\
\delta{\Phi}_{0}^{*}\end{array}\right].$ (150)

From Eqs. (143) and (150), we obtain the three gapless modes with the
dispersions,

$\displaystyle\omega_{\theta,\rho}(\bm{k},\bm{p})$ $\displaystyle=$
$\displaystyle\sqrt{\left[\varepsilon_{k}+2c_{4}\bar{\rho}-\sqrt{(2c_{4}\bar{\rho})^{2}+\left(\frac{\bm{p}\cdot\bm{k}}{m^{*}}\right)^{2}}\right]/J},$
$\displaystyle\omega_{\alpha,z}(\bm{k},\bm{p})$ $\displaystyle=$
$\displaystyle\omega_{\beta,y}(\bm{k},\bm{p})$ $\displaystyle=$
$\displaystyle\sqrt{\left[\varepsilon_{k}+2c_{4}^{\prime}\bar{\rho}-\sqrt{(2c_{4}^{\prime}\bar{\rho})^{2}+\left(\frac{\bm{p}\cdot\bm{k}}{m^{*}}\right)^{2}}\right]/J}.$

The amplitudes of the normal modes for $\omega_{\theta,\rho}$,
$\omega_{\alpha,z}$, and $\omega_{\beta,y}$ are given, respectively, as

$\displaystyle
L_{\theta,\rho}\frac{\delta\rho}{\bar{\rho}}+L_{\theta,\theta}\delta\theta,$
(153) $\displaystyle L_{\alpha,z}\frac{\delta\langle\\!\langle
F_{z}\rangle\\!\rangle}{\bar{\rho}}+L_{\alpha,\alpha}\frac{\delta\langle\\!\langle
Q_{xy}\rangle\\!\rangle}{\bar{\rho}},$ (154) $\displaystyle
L_{\beta,y}\frac{\delta\langle\\!\langle
F_{y}\rangle\\!\rangle}{\bar{\rho}}+L_{\beta,\beta}\frac{\delta\langle\\!\langle
Q_{xy}\rangle\\!\rangle}{\bar{\rho}},$ (155)

where

$\displaystyle L_{\theta,\theta}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[1-\frac{\bm{p}\cdot\bm{k}}{2c_{4}\bar{\rho}m^{*}}+\sqrt{1+\left(\frac{\bm{p}\cdot\bm{k}}{2c_{4}\bar{\rho}m^{*}}\right)^{2}}\right],$
$\displaystyle L_{\theta,\rho}$ $\displaystyle=$
$\displaystyle\frac{i}{4}\left[1+\frac{\bm{p}\cdot\bm{k}}{2c_{4}\bar{\rho}m^{*}}-\sqrt{1+\left(\frac{\bm{p}\cdot\bm{k}}{2c_{4}\bar{\rho}m^{*}}\right)^{2}}\right],$
$\displaystyle L_{\alpha,\alpha}$ $\displaystyle=$ $\displaystyle
L_{\beta,\beta}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[1-\frac{\bm{p}\cdot\bm{k}}{2c^{\prime}_{4}\bar{\rho}m^{*}}+\sqrt{1+\left(\frac{\bm{p}\cdot\bm{k}}{2c^{\prime}_{4}\bar{\rho}m^{*}}\right)^{2}}\right],$
$\displaystyle L_{\alpha,z}$ $\displaystyle=$ $\displaystyle L_{\beta,y}$
$\displaystyle=$
$\displaystyle\frac{i}{4}\left[1+\frac{\bm{p}\cdot\bm{k}}{2c^{\prime}_{4}\bar{\rho}m^{*}}-\sqrt{1+\left(\frac{\bm{p}\cdot\bm{k}}{2c^{\prime}_{4}\bar{\rho}m^{*}}\right)^{2}}\right].$

Figure 7: Dispersion relations for the collective modes in Eqs.
(LABEL:eq:omega_theta,rho) and (LABEL:eq:omega_alpha,z) in the presence of a
mass current with momentum $\bm{p}$ . We set $U_{2}/U_{0}=0.33$,
$zt/U_{0}=0.3834$, and $\mu/U_{0}=0.9831$. The solid lines in (a) and (b)
denote the real parts of Eqs. (LABEL:eq:omega_theta,rho) and
(LABEL:eq:omega_alpha,z), respectively. The dashed lines in (a) and (b)
indicate the imaginary parts of Eqs. (LABEL:eq:omega_theta,rho) and
(LABEL:eq:omega_alpha,z), respectively. Figure 8: Critical momenta for mass
currents $p_{\theta}^{\rm c}$ and $p_{\alpha}^{\rm c}$ in Eqs. (160) and (161)
as functions of $U_{0}/U_{\rm 0c}$ for the second-order phase transition.
$U_{\rm 0c}$ is the critical value of the spin-independent interaction for the
SF-MI phase transition. We set $U_{2}/U_{0}=0.33$ and $\mu/U_{0}=0.9831$.

Since the phase factor $e^{i\bm{p}\cdot\bm{r}}$ in Eq. (113) breaks the
particle-hole symmetry, a pair of gapful mode and a gapless mode in Eqs. (96),
(97), and (98) are coupled in the presence of a mass current and yields a
single gapless mode. $\omega_{\rho}$ and $\omega_{\theta}$ are coupled, for
example, and yield $\omega_{\theta,\rho}$.
Figure 7 shows the evolution of the dispersions in Eqs.
(LABEL:eq:omega_theta,rho) and (LABEL:eq:omega_alpha,z), as $p$ is increased.
Here, we set $\bm{k}$ being antiparallel to $\bm{p}$. If $p$ exceeds the
critical momenta for dynamical instabilities, the dispersion of collective
modes acquires an imaginary part. The critical momenta $p_{\theta}^{\rm c}$
and $p_{\alpha}^{\rm c}$ for the onset of the dynamical instabilities induced
by $\omega_{\theta,\rho}$ and $\omega_{\alpha,z}$ are given, respectively, by

$\displaystyle p_{\theta}^{\rm c}=\sqrt{-\frac{2}{3}C_{2}m^{*}}=\frac{p_{\rm
max}}{\sqrt{3}},$ (160) $\displaystyle p_{\alpha}^{\rm
c}=\sqrt{-\frac{2C_{2}c_{4}^{\prime}m^{*}}{2c_{4}^{\prime}+c_{4}}}=p_{\rm
max}\sqrt{\frac{c_{4}^{\prime}}{2c_{4}^{\prime}+c_{4}}}.$ (161)

The critical momentum for $\omega_{\beta,y}$ is equivalent to $p_{\alpha}^{\rm
c}$. Figure 8 shows the critical momenta $p_{\theta}^{\rm c}$ and
$p_{\alpha}^{\rm c}$ as functions of $U_{0}$ in the vicinity of the second-
order phase transition. The dynamical instability induced by
$\omega_{\theta,\rho}$ for $p>p_{\theta}^{\rm c}$ leads to exponential growth
of $\delta\theta$ and $\delta\rho$ in time.

Figure 9: Critical momentum for mass currents ${p_{\theta}^{\prime}}^{\rm c}$
in Eq. (166) as a function of $U_{0}/U_{\rm 0c}$ for the first-order phase
transition. $U_{\rm 0c}$ is the critical value of the spin-independent
interaction for the SF-MI phase transition. $U_{\rm 0c}^{\prime}$ is the value
of the spin-independent interaction at the boundary, at which the metastable
SF phase disappears. We set $U_{2}/U_{0}=0.31$ and $\mu/U_{0}=1.008$.

Meanwhile, the dynamical instabilities induced by $\omega_{\alpha,z}$ and
$\omega_{\beta,y}$ for $p>p_{\alpha}^{\rm c}$ lead to exponential growth of
$\delta\langle\\!\langle F_{y}\rangle\\!\rangle$, $\delta\langle\\!\langle
F_{z}\rangle\\!\rangle$, $\delta\langle\\!\langle Q_{xy}\rangle\\!\rangle$,
and $\delta\langle\\!\langle Q_{zx}\rangle\\!\rangle$ in time.
In the case of the first-order phase transition, using the sixth-order TDGL
equation (51), the dispersions for the gapless modes are obtained as

$\displaystyle\omega_{\theta,\rho}(\bm{k},\bm{p})$ $\displaystyle=$
$\displaystyle\left[\varepsilon_{k}+2(c_{4}+3c_{6}{\bar{\rho}}^{\prime}){\bar{\rho}}^{\prime}\right.$
$\displaystyle\left.-\sqrt{(\bm{k}\cdot\bm{p}/m^{*})^{2}+4(c_{4}+3c_{6}{\bar{\rho}}^{\prime})^{2}{\bar{\rho}}^{{}^{\prime}2}}\right]^{1/2},$
$\displaystyle\omega_{\alpha,z}(\bm{k},\bm{p})$ $\displaystyle=$
$\displaystyle\omega_{\beta,y}(\bm{k},\bm{p})$ $\displaystyle=$
$\displaystyle\left[\varepsilon_{k}+2(c_{4}^{\prime}+c_{6}^{\prime}{\bar{\rho}}^{\prime}){\bar{\rho}}^{\prime}\right.$
$\displaystyle\left.-\sqrt{(\bm{k}\cdot\bm{p}/m^{*})^{2}+4(c_{4}^{\prime}+c_{6}^{\prime}{\bar{\rho}}^{\prime})^{2}{\bar{\rho}}^{{}^{\prime}2}}\right]^{1/2},$

where the superfluid density is given by

$\displaystyle\bar{\rho}^{\prime}=\frac{-c_{4}+\sqrt{c_{4}^{2}-3c_{6}(C_{2}+\varepsilon_{p})}}{3c_{6}}.$
(164)

The maximum value of the momentum $p_{\rm max}^{\prime}$, which satisfies
$\bar{\rho}^{\prime}>0$, is given by

$\displaystyle p_{\rm
max}^{\prime}=\sqrt{2m^{*}\left(\frac{c_{4}^{2}}{3c_{6}}-C_{2}\right)}.$ (165)

From Eq. (LABEL:eq:omega'_tp), the critical momentum
${p_{\theta}^{\prime}}^{\rm c}$ for the onset of dynamical instabilities
induced by $\omega_{\theta,\rho}$ is

$\displaystyle{p_{\theta}^{\prime}}^{\rm
c}=\sqrt{\frac{m^{*}}{4}\left[\frac{c_{4}^{2}}{c_{6}}-4C_{2}+\sqrt{\frac{c_{4}^{2}}{c_{6}}\left(\frac{c_{4}^{2}}{c_{6}}-\frac{8}{3}C_{2}\right)}\right]}.$
(166)

Figure 9 shows the critical momenta $p^{\rm c}_{\theta}$ and
${p_{\theta}^{\prime}}^{\rm c}$ as functions of $U_{0}$.
In the case of the second-order phase transition, we find $p_{\theta}^{\rm
c}\rightarrow 0$ as $U_{0}\rightarrow U_{\rm 0c}$, as shown in Fig. 8. In
fact, from Eq. (160), one finds $p_{\theta}^{\rm c}=0$ due to $C_{2}=0$ at the
phase boundary. This can be understood as follows: The dynamical instability
is expected to occur if the phase gradient per healing length $\xi$ exceeds
$\pi/2$, i.e., $p\xi>\pi/2$, due to negative effective mass [39]. Since the
healing length $\xi=1/\sqrt{4|C_{2}|m^{*}}$ diverges as $|C_{2}|\rightarrow 0$
at the phase boundary, $p_{\theta}^{\rm c}$ also vanishes at $U_{0}=U_{\rm
0c}$ [53]. Meanwhile, in the case of the first-order phase transition, Fig. 9
shows that ${p_{\theta}^{\prime}}^{\rm c}>0$ at the phase boundary reflecting
the finite healing length $\xi=1/\sqrt{8m^{*}(|c_{4}|\rho^{\prime}-C_{2})}$
[53]. Furthermore, it is finite throughout the metastable SF phase and
vanishes at $U_{0}=U_{\rm 0c}^{\prime}$, at which $C_{4}^{2}=3C_{2}C_{6}$
holds from Eq. (50) and the healing length $\xi$ diverges. The finite critical
momentum means that the metastable SF state can sustain a mass current without
dissipation as far as $p<{p_{\theta}^{\prime}}^{\rm c}$ and therefore exhibits
superfluidity.
In contrast to the case of the second-order phase transition, the dispersions
of $\omega_{\alpha,z}$ and $\omega_{\beta,y}$ exhibit very little change in
the presence of a finite mass current. Thus, in the case of the first-order
phase transition, $\omega_{\alpha,z}$ and $\omega_{\beta,y}$ do not cause
dynamical instabilities.

###  IV.2 Spin currents

To study the stability of spin currents, we assume a static solution, in which
$m=\pm 1$ components flowing in the opposite direction with the same momenta
$\bm{p}$, as

$\displaystyle{\bf\Psi}^{\rm 0}_{\rm
s}=\sqrt{\frac{\bar{\rho}}{2}}\left[\begin{array}[]{c}-e^{-i\bm{p}\cdot\bm{r}}\\\
0\\\ e^{i\bm{p}\cdot\bm{r}}\end{array}\right],$ (170)

where $\bar{\rho}=-(\varepsilon_{p}+C_{2})/2c_{4}$ is the SF density. Since
$m=\pm 1$ components flow in the opposite direction with the same momenta, the
net mass current vanishes ($\bm{j}_{\rm m}={\bf 0}$). The spin current for Eq.
(170) is given as

$\displaystyle\bm{j}^{\mu}_{\rm s}=-\frac{\bm{p}}{m^{*}}\delta_{\mu,z}.$ (171)

We introduce fluctuation of the order parameter around the static solution as

$\displaystyle{\bf\Psi}={\bf\Psi}^{\rm 0}_{\rm s}+\delta{\bf\Psi},$ (172)

where

$\displaystyle\delta{\bf\Psi}$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{c}\delta\Psi_{1}e^{-i\bm{p}\cdot\bm{r}}\\\
\delta\Psi_{0}\\\ \delta\Psi_{-1}e^{i\bm{p}\cdot\bm{r}}\end{array}\right]$
(176)
$\displaystyle=\sqrt{\bar{\rho}}\left[\begin{array}[]{c}-\frac{e^{-i\bm{p}\cdot\bm{r}}}{\sqrt{2}}\left(\frac{\delta\rho}{2\bar{\rho}}+i\delta\theta-i\delta\alpha\right)\\\
-\delta\beta\\\
\frac{e^{i\bm{p}\cdot\bm{r}}}{\sqrt{2}}\left(\frac{\delta\rho}{2\bar{\rho}}+i\delta\theta+i\delta\alpha\right)\end{array}\right].$
(180)

Substituting Eq. (172) into the fourth-order TDGL equation with $K=0$ and
linearizing with respect to fluctuation $\delta{\bf\Psi}$, we obtain the
equations for the Fourier components as

$\displaystyle\left[\begin{array}[]{cccc}\varepsilon_{k}+\frac{\bm{p}\cdot\bm{k}}{m^{*}}+c_{+}\bar{\rho}&c_{+}\bar{\rho}&c_{-}\bar{\rho}&c_{-}\bar{\rho}\\\
c_{+}\bar{\rho}&\varepsilon_{k}-\frac{\bm{p}\cdot\bm{k}}{m^{*}}+c_{+}\bar{\rho}&c_{-}\bar{\rho}&c_{-}\bar{\rho}\\\
c_{-}\bar{\rho}&c_{-}\bar{\rho}&\varepsilon_{k}-\frac{\bm{p}\cdot\bm{k}}{m^{*}}+c_{+}\bar{\rho}&c_{+}\bar{\rho}\\\
c_{-}\bar{\rho}&c_{-}\bar{\rho}&c_{+}\bar{\rho}&\varepsilon_{k}+\frac{\bm{p}\cdot\bm{k}}{m^{*}}+c_{+}\bar{\rho}\\\
\end{array}\right]\left[\begin{array}[]{c}\delta{\Phi}_{1}\\\
\delta{\Phi}_{1}^{*}\\\ \delta{\Phi}_{-1}\\\ \delta{\Phi}_{-1}^{*}\\\
\end{array}\right]$ $\displaystyle=$ $\displaystyle
J\omega^{2}\left[\begin{array}[]{c}\delta{\Phi}_{1}\\\ \delta{\Phi}_{1}^{*}\\\
\delta{\Phi}_{-1}\\\ \delta{\Phi}_{-1}^{*}\\\ \end{array}\right],$ (193)
$\displaystyle\left[\begin{array}[]{cc}\varepsilon_{k}-\varepsilon_{p}+2c_{4}^{\prime}\bar{\rho}&-2c_{4}^{\prime}\bar{\rho}\\\
-2c_{4}^{\prime}\bar{\rho}&\varepsilon_{k}-\varepsilon_{p}+2c_{4}^{\prime}\bar{\rho}\end{array}\right]\left[\begin{array}[]{c}\delta{\Phi}_{0}\\\
\delta{\Phi}_{0}^{*}\end{array}\right]$ $\displaystyle=$ $\displaystyle
J\omega^{2}\left[\begin{array}[]{c}\delta{\Phi}_{0}\\\
\delta{\Phi}_{0}^{*}\end{array}\right].$ (200)

From Eq. (LABEL:eq:phi_0_spin), we find that $\omega_{\beta}$ takes the form,

$\displaystyle{\omega}_{\beta}(k,p)=\sqrt{(\varepsilon_{k}-\varepsilon_{p})/J}.$
(202)

The amplitude of the normal mode for $\omega_{\beta}$ is

$\displaystyle\delta{\Phi}_{0}+\delta{\Phi}_{0}^{*}\propto\delta\beta\
(\propto\delta\langle\\!\langle Q_{zx}\rangle\\!\rangle).$ (203)

Equation (202) shows that ${\omega}_{\beta}(k<p,p)$ has an imaginary part for
any $p\neq 0$. This indicates that even an infinitesimally small amount of
spin current induces dynamical instabilities and leads to exponential growth
of the long wave length spin-nematic fluctuation. We also find that
$\omega_{\beta}$ derived from the sixth-order TDGL equation (51) with $K=0$
has the same dispersion relation as Eq. (202). The critical momentum is thus
also zero. The zero critical velocity for spin currents has been also reported
in the study of a spin-1 BEC [54, 55].

##  V Stability of spin currents

Figure 10: (Color online) Schematic illustration of the spatial variation of
the $d$ vector in Eq. (211) and the spin current $\bm{j}_{\rm s}^{z}$ in Eq.
(215). We set $\beta=\pi/2$ and $\bm{p}=(p,0,0)^{T}$. Figure 11: (Color
online) (a) Energy functional (216) as a function of $\beta$. We set
$U_{2}/U_{0}=0.33,\ zt/U_{0}=0.3834$, and $\mu/U_{0}=0.9831$. (b), (c), and
(d) are the schematic drawings of $d$ vectors (211) for $\beta=0,\pi/2,$ and
$\pi$, respectively.

To understand the origin of the dynamical instability of spin currents, we
consider the static solution,

$\displaystyle{\bf\Psi}_{\bm{p},\beta}=\sqrt{\rho}\left[\begin{array}[]{c}-\frac{e^{-i\bm{p}\cdot\bm{r}}}{\sqrt{2}}\sin\beta\\\
\cos\beta\\\
\frac{e^{i\bm{p}\cdot\bm{r}}}{\sqrt{2}}\sin\beta\end{array}\right],$ (207)

where ${\bf\Psi}_{\bm{p},\beta=\pi/2}={\bf\Psi}_{\rm s}^{\rm 0}$. The $d$
vector for Eq. (207) is given by

$\displaystyle\bm{d}=R_{z}(\bm{p}\cdot\bm{r})\sqrt{\rho}\left[\begin{array}[]{c}\sin\beta\\\
0\\\ \cos\beta\end{array}\right],$ (211)

where $R_{z}(\theta)$ represents the rotation matrix by angle $\theta$ about
the $z$-axis.
The SF density $\rho$ and spin currents $\bm{j}^{\mu}_{\rm s}$ for Eq. (207)
are given as

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle-\frac{C_{2}+\varepsilon_{p}\sin^{2}\beta}{2c_{4}},$ (212)
$\displaystyle\bm{j}_{\rm s}^{x}$ $\displaystyle=$
$\displaystyle\frac{\bm{p}}{m^{*}}\rho\sin\beta\cos\beta\cos(\bm{p}\cdot\bm{r}),$
(213) $\displaystyle\bm{j}_{\rm s}^{y}$ $\displaystyle=$
$\displaystyle\frac{\bm{p}}{m^{*}}\rho\sin\beta\cos\beta\sin(\bm{p}\cdot\bm{r}),$
(214) $\displaystyle\bm{j}_{\rm s}^{z}$ $\displaystyle=$
$\displaystyle-\frac{\bm{p}}{m^{*}}\sin^{2}\beta.$ (215)

Equations (213)-(215) demonstrate that spatial variation of the $d$ vector
induces a spin current. When $\beta=\pi/2$, for example, as one moves in the
direction of $\bm{p}$, the $d$ vector rotating in the $xy$ plane induces a
spin current $\bm{j}^{z}_{\rm s}||\bm{p}$, as shown in Fig. 10.
We focus on the case of the second-order phase transition. Substituting Eq.
(207) into Eq. (LABEL:eq:E_func) with $c_{6}=c_{6}^{\prime}=0$, the energy
functional for Eq. (207) can be calculated as

$\displaystyle{\cal
K}(\\{{\bf\Psi}_{\bm{p},\beta}\\})=-\frac{V}{4c_{4}}(\varepsilon_{p}\sin^{2}\beta-|C_{2}|)^{2},$
(216)

where $V=N_{\rm s}a^{d}$ is the volume of the system. Figure 11 shows ${\cal
K}(\\{{\bf\Psi}_{\bm{p},\beta}\\})$ as a function of $\beta$ and the
configuration of the $d$ vector (211) for each value of $\beta$. We find that
the degeneracy of the energy functional with respect to $\beta$ due to
$S^{2}_{\alpha,\beta}$ is lifted by a spin current: $\beta=0$ and $\pi$
correspond to the energy minima and $\beta=\pi/2$ to the maximum. In the
former, since $R_{z}(\bm{p}\cdot\bm{r})$ does not change the $d$ vector that
is in parallel with the $z$-axis, no spin currents are induced and thus ${\cal
K}(\\{{\bf\Psi}_{\bm{p},\beta}\\})$ takes the minima. In the latter, since the
$d$ vector is in the $xy$ plane, spatial variation of the $d$ vector by
$R_{z}(\bm{p}\cdot\bm{r})$ is maximized and thus ${\cal
K}(\\{{\bf\Psi}_{{\bm{p}},\beta}\\})$ takes the maximum. The dynamical
instability at $\beta=\pi/2$, therefore, occurs because it corresponds to an
energetically unstable point. Even infinitesimally small fluctuations of the
order parameter around $\beta=\pi/2$ grow and drive the system away from
$\beta=\pi/2$ towards the stable states at $\beta=0$ or $\pi$.
We note that the energy functional for the first-order phase transition also
exhibits the maximum at $\beta=\pi/2$ and minima at $\beta=0$ and $\pi$. The
origin of the dynamical instability is, therefore, the same as explained
above.

Our results do not exclude the possibility of stable spin currents in general.
We expect that spin currents can be stabilized by the quadratic Zeeman effect.
In fact, stable spin currents have been observed in a spin-1 antiferromagnetic
Bose-Einstein condensate (BEC) of 23Na atoms in the presence of the quadratic
Zeeman effect [41], despite the fact that spin currents are predicted to be
unstable in this system without the quadratic Zeeman effect [54, 55].

In the above argument, we have revealed that the instability of spin currents
originates from the fact that the energy functional takes its maximum for a
superfluid state carrying a spin current (see Fig. 11) that results in the
dynamical instability of the spin nematic mode $\omega_{\beta}$. The spin
nematic mode $\omega_{\beta}$ involves small population of $m=0$ hyperfine
state as indicated in Eq. (102). In the presence of the quadratic Zeeman
effect, since populating $m=0$ state from the polar state in Eq. (55) costs
finite energy due to the quadratic Zeeman effect, it is expected that the spin
nematic mode may be stabilized and the energy functional takes minimum for a
state carrying a spin current. However, it is beyond the scope of the paper to
confirm this possibility using the TDGL equation, because it is necessary to
derive the TDGL equation taking into account the quadratic Zeeman effect from
scratch. Note that the insulating state in the presence of the quadratic
Zeeman effect is different from the spin-singlet state in Eq. (20). We examine
the possibility of stabilizing spin currents by the quadratic Zeeman effect
using the discrete Gross-Pitaevskii (DGP) equation, which is valid in the
weakly-interacting regime deep in the superfluid phase.

The spin-1 Bose-Hubbard model with the quadratic Zeeman effect [26] is given
by

$\displaystyle\hat{H}_{q}=\hat{H}-q\sum_{i,\alpha,\beta}\hat{b}_{i\alpha}^{\dagger}(F_{z}^{2})_{\alpha\beta}\hat{b}_{i\beta},$
(217)

where $q>0$. In the weakly interacting regime ($t\gg U_{0}$), introducing the
mean-field $\psi_{i\alpha}=\langle\hat{b}_{i\alpha}\rangle$, we obtain the
energy functional

$\displaystyle\langle\hat{H}_{q}\rangle$ $\displaystyle=$
$\displaystyle-q\sum_{i}(|\psi_{i1}|^{2}+|\psi_{i-1}|^{2})-t\sum_{\langle
i,j\rangle,\alpha}(\psi_{i\alpha}^{*}\psi_{j\alpha}+{\rm c.c.})$ (218)
$\displaystyle-$
$\displaystyle\mu\sum_{i,\alpha}|\psi_{i\alpha}|^{2}+\frac{U_{0}}{2}\sum_{i}(\sum_{\alpha}|\psi_{i\alpha}|^{2})^{2}$
$\displaystyle+$
$\displaystyle\frac{U_{2}}{2}\sum_{i}(\sum_{\alpha,\beta}\psi_{i\alpha}^{*}\bm{F}_{\alpha\beta}\psi_{i\beta})^{2}.$

From
$i\partial_{t}\psi_{i\alpha}=\delta\langle\hat{H}_{q}\rangle/\delta\psi_{i\alpha}^{*}$,
the discrete Gross-Pitaevskii (DGP) equation can be derived as

$\displaystyle i\partial_{t}\psi_{i\alpha}$ $\displaystyle=$
$\displaystyle-t\sum_{j\in
V(i)}\psi_{j\alpha}-q\psi_{i\alpha}(\delta_{\alpha,1}+\delta_{\alpha,-1})$
(219) $\displaystyle-$
$\displaystyle\mu\psi_{i\alpha}+U_{0}\psi_{i\alpha}\sum_{\beta}|\psi_{i\beta}|^{2}$
$\displaystyle+$ $\displaystyle
U_{2}\sum_{\beta,\gamma,\delta}\psi_{i\beta}^{*}\bm{F}_{\beta\gamma}\psi_{i\gamma}\cdot\bm{F}_{\alpha\delta}\psi_{i\delta},$

where $j\in V(i)$ denotes the summation over the nearest neighboring sites of
site $i$. We consider the static solution that corresponds to Eq. (80)

$\displaystyle\vec{\psi}_{i}^{0}(\beta)=\sqrt{n_{\rm
c}}\left[\begin{array}[]{c}-\frac{e^{-i\bm{p}\cdot\bm{r}_{i}}}{\sqrt{2}}\sin\beta\\\
\cos\beta\\\
\frac{e^{i\bm{p}\cdot\bm{r}_{i}}}{\sqrt{2}}\sin\beta\end{array}\right].$ (223)

Setting $\beta=\pi/2$, Eq. (223) reduces to

$\displaystyle\vec{\psi}_{i}^{\rm 0}(\pi/2)=\sqrt{\frac{n_{\rm
c}}{2}}\left[\begin{array}[]{c}-e^{-i\bm{p}\cdot\bm{r}_{i}}\\\ 0\\\
e^{i\bm{p}\cdot\bm{r}_{i}}\end{array}\right],$ (227)

which corresponds to Eq. (72). We introduce fluctuation of the order parameter
around $\vec{\psi}_{i}^{\rm 0}(\pi/2)$ as

$\displaystyle\vec{\psi}_{i}=\vec{\psi}_{i}^{0}(\pi/2)+\delta\vec{\psi}_{i},$
(228)

where

$\displaystyle\delta\vec{\psi}_{i}=\left[\begin{array}[]{c}\delta\psi_{i1}e^{-i\bm{p}\cdot\bm{r}_{i}}\\\
\delta\psi_{i0}\\\
\delta\psi_{i-1}e^{i\bm{p}\cdot\bm{r}_{i}}\end{array}\right].$ (232)

Substituting Eq. (228) into the DGP equation (219) and linearizing with
respect to fluctuation, we obtain

$\displaystyle\omega\left[\begin{array}[]{c}\varphi_{\bm{k},0}\\\
\varphi_{\bm{k},0}^{*}\end{array}\right]$ (235) $\displaystyle=$
$\displaystyle\left[\begin{array}[]{cc}\epsilon_{\bm{k}}-\epsilon_{\bm{p}}+U_{2}n_{\rm
c}-q&-U_{2}n_{\rm c}\\\ U_{2}n_{\rm
c}&-(\epsilon_{\bm{k}}-\epsilon_{\bm{p}}+U_{2}n_{\rm
c}-q)\end{array}\right]\left[\begin{array}[]{c}\varphi_{\bm{k},0}\\\
\varphi_{\bm{k},0}^{*}\end{array}\right],$ (240)

where $\varphi_{\bm{k},0}$ is the Fourier component of $\delta\psi_{i0}$ and
$\epsilon_{\bm{k}}$ is given in Eq. (C22). Solving Eq. (240), we obtain the
dispersion relation for the spin nematic mode corresponding to
$\omega_{\beta}$ as

$\displaystyle\tilde{\omega}_{\beta}(\bm{k},\bm{p})=\sqrt{(\epsilon_{\bm{k}}-\epsilon_{\bm{p}}+U_{2}n_{\rm
c}+q)^{2}-(U_{2}n_{\rm c})^{2}}.$ (242)

Without the quadratic Zeeman effect ($q=0$),
$\tilde{\omega}_{\beta}(\bm{k},\bm{p})$ becomes pure imaginary at $\bm{k}={\bf
0}$ for finite $\bm{p}$, which means that spin currents are dynamically
unstable. This demonstrates that the absence of spin supercurrents is not
restricted within the TDGL equation, which is valid in the vicinity of the
Mott insulating (MI) phase. Our results suggest that spin supercurrents are
unstable in the entire superfluid phase without the quadratic Zeeman effect.

In the presence of the quadratic Zeeman effect ($q>0$), one finds that
$\tilde{\omega}_{\beta}$ does not exhibit dynamical instability for $\bm{p}$
smaller than the critical momentum $\bm{p}_{\rm c}$. For small $\bm{p}$
($pa\ll 1$), the critical momentum $p_{\rm c}$ for the onset of the dynamical
instability is obtained as

$\displaystyle p_{\rm c}a=\sqrt{q/t}.$ (243)

The quadratic Zeeman effect thus stabilizes spin currents.

Substituting Eq. (223) into the energy functional (218), we obtain

$\displaystyle\frac{\langle\hat{H}_{q}\rangle}{n_{\rm c}N_{\rm
s}}=\frac{\epsilon_{\bm{p}}}{2}-zt\cos^{2}\beta-\mu-\frac{q}{2}\sin^{2}\beta.$
(244)

For small $\bm{p}$, Eq. (244) reduces to

$\displaystyle\frac{\langle\hat{H}_{q}\rangle}{n_{\rm c}N_{\rm s}}=-\mu-
zt-\frac{ta^{2}}{2}(p_{\rm c}^{2}-p^{2})\sin^{2}\beta.$ (245)

Equation (245) shows that $\beta=\pi/2$ corresponds to an energy minimum for
$p<p_{\rm c}$ and therefore the static solution (227) carrying a finite spin
current is stable as expected.

##  VI conclusion

To summarize, we have studied the stability of supercurrents in the polar
phase of antiferromagnetically interacting spin-1 bosons in an optical lattice
using the TDGL equation. We have calculated the critical momenta for
supercurrents in the vicinity of the MI phase with even filling factors. We
found that the critical momentum for mass currents is finite throughout the
metastable SF phase, which demonstrates the superfluidity of the metastable SF
state. The critical momentum for spin currents was found to be zero. We have
found that this instability of spin currents originates from the fact that the
polar state with a spin current corresponds to an energy saddle point, thereby
an infinitesimal spin current causes a dynamical instability.
Although we set $U_{2}/U_{0}$ optimal values for validity of the theory that
is not realistic in the present experimental status, the obtained results may
be applicable to ${}^{23}{\rm Na}$ with $U_{2}/U_{0}=0.04$. The first-order
phase transition dominates the larger part of the phase boundary and the
region of the metastable SF phase in the phase diagram gets larger for smaller
$U_{2}/U_{0}$ [14, 17]. The effect of the first-order phase transition,
therefore, may be remarkable and easy to study experimentally for ${}^{23}{\rm
Na}$.
Our predictions for the stability of supercurrents can be verified
experimentally with current technologies. As for mass currents, we propose a
similar setup with Ref. [33], where a mass current is induced by a moving
optical lattice. Spin currents can be induced by applying a linear magnetic
field gradient that induces counter flow of two spin components as in the
experiment for spinor gases [41].

###### Acknowledgements.

We thank R. Asaoka, I. Danshita, D. Kagamihara, Y. Kawaguchi, T. Kobayashi,
and D. Yamamoto for helpful discussions. The work of RY is supported by the
Japanese Society for the Promotion of Science Grant-in-Aid for Scientific
Research (KAKENHI Grants No. 19K14616 and No. 20H01838). This work is
supported by the Japanese Society for the Promotion of Science Grant-in-Aid
for Scientific Research (KAKENHI Grant No. 19K03691).

##  Appendix A MATRIX ELEMENTS OF CREATION AND ANNIHILATION OPERATORS

We present the details of the derivation of the matrix elements of creation
and annihilation operators summarized in Table I that are necessary for the
calculation of the ground-state energy (27). The following derivation is based
on Ref. [56].
If one operates $\hat{b}_{\alpha}^{\dagger}\ (\alpha=0,\pm 1)$ on the state
$|S,m,n\rangle$, the magnetic quantum number $m$ and particle number $n$
change as $m\rightarrow m+\alpha$ and $n\rightarrow n+1$. Meanwhile, the total
spin quantum number $S$ changes as $S\rightarrow S\pm 1$ by adding a spin-1
particle. Thus, $\hat{b}_{\alpha}^{\dagger}|S,m,n\rangle$ can be written as

$\displaystyle\hat{b}_{\alpha}^{\dagger}|S,m,n\rangle$ $\displaystyle=$
$\displaystyle M_{\alpha,S,m,n}|S+1,m+\alpha,n+1\rangle$ (246)
$\displaystyle+N_{\alpha,S,m,n}|S-1,m+\alpha,n+1\rangle.$

Analogously, $\hat{b}_{\alpha}|S,m,n\rangle$ can be written as

$\displaystyle\hat{b}_{\alpha}|S,m,n\rangle$ $\displaystyle=$ $\displaystyle
O_{\alpha,S,m,n}|S+1,m-\alpha,n-1\rangle$ (247)
$\displaystyle+P_{\alpha,S,m,n}|S-1,m-\alpha,n-1\rangle.$

We derive the recurrence formulas for $M_{\alpha,S,m,n},\ N_{\alpha,S,m,n},\
O_{\alpha,S,m,n},$ and $P_{\alpha,S,m,n}$. From $-S\leq m\leq S$, they satisfy

$\displaystyle N_{1,S,S,n}=N_{0,S,S,n}=P_{0,S,S,n}=P_{-1,S,S,n}=0.$ (248)

Then, we obtain

$\displaystyle\hat{b}_{1}^{\dagger}|S,S,n\rangle=M_{1,S,S,n}|S+1,S+1,n+1\rangle,$
(249)

where $M_{1,S,S,n}$ is given as [11]

$\displaystyle M_{1,S,S,n}=\sqrt{\frac{(S+1)(n+S+3)}{2S+3}}.$ (250)

Other matrix elements can be derived by operating the ladder operators
$\hat{S}^{\pm}=\hat{S}_{x}\pm i\hat{S}_{y}$. By simple calculations, one can
show

$\displaystyle[\hat{S}^{+},\hat{b}_{1}]=-\sqrt{2}\hat{b}_{0},\
[\hat{S}^{+},\hat{b}_{0}]=-\sqrt{2}\hat{b}_{-1},\
[\hat{S}^{+},\hat{b}_{-1}]=0,$ $\displaystyle[\hat{S}^{-},\hat{b}_{1}]=0,\
[\hat{S}^{-},\hat{b}_{0}]=-\sqrt{2}\hat{b}_{-1},\
[\hat{S}^{-},\hat{b}_{-1}]=-\sqrt{2}\hat{b}_{0},$
$\displaystyle[\hat{S}_{z},\hat{b}_{1}]=-\hat{b}_{1},\
[\hat{S}_{z},\hat{b}_{0}]=0,\ [\hat{S}_{z},\hat{b}_{-1}]=\hat{b}_{-1},$
$\displaystyle[\hat{b}_{1},\hat{\Theta}^{\dagger}]=-2\hat{b}^{\dagger}_{-1},\
[\hat{b}_{0},\hat{\Theta}^{\dagger}]=2\hat{b}_{0}^{\dagger},\
[\hat{b}_{-1},\hat{\Theta}^{\dagger}]=-2\hat{b}_{1}^{\dagger}.$

Using the above commutation relations, we obtain

$\displaystyle\hat{S}^{+}\hat{b}_{1}^{\dagger}|S,m,n\rangle$ $\displaystyle=$
$\displaystyle\sqrt{(S-m)(S+m+3)}M_{1,S,m,n}|S+1,m+2,n+1\rangle$
$\displaystyle+$
$\displaystyle\sqrt{(S-m-2)(S+m+1)}N_{1,S,m,n}|S-1,m+2,n+1\rangle$
$\displaystyle=$ $\displaystyle\hat{b}_{1}^{\dagger}\hat{S}^{+}|S,m,n\rangle$
$\displaystyle=$
$\displaystyle\sqrt{(S-m)(S+m+1)}M_{1,S,m+1,n}|S+1,m+2,n+1\rangle$
$\displaystyle+$
$\displaystyle\sqrt{(S-m)(S+m+1)}N_{1,S,m+1,n}|S-1,m+2,n+1\rangle.$

Comparing Eqs. (LABEL:eq:sb) and (LABEL:eq:bs), we obtain the recurrence
formula for $M_{1,S,m,m}$ as

$\displaystyle M_{1,S,m,n}=\sqrt{\frac{S+m+1}{S+m+3}}M_{1,S,m+1,n}.$ (253)

In the same manner, we can derive other recurrence formulas

$\displaystyle M_{0,S,m,n}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{(S-m)(S+m+1)}{(S-m+1)(S+m+2)}}M_{0,S,m+1,n}$ (254)
$\displaystyle+\sqrt{\frac{2}{(S-m+1)(S+m+2)}}M_{1,S,m,n},$ $\displaystyle
M_{-1,S,m,n}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{(S-m)(S+m+1)}{(S-m+2)(S+m+1)}}M_{-1,S,m+1,n}$
$\displaystyle+\sqrt{\frac{2}{(S-m+2)(S+m+1)}}M_{0,S,m,n},$ $\displaystyle
N_{-1,S,m,n}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{S(S+1)-m(m-1)}{S(S-1)-(m-1)(m-2)}}N_{-1,S,m-1,n},$
$\displaystyle N_{0,S,m,n}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{S(S+1)-m(m-1)}{S(S-1)-m(m-1)}}N_{-1,S,m-1,n}$ (257)
$\displaystyle+\sqrt{\frac{2}{S(S-1)-m(m-1)}}N_{-1,S,m,n},$ $\displaystyle
N_{1,S,m,n}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{S(S+1)-m(m-1)}{S(S-1)-m(m+1)}}N_{1,S,m-1,n}$ (258)
$\displaystyle+\sqrt{\frac{2}{S(S+1)-m(m-1)}}N_{0,S,m,n}.$

Here, $M_{\alpha,S,S,n}$ and $N_{\alpha,S,S,n}$ satisfy

$\displaystyle n$ $\displaystyle=$ $\displaystyle\sum_{\alpha}\langle
S,S,n|\hat{b}_{\alpha}^{\dagger}\hat{b}_{\alpha}|S,S,n\rangle$ (259)
$\displaystyle=$
$\displaystyle\sum_{\alpha}[M_{\alpha,S,S,n}^{2}+N_{\alpha,S,S,n}^{2}]-3.$

From Eq. (259) and $N_{1,S,S,n}=N_{0,S,S,n}=0$, the initial value of
$N_{\alpha,S,m,n}$ is given by

$\displaystyle N_{-1,S,S,n}=-\sqrt{3+n-\sum_{\alpha}M^{2}_{\alpha,S,S,n}}.$
(260)

The sign of $N_{-1,S,S,n}$ is set to be consistent with Ref. [11]. Solving the
recurrence formulas (254)-(258) and (260), one obtains Table 1 [57].
One can calculate $O_{\alpha,S,m,n}$ and $P_{\alpha,S,m,n}$ as [56]

$\displaystyle\langle S,m,n|\hat{b}^{\dagger}_{\alpha}|S-1,m-\alpha,n\rangle$
$\displaystyle=$ $\displaystyle M_{\alpha,S-1,m-\alpha,n-1}$ (261)
$\displaystyle=$ $\displaystyle P_{\alpha,S,m,n},$ $\displaystyle\langle
S,m,n|\hat{b}^{\dagger}_{\alpha}|S+1,m-\alpha,n\rangle$ $\displaystyle=$
$\displaystyle M_{\alpha,S+1,m-\alpha,n-1}$ (262) $\displaystyle=$
$\displaystyle O_{\alpha,S,m,n}.$

##  Appendix B PERTURBATIVE MEAN-FIELD CALCULATION

In this appendix, we present the details of the perturbative mean-field
calculation in Sec. II.2. Since the ground-state energy should be invariant
under the ${\rm U}(1)$ gauge transformation and spin rotations, it is
convenient to expand the ground-state energy by the $d$ vector that transforms
as a vector under spin rotations as

$\displaystyle E$ $\displaystyle=$ $\displaystyle E^{(2)}+E^{(4)}+E^{(6)},$
(263) $\displaystyle E^{(2)}$ $\displaystyle=$ $\displaystyle
C_{2}(\bm{d}^{\dagger}\bm{d})=C_{2}n_{\rm c},$ (264) $\displaystyle E^{(4)}$
$\displaystyle=$ $\displaystyle
D_{4}(\bm{d}^{\dagger}\bm{d})^{2}+D_{4}^{\prime}|\bm{d}\cdot\bm{d}|^{2}$ (265)
$\displaystyle=$ $\displaystyle C_{4}n_{\rm c}^{2}+C_{4}^{\prime}n_{\rm
c}^{2}\langle\bm{F}\rangle^{2},$ $\displaystyle E^{(6)}$ $\displaystyle=$
$\displaystyle
D_{6}(\bm{d}^{\dagger}\bm{d})^{3}+D_{6}^{\prime}|\bm{d}\cdot\bm{d}|^{2}(\bm{d}^{\dagger}\bm{d})$
(266) $\displaystyle=$ $\displaystyle C_{6}n_{\rm c}^{3}+C_{6}^{\prime}n_{\rm
c}^{3}\langle\bm{F}\rangle^{2},$

where $C_{m}=D_{m}+D_{m}^{\prime}$, $C_{m}^{\prime}=-D_{m}^{\prime}\ (m=4,6)$.
Note that $\bm{d}^{\dagger}\bm{d}=n_{\rm c},\ \bm{d}\cdot\bm{d}=n_{\rm
c}(\zeta_{0}^{2}-2\zeta_{1}\zeta_{-1}),$ and $|\bm{d}\cdot\bm{d}|^{2}=n_{\rm
c}^{2}(1-\langle\bm{F}\rangle^{2})$. To reduce the number of terms that need
to be calculated, we set $\zeta_{0}=0$ in Eqs. (264), (265), (266), and (23).
We thus obtain

$\displaystyle E^{(2)}$ $\displaystyle=$ $\displaystyle C_{2}n_{\rm
c}(|\zeta_{1}|^{2}+|\zeta_{-1}|^{2}),$ (267) $\displaystyle E^{(4)}$
$\displaystyle=$ $\displaystyle D_{4}n_{\rm
c}^{2}(|\zeta_{1}|^{4}+|\zeta_{-1}|^{4})$ (268)
$\displaystyle+(2D_{4}+4D_{4}^{\prime})n_{\rm
c}^{2}|\zeta_{1}|^{2}|\zeta_{-1}|^{2},$ $\displaystyle E^{(6)}$
$\displaystyle=$ $\displaystyle D_{6}n_{\rm
c}^{3}(|\zeta_{1}|^{6}+|\zeta_{-1}|^{6})$
$\displaystyle+(3D_{6}+4D_{6}^{\prime})n_{\rm
c}^{3}|\zeta_{1}|^{2}|\zeta_{-1}|^{2}(|\zeta_{1}|^{2}+|\zeta_{-1}|^{2}),$
$\displaystyle\hat{V}$ $\displaystyle=$ $\displaystyle-
zt(\psi_{1}\hat{b}^{\dagger}_{1}+\psi_{-1}\hat{b}^{\dagger}_{-1}+\mathrm{h.c.}).$
(270)

We employ the formulas of the standard perturbation theory [58]:

$\displaystyle E^{(2)}$ $\displaystyle=$ $\displaystyle-\sum_{\alpha\neq
g}\frac{|\langle\alpha|\hat{V}^{\mathrm{pert}}|g\rangle|^{2}}{\Delta
E_{\alpha}},$ (271) $\displaystyle E^{(4)}$ $\displaystyle=$
$\displaystyle-\sum_{\alpha\beta\gamma\neq g}\langle
g|\hat{V}^{\mathrm{pert}}|\alpha\rangle\frac{\langle\alpha|\hat{V}^{\mathrm{pert}}|\beta\rangle}{\Delta
E_{\alpha}}\frac{\langle\beta|\hat{V}^{\mathrm{pert}}|\gamma\rangle}{\Delta
E_{\beta}}\frac{\langle\gamma|\hat{V}^{\mathrm{pert}}|g\rangle}{\Delta
E_{\gamma}}-E^{(2)}\sum_{\alpha\neq
g}\frac{|\langle\alpha|\hat{V}^{\mathrm{pert}}|g\rangle|^{2}}{\Delta
E_{\alpha}^{2}},$ (272) $\displaystyle E^{(6)}$ $\displaystyle=$
$\displaystyle-\sum_{\alpha\beta\gamma\delta\epsilon\neq g}\langle
g|\hat{V}^{\mathrm{pert}}|\alpha\rangle\frac{\langle\alpha|\hat{V}^{\mathrm{pert}}|\beta\rangle}{\Delta
E_{\alpha}}\frac{\langle\beta|\hat{V}^{\mathrm{pert}}|\gamma\rangle}{\Delta
E_{\beta}}\frac{\langle\gamma|\hat{V}^{\mathrm{pert}}|\delta\rangle}{\Delta
E_{\gamma}}\frac{\langle\delta|\hat{V}^{\mathrm{pert}}|\epsilon\rangle}{\Delta
E_{\delta}}\frac{\langle\epsilon|\hat{V}^{\mathrm{pert}}|g\rangle}{\Delta
E_{\epsilon}}$ (273) $\displaystyle+$
$\displaystyle\sum_{\alpha\beta\gamma\delta\neq g}\langle
g|\hat{V}^{\mathrm{pert}}|\alpha\rangle\frac{\langle\alpha|\hat{V}^{\mathrm{pert}}|\beta\rangle}{\Delta
E_{\alpha}}\frac{\langle\beta|\hat{V}^{\mathrm{pert}}|\gamma\rangle}{\Delta
E_{\beta}}\frac{\langle\gamma|\hat{V}^{\mathrm{pert}}|\delta\rangle}{\Delta
E_{\gamma}}\frac{\langle\delta|\hat{V}^{\mathrm{pert}}|g\rangle}{\Delta
E_{\delta}}$ $\displaystyle\left(\frac{1}{\Delta E_{\alpha}}+\frac{1}{\Delta
E_{\beta}}+\frac{1}{\Delta E_{\gamma}}+\frac{1}{\Delta
E_{\delta}}\right)-\sum_{\alpha\neq
g}\frac{|\langle\alpha|\hat{V}^{\mathrm{pert}}|g\rangle|^{6}}{\Delta
E_{\alpha}^{5}},$

where $|g\rangle$ is an unperturbative state, $\hat{V}^{\rm pert}$ represents
the perturbation,
$|\alpha\rangle,\,|\beta\rangle,\,|\gamma\rangle,\,|\delta\rangle,$ and
$|\epsilon\rangle$ are the intermediate states, and $\Delta E_{\alpha},\
\Delta E_{\beta},\ \Delta E_{\gamma},\ \Delta E_{\delta}$, and $\Delta
E_{\epsilon}$ denote the excitation energy for $|\alpha\rangle,\
|\beta\rangle,\ |\gamma\rangle,\ |\delta\rangle,$ and $|\epsilon\rangle$ from
$|g\rangle$. Using the formulas (271), (272), (273), and Table 1, the
coefficients $C_{2},\ D_{4},\ D_{4}^{\prime},\ D_{6},$ and $D_{6}^{\prime}$
can be calculated as

$\displaystyle C_{2}$ $\displaystyle=$
$\displaystyle(zt)^{2}-\frac{(zt)^{2}}{3}\left[\frac{n+3}{\Delta
E_{1,n+1}}+\frac{n}{\Delta E_{1,n-1}}\right],$ (274) $\displaystyle D_{4}$
$\displaystyle=$
$\displaystyle\frac{(zt)^{4}}{9}\left[-\frac{6}{5}\frac{n(n-2)}{\Delta
E_{1,n-1}^{2}\Delta E_{2,n-2}}-\frac{1}{5}\frac{n(n+3)}{\Delta
E_{2,n}}\left(\frac{1}{\Delta E_{1,n+1}}+\frac{1}{\Delta
E_{1,n-1}}\right)^{2}\right.$ (275) $\displaystyle-$
$\displaystyle\left.\frac{6}{5}\frac{(n+3)(n+5)}{\Delta E_{1,n+1}^{2}\Delta
E_{2,n+2}}+\left(\frac{n+3}{\Delta E_{1,n+1}}+\frac{n}{\Delta
E_{1,n-1}}\right)\left(\frac{n+3}{\Delta E_{1,n+1}^{2}}+\frac{n}{\Delta
E_{1,n-1}^{2}}\right)\right],$ $\displaystyle D_{4}^{\prime}$ $\displaystyle=$
$\displaystyle-\frac{(zt)^{4}}{9}\left[\frac{n(n+1)}{\Delta
E_{1,n-1}^{2}\Delta E_{0,n-2}}+\frac{(n+2)(n+3)}{\Delta E_{1,n+1}^{2}\Delta
E_{0,n+2}}-\frac{2}{5}\frac{n(n-2)}{\Delta E_{1,n-1}^{2}\Delta
E_{2,n-2}}\right.$ (276) $\displaystyle+$
$\displaystyle\frac{3}{5}\frac{n(n+3)}{\Delta E_{2,n}}\left(\frac{1}{\Delta
E_{1,n+1}}+\frac{1}{\Delta
E_{1,n-1}}\right)^{2}-\left.\frac{2}{5}\frac{(n+3)(n+5)}{\Delta
E_{1,n+1}^{2}\Delta E_{2,n+2}}\right],$ $\displaystyle D_{6}$ $\displaystyle=$
$\displaystyle\left[-\frac{2}{35}\frac{(n+7)(n+5)(n+3)}{\Delta E_{3,n+3}\Delta
E_{2,n+2}^{2}\Delta
E_{1,n+1}^{2}}\right.-\frac{2}{35}\frac{(n-4)(n-2)n}{\Delta E_{3,n-3}\Delta
E_{2,n-2}^{2}\Delta E_{1,n-1}^{2}}$ $\displaystyle-$
$\displaystyle\frac{2}{525}\frac{n(n+3)(n+5)}{\Delta
E_{3,n+1}}\left(\frac{1}{\Delta E_{2,n+2}\Delta E_{1,n+1}}+\frac{1}{\Delta
E_{2,n}\Delta E_{1,n+1}}+\frac{1}{\Delta E_{2,n}\Delta E_{1,n-1}}\right)^{2}$
$\displaystyle-$ $\displaystyle\frac{2}{525}\frac{n(n-2)(n+3)}{\Delta
E_{3,n-1}}\left(\frac{1}{\Delta E_{2,n-2}\Delta E_{1,n-1}}+\frac{1}{\Delta
E_{2,n}\Delta E_{1,n+1}}+\frac{1}{\Delta E_{2,n}\Delta E_{1,n-1}}\right)^{2}$
$\displaystyle-$ $\displaystyle\frac{1}{3}\frac{(n+3)}{\Delta
E_{1,n+1}}\left(\frac{2(n+5)}{5}\frac{1}{\Delta E_{2,n+2}\Delta
E_{1,n+1}}+\frac{n}{15}\frac{1}{\Delta E_{2,n}\Delta
E_{1,n+1}}+\frac{n}{15}\frac{1}{\Delta E_{2,n}\Delta E_{1,n-1}}\right)^{2}$
$\displaystyle-$ $\displaystyle\frac{1}{3}\frac{n}{\Delta
E_{1,n-1}}\left(\frac{(n+3)}{15}\frac{1}{\Delta E_{2,n}\Delta
E_{1,n+1}}+\frac{(n+3)}{15}\frac{1}{\Delta E_{2,n}\Delta
E_{1,n-1}}+\frac{2(n-2)}{5}\frac{1}{\Delta E_{2,n-2}\Delta
E_{1,n-1}}\right)^{2}$ $\displaystyle+$
$\displaystyle\left\\{\frac{2}{15}\frac{(n+3)(n+5)}{\Delta E_{2,n+2}\Delta
E_{1,n+1}^{2}}\right.+\frac{1}{45}\frac{n(n+3)}{\Delta
E_{2,n}}\left(\frac{1}{\Delta E_{1,n+1}}+\frac{1}{\Delta
E_{1,n-1}}\right)^{2}$ $\displaystyle+$
$\displaystyle\left.\frac{2}{15}\frac{(n-2)n}{\Delta E_{2,n-2}\Delta
E_{1,n-1}^{2}}\right\\}\left(\frac{(n+3)}{3}\frac{1}{\Delta
E_{1,n+1}^{2}}+\frac{n}{3}\frac{1}{\Delta E_{1,n-1}^{2}}\right)$
$\displaystyle+$ $\displaystyle\left\\{\frac{4}{15}\frac{(n+3)(n+5)}{\Delta
E_{2,n+2}\Delta E_{1,n+1}^{3}}\right.+\frac{2}{45}\frac{n(n+3)}{\Delta
E_{2,n}}\left(\frac{1}{\Delta E_{1,n+1}^{2}}+\frac{1}{\Delta
E_{1,n-1}^{2}}\right)\left(\frac{1}{\Delta E_{1,n+1}}+\frac{1}{\Delta
E_{1,n-1}}\right)$ $\displaystyle+$
$\displaystyle\left.\frac{4}{15}\frac{(n-2)n}{\Delta E_{2,n-2}\Delta
E_{1,n-1}^{3}}\right\\}\left(\frac{(n+3)}{3\Delta E_{1,n+1}}+\frac{n}{3\Delta
E_{1,n-1}}\right)$ $\displaystyle+$
$\displaystyle\left\\{\frac{2}{15}\frac{(n+3)(n+5)}{\Delta E_{2,n+2}^{2}\Delta
E_{1,n+1}^{2}}+\frac{1}{45}\frac{n(n+3)}{\Delta
E_{2,n}^{2}}\left(\frac{1}{\Delta E_{1,n+1}}+\frac{1}{\Delta
E_{1,n-1}}\right)^{2}\right.$ $\displaystyle+$
$\displaystyle\left.\frac{2}{15}\frac{(n-2)n}{\Delta E_{1,n-1}^{2}\Delta
E_{2,n-2}^{2}}\right\\}\left(\frac{(n+3)}{3\Delta E_{1,n+1}}+\frac{n}{3\Delta
E_{1,n-1}}\right)$ $\displaystyle-$ $\displaystyle\left(\frac{(n+3)}{3\Delta
E_{1,n+1}^{3}}+\frac{n}{3\Delta
E_{1,n-1}^{3}}\right)\left(\frac{(n+3)}{3\Delta E_{1,n+1}}+\frac{n}{3\Delta
E_{1,n-1}}\right)^{2}$ $\displaystyle-$
$\displaystyle\left.\left(\frac{(n+3)}{3\Delta E_{1,n+1}}+\frac{n}{3\Delta
E_{1,n-1}}\right)\left(\frac{n+3}{3\Delta E_{1,n+1}^{2}}+\frac{n}{3\Delta
E_{1,n-1}^{2}}\right)^{2}\right](zt)^{6},$ $\displaystyle D_{6}^{\prime}$
$\displaystyle=$
$\displaystyle\frac{(zt)^{6}}{4}\left[-\frac{46}{525}\frac{n(n+3)(n+5)}{\Delta
E_{3,n+1}}\left(\frac{1}{\Delta E_{1,n-1}\Delta E_{2,n}}+\frac{1}{\Delta
E_{1,n+1}\Delta E_{2,n}}+\frac{1}{\Delta E_{1,n+1}\Delta
E_{2,n+2}}\right)^{2}\right.$ $\displaystyle-$
$\displaystyle\frac{46}{525}\frac{(n-2)n(n+3)}{\Delta
E_{3,n-1}}\left(\frac{1}{\Delta E_{1,n-1}\Delta E_{2,n}}+\frac{1}{\Delta
E_{1,n+1}\Delta E_{2,n}}+\frac{1}{\Delta E_{1,n-1}\Delta
E_{2,n-2}}\right)^{2}$ $\displaystyle-$
$\displaystyle\frac{18}{525}\frac{(n+3)(n+5)(n+7)}{\Delta E_{3,n+3}\Delta
E_{1,n+1}^{2}\Delta E_{2,n+2}^{2}}-\frac{18}{525}\frac{(n-4)(n-2)n}{\Delta
E_{3,n-3}\Delta E_{1,n-1}^{2}\Delta E_{2,n-2}^{2}}$ $\displaystyle-$
$\displaystyle\frac{1}{3}\frac{(n+2)(n+3)(n+5)}{\Delta
E_{1,n+3}}\left(\frac{8}{15\Delta E_{2,n+2}\Delta E_{1,n+1}}+\frac{2}{3\Delta
E_{0,n+2}\Delta E_{1,n+1}}\right)^{2}$ $\displaystyle-$
$\displaystyle\frac{2}{3}\frac{(n+3)}{\Delta
E_{1,n+1}}\left(\frac{7n}{15\Delta E_{2,n}\Delta E_{1,n+1}}+\frac{7n}{15\Delta
E_{2,n}\Delta E_{1,n-1}}+\frac{2(n+5)}{15\Delta E_{2,n+2}\Delta
E_{1,n+1}}+\frac{2(n+2)}{3\Delta E_{0,n+2}\Delta E_{1,n+1}}\right)$
$\displaystyle\left(\frac{2(n+5)}{5\Delta E_{2,n+2}\Delta
E_{1,n+1}}+\frac{n}{15\Delta E_{2,n}\Delta E_{1,n+1}}+\frac{n}{15\Delta
E_{2,n}\Delta E_{1,n-1}}\right)$ $\displaystyle-$
$\displaystyle\frac{1}{3}\frac{n}{\Delta
E_{1,n-1}}\left(\frac{7(n+3)}{15\Delta E_{2,n}\Delta
E_{1,n+1}}+\frac{7(n+3)}{15\Delta E_{2,n}\Delta
E_{1,n-1}}+\frac{2(n-2)}{15\Delta E_{2,n-2}\Delta
E_{1,n-1}}+\frac{2(n+1)}{3\Delta E_{0,n-2}\Delta E_{1,n-1}}\right)^{2}$
$\displaystyle-$ $\displaystyle\frac{1}{3}\frac{(n+3)}{\Delta
E_{1,n+1}}\left(\frac{2(n+5)}{15\Delta E_{2,n+2}\Delta
E_{1,n+1}}+\frac{2(n+2)}{3\Delta E_{0,n+2}\Delta E_{1,n+1}}+\frac{7n}{15\Delta
E_{2,n}\Delta E_{1,n-1}}+\frac{7n}{15\Delta E_{2,n}\Delta
E_{1,n+1}}\right)^{2}$ $\displaystyle-$
$\displaystyle\frac{2}{3}\frac{n}{\Delta E_{1,n-1}}\left(\frac{(n+3)}{15\Delta
E_{2,n}\Delta E_{1,n+1}}+\frac{(n+3)}{15\Delta E_{2,n}\Delta
E_{1,n-1}}+\frac{2(n-2)}{5\Delta E_{2,n-2}\Delta E_{1,n-1}}\right)$
$\displaystyle\left(\frac{7(n+3)}{15\Delta E_{2,n}\Delta
E_{1,n+1}}+\frac{7(n+3)}{15\Delta E_{2,n}\Delta
E_{1,n-1}}+\frac{2(n-2)}{15\Delta E_{2,n-2}\Delta
E_{1,n-1}}+\frac{2(n+1)}{3\Delta E_{0,n-2}\Delta E_{1,n-1}}\right)$
$\displaystyle-$ $\displaystyle\frac{1}{3}\frac{(n-2)n(n+1)}{\Delta
E_{1,n-3}}\left(\frac{8}{15\Delta E_{2,n-2}\Delta E_{1,n-1}}+\frac{2}{3\Delta
E_{0,n-2}\Delta E_{1,n-1}}\right)^{2}$ $\displaystyle+$
$\displaystyle\left\\{\frac{4}{9}\frac{n(n+1)}{\Delta E_{1,n-1}^{2}\Delta
E_{0,n-2}}+\frac{4}{9}\frac{(n+2)(n+3)}{\Delta E_{1,n+1}^{2}\Delta
E_{0,n+2}}+\frac{2}{9}\frac{n(n-2)}{\Delta E_{1,n-1}^{2}\Delta
E_{2,n-2}}\right.$ $\displaystyle+$
$\displaystyle\frac{1}{3}\frac{n(n+3)}{\Delta E_{2,n}}\left(\frac{1}{\Delta
E_{1,n+1}}+\frac{1}{\Delta
E_{1,n-1}}\right)^{2}\left.+\frac{2}{9}\frac{(n+3)(n+5)}{\Delta
E_{1,n+1}^{2}\Delta E_{2,n+2}}\right\\}\left(\frac{(n+3)}{3\Delta
E_{1,n+1}^{2}}+\frac{n}{3\Delta E_{1,n-1}^{2}}\right)$ $\displaystyle+$
$\displaystyle 2\left\\{\frac{4}{9}\frac{n(n+1)}{\Delta E_{1,n-1}^{3}\Delta
E_{0,n-2}}+\frac{4}{9}\frac{(n+2)(n+3)}{\Delta E_{1,n+1}^{3}\Delta
E_{0,n+2}}+\frac{2}{9}\frac{n(n-2)}{\Delta E_{1,n-1}^{3}\Delta
E_{2,n-2}}+\frac{2}{9}\frac{(n+3)(n+5)}{\Delta E_{1,n+1}^{3}\Delta
E_{2,n+2}}\right.$ $\displaystyle+$
$\displaystyle\frac{1}{3}\frac{n(n+3)}{\Delta E_{2,n}}\left(\frac{1}{\Delta
E_{1,n+1}}+\frac{1}{\Delta E_{1,n-1}}\right)\left.\left(\frac{1}{\Delta
E_{1,n+1}^{2}}+\frac{1}{\Delta
E_{1,n-1}^{2}}\right)\right\\}\left(\frac{(n+3)}{3\Delta
E_{1,n+1}}+\frac{n}{3\Delta E_{1,n-1}}\right)$ $\displaystyle+$
$\displaystyle\left\\{\frac{4}{9}\frac{n(n+1)}{\Delta E_{1,n-1}^{2}\Delta
E_{0,n-2}^{2}}+\frac{4}{9}\frac{(n+2)(n+3)}{\Delta E_{1,n+1}^{2}\Delta
E_{0,n+2}^{2}}+\frac{2}{9}\frac{n(n-2)}{\Delta E_{1,n-1}^{2}\Delta
E_{2,n-2}^{2}}\right.$ $\displaystyle+$
$\displaystyle\frac{1}{3}\frac{n(n+3)}{\Delta
E_{2,n}^{2}}\left.\left(\frac{1}{\Delta E_{1,n+1}}+\frac{1}{\Delta
E_{1,n-1}}\right)^{2}+\frac{2}{9}\frac{(n+3)(n+5)}{\Delta E_{1,n+1}^{2}\Delta
E_{2,n+2}^{2}}\right\\}\left(\frac{(n+3)}{3\Delta E_{1,n+1}}+\frac{n}{3\Delta
E_{1,n-1}}\right)$ $\displaystyle-$ $\displaystyle 3\left(\frac{(n+3)}{3\Delta
E_{1,n+1}^{3}}+\frac{n}{3\Delta
E_{1,n-1}^{3}}\right)\left(\frac{(n+3)}{3\Delta E_{1,n+1}}+\frac{n}{3\Delta
E_{1,n-1}}\right)^{2}$ $\displaystyle-$ $\displaystyle
3\left.\left(\frac{(n+3)}{3\Delta E_{1,n+1}}+\frac{n}{3\Delta
E_{1,n-1}}\right)\left(\frac{n+3}{3\Delta E_{1,n+1}^{2}}+\frac{n}{3\Delta
E_{1,n-1}^{2}}\right)^{2}\right]-\frac{3}{4}D_{6}.$

Here, $\Delta E_{S,n+a}=E^{0}(S,n+a)-E^{0}(0,n)$.

##  Appendix C DERIVATION OF THE TDGL EQUATION

In this Appendix, we summarize the derivation of the TDGL equation (303).
Using the coherent-state path integral [59], the grand partition function for
the spin-1 Bose-Hubbard model can be written as

$\displaystyle\Xi=\int\prod_{\alpha}\mathcal{D}b^{*}_{i\alpha}\mathcal{D}b_{i\alpha}e^{-S(\\{b_{i\alpha}\\})},$
(277)

where the action $S(\\{b_{i\alpha}\\})$ is given by

$\displaystyle S(\\{b_{i\alpha}\\})$ $\displaystyle=$
$\displaystyle\int^{\beta}_{0}{\rm
d}\tau\left[\sum_{i,\alpha}(b_{i\alpha}^{*}\partial_{\tau}b_{i\alpha}-\mu
b_{i\alpha}^{*}b_{i\alpha})\right.$ (278)
$\displaystyle-\left.\sum_{i,j,\alpha}t_{ij}b^{*}_{i\alpha}b_{j\alpha}+\frac{U_{0}}{2}\sum_{i,\alpha,\beta}b^{*}_{i\alpha}b^{*}_{i\beta}b_{i\beta}b_{i\alpha}\right.$
$\displaystyle+\left.\frac{U_{2}}{2}\sum_{i,\alpha,\beta,\gamma,\delta}b_{i\alpha}^{*}b_{i\gamma}^{*}\bm{F}_{\alpha\beta}\cdot\bm{F}_{\gamma\delta}b_{i\beta}b_{i\delta}\right].$

Here, the creation and annihilation operators in Eq. (1) are replaced by the
$c$-number field $b_{i\alpha}(\tau)$ that depends on imaginary time $\tau$ and
site $i$. $\beta=1/T$ is the inverse temperature and $t_{ij}$ is an element of
the hopping matrix $\hat{t}$. For the nearest-neighbor hopping term, $t_{ij}$
is given as

$\displaystyle t_{ij}=\left\\{\begin{array}[]{l}t,\ \ \mathrm{if\ {\it i}\
and\ {\it j}\ are\ nearest\ neighboring},\\\ 0,\ \
\mathrm{otherwise}.\end{array}\right.$ (281)

We introduce the auxiliary fields $\phi_{i\alpha}$ and $\phi_{i\alpha}^{*}$
via the Hubbard-Stratonovich transformation [60]:

$\displaystyle\int\mathcal{D}\phi^{*}_{i\alpha}\mathcal{D}\phi_{i\alpha}\exp\left[-\int{\rm
d}\tau(\vec{\phi}_{\alpha}^{\dagger}-\vec{b}_{\alpha}^{\dagger}\hat{t})(\hat{t})^{-1}(\vec{\phi}_{\alpha}-\hat{t}\vec{b}_{\alpha})\right]$
$\displaystyle=\mathrm{const},$ (282)

where $\vec{\phi}_{\alpha}=(\\{\phi_{i\alpha}\\})^{T}$ and
$\vec{b}_{\alpha}=(\\{b_{i\alpha}\\})^{T}$. Multiplying Eq. (282) to Eq.
(278), we obtain

$\displaystyle\Xi=\int\prod_{\alpha}\mathcal{D}\phi_{i\alpha}^{*}\mathcal{D}\phi_{i\alpha}\prod_{\alpha}\mathcal{D}b^{*}_{i\alpha}\mathcal{D}b_{i\alpha}e^{-S(\\{b_{i\alpha}\\},\\{\phi_{i\alpha}\\})},$

where

$\displaystyle S$ $\displaystyle=$ $\displaystyle\int^{\beta}_{0}{\rm
d}\tau\sum_{i,j,\alpha}\phi_{i\alpha}^{*}(\hat{t}^{-1})_{ij}\phi_{j\alpha}+S^{0}+S^{\mathrm{pert}},$
(284) $\displaystyle S^{0}$ $\displaystyle=$
$\displaystyle\int^{\beta}_{0}{\rm
d}\tau\left[\sum_{i,\alpha}(b^{*}_{i\alpha}\partial_{\tau}b_{i\alpha}-\mu
b^{*}_{i\alpha}b_{i\alpha})\right.$ (285)
$\displaystyle+\frac{U_{0}}{2}\sum_{i,\alpha,\beta}b^{*}_{i\alpha}b^{*}_{i\beta}b_{i\beta}b_{i\alpha}$
$\displaystyle\left.+\frac{U_{2}}{2}\sum_{i,\alpha,\beta,\gamma,\delta}b^{*}_{i\alpha}b^{*}_{i\gamma}\bm{F}_{\alpha\beta}\cdot\bm{F}_{\gamma\delta}b_{i\beta}b_{i\delta}\right],$
$\displaystyle S^{\mathrm{pert}}$ $\displaystyle=$
$\displaystyle-\int^{\beta}_{0}{\rm
d}\tau\sum_{i,\alpha}(\phi_{i\alpha}b^{*}_{i\alpha}+\mathrm{c.c.}).$ (286)

Integrating out $b_{i\alpha}$ and $b^{*}_{i\alpha}$, we obtain

$\displaystyle\Xi=\int\prod_{\alpha}\mathcal{D}\phi_{i\alpha}^{*}\mathcal{D}\phi_{i\alpha}e^{-S^{\mathrm{eff}}(\\{\phi_{i\alpha}\\})},$
(287)

where the effective action $S^{\mathrm{eff}}$ is given by

$\displaystyle S^{\mathrm{eff}}$ $\displaystyle=$ $\displaystyle
E^{0}(0,n)N_{\rm s}\beta-\log\langle e^{-S^{\mathrm{pert}}}\rangle_{0}$ (288)
$\displaystyle+\int^{\beta}_{0}{\rm
d}\tau\sum_{i,j,\alpha}(\hat{t}^{-1})_{ij}\phi^{*}_{i\alpha}\phi_{j\alpha}.$

Here, $\langle\hat{\cal O}\rangle_{0}$ denotes the expectation value of the
operator $\hat{\cal O}$ with respect to $S^{0}$:

$\displaystyle\langle\hat{\mathcal{O}}\rangle_{0}\equiv\frac{\int\prod_{\alpha}\mathcal{D}b_{i\alpha}\mathcal{D}b^{*}_{i\alpha}\hat{\mathcal{O}}e^{-S^{0}}}{\int\prod_{\alpha}\mathcal{D}b_{i\alpha}\mathcal{D}b^{*}_{i\alpha}e^{-S^{0}}}.$
(289)

We perform the cumulant expansion in the second term in Eq. (288) as [61]

$\displaystyle\log\langle e^{-S^{\mathrm{pert}}}\rangle_{0}$
$\displaystyle\sim$
$\displaystyle\frac{1}{2!}\langle(S^{\mathrm{pert}})^{2}\rangle_{\rm
c}+\frac{1}{4!}\langle(S^{\mathrm{pert}})^{4}\rangle_{\rm c}$ (290)
$\displaystyle+\frac{1}{6!}\langle(S^{\mathrm{pert}})^{6}\rangle_{\rm c},$

where $\langle\hat{\mathcal{O}}\rangle_{\rm c}$ denotes the cumulant average
of $\hat{\mathcal{O}}$ defined as

$\displaystyle\log\langle
e^{\hat{\mathcal{O}}}\rangle_{0}=\sum_{n=1}^{\infty}\frac{\langle(\hat{\mathcal{O}})^{n}\rangle_{\rm
c}}{n!}.$ (291)

We note that, since Eq. (284) is invariant under the U(1) gauge transformation
$b_{i\alpha}\rightarrow e^{i\theta}b_{i\alpha}$ and $\phi_{i\alpha}\rightarrow
e^{i\theta}\phi_{i\alpha}$, all the odd-order terms in the cumulant expansion
(290) vanish. Each term in the right-hand side of Eq. (290) can be written as

$\displaystyle\langle(S^{\mathrm{pert}})^{2}\rangle_{\rm c}$ $\displaystyle=$
$\displaystyle\langle(S^{\mathrm{pert}})^{2}\rangle_{0},$ (292)
$\displaystyle\langle(S^{\mathrm{pert}})^{4}\rangle_{\rm c}$ $\displaystyle=$
$\displaystyle\langle(S^{\mathrm{pert}})^{4}\rangle_{0}-3\langle(S^{\mathrm{pert}})^{2}\rangle_{0}^{2},$
(293) $\displaystyle\langle(S^{\mathrm{pert}})^{6}\rangle_{\rm c}$
$\displaystyle=$
$\displaystyle\langle(S^{\mathrm{pert}})^{6}\rangle_{0}-15\langle(S^{\mathrm{pert}})^{2}\rangle_{0}\langle(S^{\mathrm{pert}})^{4}\rangle_{0}$
(294) $\displaystyle+30\langle(S^{\mathrm{pert}})^{2}\rangle^{3}_{0}.$

To describe the critical behavior of the SF phase, it is sufficient to retain
the first and second-order time-derivative terms [59]. Using Eqs. (292)-(294)
and taking the zero-temperature limit $\beta\rightarrow\infty$, we obtain

$\displaystyle S^{\mathrm{eff}}$ $\displaystyle\sim$
$\displaystyle\int^{\infty}_{0}{\rm d}\tau\left\\{E^{0}(0,n)N_{\rm
s}+\sum_{i,j,\alpha}(\hat{t}^{-1})_{ij}\phi_{i\alpha}^{*}\phi_{j\alpha}\right.$
(295) $\displaystyle+$
$\displaystyle\sum_{i,\alpha}\left[K\frac{\phi_{i\alpha}^{*}}{zt}\frac{\partial_{\tau}\phi_{i\alpha}}{zt}+J\left|\frac{\partial_{\tau}\phi_{i\alpha}^{*}}{zt}\right|^{2}+(C_{2}-zt)\left|\frac{\phi_{i\alpha}}{zt}\right|^{2}+C_{4}\left|\frac{\phi_{i\alpha}}{zt}\right|^{4}+C_{6}\left|\frac{\phi_{i\alpha}}{zt}\right|^{6}\right]$
$\displaystyle+$
$\displaystyle\left.\sum_{i}\left[C_{4}^{\prime}\left(\sum_{\alpha,\beta}\frac{\phi_{i\alpha}^{*}}{zt}\bm{F}_{\alpha\beta}\frac{\phi_{i\beta}}{zt}\right)^{2}+C_{6}^{\prime}\left(\sum_{\alpha,\beta}\frac{\phi_{i\alpha}^{*}}{zt}\bm{F}_{\alpha\beta}\frac{\phi_{i\beta}}{zt}\right)^{2}\sum_{\alpha}\left|\frac{\phi_{i\alpha}}{zt}\right|^{2}\right]\right\\},$

where $K$ and $J$ are given by

$\displaystyle K=\frac{(zt)^{2}}{3}\left[\frac{n+3}{\Delta
E_{1,n+1}^{2}}-\frac{n}{\Delta E_{1,n-1}^{2}}\right],$ (296) $\displaystyle
J=\frac{(zt)^{2}}{3}\left[\frac{n+3}{\Delta E_{1,n+1}^{3}}+\frac{n}{\Delta
E_{1,n-1}^{3}}\right].$ (297)

By the Fourier transform, the second term in Eq. (295) can be written as

$\displaystyle\int^{\infty}_{0}{\rm
d}\tau\sum_{i,j,\alpha}(\hat{t}^{-1})_{ij}\phi_{i\alpha}^{*}\phi_{j\alpha}$
$\displaystyle=$
$\displaystyle-\sum_{\bm{k},\omega,\alpha}\frac{|\bar{\phi}_{\alpha}(\bm{k},\omega)|^{2}}{\epsilon_{\bm{k}}}$
$\displaystyle\sim$
$\displaystyle\sum_{\bm{k},\omega,\alpha}\left|\frac{\bar{\phi}_{\alpha}(\bm{k},\omega)}{zt}\right|^{2}[zt+(ka)^{2}],$

where

$\displaystyle\phi_{i\alpha}(\tau)$ $\displaystyle=$
$\displaystyle\sum_{\bm{k},\omega}\bar{\phi}_{\alpha}(\bm{k},\omega)e^{i(\bm{k}\cdot\bm{r}_{i}-\omega\tau)},$
(299) $\displaystyle\epsilon_{\bm{k}}$ $\displaystyle=$
$\displaystyle-2t\sum_{l=1}^{d}\cos(k_{l}a).$ (300)

Here, $\bm{r}_{i}\equiv(x_{i}a,y_{i}a,z_{i}a)^{T}$ and $a$ is the lattice
constant. We have taken the long-wavelength limit, i.e., $ka\ll 1$ in Eq.
(LABEL:eq:kinetic). In the continuum limit that is effective in the vicinity
of the phase boundary, making replacements
$\phi_{i\alpha}(\tau)/(a^{d/2}zt)\rightarrow\Psi_{\alpha}(\tau,\bm{r})$,
$S^{\rm eff}$ reduces to the Ginzburg-Landau (GL) action $S^{\mathrm{GL}}$:

$\displaystyle S^{\mathrm{GL}}$ $\displaystyle=$
$\displaystyle\int^{\infty}_{0}{\rm d}\tau\int{\rm
d}^{d}x\left[E^{0}(0,n)+K{\bf\Psi}^{\dagger}\partial_{\tau}{\bf\Psi}+J(\partial_{\tau}{\bf\Psi}^{\dagger})(\partial_{\tau}{\bf\Psi})+\frac{1}{2m^{*}}(\nabla{\bf\Psi}^{\dagger})\cdot(\nabla{\bf\Psi})\right.$
(301)
$\displaystyle\left.+C_{2}({\bf\Psi}^{\dagger}{\bf\Psi})+c_{4}({\bf\Psi}^{\dagger}{\bf\Psi})^{2}+c_{6}({\bf\Psi}^{\dagger}{\bf\Psi})^{3}+c_{4}^{\prime}\langle\\!\langle\bm{F}\rangle\\!\rangle^{2}+c_{6}^{\prime}\langle\\!\langle\bm{F}\rangle\\!\rangle^{2}({\bf\Psi}^{\dagger}{\bf\Psi})\right].$

Here, $c_{4}=a^{d}C_{4},\
c_{4}^{\prime}=a^{d}C_{4}^{\prime},c_{6}=a^{2d}C_{6},\
c_{6}^{\prime}=a^{2d}C_{6}^{\prime},\
{\bf\Psi}=(\Psi_{1},\Psi_{0},\Psi_{-1})^{T}\ {\rm is\ the\ SF\ order\
parameter},\
\langle\\!\langle\bm{F}\rangle\\!\rangle=\sum_{\alpha,\beta}\Psi_{\alpha}^{*}\bm{F}_{\alpha\beta}\Psi_{\beta}\
{\rm is\ the\ spin\ average},$ and $m^{*}=1/(2ta^{2})$ is the effective mass.
Setting $\partial_{\tau}\Psi_{\alpha}=0$ in Eq. (301), the energy functional
for a static solution ${\bf\Psi}$ can be written as

$\displaystyle{\cal K}(\\{{\bf\Psi}\\})=\int{\rm d}^{d}x$
$\displaystyle\left[\frac{1}{2m^{*}}(\nabla{\bf\Psi}^{\dagger})\cdot(\nabla{\bf\Psi})+C_{2}({\bf\Psi}^{\dagger}{\bf\Psi})+c_{4}({\bf\Psi}^{\dagger}{\bf\Psi})^{2}+c_{6}({\bf\Psi}^{\dagger}{\bf\Psi})^{3}+c_{4}^{\prime}\langle\\!\langle{\bm{F}}\rangle\\!\rangle^{2}+c_{6}^{\prime}\langle\\!\langle{\bm{F}}\rangle\\!\rangle^{2}({\bf\Psi}^{\dagger}{\bf\Psi})\right].$

From $\frac{\delta S^{\rm GL}}{\delta\Psi_{\alpha}^{*}}=0$, we finally obtain
the TDGL equation

$\displaystyle iK\partial_{t}\Psi_{\alpha}-J\partial_{t}^{2}\Psi_{\alpha}$
$\displaystyle=$
$\displaystyle-\frac{\nabla^{2}}{2m^{*}}\Psi_{\alpha}+C_{2}\Psi_{\alpha}+2c_{4}({\bf\Psi}^{\dagger}{\bf\Psi})\Psi_{\alpha}+3c_{6}({\bf\Psi}^{\dagger}{\bf\Psi})^{2}\Psi_{\alpha}+c_{6}^{\prime}\langle\\!\langle\bm{F}\rangle\\!\rangle^{2}\Psi_{\alpha}$
(303)
$\displaystyle+2c_{4}^{\prime}\langle\\!\langle\bm{F}\rangle\\!\rangle\cdot\sum_{\beta}(\bm{F}_{\alpha\beta}\Psi_{\beta})+2c_{6}^{\prime}({\bf\Psi}^{\dagger}{\bf\Psi})\langle\\!\langle\bm{F}\rangle\\!\rangle\cdot\sum_{\beta}(\bm{F}_{\alpha\beta}\Psi_{\beta}).$

Note that Eq. (303) is valid for low energy dynamics, in which slow temporal
variation of ${\bf\Psi}$ is allowed. If $K\neq 0$, the first-order time-
derivative term is dominant and the second-order one should be neglected.

Table 1: Matrix elements $S$ | $m$ | $M_{1,S,m,n}$ | $M_{0,S,m,n}$ | $M_{-1,S,m,n}$ | $N_{1,S,m,n}$ | $N_{0,S,m,n}$ | $N_{-1,S,m,n}$
---|---|---|---|---|---|---|---
0 | 0 | $\sqrt{\frac{n+3}{3}}$ | $\sqrt{\frac{n+3}{3}}$ | $\sqrt{\frac{n+3}{3}}$ | 0 | 0 | 0
1 | 1 | $\sqrt{\frac{2(n+4)}{5}}$ | $\sqrt{\frac{n+4}{5}}$ | $\sqrt{\frac{n+4}{15}}$ | 0 | 0 | $-\sqrt{\frac{n+1}{3}}$
1 | 0 | $\sqrt{\frac{n+4}{5}}$ | $2\sqrt{\frac{n+4}{15}}$ | $\sqrt{\frac{n+4}{5}}$ | 0 | $\sqrt{\frac{n+1}{3}}$ | 0
1 | -1 | $\sqrt{\frac{n+4}{15}}$ | $\sqrt{\frac{n+4}{5}}$ | $\sqrt{\frac{2(n+4)}{5}}$ | $-\sqrt{\frac{n+1}{3}}$ | 0 | 0
2 | 2 | $\sqrt{\frac{3(n+5)}{7}}$ | $\sqrt{\frac{n+5}{7}}$ | $\sqrt{\frac{3(n+5)}{105}}$ | 0 | 0 | $-\sqrt{\frac{2n}{5}}$
2 | 1 | $\sqrt{\frac{2(n+5)}{7}}$ | $2\sqrt{\frac{6(n+5)}{105}}$ | $3\sqrt{\frac{n+5}{105}}$ | 0 | $\sqrt{\frac{n}{5}}$ | $-\sqrt{\frac{n}{5}}$
2 | 0 | $3\sqrt{\frac{2(n+5)}{105}}$ | $3\sqrt{\frac{3(n+5)}{105}}$ | $3\sqrt{\frac{2(n+5)}{105}}$ | $-\sqrt{\frac{n}{15}}$ | $2\sqrt{\frac{n}{15}}$ | $-\sqrt{\frac{n}{15}}$
2 | -1 | $3\sqrt{\frac{n+5}{105}}$ | $2\sqrt{\frac{6(n+5)}{105}}$ | $\sqrt{\frac{2(n+5)}{7}}$ | $-\sqrt{\frac{n}{5}}$ | $\sqrt{\frac{n}{5}}$ | 0
2 | -2 | $\sqrt{\frac{3(n+5)}{105}}$ | $\sqrt{\frac{n+5}{7}}$ | $\sqrt{\frac{3(n+5)}{7}}$ | $-\sqrt{\frac{2n}{5}}$ | 0 | 0
3 | 3 | $\frac{2}{3}\sqrt{n+6}$ | $\frac{1}{3}\sqrt{n+6}$ | $\frac{1}{3}\sqrt{\frac{n+6}{7}}$ | 0 | 0 | $-\sqrt{\frac{3(n-1)}{7}}$
3 | 2 | $\sqrt{\frac{n+6}{3}}$ | $2\sqrt{\frac{n+6}{21}}$ | $\sqrt{\frac{n+6}{21}}$ | 0 | $\sqrt{\frac{n-1}{7}}$ | $-\sqrt{\frac{2(n-1)}{7}}$
3 | 1 | $\sqrt{\frac{5(n+6)}{21}}$ | $\sqrt{\frac{5(n+6)}{21}}$ | $\sqrt{\frac{2(n+6)}{21}}$ | $-\sqrt{\frac{n-1}{35}}$ | $2\sqrt{\frac{2(n-1)}{35}}$ | $-\sqrt{\frac{6(n-1)}{35}}$
3 | 0 | $\sqrt{\frac{10(n+6)}{63}}$ | $\sqrt{\frac{16(n+6)}{63}}$ | $\sqrt{\frac{10(n+6)}{63}}$ | $-\sqrt{\frac{3(n-1)}{35}}$ | $3\sqrt{\frac{n-1}{35}}$ | $-\sqrt{\frac{3(n-1)}{35}}$
3 | -1 | $\sqrt{\frac{2(n+6)}{21}}$ | $\sqrt{\frac{5(n+6)}{21}}$ | $\sqrt{\frac{5(n+6)}{21}}$ | $-\sqrt{\frac{6(n-1)}{35}}$ | $2\sqrt{\frac{2(n-1)}{35}}$ | $-\sqrt{\frac{n-1}{35}}$
3 | -2 | $\sqrt{\frac{n+6}{21}}$ | $2\sqrt{\frac{n+6}{21}}$ | $\sqrt{\frac{n+6}{3}}$ | $-\sqrt{\frac{2(n-1)}{7}}$ | $\sqrt{\frac{n-1}{7}}$ | 0
3 | -3 | $\frac{1}{3}\sqrt{\frac{n+6}{7}}$ | $\frac{1}{3}\sqrt{n+6}$ | $\frac{2}{3}\sqrt{n+6}$ | $-\sqrt{\frac{3(n-1)}{7}}$ | 0 | 0

##  Appendix D CONSERVATION LAWS

In this Appendix, we derive the conservation laws associated with the ${\rm
U}(1)$ gauge and spin rotational symmetries. We also derive the expressions
for mass current (116) and spin current (117) from the GL action (301).
We first study the conservation law associated with the ${\rm U}(1)$ gauge
symmetry. The following argument is based on Refs. [62, 38]. Under an
infinitesimal ${\rm U}(1)$ gauge transformation $\Psi_{\alpha}\rightarrow
e^{i\delta\theta(\bm{r})}\Psi_{\alpha}\sim\Psi_{\alpha}+i\delta\theta(\bm{r})\Psi_{\alpha}$,
the variation of the GL action (301) is given as

$\displaystyle\delta S^{\mathrm{GL}}$ $\displaystyle=$ $\displaystyle
i\delta\theta\int^{\infty}_{0}{\rm d}\tau\int{\rm
d}^{d}x\left\\{K\partial_{\tau}({\bf\Psi}^{\dagger}{\bf\Psi})\right.$ (304)
$\displaystyle+J\partial_{\tau}[(\partial_{\tau}{\bf\Psi}^{\dagger}){\bf\Psi}-{\bf\Psi}^{\dagger}\partial_{\tau}{\bf\Psi}]$
$\displaystyle+\left.\frac{1}{2m^{*}}\nabla\cdot[(\nabla{\bf\Psi}^{\dagger}){\bf\Psi}-{\bf\Psi}^{\dagger}\nabla{\bf\Psi}]\right\\}.$

From $\delta S^{\mathrm{GL}}=0$ for arbitrary $\delta\theta$, we obtain the
continuity equation

$\displaystyle\partial_{t}(K\rho+Jq)+\nabla\cdot\bm{j}_{\rm m}=0,$ (305)

where $\rho$ and $\bm{j}_{\rm m}$ denote the superfluid density and superfluid
mass current, respectively, defined as

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle{\bf\Psi}^{\dagger}{\bf\Psi},$ (306) $\displaystyle q$
$\displaystyle=$
$\displaystyle\frac{1}{i}\left[(\partial_{t}{\bf\Psi}^{\dagger}){\bf\Psi}-{\bf\Psi}^{\dagger}(\partial_{t}{\bf\Psi})\right],$
(307) $\displaystyle\bm{j}_{\rm m}$ $\displaystyle=$
$\displaystyle\frac{1}{2im^{*}}\left[{\bf\Psi}^{\dagger}\nabla{\bf\Psi}-(\nabla{\bf\Psi}^{\dagger}){\bf\Psi}\right].$
(308)

We clarify the physical meaning of $q$ below. If $K\neq 0$, neglecting $q$,
Eq. (305) reduces to

$\displaystyle\partial_{t}\rho+\frac{1}{K}\nabla\cdot\bm{j}_{\rm m}=0.$ (309)

Equation (309) represents the conservation of superfluid density.
If $K=0$, keeping $q$, Eq. (309) reduces to

$\displaystyle\partial_{t}q+\frac{1}{J}\nabla\cdot\bm{j}_{\rm m}=0.$ (310)

Before discussing the interpretation of $q$, we show that $K=0$ holds for a
commensurate filling in the SF phase. We first note that Eq. (284) is
invariant under the following local gauge transformation

$\displaystyle\begin{split}b_{i\alpha}&\rightarrow
e^{i\eta(\tau)}b_{i\alpha},\\\ \phi_{i\alpha}&\rightarrow
e^{i\eta(\tau)}\phi_{i\alpha},\\\
\mu&\rightarrow\mu+i\partial_{\tau}\eta(\tau).\end{split}$ (311)

The GL action (301) should be invariant under the corresponding transformation
$\Psi_{\alpha}\rightarrow e^{i\eta(\tau)}\Psi_{\alpha}$, which requires the
condition

$\displaystyle K=-\frac{\partial C_{2}}{\partial\mu}.$ (312)

We can directly confirm Eq. (312) from Eqs. (274) and (296). On the other
hand, in the vicinity of the phase boundary, the ground-state energy can be
expanded by the superfluid density $n_{\rm c}$ as

$\displaystyle E=E^{0}(0,n)+C_{2}n_{\rm c}+O(n_{\rm c}^{2}).$ (313)

The expectation value of the filling number is given by

$\displaystyle\langle\hat{n}\rangle=-\frac{\partial E}{\partial\mu}\sim
n-\frac{\partial C_{2}}{\partial\mu}n_{\rm c}.$ (314)

In the MI phase, $\langle\hat{n}\rangle$ is an integer because $n_{\rm c}=0$
as expected. Meanwhile, since $n_{\rm c}\neq 0$, $\langle\hat{n}\rangle$ takes
an integer value only if $-\partial C_{2}/\partial\mu=0$ in the SF phase.
Thus, $K=0$ holds for a commensurate filling in the SF phase. Furthermore,
when $K=0$, the TDGL equation (303) is invariant under the transformation
$\Psi_{\alpha}\leftrightarrow\Psi_{\alpha}^{*}$, i.e., particle-hole symmetric
[38, 59, 63].
Equation (311) implies that the transformation $b_{i\alpha}\rightarrow
e^{i\eta}b_{i\alpha}$ and $\phi_{i\alpha}\rightarrow e^{i\eta}\phi_{i\alpha}$
is equivalent to the shift of $\mu$: $\mu\rightarrow\mu-i\partial_{\tau}\eta$
[59]. The infinitesimal shift of the chemical potential
$\mu\rightarrow\mu+i\partial\delta\eta$ is, therefore, equivalent to the
transformation of $\Psi_{\alpha}$:

$\displaystyle\Psi_{\alpha}\rightarrow
e^{-i\delta\eta(\tau)}\Psi_{\alpha}\sim\Psi_{\alpha}-i\delta\eta(\tau)\Psi_{\alpha}.$
(315)

The variation of $S^{\rm GL}$ under Eq. (315) is given by

$\displaystyle\delta S^{\mathrm{GL}}=-\delta\mu\int{\rm d}^{d}x(K\rho+Jq),$
(316)

where $\delta\mu=i\partial_{\tau}\eta(\tau)$. Comparing Eq. (316) with $\delta
S^{\mathrm{GL}}=(\delta S^{\mathrm{GL}}/\delta\mu(\tau))\delta\mu(\tau)$, we
obtain

$\displaystyle-\frac{\delta S^{\mathrm{GL}}}{\delta\mu(\tau)}=\int{\rm
d}^{d}x(K\rho+Jq).$ (317)

The left-hand side of Eq. (317) represents the deviation of particle number
from the static value. For a commensurate filling, setting $K=0$, Eq. (317)
reduces to

$\displaystyle-\frac{\delta S^{\mathrm{GL}}}{\delta\mu(\tau)}=\int{\rm
d}^{d}xJq.$ (318)

We thus find that $q$ means the deviation of particle density from a
commensurate filling.
We next derive the conservation law associated with the spin rotational
symmetry. We introduce an infinitesimal spin rotation about a unit vector
$\bm{n}$:

$\displaystyle{\bf\Psi}\rightarrow
e^{-i(\bm{n}\cdot\bm{F})\delta\lambda}{\bf\Psi}\sim{\bf\Psi}-i(\bm{n}\cdot\bm{F})\delta\lambda{\bf\Psi}.$
(319)

The variation of the GL action (301) under Eq. (319) is given by

$\displaystyle\delta S^{\mathrm{GL}}$ $\displaystyle=$
$\displaystyle-i\delta\lambda\int^{\infty}_{0}{\rm d}\tau\int{\rm
d}^{d}x\left\\{K\partial_{\tau}\left[{\bf\Psi}^{\dagger}(\bm{n}\cdot\bm{F}){\bf\Psi}\right]\right.$
$\displaystyle+$ $\displaystyle
J\partial_{\tau}\left[(\partial_{\tau}{\bf\Psi}^{\dagger})(\bm{n}\cdot\bm{F}){\bf\Psi}-{\bf\Psi}^{\dagger}(\bm{n}\cdot\bm{F})\partial_{\tau}{\bf\Psi}\right]$
$\displaystyle+$
$\displaystyle\left.\frac{1}{2m^{*}}\nabla\cdot\left[(\nabla{\bf\Psi}^{\dagger})(\bm{n}\cdot\bm{F}){\bf\Psi}-{\bf\Psi}^{\dagger}(\bm{n}\cdot\bm{F})\nabla{\bf\Psi}\right]\right\\}.$

$\delta S^{\mathrm{GL}}=0$ for arbitrary $\delta\lambda$ yields the continuity
equation

$\displaystyle\partial_{t}(K\rho_{\rm s}^{\bm{n}}+Jq_{\rm
s}^{\bm{n}})+\nabla\cdot\bm{j}_{\rm s}^{\bm{n}}=0,$ (321)

where $\rho_{\rm s}^{\bm{n}}$ denotes the magnetization density for the
component along $\bm{n}$ and $\bm{j}_{\rm s}^{\bm{n}}$ the spin current for
$\rho_{\rm s}^{\bm{n}}$. They are defined as

$\displaystyle\rho_{\rm s}^{\bm{n}}$ $\displaystyle=$
$\displaystyle{\bf\Psi}^{\dagger}(\bm{n}\cdot\bm{F}){\bf\Psi},$ (322)
$\displaystyle q_{\rm s}^{\bm{n}}$ $\displaystyle=$
$\displaystyle\frac{1}{i}\left[(\partial_{t}{\bf\Psi}^{\dagger})(\bm{n}\cdot\bm{F}){\bf\Psi}-{\bf\Psi}^{\dagger}(\bm{n}\cdot\bm{F})\partial_{t}{\bf\Psi})\right],$
(323) $\displaystyle\bm{j}_{\rm s}^{\bm{n}}$ $\displaystyle=$
$\displaystyle\frac{1}{2im^{*}}\left[{\bf\Psi}^{\dagger}(\bm{n}\cdot\bm{F})\nabla{\bf\Psi}-(\nabla{\bf\Psi}^{\dagger})(\bm{n}\cdot\bm{F}){\bf\Psi}\right].$

We clarify the physical meaning of $q^{\bm{n}}_{\rm s}$ below. If $K\neq 0$,
neglecting $q^{\bm{n}}_{\rm s}$, Eq. (321) represents the conservation of
magnetization for the component along $\bm{n}$.
Analogous to Eq. (311), we consider the following trans formation under a spin
rotation

$\displaystyle\begin{split}\vec{b}_{i}&\rightarrow
e^{-i(\bm{n}\cdot\bm{F})\nu(\tau)}\vec{b}_{i},\\\ \vec{\phi}_{i}&\rightarrow
e^{-i(\bm{n}\cdot\bm{F})\nu(\tau)}\vec{\phi}_{i},\end{split}$ (325)

where $\vec{b}_{i}=(b_{i1},b_{i0},b_{i-1})^{T}$ and
$\vec{\phi}_{i}=(\phi_{i1},\phi_{i0},\phi_{i-1})^{T}$. Under the
transformation (325), Eq. (284) is transformed as

$\displaystyle S\rightarrow S-\int^{\beta}_{0}{\rm
d}\tau\bm{H}(\tau)\cdot\int{\rm d}^{d}x\langle\bm{F}\rangle_{b}.$ (326)

Here, $\bm{H}=i\partial_{\tau}\nu\bm{n}$ and
$\langle\bm{F}\rangle_{b}=\sum_{\alpha,\beta}b_{i\alpha}^{*}(e^{i\bm{F}\cdot\bm{n}\nu}\bm{F}e^{-i\bm{F}\cdot\bm{n}\nu})_{\alpha\beta}b_{i\beta}$.
Equation (326) shows that the transformation (325) is equivalent to
application of magnetic field $\bm{H}$ in the direction of $\bm{n}$.

The transformation of ${\bf\Psi}$ that corresponds to Eq. (325) is
${\bf\Psi}\rightarrow e^{-i(\bm{n}\cdot\bm{F})\nu(\tau)}{\bf\Psi}$. Under an
infinitesimal spin rotation ${\bf\Psi}\rightarrow
e^{-i(\bm{n}\cdot\bm{F})\delta\nu}{\bf\Psi}\sim{\bf\Psi}-i(\bm{n}\cdot\bm{F})\delta\nu{\bf\Psi}$,
the variation of Eq. (301) is given by

$\displaystyle\delta S^{\mathrm{GL}}=-\delta\bm{H}(\tau)\cdot\int{\rm
d}^{d}x(K\langle\\!\langle\bm{F}\rangle\\!\rangle+J\bm{q}_{\rm s}),$ (327)

where $\delta\bm{H}(\tau)=i\delta\partial_{\tau}\nu(\tau)\bm{n}$. Comparing
Eq. (327) and

$\displaystyle\delta S^{\mathrm{GL}}=\frac{\delta
S^{\mathrm{GL}}}{\delta\bm{H}}\cdot\delta\bm{H},$ (328)

we obtain

$\displaystyle-\frac{\delta S^{\mathrm{GL}}}{\delta\bm{H}}=\int{\rm
d}^{d}x(K\langle\\!\langle\bm{F}\rangle\\!\rangle+J\bm{q}_{\rm s}),$ (329)

where $\bm{n}\cdot\bm{q}_{\rm s}=q_{\rm s}^{\bm{n}}$. The left-hand side of
Eq. (329) represents the magnetization along $\bm{H}$. In the case of a
commensurate filling, i.e., $K=0$, we obtain

$\displaystyle-\frac{\delta S^{\mathrm{GL}}}{\delta\bm{H}}=\int{\rm
d}^{d}xJ\bm{q}_{\rm s}.$ (330)

We thus find that the $q_{\rm s}^{\bm{n}}$ represents the deviation of
magnetization density from the value for a commensurate filling.

## References

  * Lewenstein  _et al._ [2007] M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, A. Sen(De), and U. Sen, Advances in Physics 56, 243 (2007).
  * Lewenstein  _et al._ [2012] M. Lewenstein, A. Sanpera, and V. Ahufinger,  _Ultracold Atoms in Optical Lattices: Simulating quantum many-body systems_ (OUP Oxford, 2012).
  * García-Ripoll  _et al._ [2004] J. J. García-Ripoll, M. A. Martin-Delgado, and J. I. Cirac, Phys. Rev. Lett. 93, 250405 (2004).
  * Imambekov  _et al._ [2004] A. Imambekov, M. Lukin, and E. Demler, Phys. Rev. Lett. 93, 120405 (2004).
  * Widera  _et al._ [2005] A. Widera, F. Gerbier, S. Fölling, T. Gericke, O. Mandel, and I. Bloch, Phys. Rev. Lett. 95, 190405 (2005).
  * Yamashita and Jack [2007] M. Yamashita and M. W. Jack, Phys. Rev. A 76, 023606 (2007).
  * Shinozaki  _et al._ [2013] M. Shinozaki, S. Tsuchiya, S. Abe, T. Ozaki, and T. Nikuni, Journal of Low Temperature Physics 175, 236 (2013).
  * Stamper-Kurn and Ueda [2013] D. M. Stamper-Kurn and M. Ueda, Rev. Mod. Phys. 85, 1191 (2013).
  * Fujimoto  _et al._ [2019] K. Fujimoto, R. Hamazaki, and M. Ueda, Phys. Rev. Lett. 122, 173001 (2019).
  * Demler and Zhou [2002] E. Demler and F. Zhou, Phys. Rev. Lett. 88, 163001 (2002).
  * Tsuchiya  _et al._ [2004] S. Tsuchiya, S. Kurihara, and T. Kimura, Phys. Rev. A 70, 043628 (2004).
  * Rizzi  _et al._ [2005] M. Rizzi, D. Rossini, G. De Chiara, S. Montangero, and R. Fazio, Phys. Rev. Lett. 95, 240404 (2005).
  * Łącki  _et al._ [2011] M. Łącki, S. Paganelli, V. Ahufinger, A. Sanpera, and J. Zakrzewski, Phys. Rev. A 83, 013605 (2011).
  * Kimura  _et al._ [2005] T. Kimura, S. Tsuchiya, and S. Kurihara, Phys. Rev. Lett. 94, 110403 (2005).
  * Krutitsky  _et al._ [2005] K. V. Krutitsky, M. Timmer, and R. Graham, Phys. Rev. A 71, 033623 (2005).
  * Batrouni  _et al._ [2005] G. G. Batrouni, F. F. Assaad, R. T. Scalettar, and P. J. H. Denteneer, Phys. Rev. A 72, 031601 (2005).
  * Yamamoto  _et al._ [2013] D. Yamamoto, T. Ozaki, C. A. R. Sá de Melo, and I. Danshita, Phys. Rev. A 88, 033624 (2013).
  * Yip [2003] S. K. Yip, Phys. Rev. Lett. 90, 250402 (2003).
  * Imambekov  _et al._ [2003] A. Imambekov, M. Lukin, and E. Demler, Phys. Rev. A 68, 063602 (2003).
  * Zhou and Snoek [2003] F. Zhou and M. Snoek, Annals of Physics 308, 692 (2003).
  * Snoek and Zhou [2004] M. Snoek and F. Zhou, Phys. Rev. B 69, 094410 (2004).
  * Apaja and Syljuåsen [2006] V. Apaja and O. F. Syljuåsen, Phys. Rev. A 74, 035601 (2006).
  * Bernier  _et al._ [2006] J.-S. Bernier, K. Sengupta, and Y. B. Kim, Phys. Rev. B 74, 155124 (2006).
  * de Forges de Parny  _et al._ [2013] L. de Forges de Parny, F. Hébert, V. G. Rousseau, and G. G. Batrouni, Phys. Rev. B 88, 104509 (2013).
  * Alavani  _et al._ [2018] B. K. Alavani, A. Das, and R. V. Pai, Journal of Physics B: Atomic, Molecular and Optical Physics 51, 145302 (2018).
  * de Forges de Parny and Rousseau [2018] L. de Forges de Parny and V. G. Rousseau, Phys. Rev. A 97, 023628 (2018).
  * Jiang  _et al._ [2016] J. Jiang, L. Zhao, S.-T. Wang, Z. Chen, T. Tang, L.-M. Duan, and Y. Liu, Phys. Rev. A 93, 063607 (2016).
  * Austin  _et al._ [2021] J. O. Austin, Z. Chen, Z. N. Shaw, K. W. Mahmud, and Y. Liu, Communications Physics 4 (2021).
  * Burger  _et al._ [2001] S. Burger, F. S. Cataliotti, C. Fort, F. Minardi, M. Inguscio, M. L. Chiofalo, and M. P. Tosi, Phys. Rev. Lett. 86, 4447 (2001).
  * Cataliotti  _et al._ [2003] F. S. Cataliotti, L. Fallani, F. Ferlaino, C. Fort, P. Maddaloni, and M. Inguscio, New Journal of Physics 5, 71 (2003).
  * Fallani  _et al._ [2004] L. Fallani, L. De Sarlo, J. E. Lye, M. Modugno, R. Saers, C. Fort, and M. Inguscio, Phys. Rev. Lett. 93, 140406 (2004).
  * De Sarlo  _et al._ [2005] L. De Sarlo, L. Fallani, J. E. Lye, M. Modugno, R. Saers, C. Fort, and M. Inguscio, Phys. Rev. A 72, 013603 (2005).
  * Mun  _et al._ [2007] J. Mun, P. Medley, G. K. Campbell, L. G. Marcassa, D. E. Pritchard, and W. Ketterle, Phys. Rev. Lett. 99, 150604 (2007).
  * Ferris  _et al._ [2008] A. J. Ferris, M. J. Davis, R. W. Geursen, P. B. Blakie, and A. C. Wilson, Phys. Rev. A 77, 012712 (2008).
  * Wu and Niu [2001] B. Wu and Q. Niu, Phys. Rev. A 64, 061603 (2001).
  * Altman and Auerbach [2002] E. Altman and A. Auerbach, Phys. Rev. Lett. 89, 250404 (2002).
  * Biao and Qian [2003] W. Biao and N. Qian, New Journal of Physics 5, 104 (2003).
  * Polkovnikov  _et al._ [2005] A. Polkovnikov, E. Altman, E. Demler, B. Halperin, and M. D. Lukin, Phys. Rev. A 71, 063613 (2005).
  * Altman  _et al._ [2005] E. Altman, A. Polkovnikov, E. Demler, B. I. Halperin, and M. D. Lukin, Phys. Rev. Lett. 95, 020402 (2005).
  * Asaoka  _et al._ [2016] R. Asaoka, H. Tsuchiura, M. Yamashita, and Y. Toga, Phys. Rev. A 93, 013628 (2016).
  * Kim  _et al._ [2017] J. H. Kim, S. W. Seo, and Y. Shin, Phys. Rev. Lett. 119, 185302 (2017).
  * Wu [1996] Y. Wu, Phys. Rev. A 54, 4534 (1996).
  * Ho and Yip [2000] T.-L. Ho and S. K. Yip, Phys. Rev. Lett. 84, 4031 (2000).
  * Ho [1998] T.-L. Ho, Phys. Rev. Lett. 81, 742 (1998).
  * Zhou [2001] F. Zhou, Phys. Rev. Lett. 87, 080401 (2001).
  * Ohmi and Machida [1998] T. Ohmi and K. Machida, Journal of the Physical Society of Japan 67, 1822 (1998).
  * Landau and Lifshitz [2013] L. Landau and E. Lifshitz,  _Statistical Physics: Volume 5_ , 5 (Elsevier Science, 2013).
  * Goldstone  _et al._ [1962] J. Goldstone, A. Salam, and S. Weinberg, Phys. Rev. 127, 965 (1962).
  * de Gennes and Prost [1993] P. de Gennes and J. Prost,  _The Physics of Liquid Crystals_ , International Series of Monogr (Clarendon Press, 1993).
  * Barnett  _et al._ [2006] R. Barnett, A. Turner, and E. Demler, Phys. Rev. Lett. 97, 180412 (2006).
  * Yukawa and Ueda [2012] E. Yukawa and M. Ueda, Phys. Rev. A 86, 063614 (2012).
  * Patrick  _et al._ [2014] N. Patrick, Q. Friedemann, and S. Ralf, Journal of Physics A: Mathematical and Theoretical 47, 225004 (2014).
  * [53] The healing length $\xi$ can be derived as follows: substituting ${\bf\Psi}=\sqrt{\rho^{\rm 0}+\delta\rho(x)}(0,1,0)^{T}$ into the TDGL equation (51) and linearizing by $\delta\rho(x)$, one obtains $\delta\rho(x)\propto e^{-x/\xi}$, where $\rho^{\rm 0}$ is the equilibrium SF density and $\delta\rho(\infty)=0$ is assumed.
  * Fujimoto and Tsubota [2012] K. Fujimoto and M. Tsubota, Phys. Rev. A 85, 033642 (2012).
  * Zhu  _et al._ [2015] Q. Zhu, Q.-f. Sun, and B. Wu, Phys. Rev. A 91, 023633 (2015).
  * Mohamed and Axel [2013a] M. Mohamed and P. Axel, Tuning superfluid phases of spin-1 bosons in cubic optical lattice with linear zeeman effect (2013a), arXiv:1310.0600 [cond-mat.quant-gas] .
  * Mohamed and Axel [2013b] M. Mohamed and P. Axel, Laser Physics Letters 10, 115501 (2013b).
  * Messiah [1999] A. Messiah,  _Quantum Mechanics_ , Dover books on physics No. 2 (Dover Publications, 1999).
  * Fisher  _et al._ [1989] M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, Phys. Rev. B 40, 546 (1989).
  * Negele and Orland [1998] J. W. Negele and H. Orland,  _Quantum Many-particle Systems_ (Westview Press, 1998).
  * Kubo [1962] R. Kubo, Journal of the Physical Society of Japan 17, 1100 (1962).
  * Sachdev [2011] S. Sachdev,  _Quantum Phase Transitions_, 2nd ed. (Cambridge University Press, 2011).
  * Nakayama  _et al._ [2015] T. Nakayama, I. Danshita, T. Nikuni, and S. Tsuchiya, Phys. Rev. A 92, 043610 (2015).