#  Tachyonic and Plasma Instabilities of $\eta$-Pairing States Coupled to
Electromagnetic Fields

Naoto Tsuji  Department of Physics, University of Tokyo, Hongo, Tokyo
113-0033, Japan  RIKEN Center for Emergent Matter Science (CEMS), Wako
351-0198, Japan  Masaya Nakagawa  Department of Physics, University of Tokyo,
Hongo, Tokyo 113-0033, Japan  Masahito Ueda  Department of Physics, University
of Tokyo, Hongo, Tokyo 113-0033, Japan  RIKEN Center for Emergent Matter
Science (CEMS), Wako 351-0198, Japan  Institute for Physics of Intelligence,
University of Tokyo, Hongo, Tokyo 113-0033, Japan

###### Abstract

Cooper pairs featuring a nonzero center-of-mass crystal momentum
$\boldsymbol{Q}=(\pi,\pi,\dots)$ and an off-diagonal long-range order
($\eta$-pairing states) constitute exact eigenstates of a Hubbard model [C. N.
Yang, Phys. Rev. Lett. 63, 2144 (1989)]. Here we show that the $\eta$-pairing
states are rendered dynamically unstable via coupling to dynamical
electromagnetic fields. The instability is caused by “tachyonic”
electromagnetic fields and unstable plasma modes for attractive and repulsive
interactions, respectively. The typical time scale of the growth of the
instability is of the order of femtoseconds for electron systems in solids,
which places a strict bound on the lifetime of the $\eta$-pairing states. The
decay of the $\eta$-pairing states leads to enhanced light emission with
frequencies shifted from the Hubbard interaction strength for the repulsive
case and unattenuated electromagnetic penetration for the attractive case.

Introduction.— A search for long-lived non-thermal excited states that support
macroscopic long-range order is a modern challenge in nonequilibrium condensed
matter physics. If such a state exists, a rich variety of possibilities arises
that could extend the landscape of long-range order in solid-state materials.
In fact, there have been a number of experimental reports on the possible
existence of nonequilibrium long-range order, including light-induced Fausti
_et al._ (2011); Kaiser  _et al._ (2014); Hu  _et al._ (2014); Mitrano  _et
al._ (2016); Cantaluppi  _et al._ (2018); Buzzi  _et al._ (2020); Bud  or
quench-induced Oike  _et al._ (2018) superconductivity (see also Niwa  _et
al._ (2019); Zhang  _et al._ (2020)), photoinduced ferromagnetism Matsubara
_et al._ (2007); McLeod  _et al._ (2020) and charge density waves Stojchevska
_et al._ (2014); Vaskivskyi  _et al._ (2015); Sun  _et al._ (2018). While the
experimental progress is underway, the theoretical understanding is yet to be
made. A major challenge is that an analysis of excited states in quantum many-
body systems often requires approximations that render the conclusion on the
existence of such a state highly nontrivial.

A very exception to this situation is the $\eta$-pairing states, which are
known to be exact eigenstates of a Hubbard model as revealed by C. N. Yang
Yang (1989). The $\eta$-pairing states exhibit a number of remarkable
features. In particular, they have an off-diagonal long-range order (ODLRO) in
arbitrary dimensions even though their eigenenergies lie much higher than the
ground state. This is to be contrasted with finite-temperature thermal states,
which cannot show ODLROs in one and two dimensions due to the Mermin-Wagner
theorem. The non-thermal nature of the $\eta$-pairing states has also been
discussed recently in the context of quantum many-body scars Vafek  _et al._
(2017); Mark and Motrunich (2020); Moudgalya  _et al._ (2020).

The presence of such a non-thermal state suggests that the $\eta$-pairing
states with an ODLRO (and hence superconductivity) might be realized in
nonequilibrium situations. Recent theoretical studies have demonstrated that
this is indeed possible in several different setups, including periodic
Kitamura and Aoki (2016); Peronaci  _et al._ (2020); Cook and Clark (2020);
Tindall  _et al._ (2021) and pulsed Kaneko  _et al._ (2019, 2020); Werner  _et
al._ (2019); Li  _et al._ (2020) electric-field drives, dissipation
engineering Diehl  _et al._ (2008); Kraus  _et al._ (2008), spin-dependent
dephasing Bernier  _et al._ (2013); Tindall  _et al._ (2019), and spontaneous
light emission Nak . These mechanisms will work for the Hubbard model with or
without coupling to an external bath, which may be realized in electrically
neutral ultracold atoms trapped in an optical lattice.

In view of applications to real materials, one cannot ignore the coupling of
electrons to dynamical electromagnetic fields, since electrons have electric
charges. This point is crucial for the stability of the $\eta$-pairing states
supported by the long lifetime of doublons. If doublons decay into single
particles or lose their momenta, they induce local electric currents due to
charge transfer, which then generate dynamical electromagnetic fields. The
effect of the latter feedbacks to electrons, and causes collective modes of
electromagnetic fields, which accelerate the relaxation of doublons. Such a
dynamical instability deserves careful scrutiny in view of growing attention
in nonequilibrium superconductivity.

In this Letter, we study the dynamics of the $\eta$-pairing states in the
Hubbard model coupled to dynamical electromagnetic fields. Our approach is
based on the exact solution of the electromagnetic response function (or the
Meissner kernel) $K^{\mu\nu}(\boldsymbol{q},\omega)$ with full momentum
($\boldsymbol{q}$) and frequency ($\omega$) dependences. In contrast, previous
studies have focused on the static and uniform limit (i.e.,
$\boldsymbol{q}=\omega=0$) Su  _et al._ (1991, 1992); Kaneko  _et al._ (2020).
As we will see, the momentum and frequency dependences play a pivotal role in
dynamical instabilities of $\eta$-pairing states. Combining the obtained
results with the Maxwell equations, we rigorously prove the existence of the
“tachyonic” and plasma instabilities for attractively and repulsively
interacting systems, respectively. The time scale of the growth of the
instability is surprisingly short, being of the order of femtoseconds or even
shorter than that for ordinary materials. This puts a severe constraint on the
lifetime of the $\eta$-pairing states in electron systems. Finally, we discuss
that the decay of the $\eta$-pairing states leads to intense light emission
with frequencies shifted from the interaction strength in the repulsive case,
and unattenuated penetration of electromagnetic fields in the attractive case.

$\eta$ pairing in the Hubbard model.— We consider the Hubbard model on a
$d$-dimensional cubic lattice subject to the periodic boundary condition with
the Hamiltonian,

$\displaystyle H$ $\displaystyle=-t_{h}\sum_{\langle
ij\rangle,\sigma}(c_{i\sigma}^{\dagger}c_{j\sigma}+\mbox{H.c.})+U\sum_{i}n_{i\uparrow}n_{i\downarrow}-\frac{U}{2}\sum_{i\sigma}n_{i\sigma},$
(1)

where $t_{h}$ ($>0$) is the hopping amplitude, $c_{i\sigma}^{\dagger}$ is a
creation operator of an electron at site $i$ with spin
$\sigma=\uparrow,\downarrow$, $\langle ij\rangle$ represents a pair of
nearest-neighbor lattice sites, $U$ is the on-site interaction strength, and
$n_{i\sigma}=c_{i\sigma}^{\dagger}c_{i\sigma}$ is the particle-number
operator. Since we fix the total number of electrons throughout this Letter,
the last term in Eq. (1) is a constant. We set the lattice constant $a=1$ and
the Planck constant $\hbar=1$ unless otherwise noted.

The Hubbard model (1) has the spin SU(2) symmetry together with the “hidden”
$\eta$ SU(2) symmetry, which altogether form the symmetry of ${\rm
SU(2)}\times{\rm SU(2)}/\mathbb{Z}_{2}\simeq{\rm SO(4)}$ Yang and Zhang
(1990). The existence of $\eta$-pairing states as the exact eigenstates of the
Hubbard model essentially relies on this fact. To see the $\eta$ symmetry, we
define the $\eta$ operators,
$\eta^{+}:=\sum_{j}e^{i\boldsymbol{Q}\cdot\boldsymbol{R}_{j}}c_{j\uparrow}^{\dagger}c_{j\downarrow}^{\dagger}$,
$\eta^{-}:=(\eta^{+})^{\dagger}$, and
$\eta^{z}:=\frac{1}{2}\sum_{j}(n_{j\uparrow}+n_{j\downarrow}-1)$, where
$\boldsymbol{Q}=(\pi,\pi,\dots)$ is the momentum at the Brillouin-zone corner,
and $\boldsymbol{R}_{j}$ is the position vector of lattice site $j$. The
$\eta$ operators satisfy the ordinary su(2) algebra, i.e.,
$[\eta^{+},\eta^{-}]=2\eta^{z}$ and $[\eta^{z},\eta^{\pm}]=\pm\eta^{\pm}$.
From direct calculations, one can confirm that they all commute with the
Hamiltonian (1): $[H,\eta^{\alpha}]=0$ ($\alpha=\pm,z$).

Using the $\eta$ operators, one can construct Yang’s $\eta$-pairing states.
The simplest one is

$\displaystyle|\psi_{N}\rangle$
$\displaystyle=\frac{1}{\sqrt{\mathcal{N}_{N}}}(\eta^{+})^{\frac{N}{2}}|0\rangle,$
(2)

where $N$ is the number of electrons which is assumed to be an even integer,
$\mathcal{N}_{N}$ is the normalization constant (such that
$\langle\psi_{N}|\psi_{N}\rangle=1$), and $|0\rangle$ is the vacuum state. The
$\eta$-pairing state $|\psi_{N}\rangle$ consists of $\frac{N}{2}$ doublons
having momentum $\boldsymbol{Q}$. Since $\eta^{+}$ commutes with $H$ (1),
$|\psi_{N}\rangle$ (2) is indeed the exact eigenstate of $H$ with the
eigenenergy $E_{N}=0$ in arbitrary dimensions. The state $|\psi_{N}\rangle$
(2) has the ODLRO
$\frac{1}{2}\langle\psi_{N}|(c_{i\uparrow}^{\dagger}c_{i\downarrow}^{\dagger}c_{j\downarrow}c_{j\uparrow}+\mbox{H.c.})|\psi_{N}\rangle=e^{i\boldsymbol{Q}\cdot(\boldsymbol{R}_{i}-\boldsymbol{R}_{j})}C_{M,N}$
Yang (1989) ($C_{M,N}:=\frac{\frac{N}{2}(M-\frac{N}{2})}{M(M-1)}$ and $M$ is
the number of lattice sites), which saturates the upper bound of ODLRO Yang
(1962); Nak . Physically, $|\psi_{N}\rangle$ (2) corresponds to the condensate
of spin-singlet Cooper pairs with the center-of-mass momentum
$\boldsymbol{Q}$.

Electromagnetic response of $\eta$-pairing states.— We study the
electromagnetic response of the $\eta$-pairing state $|\psi_{N}\rangle$ (2)
within the linear-response regime. We focus on the three-dimensional case
($d=3$). However, most of the results in the present Letter can
straightforwardly be extended to other dimensions. The response of the current
against an external electromagnetic field with momentum $\boldsymbol{q}$ and
frequency $\omega$ is given by
$j^{\mu}(\boldsymbol{q},\omega)=-K^{\mu\nu}(\boldsymbol{q},\omega)A_{\nu}(\boldsymbol{q},\omega)$
($\mu,\nu=x,y,z$), where $K^{\mu\nu}(\boldsymbol{q},\omega)$ is the Meissner
kernel Schrieffer (1983) and $A_{\nu}(\boldsymbol{q},\omega)$ is the vector
potential.

In general, the kernel $K^{\mu\nu}(\boldsymbol{q},\omega)$ consists of the
paramagnetic and diamagnetic components Schrieffer (1983). In the case of
$\eta$-pairing states, the diamagnetic component vanishes exactly, since it is
proportional to the kinetic energy kin , which vanishes for the $\eta$-pairing
states. This is in stark contrast to ordinary superconductors, in which
perfect diamagnetism arises from the diamagnetic component of the Meissner
kernel. In the $\eta$-pairing states, the paramagnetic component takes over
the role of the diamagnetic one in ordinary superconductors.

The paramagnetic component is given by the Kubo formula,

$\displaystyle K_{\rm para}^{\mu\nu}(\boldsymbol{R}_{j},t)$
$\displaystyle=-i\theta(t)\langle\psi_{N}|[J^{\mu}(\boldsymbol{R}_{j},t),J^{\nu}(0,0)]|\psi_{N}\rangle,$
(3)

where $\theta(t)$ is the unit-step function ($\theta(t)=1$ for $t\geq 0$ and
$\theta(t)=0$ otherwise), and $J^{\mu}(\boldsymbol{R}_{j},t)$ is the local
current operator at site $j$ and time $t$ in the Heisenberg picture. The local
current $J^{\mu}(\boldsymbol{R}_{j},0)=J^{\mu}(\boldsymbol{R}_{j})$ is
expressed explicitly as
$J^{\mu}(\boldsymbol{R}_{j})=-iet_{h}\sum_{\sigma}(c_{j+\mu,\sigma}^{\dagger}c_{j\sigma}-c_{j\sigma}^{\dagger}c_{j+\mu,\sigma})$,
where $e$ is the electric charge, and $j+\mu$ represents the nearest-neighbor
site of $j$ in the $\mu$ direction.

We can evaluate Eq. (3) exactly for arbitrary $\boldsymbol{R}_{j}$ and $t$
using the following algebraic relations:
$[J^{\mu}(\boldsymbol{R}_{j}),\eta^{\pm}]=:\pm
2J_{\eta}^{\mu\pm}(\boldsymbol{R}_{j})$,
$[J_{\eta}^{\mu\pm}(\boldsymbol{R}_{j}),\eta^{\pm}]=0$, and
$[J_{\eta}^{\mu\pm}(\boldsymbol{R}_{j}),\eta^{\mp}]=\pm
J^{\mu}(\boldsymbol{R}_{j})$. They allow us to reduce the $N$-particle
correlation function (3) to that of the vacuum state sup ,

$\displaystyle K_{\rm para}^{\mu\nu}(\boldsymbol{R}_{j},t)$
$\displaystyle=-4i\theta(t)C_{M,N}[\langle
0|J_{\eta}^{\mu-}(\boldsymbol{R}_{j},t)J_{\eta}^{\nu+}(0,0)|0\rangle$
$\displaystyle\quad-\langle
0|J_{\eta}^{\nu-}(0,0)J_{\eta}^{\mu+}(\boldsymbol{R}_{j},t)|0\rangle].$ (4)

In this way, the $N$-particle problem reduces to the two-particle problem,
which is exactly solvable Essler  _et al._ (2005).

We further decompose the kernel into the transverse and longitudinal
components. Without loss of generality, we assume that the momentum
$\boldsymbol{q}$ of the vector potential points in the $z$ direction. The
transverse component is defined as
$K^{\perp}(\boldsymbol{q},\omega):=K^{\mu\mu}(\boldsymbol{q},\omega)$
($\mu=x,y$), while the longitudinal one is
$K^{\parallel}(\boldsymbol{q},\omega):=K^{zz}(\boldsymbol{q},\omega)$. The
kernel does not have the off-diagonal components
($K^{\mu\nu}(\boldsymbol{q},\omega)=0$ for $\mu\neq\nu$) sup . In the two-
particle dynamics involved in Eq. (4), the center-of-mass momentum and the
relative coordinates in the $x$ and $y$ directions of the two particles are
conserved. In addition, for the transverse components the two particles never
sit at the same site, making the dynamics effectively noninteracting. These
observations lead us to an analytical solution for the transverse component
sup ,

$\displaystyle K^{\perp}(\boldsymbol{q},\omega)$
$\displaystyle=-8ie^{2}t_{h}C_{M,N}\left[f_{q}\big{(}\tfrac{\omega+U}{t_{h}}\big{)}-f_{q}\big{(}\tfrac{\omega-U}{t_{h}}\big{)}\right],$
$\displaystyle f_{q}(x)$ $\displaystyle:=\begin{cases}\frac{i{\rm
sgn}(x)}{\sqrt{x^{2}-16\sin^{2}\frac{q}{2}}}&{\rm
for}\;|x|>4\sin\frac{q}{2};\\\
\frac{1}{\sqrt{16\sin^{2}\frac{q}{2}-x^{2}}}&{\rm
for}\;|x|<4\sin\frac{q}{2},\end{cases}$ (5)

where $q$ $(>0)$ is the $z$ component of $\boldsymbol{q}$. The longitudinal
component does not have such a compact expression, but can be evaluated
exactly in a similar manner sup .

The exact solution (5) for the electromagnetic response function reveals a
number of important properties of the $\eta$-pairing states. By taking the
limit $\lim_{\boldsymbol{q}\to 0}\lim_{\omega\to
0}K^{\perp}(\boldsymbol{q},\omega)=\frac{16e^{2}t_{h}^{2}}{U}C_{M,N}=:\frac{1}{\pi}D_{s}$,
one can recover the Meissner weight or superfluid stiffness $D_{s}$ in Refs.
Su  _et al._ (1991, 1992). By taking another limit $\lim_{\boldsymbol{q}\to
0}K^{\perp}(\boldsymbol{q},\omega)=:-i\omega\sigma(\omega)$, one obtains the
optical conductivity,
$\sigma(\omega)=\frac{8ie^{2}t_{h}^{2}C_{M,N}}{\omega+i\delta}(\frac{1}{\omega+U+i\delta}-\frac{1}{\omega-U+i\delta})$
($\delta$ is a positive infinitesimal constant), the real part of which shows
delta-function-like peaks at $\omega=0$ and $\pm U$. If one split the optical
conductivity into the singular part at $\omega=0$ and the regular part as
${\rm Re}\,\sigma(\omega)=D\delta(\omega)+\sigma_{\rm reg}(\omega)$, one
obtains the Drude weight or charge stiffness $D=\frac{16\pi
e^{2}t_{h}^{2}}{U}C_{M,N}$ Kaneko  _et al._ (2020). For the distinction
between $D$ and $D_{s}$, we refer to Ref. Scalapino  _et al._ (1993).

Figure 1: (a) Schematic dispersion of electromagnetic fields and their
effective potentials (3D images) for ordinary superconductors and tachyonic
superconductors. The dashed lines show the dispersion relation of
electromagnetic fields in the vacuum. (b) Phase diagram of the $\eta$-pairing
states in the Hubbard model.

For $U>0$, $D_{s}$ becomes positive, where the electromagnetic field acquires
a mass due to the Anderson-Higgs mechanism Anderson (1963); Shimano and Tsuji
(2020) [$\omega^{2}\simeq c^{2}q^{2}+m^{2}c^{4}$ with $c$ the speed of light
and $m^{2}\propto D_{s}>0$, see Fig. 1(a)]. On the other hand, for $U<0$,
$D_{s}$ takes a negative value, implying that the electromagnetic field has a
negative squared mass [$m^{2}\propto D_{s}<0$, Fig. 1(a)]. Thus, the system is
a “tachyonic” superconductor tac  [see the phase diagram in Fig. 1(b)], in
which the vacuum of the electromagnetic field lies at the local maximum of the
effective potential $V_{\rm eff}(\boldsymbol{A})\propto
m^{2}\boldsymbol{A}^{2}$ [Fig. 1(a)]. The electromagnetic field in the
tachyonic superconductor becomes unstable, and starts to grow exponentially in
time. The repulsive case ($U>0$) does not have such a tachyonic instability,
but shows a different type of instability, as discussed below.

If we look at the electromagnetic response closer, we find that there is a
phase transition at $U=\pm 4t_{h}$ [Fig. 1(b)]. This can be seen from the
behavior of the Meissner kernel represented in real space
$K^{\perp}(j,\omega=0)$ sup , which asymptotically decays exponentially as
$\sim\exp(-j/\xi)$ with $\xi=1/\cosh^{-1}(U^{2}/8t_{h}^{2}-1)$ for
$|U|>4t_{h}$ sup . Here $\xi$ is Pippard’s coherence length Schrieffer (1983),
which diverges at $U_{c}=\pm 4t_{h}$ as $\xi\sim|U-U_{c}|^{-1/2}$. For
$|U|<4t_{h}$, the kernel shows a power-law decay as $\sim j^{-1/2}$ sup .

Let us compare the coherence length $\xi$ with London’s penetration depth
$\lambda$ defined by $\frac{1}{\lambda^{2}}=\frac{\mu_{0}}{\pi}D_{s}$
($\mu_{0}$ is the vacuum permeability). For $U>0$, we have
$\lambda=\sqrt{\frac{U}{16\mu_{0}e^{2}t_{h}^{2}C_{M,N}}}$, which grows
smoothly as $U$ increases. In the region of $0<U<4t_{h}$, the penetration
depth is smaller than the coherence length ($\lambda<\xi$), and the system
belongs to type-I superconductors (Fig. 1). For $U>4t_{h}$, on the other hand,
the relation becomes opposite ($\lambda>\xi$), and the system turns to a type-
II superconductor (Fig. 1) typ . In an analogous way, we call the region
$-4t_{h}<U<0$ ($U<-4t_{h}$) a type-I (type-II) tachyonic superconductor (Fig.
1). They have different magnetic properties (for details, see sup ).

Figure 2: (a), (b): Energy dispersion of the transverse electromagnetic field
coupled to the $\eta$-pairing states in the Hubbard model. The solid curves
show the real frequencies, while the dashed ones represent the real part of
the complex frequencies. The sloped dashed line shows the dispersion relation
in the vacuum ($\omega=cq$), and the horizontal dashed lines correspond to
$\omega=|U|$. (c) Real part of $\omega$ at $q=q_{\ast}$, where the imaginary
part takes the maximal value. The dashed line shows $\omega=U$. (d), (e):
Inverse of the imaginary part of $\omega$ corresponding to the time scale of
the growth of the dynamical instability. (f): Inverse of the imaginary part of
$\omega$ at $q=q_{\ast}$. We set $\rho=0.5$ in (a), (b), (d), (e) and
$t_{h}=1$ [eV] and $a=1$ [Å] in (a)-(f).

Dynamical instability of $\eta$-pairing states.— Now, let us study the
dynamics of electromagnetic fields coupled to the $\eta$-pairing states for
$U>0$. To this end, we consider the Maxwell equation in the Lorenz gauge,
$-\frac{\omega^{2}}{c^{2}}\boldsymbol{A}+q^{2}\boldsymbol{A}=\mu_{0}\boldsymbol{j}$,
combined with the response of the $\eta$-pairing states,
$j^{\mu}=-K^{\mu\nu}(\boldsymbol{q},\omega)A_{\nu}$. The equation of motion
determines the energy dispersion of collective modes of electromagnetic fields
coupled with the $\eta$-pairing states. We focus on the transverse mode, whose
energy dispersion is given by

$\displaystyle\frac{\omega^{2}}{c^{2}}-q^{2}$
$\displaystyle=\mu_{0}K^{\perp}(\boldsymbol{q},\omega).$ (6)

At $\boldsymbol{q}=0$, Eq. (6) becomes
$\frac{\omega^{2}}{c^{2}}=-\frac{16e^{2}t_{h}^{2}C_{M,N}U}{\omega^{2}-U^{2}}$,
which has imaginary-frequency solutions when
$U^{3}<64\mu_{0}e^{2}c^{2}t_{h}^{2}C_{M,N}$. If we input $t_{h}=1$ [eV] and
$a=1$ [Å] for ordinary materials, the condition reads
$U/t_{h}<22.6\cdot(\rho(1-\rho))^{1/3}$, where $\rho:=(N/2)/M$ is the number
of doublons per site ($0\leq\rho\leq 1$). Surprisingly, the $\eta$-pairing
states coupled to electromagnetic fields are dynamically unstable over a wide
range of the parameter space against $\boldsymbol{q}=0$ modes. More generally,
we find that the $\eta$-pairing states are unstable for all the parameters if
we take into account arbitrary $\boldsymbol{q}$ modes.

In Fig. 2, we plot the numerical solutions of Eq. (6) for various parameters.
In Figs. 2(a) ($U>0$) and (b) ($U<0$), the solid curves show the real-
frequency solutions, while the dashed curves represent the real part of the
complex-frequency solutions. When $U$ is positive and sufficiently large,
there are two branches of the real solutions with the gaps near $q=0$. As $q$
increases, the two branches merge at some point, and turn into a conjugate
pair of complex frequencies non . After going across the vacuum dispersion
($\omega=cq$), the solutions become real and split into two branches again. At
high momentum, the two branches approach $\omega=cq$ and $\omega=U$. For
$0<U/t_{h}<14.3$ ($\rho=0.5$), the real branches near $q=0$ vanish as
discussed above, and only complex solutions exist at low momentum. For $U<0$
[Fig. 2(b)], the dispersion shows a tachyonic spectrum. In general, we can
prove that complex frequencies appear for all $U$ and $\rho$ sup , indicating
that the electromagnetic field (and hence the $\eta$-pairing state) is always
dynamically unstable.

In Figs. 2(d) and (e), we plot $[{\rm Im}\,\omega]^{-1}$, i.e., the time scale
of the growth of the instability. One can see that the shortest time scale
among the $\boldsymbol{q}$ modes (whose momentum is denoted by $q_{\ast}$) is
of the order of $\hbar/t_{h}$, which is in the femtosecond regime. In the
decaying process, the energy of the electromagnetic field is transferred from
the binding energy of doublons for $U>0$, and from the kinetic energy of
doublons for $U<0$. In the former, the doublons break up into two particles,
while in the latter the doublons lose their momentum $\boldsymbol{Q}$. In both
cases, the $\eta$-pairing states will eventually disappear.

For $U>0$, the complex frequencies have nonzero real parts [see Fig. 2(a)], so
that the exponential growth of the electromagnetic field is accompanied by
plasma oscillations. They induce intense light emission, where doublons’
binding energies are released collectively. In Figs. 2(c) and (f), we plot
${\rm Re}\,\omega$ and $[{\rm Im}\,\omega]^{-1}$ at momentum $q_{\ast}$,
corresponding to the characteristic frequency and the growing time scale of
the dominant emitted light waves, respectively. The characteristic frequency
is shifted from $U$. In particular, at $U\ll t_{h}$ it is proportional to
$U^{1/4}$ with $q_{\ast}=0$ sup . As $U$ increases, $q_{\ast}$ starts to take
a nonzero value around $U/t_{h}\approx 9$ [Fig. 2(d)], making a kink-like
structure in Fig. 2(c). The time scale of the growth increases as the density
decreases [Fig. 2(f)], but stays within the femtosecond regime even at
$\rho=0.01$. For $U<0$, the decay of the $\eta$-pairing states is accompanied
by unattenuated penetration of electromagnetic fields sup  in such a way that
a tachyonic field grows exponentially as in order-parameter dynamics near
critical points Kibble (1976); Zurek (1985); Polkovnikov  _et al._ (2011).

We emphasize that the mechanism of the plasma instability at $U>0$ is
different from that of spontaneous light emission, the latter of which is
caused by a quantum-mechanical effect of electromagnetic fields and has the
frequency $\omega=U$. The decay width of spontaneous emission is evaluated by
$\Gamma=\frac{\omega}{3\pi\varepsilon_{0}\hbar
c^{3}}\sum_{n}|\langle\psi_{N}|J^{\mu}|\phi_{n}\rangle|^{2}$ Loudon (2000),
where $\varepsilon_{0}$ is the vacuum permittivity and the sum runs over all
the eigenstates of the Hubbard model. For $\rho=0.5$, $t_{h}=U=1$ [eV], and
$a=1$ [Å], we have $\Gamma/M=2.3\times 10^{7}$ [s-1] sup . Thus, spontaneous
emission takes place for each site in the time scale of $10^{2}$ [ns], which
is much slower than the plasma instability.

Summary and outlook.— We have shown that Yang’s $\eta$-pairing states have the
intrinsic plasma instability for $U>0$ and the tachyonic instability for $U<0$
when the system is coupled to electromagnetic fields. The time scales of both
of these instabilities are of the order of femtoseconds, which puts a strong
constraint on the realization of the $\eta$-pairing states in real materials.
The decay of the $\eta$-pairing states leads to enhanced light emission with
characteristic frequencies shifted from the Hubbard interaction $U$ for the
repulsive case, and unattenuated penetration of electromagnetic fields for the
attractive case. While we have focused on the simplest form of the
$\eta$-pairing eigenstates (2), we expect that similar instabilities might
exist for more general states having unpaired particles (at least if they are
dilute enough). Stabilizing the $\eta$-pairing states coupled to
electromagnetic fields is an interesting open problem, which merits further
studies.

###### Acknowledgements.

N.T. acknowledges support by KAKENHI Grant No. JP20K03811. M.N. acknowledges
support by KAKENHI Grant No. JP20K14383. M.U. acknowledges support by KAKENHI
Grant No. JP18H01145.

## References

  * Fausti  _et al._ (2011) D. Fausti, R. I. Tobey, N. Dean, S. Kaiser, A. Dienst, M. C. Hoffmann, S. Pyon, T. Takayama, H. Takagi, and A. Cavalleri, Science 331, 189 (2011).
  * Kaiser  _et al._ (2014) S. Kaiser, C. R. Hunt, D. Nicoletti, W. Hu, I. Gierz, H. Y. Liu, M. Le Tacon, T. Loew, D. Haug, B. Keimer, and A. Cavalleri, Phys. Rev. B 89, 184516 (2014).
  * Hu  _et al._ (2014) W. Hu, S. Kaiser, D. Nicoletti, C. R. Hunt, I. Gierz, M. C. Hoffmann, M. Le Tacon, T. Loew, B. Keimer, and A. Cavalleri, Nat. Mater. 13, 705 (2014).
  * Mitrano  _et al._ (2016) M. Mitrano, A. Cantaluppi, D. Nicoletti, S. Kaiser, A. Perucchi, S. Lupi, P. Di Pietro, D. Pontiroli, M. Riccò, S. R. Clark, D. Jaksch, and A. Cavalleri, Nature 530, 461 (2016).
  * Cantaluppi  _et al._ (2018) A. Cantaluppi, M. Buzzi, G. Jotzu, D. Nicoletti, M. Mitrano, D. Pontiroli, M. Riccò, A. Perucchi, P. Di Pietro, and A. Cavalleri, Nat. Phys. 14, 837 (2018).
  * Buzzi  _et al._ (2020) M. Buzzi, D. Nicoletti, M. Fechner, N. Tancogne-Dejean, M. A. Sentef, A. Georges, T. Biesner, E. Uykur, M. Dressel, A. Henderson, T. Siegrist, J. A. Schlueter, K. Miyagawa, K. Kanoda, M.-S. Nam, A. Ardavan, J. Coulthard, J. Tindall, F. Schlawin, D. Jaksch, and A. Cavalleri, Phys. Rev. X 10, 031028 (2020).
  * (7) M. Budden, T. Gebert, M. Buzzi, G. Jotzu, E. Wang, T. Matsuyama, G. Meier, Y. Laplace, D. Pontiroli, M. Riccò, F. Schlawin, D. Jaksch, A. Cavalleri, “Evidence for metastable photo-induced superconductivity in K3C60”, arXiv:2002.12835.
  * Oike  _et al._ (2018) H. Oike, M. Kamitani, Y. Tokura, and F. Kagawa, Sci. Adv. 4, eaau3489 (2018).
  * Niwa  _et al._ (2019) H. Niwa, N. Yoshikawa, K. Tomari, R. Matsunaga, D. Song, H. Eisaki, and R. Shimano, Phys. Rev. B 100, 104507 (2019).
  * Zhang  _et al._ (2020) S. J. Zhang, Z. X. Wang, H. Xiang, X. Yao, Q. M. Liu, L. Y. Shi, T. Lin, T. Dong, D. Wu, and N. L. Wang, Phys. Rev. X 10, 011056 (2020).
  * Matsubara  _et al._ (2007) M. Matsubara, Y. Okimoto, T. Ogasawara, Y. Tomioka, H. Okamoto, and Y. Tokura, Phys. Rev. Lett. 99, 207401 (2007).
  * McLeod  _et al._ (2020) A. S. McLeod, J. Zhang, M. Q. Gu, F. Jin, G. Zhang, K. W. Post, X. G. Zhao, A. J. Millis, W. B. Wu, J. M. Rondinelli, R. D. Averitt, and D. N. Basov, Nat. Mater. 19, 397 (2020).
  * Stojchevska  _et al._ (2014) L. Stojchevska, I. Vaskivskyi, T. Mertelj, P. Kusar, D. Svetin, S. Brazovskii, and D. Mihailovic, Science 344, 177 (2014).
  * Vaskivskyi  _et al._ (2015) I. Vaskivskyi, J. Gospodaric, S. Brazovskii, D. Svetin, P. Sutar, E. Goreshnik, I. A. Mihailovic, T. Mertelj, and D. Mihailovic, Sci. Adv. 1, e1500168 (2015).
  * Sun  _et al._ (2018) K. Sun, S. Sun, C. Zhu, H. Tian, H. Yang, and J. Li, Sci. Adv. 4, eaas9660 (2018).
  * Yang (1989) C. N. Yang, Phys. Rev. Lett. 63, 2144 (1989).
  * Vafek  _et al._ (2017) O. Vafek, N. Regnault, and B. A. Bernevig, SciPost Phys. 3, 043 (2017).
  * Mark and Motrunich (2020) D. K. Mark and O. I. Motrunich, Phys. Rev. B 102, 075132 (2020).
  * Moudgalya  _et al._ (2020) S. Moudgalya, N. Regnault, and B. A. Bernevig, Phys. Rev. B 102, 085140 (2020).
  * Kitamura and Aoki (2016) S. Kitamura and H. Aoki, Phys. Rev. B 94, 174503 (2016).
  * Peronaci  _et al._ (2020) F. Peronaci, O. Parcollet, and M. Schiró, Phys. Rev. B 101, 161101 (2020).
  * Cook and Clark (2020) M. W. Cook and S. R. Clark, Phys. Rev. A 101, 033604 (2020).
  * Tindall  _et al._ (2021) J. Tindall, F. Schlawin, M. A. Sentef, and D. Jaksch, Phys. Rev. B 103, 035146 (2021).
  * Kaneko  _et al._ (2019) T. Kaneko, T. Shirakawa, S. Sorella, and S. Yunoki, Phys. Rev. Lett. 122, 077002 (2019).
  * Kaneko  _et al._ (2020) T. Kaneko, S. Yunoki, and A. J. Millis, Phys. Rev. Research 2, 032027 (2020).
  * Werner  _et al._ (2019) P. Werner, J. Li, D. Golež, and M. Eckstein, Phys. Rev. B 100, 155130 (2019).
  * Li  _et al._ (2020) J. Li, D. Golez, P. Werner, and M. Eckstein, Phys. Rev. B 102, 165136 (2020).
  * Diehl  _et al._ (2008) S. Diehl, A. Micheli, A. Kantian, B. Kraus, H. P. Büchler, and P. Zoller, Nat. Phys. 4, 878 (2008).
  * Kraus  _et al._ (2008) B. Kraus, H. P. Büchler, S. Diehl, A. Kantian, A. Micheli, and P. Zoller, Phys. Rev. A 78, 042307 (2008).
  * Bernier  _et al._ (2013) J.-S. Bernier, P. Barmettler, D. Poletti, and C. Kollath, Phys. Rev. A 87, 063608 (2013).
  * Tindall  _et al._ (2019) J. Tindall, B. Buča, J. R. Coulthard, and D. Jaksch, Phys. Rev. Lett. 123, 030603 (2019).
  * (32) M. Nakagawa, N. Tsuji, N. Kawakami, and M. Ueda, to be published.
  * Su  _et al._ (1991) G. Su, M.-L. Ge, and B.-H. Zhao, Phys. Lett. A 155, 425 (1991).
  * Su  _et al._ (1992) G. Su, M.-L. Ge, and B.-H. Zhao, Phys. Stat. Sol. B 173, 709 (1992).
  * Yang and Zhang (1990) C. N. Yang and S. C. Zhang, Mod. Phys. Lett. B 4, 759 (1990).
  * Yang (1962) C. N. Yang, Rev. Mod. Phys. 34, 694 (1962).
  * Schrieffer (1983) J. R. Schrieffer,  _Theory of Superconductivity_ (Perseus, 1983).
  * (38) The diamagnetic component of the electromagnetic response function for the $\eta$-pairing states in the Hubbard model is given by $K_{\rm dia}^{\mu\nu}(\boldsymbol{q},\omega)$ $\propto\delta_{\mu\nu}\sum_{\boldsymbol{k}\sigma}\frac{\partial^{2}\varepsilon_{\boldsymbol{k}}}{\partial k_{\mu}^{2}}\langle\psi_{N}|c_{\boldsymbol{k}\sigma}^{\dagger}c_{\boldsymbol{k}\sigma}|\psi_{N}\rangle\propto\delta_{\mu\nu}\sum_{\boldsymbol{k}\sigma}\varepsilon_{\boldsymbol{k}}\langle\psi_{N}|c_{\boldsymbol{k}\sigma}^{\dagger}c_{\boldsymbol{k}\sigma}|\psi_{N}\rangle=0$, where $\varepsilon_{\boldsymbol{k}}=-2t_{h}\sum_{\mu=1}^{3}\cos k_{\mu}$ is the single-particle band dispersion and $c_{\boldsymbol{k}\sigma}=\frac{1}{\sqrt{M}}\sum_{j}e^{-i\boldsymbol{k}\cdot\boldsymbol{R}_{j}}c_{j\sigma}$.
  * (39) See Supplementary Material for the details of the derivation of the electromagnetic response function for the $\eta$-pairing states and its properties including the symmetry constraint, the solution of the Maxwell equations, magnetic properties of tachyonic superconductors, and spontaneous light emission.
  * Essler  _et al._ (2005) F. H. L. Essler, H. Frahm, F. Göhmann, A. Klümper, and V. E. Korepin,  _The One-Dimensional Hubbard Model_ (Cambridge University Press, 2005).
  * Scalapino  _et al._ (1993) D. J. Scalapino, S. R. White, and S. Zhang, Phys. Rev. B 47, 7995 (1993).
  * Anderson (1963) P. W. Anderson, Phys. Rev. 130, 439 (1963).
  * Shimano and Tsuji (2020) R. Shimano and N. Tsuji, Annu. Rev. Condens. Matter Phys. 11, 103 (2020).
  * (44) Here, by “tachyonic”, we mean that the vacuum state is unstable since the system sits on a local maximum of a potential, and do not mean a hypothetical superluminal particle. When an ODLRO is present, a uniform magnetic field cannot exist Sewell (1990); Nieh  _et al._ (1995). In other words, the energy dispersion of the electromagnetic field cannot cross the origin ($\omega=\boldsymbol{q}=0$). To satisfy this condition, the system must be either a superconductor (with the dispersion of the electromagnetic field $\omega^{2}\simeq c^{2}q^{2}+m^{2}c^{4}$ with $m^{2}>0$) or a tachyonic superconductor ($\omega^{2}\simeq c^{2}q^{2}+m^{2}c^{4}$ with $m^{2}<0$).
  * (45) The classification is rather formal here. To identify the nature of vortices in $\eta$ pairing states, one has to go beyond the linear-response theory of the present analysis. The point at which $\lambda$ exceeds $\xi$ is, precisely speaking, not exactly at $|U|=4t_{h}$ but very close to it sup .
  * (46) This behavior is reminiscent of exceptional points in non-Hermitian systems Ash .
  * Kibble (1976) T. W. B. Kibble, J. Phys. A: Math. Gen. 9, 1387 (1976).
  * Zurek (1985) W. H. Zurek, Nature 317, 505 (1985).
  * Polkovnikov  _et al._ (2011) A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Rev. Mod. Phys. 83, 863 (2011).
  * Loudon (2000) R. Loudon,  _The Quantum Theory of Light_ , 3rd ed. (Oxford University Press, New York, 2000).
  * Sewell (1990) G. L. Sewell, J. Stat. Phys. 61, 415 (1990).
  * Nieh  _et al._ (1995) H. T. Nieh, G. Su, and B.-H. Zhao, Phys. Rev. B 51, 3760 (1995).
  * (53) Y. Ashida, Z. Gong, and M. Ueda, “Non-Hermitian Physics”, arXiv:2006.01837.
  * Gradshteyn and Ryzhik (1995) I. S. Gradshteyn and I. M. Ryzhik,  _Table of Integrals, Series, and Products_ (Academic Press, New York, 1995).

Supplemental Material for “Tachyonic and Plasma Instabilities of
$\eta$-Pairing States Coupled to Electromagnetic Fields”
Naoto Tsuji1,2, Masaya Nakagawa1, and Masahito Ueda1,2,3

1Department of Physics, University of Tokyo, Hongo, Tokyo 113-0033, Japan
2RIKEN Center for Emergent Matter Science (CEMS), Wako 351-0198, Japan
3Institute for Physics of Intelligence, University of Tokyo, Hongo, Tokyo
113-0033, Japan
(Dated: )

##  I I. Derivation of the electromagnetic response function

In this section, we describe how to analytically evaluate the electromagnetic
response function (Meissner kernel) [Eq. (3) in the main text],

$\displaystyle K^{\mu\nu}(\boldsymbol{R}_{j},t)$
$\displaystyle=-i\theta(t)\langle\psi_{N}|[J^{\mu}(\boldsymbol{R}_{j},t),J^{\nu}(0,0)]|\psi_{N}\rangle,$
(S1)

at arbitrary lattice coordinate $\boldsymbol{R}_{j}$ and time $t$ for Yang’s
$\eta$-pairing state $|\psi_{N}\rangle$ [Eq. (4)] in the Hubbard model.

###  I.1 A. Reduction to the two-particle correlation function

The first step is to reduce the correlation function of $N$ particles (S1) to
that of two particles by shifting all the $\eta$ operators in
$|\psi_{N}\rangle$ to the left of the current operators using the commutation
relation

$\displaystyle[J^{\mu}(\boldsymbol{R}_{j}),\eta^{+}]$
$\displaystyle=-2iet_{h}e^{i\boldsymbol{Q}\cdot\boldsymbol{R}_{j}}(c_{j+\mu\uparrow}^{\dagger}c_{j\downarrow}^{\dagger}+c_{j\uparrow}^{\dagger}c_{j+\mu\downarrow}^{\dagger}).$
(S2)

For convenience, we define the operator that appeared on the right-hand side
of Eq. (S2) as

$\displaystyle J_{\eta}^{\mu+}(\boldsymbol{R}_{j})$
$\displaystyle:=-iet_{h}e^{i\boldsymbol{Q}\cdot\boldsymbol{R}_{j}}(c_{j+\mu\uparrow}^{\dagger}c_{j\downarrow}^{\dagger}+c_{j\uparrow}^{\dagger}c_{j+\mu\downarrow}^{\dagger}).$
(S3)

With this definition, we can write

$\displaystyle[J^{\mu}(\boldsymbol{R}_{j}),\eta^{+}]$
$\displaystyle=2J_{\eta}^{\mu+}(\boldsymbol{R}_{j}).$ (S4)

Since $J_{\eta}^{\mu+}(\boldsymbol{R}_{j})$ only involves creation operators,
$J_{\eta}^{\mu+}(\boldsymbol{R}_{j})$ commutes with $\eta^{+}$. Therefore, we
can repeatedly use the relation (S4) to obtain

$\displaystyle[J^{\mu}(\boldsymbol{R}_{j}),(\eta^{+})^{n}]$
$\displaystyle=2n(\eta^{+})^{n-1}J_{\eta}^{\mu+}(\boldsymbol{R}_{j})$ (S5)

for $n=1,2,\dots$. By using Eq. (S5), we can evaluate the current-current
correlation function as

$\displaystyle\langle\psi_{N}|J^{\mu}(\boldsymbol{R}_{j},t)J^{\nu}(0,0)|\psi_{N}\rangle$
$\displaystyle=\frac{1}{\sqrt{\mathcal{N}_{N}}}\langle\psi_{N}|J^{\mu}(\boldsymbol{R}_{j},t)J^{\nu}(0,0)(\eta^{+})^{\frac{N}{2}}|0\rangle$
$\displaystyle=\frac{1}{\sqrt{\mathcal{N}_{N}}}N\langle\psi_{N}|J^{\mu}(\boldsymbol{R}_{j},t)(\eta^{+})^{\frac{N}{2}-1}J_{\eta}^{\nu+}(0)|0\rangle$
$\displaystyle=\frac{1}{\sqrt{\mathcal{N}_{N}}}N\langle\psi_{N}|(\eta^{+})^{\frac{N}{2}-1}J^{\mu}(\boldsymbol{R}_{j},t)J_{\eta}^{\nu+}(0)|0\rangle$
$\displaystyle\quad+\frac{1}{\sqrt{\mathcal{N}_{N}}}N(N-2)\langle\psi_{N}|(\eta^{+})^{\frac{N}{2}-2}J_{\eta}^{\mu+}(\boldsymbol{R}_{j},t)J_{\eta}^{\nu+}(0)|0\rangle.$
(S6)

Then, straightforward calculations show that

$\displaystyle\langle\psi_{N}|(\eta^{+})^{\frac{N}{2}-1}$
$\displaystyle=\frac{1}{\sqrt{\mathcal{N}_{N}}}\frac{(\frac{N}{2})!(M-1)!}{(M-\frac{N}{2})!}\langle
0|\eta^{-},$ (S7) $\displaystyle\langle\psi_{N}|(\eta^{+})^{\frac{N}{2}-2}$
$\displaystyle=\frac{1}{\sqrt{\mathcal{N}_{N}}}\frac{(\frac{N}{2})!(M-2)!}{2(M-\frac{N}{2})!}\langle
0|(\eta^{-})^{2},$ (S8)

which can be used to rewrite the correlation function (S6) as

$\displaystyle\langle\psi_{N}|J^{\mu}(\boldsymbol{R}_{j},t)J^{\nu}(0,0)|\psi_{N}\rangle$
$\displaystyle=\frac{N}{\mathcal{N}_{N}}\frac{(\frac{N}{2})!(M-1)!}{(M-\frac{N}{2})!}\langle
0|\eta^{-}J^{\mu}(\boldsymbol{R}_{j},t)J_{\eta}^{\nu+}(0)|0\rangle$
$\displaystyle+\frac{N(N-2)}{\mathcal{N}_{N}}\frac{(\frac{N}{2})!(M-2)!}{2(M-\frac{N}{2})!}\langle
0|(\eta^{-})^{2}J_{\eta}^{\mu+}(\boldsymbol{R}_{j},t)J_{\eta}^{\nu+}(0)|0\rangle.$
(S9)

To further simplify the expression, we use the commutation relations,

$\displaystyle[\eta^{-},J^{\mu}(\boldsymbol{R}_{j})]$
$\displaystyle=2J_{\eta}^{\mu-}(\boldsymbol{R}_{j}),$ (S10)
$\displaystyle[\eta^{-},J_{\eta}^{\mu+}(\boldsymbol{R}_{j})]$
$\displaystyle=-J^{\mu}(\boldsymbol{R}_{j}),$ (S11)

where
$J_{\eta}^{\mu-}(\boldsymbol{R}_{j}):=[J_{\eta}^{\mu+}(\boldsymbol{R}_{j})]^{\dagger}$.
Thus, the correlation function becomes

$\displaystyle\langle\psi_{N}|J^{\mu}(\boldsymbol{R}_{j},t)J^{\nu}(0,0)|\psi_{N}\rangle$
$\displaystyle=\frac{2N}{\mathcal{N}_{N}}\frac{(\frac{N}{2})!(M-1)!}{(M-\frac{N}{2})!}\langle
0|J_{\eta}^{\mu-}(\boldsymbol{R}_{j},t)J_{\eta}^{\nu+}(0)|0\rangle$
$\displaystyle\quad-\frac{N(N-2)}{\mathcal{N}_{N}}\frac{(\frac{N}{2})!(M-2)!}{2(M-\frac{N}{2})!}\langle
0|\eta^{-}J^{\mu}(\boldsymbol{R}_{j},t)J_{\eta}^{\nu+}(0)|0\rangle$
$\displaystyle=\frac{2N(\frac{N}{2})!(M-2)!}{(M-\frac{N}{2}-1)!\mathcal{N}_{N}}\langle
0|J_{\eta}^{\mu-}(\boldsymbol{R}_{j},t)J_{\eta}^{\nu+}(0)|0\rangle.$ (S12)

Here, let us recall that the normalization constant for the $\eta$-pairing
state is explicitly given by

$\displaystyle\mathcal{N}_{N}$
$\displaystyle=\frac{(\frac{N}{2})!M!}{(M-\frac{N}{2})!}.$ (S13)

We use Eq. (S12) to reduce the current-current correlation function of $N$
particles to that of two particles,

$\displaystyle\langle\psi_{N}|J^{\mu}(\boldsymbol{R}_{j},t)J^{\nu}(0,0)|\psi_{N}\rangle$
$\displaystyle=4C_{M,N}\langle
0|J_{\eta}^{\mu-}(\boldsymbol{R}_{j},t)J_{\eta}^{\nu+}(0)|0\rangle,$ (S14)

where

$\displaystyle C_{M,N}$
$\displaystyle:=\frac{\frac{N}{2}(M-\frac{N}{2})}{M(M-1)}.$ (S15)

Equation (S14) leads to Eq. (4) in the main text.

The technique used here (i.e., reduction of $N$-particle to few-particle
correlation functions) can be applied not only to the electromagnetic response
function (S1) but also to arbitrary correlation functions constructed from
few-body operators.

###  I.2 B. Evaluation of the two-particle dynamics

In the previous subsection, we have shown that the $N$-particle correlation
function (S1) can be reduced to the two-particle correlation function (S14).
Since the two-particle problem in the Hubbard model is exactly solvable, we
can evaluate the two-particle correlation function exactly. Here we describe
the details of the evaluation.

First, we Fourier transform Eq. (S14) to obtain

$\displaystyle\sum_{j}e^{i\boldsymbol{q}\cdot\boldsymbol{R}_{j}}\langle\psi_{N}|J^{\mu}(\boldsymbol{R}_{j},t)J^{\nu}(0,0)|\psi_{N}\rangle$
$\displaystyle=\frac{4}{M}C_{M,N}\langle
0|J_{\eta}^{\mu-}(\boldsymbol{q},t)J_{\eta}^{\nu+}(\boldsymbol{q})|0\rangle,$
(S16)

where

$\displaystyle J_{\eta}^{\mu+}(\boldsymbol{q})$
$\displaystyle:=\sum_{j}e^{-i\boldsymbol{q}\cdot\boldsymbol{R}_{j}}J_{\eta}^{\mu+}(\boldsymbol{R}_{j}),$
(S17) $\displaystyle J_{\eta}^{\mu-}(\boldsymbol{q})$
$\displaystyle:=[J_{\eta}^{\mu+}(\boldsymbol{q})]^{\dagger}.$ (S18)

Acting $J_{\eta}^{\nu+}(\boldsymbol{q})$ on the vacuum state, we obtain

$\displaystyle
J_{\eta}^{\nu+}(\boldsymbol{q})|0\rangle=-iet_{h}\sum_{j}e^{i(\boldsymbol{Q}-\boldsymbol{q})\cdot\boldsymbol{R}_{j}}(c_{j+\nu\uparrow}^{\dagger}c_{j\downarrow}^{\dagger}+c_{j\uparrow}^{\dagger}c_{j+\nu\downarrow}^{\dagger})|0\rangle$
$\displaystyle=-iet_{h}\sum_{j}e^{i(\boldsymbol{Q}-\boldsymbol{q})\cdot\boldsymbol{R}_{j}}(|\boldsymbol{R}_{j}+\boldsymbol{e}_{\nu},\boldsymbol{R}_{j}\rangle+|\boldsymbol{R}_{j},\boldsymbol{R}_{j}+\boldsymbol{e}_{\nu}\rangle),$
(S19)

where we have introduced the notation
$|\boldsymbol{R}_{j},\boldsymbol{R}_{k}\rangle:=c_{j\uparrow}^{\dagger}c_{k\downarrow}^{\dagger}|0\rangle$
to represent a two-particle state with $\uparrow$ spin at site $j$ and
$\downarrow$ spin at site $k$, and $\boldsymbol{e}_{\nu}$ is the unit vector
in the $\nu$ direction.

Since the system has the (discrete) translation symmetry, the center-of-mass
(crystal) momentum of two particles under consideration is conserved. Let us
define the translation operator $T_{\mu}$ that shifts two particles by one
lattice site in the $\mu$ direction. We express the eigenstates of $T_{\mu}$
in terms of the center-of-mass momentum $\boldsymbol{K}$ and the relative
coordinate $\boldsymbol{r}$ of the two particles as

$\displaystyle|\boldsymbol{K},\boldsymbol{r}\rangle$
$\displaystyle=\frac{1}{\sqrt{M}}\sum_{j}e^{i\boldsymbol{K}\cdot\boldsymbol{R}_{j}}|\boldsymbol{R}_{j}+\boldsymbol{r},\boldsymbol{R}_{j}\rangle,$
(S20)

which satisfies
$T_{\mu}|\boldsymbol{K},\boldsymbol{r}\rangle=e^{-i\boldsymbol{K}\cdot\boldsymbol{e}_{\mu}}|\boldsymbol{K},\boldsymbol{r}\rangle$.
Using the eigenstates (S20), $J_{\eta}^{\nu+}(\boldsymbol{q})|0\rangle$ (S19)
can be written as

$\displaystyle J_{\eta}^{\nu+}(\boldsymbol{q})|0\rangle$
$\displaystyle=-iet_{h}\sqrt{M}(|\boldsymbol{Q}-\boldsymbol{q},+\boldsymbol{e}_{\nu}\rangle$
$\displaystyle\quad+e^{-i(\boldsymbol{Q}-\boldsymbol{q})\cdot\boldsymbol{e}_{\nu}}|\boldsymbol{Q}-\boldsymbol{q},-\boldsymbol{e}_{\nu}\rangle),$
(S21)

which has the center-of-mass momentum $\boldsymbol{Q}-\boldsymbol{q}$. The
action of the Hamiltonian in Eq. (1) in the main text on the state
$|\boldsymbol{K},\boldsymbol{r}\rangle$ (S20) is given by

$\displaystyle H|\boldsymbol{K},\boldsymbol{r}\rangle$
$\displaystyle=-t_{h}\sum_{\mu}[(1+e^{i\boldsymbol{K}\cdot\boldsymbol{e}_{\mu}})|\boldsymbol{K},\boldsymbol{r}+\boldsymbol{e}_{\mu}\rangle$
$\displaystyle\quad+(1+e^{-i\boldsymbol{K}\cdot\boldsymbol{e}_{\mu}})|\boldsymbol{K},\boldsymbol{r}-\boldsymbol{e}_{\mu}\rangle]$
$\displaystyle\quad+U(\delta_{\boldsymbol{r},0}-1)|\boldsymbol{K},\boldsymbol{r}\rangle.$
(S22)

If we define operators $\Delta_{\mu}^{\pm}$ that shift the relative coordinate
of two particles by $\pm\boldsymbol{e}_{\mu}$, the Hamiltonian can be
represented in the Hilbert subspace of two particles with the center-of-mass
momentum $\boldsymbol{K}$ as

$\displaystyle H(\boldsymbol{K})$
$\displaystyle=-t_{h}\sum_{\mu}[(1+e^{i\boldsymbol{K}\cdot\boldsymbol{e}_{\mu}})\Delta_{\mu}^{+}+(1+e^{-i\boldsymbol{K}\cdot\boldsymbol{e}_{\mu}})\Delta_{\mu}^{-}]$
$\displaystyle\quad+U(\delta_{\boldsymbol{r},0}-1).$ (S23)

Using the representation (S23), the correlation function (S16) can be written
as

$\displaystyle\sum_{j}e^{i\boldsymbol{q}\cdot\boldsymbol{R}_{j}}\langle\psi_{N}|J^{\mu}(\boldsymbol{R}_{j},t)J^{\nu}(0,0)|\psi_{N}\rangle$
$\displaystyle=\frac{4}{M}C_{M,N}\langle
0|J_{\eta}^{\mu-}(\boldsymbol{q})e^{-iH(\boldsymbol{Q}-\boldsymbol{q})t}J_{\eta}^{\nu+}(\boldsymbol{q})|0\rangle.$
(S24)

Without loss of generality, we assume that $\boldsymbol{q}$ is parallel to the
$z$ direction [$\boldsymbol{q}=(0,0,q)$]. Then,

$\displaystyle H(\boldsymbol{Q}-\boldsymbol{q})$
$\displaystyle=-t_{h}[(1-e^{-iq})\Delta_{z}^{+}+(1-e^{iq})\Delta_{z}^{-}]$
$\displaystyle\quad+U(\delta_{\boldsymbol{r},0}-1).$ (S25)

Hence, the two-particle dynamics that we have to consider is essentially a
one-dimensional problem, in which the relative coordinates $r_{x}$ and $r_{y}$
are conserved. Below, we decompose the electromagnetic response function into
the transverse component
$K^{\perp}(\boldsymbol{q},\omega):=K^{\mu\mu}(\boldsymbol{q},\omega)$
($\mu=x,y$) and the longitudinal one
$K^{\parallel}(\boldsymbol{q},\omega):=K^{zz}(\boldsymbol{q},\omega)$. The
off-diagonal components are absent, i.e.,
$K^{\mu\nu}(\boldsymbol{q},\omega)=0$ for $\mu\neq\nu$, since the state
$e^{-iH(\boldsymbol{Q}-\boldsymbol{q})t}J_{\eta}^{\nu+}(\boldsymbol{q})|0\rangle$
never has an overlap with the state $J_{\eta}^{\mu+}|0\rangle$ for
$\mu\neq\nu$.

###  I.3 C. Transverse component

For the transverse component, the two-particle state
$J_{\eta}^{\nu+}(\boldsymbol{q})|0\rangle$ (S21) has the relative coordinates,
$(r_{x},r_{y})\neq(0,0)$. Since $r_{x}$ and $r_{y}$ are conserved during the
time evolution, the two particles do not sit on the same site. Therefore, they
do not interact with each other, and the dynamics becomes effectively
noninteracting. The correlation function (S16) now reads

$\displaystyle\sum_{j}e^{i\boldsymbol{q}\cdot\boldsymbol{R}_{j}}\langle\psi_{N}|J^{\mu}(\boldsymbol{R}_{j},t)J^{\mu}(0,0)|\psi_{N}\rangle\quad(\mu=x,y)$
$\displaystyle=4e^{2}t_{h}^{2}C_{M,N}\sum_{s=\pm}\langle\boldsymbol{Q}-\boldsymbol{q},s\boldsymbol{e}_{\mu}|e^{-iH(\boldsymbol{Q}-\boldsymbol{q})t}|\boldsymbol{Q}-\boldsymbol{q},s\boldsymbol{e}_{\mu}\rangle,$
(S26)

where

$\displaystyle H(\boldsymbol{Q}-\boldsymbol{q})$
$\displaystyle=-t_{h}[(1-e^{-iq})\Delta_{z}^{+}+(1-e^{iq})\Delta_{z}^{-}]-U$
(S27)

is the noninteracting Hamiltonian which can be diagonalized by Fourier
transformation with respect to the relative coordinate $r_{z}$. The result is

$\displaystyle\sum_{j}e^{i\boldsymbol{q}\cdot\boldsymbol{R}_{j}}\langle\psi_{N}|J^{\mu}(\boldsymbol{R}_{j},t)J^{\mu}(0,0)|\psi_{N}\rangle\quad(\mu=x,y)$
$\displaystyle=8e^{2}t_{h}^{2}C_{M,N}\frac{e^{iUt}}{M}\sum_{\boldsymbol{k}}e^{2it_{h}(\cos
k_{z}-\cos(k_{z}+q))t}$
$\displaystyle=8e^{2}t_{h}^{2}C_{M,N}e^{iUt}\int_{-\pi}^{\pi}\frac{dk_{z}}{2\pi}e^{2it_{h}(\cos
k_{z}-\cos(k_{z}+q))t}$
$\displaystyle=8e^{2}t_{h}^{2}C_{M,N}e^{iUt}J_{0}\left(4t_{h}\sin\frac{q}{2}t\right),$
(S28)

where $J_{0}(x)$ is the zeroth-order Bessel function of the first kind. We
thus obtain the transverse component of the electromagnetic response function
as

$\displaystyle K^{\perp}(\boldsymbol{q},t)$
$\displaystyle=-8i\theta(t)e^{2}t_{h}^{2}C_{M,N}(e^{iUt}-e^{-iUt})$
$\displaystyle\quad\times J_{0}\left(4t_{h}\sin\frac{q}{2}t\right).$ (S29)

Using the integral formula for the Bessel function ($a>0$),

$\displaystyle\int_{0}^{\infty}dte^{i\omega t}J_{0}(at)$
$\displaystyle=\begin{cases}\frac{i{\rm
sgn}(\omega)}{\sqrt{\omega^{2}-a^{2}}}&\mbox{for}\;|\omega|>a;\\\
\frac{1}{\sqrt{a^{2}-\omega^{2}}}&\mbox{for}\;|\omega|<a,\end{cases}$ (S30)

we obtain

$\displaystyle
K^{\perp}(\boldsymbol{q},\omega)=-8ie^{2}t_{h}C_{M,N}\left[f_{q}\left(\frac{\omega+U}{t_{h}}\right)-f_{q}\left(\frac{\omega-U}{t_{h}}\right)\right],$
(S31)

where

$\displaystyle f_{q}(x)$ $\displaystyle:=\begin{cases}\frac{i{\rm
sgn}(x)}{\sqrt{x^{2}-16\sin^{2}\frac{q}{2}}}&\mbox{for}\;|x|>4\sin\frac{q}{2};\\\
\frac{1}{\sqrt{16\sin^{2}\frac{q}{2}-x^{2}}}&\mbox{for}\;|x|<4\sin\frac{q}{2}.\end{cases}$
(S32)

This is the final result for the transverse component in Eq. (5) in the main
text.

Figure S1: Transverse component $K^{\perp}(\boldsymbol{q},\omega)$ of the
electromagnetic response function for the $\eta$-pairing states in the Hubbard
model in units of $4e^{2}t_{h}C_{M,N}$. (a), (b): Real (a) and imaginary (b)
parts of $K^{\perp}(\boldsymbol{q},\omega)$ for $U/t_{h}=2$. (c), (d): Real
(c) and imaginary (d) parts of $K^{\perp}(\boldsymbol{q},\omega)$ for
$U/t_{h}=6$.

In Fig. S1, we plot $K^{\perp}(\boldsymbol{q},\omega)$ for $U/t_{h}=2$ and
$6$. The imaginary part of $K^{\perp}(\boldsymbol{q},\omega)$ represents the
absorption (emission) of light for the $\eta$-pairing states. When a doublon
decays by emitting light with momentum $\boldsymbol{q}$, quasiparticles with
momentum $\boldsymbol{k}$ and $-\boldsymbol{k}+\boldsymbol{Q}-\boldsymbol{q}$
are created. As discussed above, these particles never interact with each
other, so that the total energy of the two particles is given by
$-2t_{h}\sum_{\mu}[(\cos
k_{\mu})+\cos(-k_{\mu}+Q_{\mu}-q_{\mu})]=-4t_{h}\sin\frac{q}{2}\sin(k_{z}+\frac{q}{2})$.
The energy of quasiparticles ranges from $-4t_{h}\sin\frac{q}{2}$ to
$4t_{h}\sin\frac{q}{2}$. Since the energy of a single doublon is $U$, the
condition for light emission to take place is
$|\omega-U|<4t_{h}\sin\frac{q}{2}$. This is exactly the condition of ${\rm
Im}\,K^{\perp}(\boldsymbol{q},\omega)>0$ (for $U/t_{h}>4$). The gap for the
electromagnetic response function closes at $q=\pi$ when $U/t_{h}=4$.

###  I.4 D. Asymptotic behavior at long distance

Here we derive the asymptotic behavior of the transverse electromagnetic
response function (S31) at long distance and low frequency, which is related
to Pippard’s coherence length. In the low-frequency limit, the kernel is given
by

$\displaystyle K^{\perp}(\boldsymbol{q},\omega=0)$
$\displaystyle=16e^{2}t_{h}^{2}C_{M,N}\begin{cases}\frac{{\rm
sgn}(U)}{\sqrt{U^{2}-16t_{h}^{2}\sin^{2}\frac{q}{2}}}&\mbox{for}\;|U|>4t_{h}\sin\frac{q}{2};\\\
0&\mbox{for}\;|U|<4t_{h}\sin\frac{q}{2}.\end{cases}$ (S33)

We consider the problem in two distinct regimes: $|U|>4t_{h}$ and
$|U|<4t_{h}$.

In the first case ($|U|>4t_{h}$), the Fourier transform of
$K^{\perp}(\boldsymbol{q},\omega=0)$ (S33) to real space in the $z$ direction
reads

$\displaystyle K^{\perp}(j,\omega=0)$
$\displaystyle=16e^{2}t_{h}^{2}C_{M,N}\int_{-\pi}^{\pi}\frac{dq}{2\pi}e^{iqj}\frac{{\rm
sgn}(U)}{\sqrt{U^{2}-16t_{h}^{2}\sin^{2}\frac{q}{2}}}.$ (S34)

Precisely speaking, Eq. (S34) shows the electromagnetic response function for
$q_{x}=q_{y}=0$, $R_{z}=j$, and $\omega=0$. The asymptotic behavior of
$K^{\perp}(j,\omega=0)$ for $j\gg 1$ is qualitatively similar to that of
$K^{\perp}(\boldsymbol{R}_{j},\omega=0)$ (i.e., the full Fourier transform of
$K^{\perp}(\boldsymbol{q},\omega=0$) in all directions) for
$|\boldsymbol{R}_{j}|\gg 1$, since the dominant contribution in
$K^{\perp}(j,\omega=0)=\sum_{R_{x},R_{y}}K^{\perp}(R_{x},R_{y},R_{z}=j,\omega=0)$
for $j\gg 1$ arises near $R_{x}\sim R_{y}\sim 0$. In particular, if
$K^{\perp}(\boldsymbol{R}_{j},\omega=0)$ decays exponentially for
$|\boldsymbol{R}_{j}|\gg 1$, then $K^{\perp}(j,\omega=0)$ (S34) also decays
exponentially for $j\gg 1$ with the same correlation length.

The integral in Eq. (S34) can be evaluated analytically as

$\displaystyle K^{\perp}(j,\omega=0)$
$\displaystyle=16e^{2}t_{h}^{2}C_{M,N}\frac{{\rm
sgn}(U)}{\sqrt{U^{2}-16t_{h}^{2}}}$
$\displaystyle\quad\times\frac{{}_{3}F_{2}(\tfrac{1}{2},\tfrac{1}{2},1;1-j,1+j;\tfrac{1}{1-U^{2}/16t_{h}^{2}})}{\Gamma(1-j)\Gamma(1+j)},$
(S35)

where ${}_{3}F_{2}(a_{1},a_{2},a_{3};b_{1},b_{2};z)$ is the generalized
hypergeometric function Gradshteyn and Ryzhik (1995), and $\Gamma(z)$ is the
gamma function. The explicit expression (S35), however, does not directly tell
us about the long-distance behavior. Hence we take a different approach.

Figure S2: Contour for the integral (S39) in the complex plane depicted by the
dashed circle, which can continuously be deformed to the solid curve without
crossing the branch cuts shown by the wavy lines.

The kernel (S34) can be written as

$\displaystyle K^{\perp}(j,\omega=0)$ $\displaystyle=4e^{2}t_{h}C_{M,N}{\rm
sgn}(U)I_{j}(u),$ (S36)

where

$\displaystyle u$ $\displaystyle=\frac{|U|}{4t_{h}}$ (S37)

is the normalized interaction strength, and

$\displaystyle I_{j}(u)$
$\displaystyle:=\int_{-\pi}^{\pi}\frac{dq}{2\pi}e^{-iqj}\frac{1}{\sqrt{u^{2}-\sin^{2}\frac{q}{2}}}$
(S38)

for $u>1$. By putting $z=e^{-iq}$, we transform the integral (S38) to a
complex contour integral,

$\displaystyle I_{j}(u)$ $\displaystyle=\oint\frac{dz}{2\pi
i}z^{j-1}\frac{1}{\sqrt{u^{2}-\frac{1}{2}(1-\frac{1}{2}(z+z^{-1}))}}$
$\displaystyle=2\oint\frac{dz}{2\pi
i}\frac{z^{j-\frac{1}{2}}}{\sqrt{z^{2}+2(2u^{2}-1)z+1}},$ (S39)

where the contour is taken to be the circle around the origin with unit radius
$|z|=1$ (dashed curve in Fig. S2). In the case of $|U|>4t_{h}$, we have $u>1$.
The roots of the quadratic polynomial in the denominator of (S39) are given by

$\displaystyle\alpha_{\pm}$
$\displaystyle=1-2u^{2}\pm\sqrt{(2u^{2}-1)^{2}-1},$ (S40)

which are real numbers. These roots satisfy $\alpha_{-}<-1<\alpha_{+}<0$ and
$\alpha_{+}\alpha_{-}=1$. Using $\alpha_{\pm}$, we can write the contour
integral (S39) as

$\displaystyle I_{j}(u)$ $\displaystyle=2\oint\frac{dz}{2\pi
i}\frac{z^{j-\frac{1}{2}}}{\sqrt{(z-\alpha_{+})(z-\alpha_{-})}}.$ (S41)

In Fig. S2, we show the branch cuts that we adopt in evaluating Eq. (S41) by
wavy lines. With this configuration, the contour can be smoothly deformed to
the solid curve in Fig. S2 without crossing the branch cuts, where the
integral is evaluated as

$\displaystyle I_{j}(u)$
$\displaystyle=2(e^{i(j-\frac{1}{2})\pi+i\pi}-e^{i(j-\frac{1}{2})\pi})$
$\displaystyle\quad\times\int_{\alpha_{+}}^{0}\frac{dx}{2\pi
i}\frac{(-x)^{j-\frac{1}{2}}}{\sqrt{(x-\alpha_{+})(x-\alpha_{-})}}$
$\displaystyle=-4e^{i(j-\frac{1}{2})\pi}\int_{\alpha_{+}}^{0}\frac{dx}{2\pi
i}\frac{(-x)^{j-\frac{1}{2}}}{\sqrt{(x-\alpha_{+})(x-\alpha_{-})}}.$ (S42)

The asymptotic behavior of $I_{j}(u)$ can be read off as follows: for $j\gg
1$, the function $\frac{(-x)^{j-\frac{1}{2}}}{\sqrt{x-\alpha_{+}}}$ in the
integrand has a concentrated contribution near $x=\alpha_{+}$, while
$\frac{1}{\sqrt{x-\alpha_{-}}}$ is a smooth function in $(\alpha_{+},0]$.
Therefore, one can replace $\frac{1}{\sqrt{x-\alpha_{-}}}$ by
$\frac{1}{\sqrt{\alpha_{+}-\alpha_{-}}}$ in the integrand of Eq. (S42),
obtaining

$\displaystyle I_{j}(u)$
$\displaystyle\approx-4e^{i(j-\frac{1}{2})\pi}\frac{1}{\sqrt{\alpha_{+}-\alpha_{-}}}\int_{\alpha_{+}}^{0}\frac{dx}{2\pi
i}\frac{(-x)^{j-\frac{1}{2}}}{\sqrt{x-\alpha_{+}}}$ (S43)

for $j\gg 1$. The rest of the integral can be evaluated as

$\displaystyle I_{j}(u)$
$\displaystyle\approx-4e^{i(j-\frac{1}{2})\pi}\frac{1}{\sqrt{\alpha_{+}-\alpha_{-}}}\frac{1}{2\pi
i}\frac{\sqrt{\pi}\Gamma(j+\frac{1}{2})}{\Gamma(j+1)}(-\alpha_{+})^{j}$
$\displaystyle=\frac{2}{\sqrt{\pi}}\frac{1}{\sqrt{\alpha_{+}-\alpha_{-}}}\frac{\Gamma(j+\frac{1}{2})}{\Gamma(j+1)}\alpha_{+}^{j}.$
(S44)

Using Stirling’s formula, we obtain the asymptotic form of $I_{j}(u)$ as

$\displaystyle I_{j}(u)$
$\displaystyle\approx\frac{2}{\sqrt{\pi}}\frac{1}{\sqrt{\alpha_{+}-\alpha_{-}}}\frac{1}{\sqrt{j}}\alpha_{+}^{j}.$
(S45)

Thus, the transverse electromagnetic response function behaves in the long
distance ($j\to\infty$) as

$\displaystyle K^{\perp}(j,\omega=0)$
$\displaystyle\approx\mbox{cosnt.}\times\frac{\alpha_{+}^{j}}{\sqrt{j}}.$
(S46)

Since $|\alpha_{+}|<1$, $K^{\perp}(j,\omega=0)$ decays exponentially in space
with Pippard’s coherence length $\xi$ defined by

$\displaystyle K^{\perp}(j,\omega=0)$
$\displaystyle\approx\mbox{const.}\times\exp\left(-\frac{j}{\xi}\right).$
(S47)

Physically, $\xi$ represents the length scale over which a response against a
local perturbation of electromagnetic fields propagates in space (Pippard’s
nonlocal electrodynamics Schrieffer (1983)).

From the result (S46), $\xi$ is identified as

$\displaystyle\xi$
$\displaystyle=-\frac{a}{\ln|\alpha_{+}|}=\frac{a}{\cosh^{-1}(2u^{2}-1)},$
(S48)

which does not depend on the doublon density $\rho$. At $|U|=4t_{h}=:U_{c}$,
the coherence length diverges as

$\displaystyle\xi$ $\displaystyle\propto\frac{1}{|U-U_{c}|^{\frac{1}{2}}}.$
(S49)

This is exactly the point where the electromagnetic gap closes at $q=\pi$.

Figure S3: Contour $C$ for the integral (S52) in the complex plane depicted by
the dashed arc, which can continuously be deformed to the solid lines without
crossing the branch cuts as shown by the wavy lines.

In the second case ($|U|<4t_{h}$), the Meissner kernel in the low-frequency
limit is expressed as

$\displaystyle K^{\perp}(j,\omega=0)$ $\displaystyle=4e^{2}t_{h}C_{M,N}{\rm
sgn}(U)I_{j}(u)$ (S50)

with

$\displaystyle I_{j}(u)$
$\displaystyle=\int_{u\geq|\sin\frac{q}{2}|}\frac{dq}{2\pi}e^{-iqj}\frac{1}{\sqrt{u^{2}-\sin^{2}\frac{q}{2}}}$
(S51)

for $u<1$. Similarly to the first case, we put $z=e^{-iq}$ to rewrite the
integral (S51) as

$\displaystyle I_{j}(u)$ $\displaystyle=2\int_{C}\frac{dz}{2\pi
i}\frac{z^{j-\frac{1}{2}}}{\sqrt{(z-\alpha_{+})(z-\alpha_{-})}},$ (S52)

where the roots in the denominator are given by

$\displaystyle\alpha_{\pm}$ $\displaystyle=1-2u^{2}\pm
i\sqrt{1-(2u^{2}-1)^{2}},$ (S53)

and the contour $C$ is taken to be the arc of the circle with the unit radius
connecting $\alpha_{-}$ and $\alpha_{+}$ (dashed curve in Fig. S3). We choose
the branch cuts in the integrand of Eq. (S52) as shown by wavy lines in Fig.
S3.

Following the steepest descent method, we deform the contour from $C$ to the
solid lines in Fig. S3 without crossing the branch cuts, where we put
$z=r\alpha_{\pm}$. Now the integral (S52) can be evaluated as

$\displaystyle I_{j}(u)$ $\displaystyle=2\int_{0}^{1}\frac{dr}{2\pi
i}\alpha_{+}\frac{(r\alpha_{+})^{j-\frac{1}{2}}}{\sqrt{(r\alpha_{+}-\alpha_{+})(r\alpha_{+}-\alpha_{-})}}$
$\displaystyle\quad-2\int_{0}^{1}\frac{dr}{2\pi
i}\alpha_{-}\frac{(r\alpha_{-})^{j-\frac{1}{2}}}{\sqrt{(r\alpha_{-}-\alpha_{+})(r\alpha_{-}-\alpha_{-})}}.$
(S54)

For $j\gg 1$, the function
$\frac{r^{j-\frac{1}{2}}}{\sqrt{r\alpha_{+}-\alpha_{+}}}$ in the first
integral is dominantly contributed from a region near $r=1$, which allows us
to replace $\frac{1}{\sqrt{r\alpha_{+}-\alpha_{-}}}$ by
$\frac{1}{\sqrt{\alpha_{+}-\alpha_{-}}}$ in the integrand. A similar
approximation can be applied to the second term. Taking care of the branch
cuts, we obtain

$\displaystyle I_{j}(u)$
$\displaystyle\approx\frac{2(\alpha_{+}^{j}e^{\frac{i\pi}{4}}-\alpha_{-}^{j}e^{-\frac{i\pi}{4}})}{\sqrt{|\alpha_{+}-\alpha_{-}|}}\int_{0}^{1}\frac{dr}{2\pi
i}\frac{r^{j-\frac{1}{2}}}{\sqrt{1-r}}.$ (S55)

If we define $\alpha_{\pm}=:e^{\pm i\varphi}$, the integral (S51) can be
approximated as

$\displaystyle I_{j}(u)$
$\displaystyle\approx\frac{2\sin(j\varphi+\frac{\pi}{4})}{\sqrt{\pi}\sqrt{|\alpha_{+}-\alpha_{-}|}}\frac{\Gamma(j+\frac{1}{2})}{\Gamma(j+1)}.$
(S56)

Using Stirling’s formula, the asymptotic form of $I_{j}(u)$ for $j\gg 1$ is
given by

$\displaystyle I_{j}(u)$
$\displaystyle\approx\frac{2\sin(j\varphi+\frac{\pi}{4})}{\sqrt{\pi}\sqrt{|\alpha_{+}-\alpha_{-}|}}\frac{1}{\sqrt{j}}.$
(S57)

Therefore, the kernel decays in a long distance according to a power law as

$\displaystyle K^{\perp}(j,\omega=0)$
$\displaystyle\approx\mbox{const.}\times\frac{1}{\sqrt{j}}.$ (S58)

This means that the coherence length diverges ($\xi=\infty$) for $|U|<4t_{h}$.

Figure S4: Comparison between London’s penetration depth $\lambda$ and
Pippard’s coherence length $\xi$ for the $\eta$-pairing states in the Hubbard
model with $t_{h}=1$ [eV], $a=1$ [Å], and $\rho=0.5$.

In Fig. S4, we plot the coherence length $\xi$ for the $\eta$-pairing states
with $t_{h}=1$ [eV], $a=1$ [Å], and $\rho=0.5$ in comparison with London’s
penetration depth defined by
$\lambda=\sqrt{\frac{\hbar^{2}aU}{16\mu_{0}e^{2}t_{h}^{2}C_{M,N}}}$. At $0\leq
U<4t_{h}$, $\lambda$ grows smoothly as a function of $U$, and satisfies
$\lambda<\xi=\infty$. When $U$ exceeds $4t_{h}$, $\xi$ immediately decays to
the order of 1 [Å], whereas $\lambda$ stays on the order of 10 [nm]. The point
at which $\lambda$ becomes equal to $\xi$ is very close to $U=4t_{h}$, beyond
which $\lambda$ becomes larger than $\xi$. Thus, for $0\leq U<4t_{h}$ the
$\eta$-pairing state is a type-I superconductor, whereas for $U>4t_{h}$ the
$\eta$-pairing state is classified to a type-II superconductor.

For $U<0$, we analytically continue $\lambda$ to complex values, which has a
physical meaning as discussed in Sec. IV. We will see that the $\eta$-pairing
state has different magnetic properties depending on whether $|\lambda|$ is
larger than $\xi$ or not. For $-4t_{h}<U<0$, we have the relation
$|\lambda|<\xi$, where the $\eta$-pairing state is called a type-I tachyonic
superconductor (see the main text). For $U<-4t_{h}$, we have $|\lambda|>\xi$,
where the $\eta$-pairing state is called a type-II tachyonic superconductor
(see the phase diagram in Fig. 1(b) in the main text).

###  I.5 E. Longitudinal component

The longitudinal component of the electromagnetic response function is defined
by $K^{\parallel}(\boldsymbol{q},\omega):=K^{zz}(\boldsymbol{q},\omega)$,
where we take $\boldsymbol{q}\parallel\boldsymbol{e}_{z}$. The $zz$ component
of the current-current correlation function (S14) reads

$\displaystyle\sum_{j}e^{i\boldsymbol{q}\cdot\boldsymbol{R}_{j}}\langle\psi_{N}|J^{z}(\boldsymbol{R}_{j},t)J^{z}(0,0)|\psi_{N}\rangle$
$\displaystyle=4e^{2}t_{h}^{2}C_{M,N}(\langle\boldsymbol{Q}-\boldsymbol{q},+\boldsymbol{e}_{z}|-e^{-iq}\langle\boldsymbol{Q}-\boldsymbol{q},-\boldsymbol{e}_{z}|)$
$\displaystyle\quad\times
e^{-iH(\boldsymbol{Q}-\boldsymbol{q})t}(|\boldsymbol{Q}-\boldsymbol{q},+\boldsymbol{e}_{z}\rangle-e^{iq}|\boldsymbol{Q}-\boldsymbol{q},-\boldsymbol{e}_{z}\rangle),$
(S59)

where $H(\boldsymbol{Q}-\boldsymbol{q})$ is given in Eq. (S25). During the
time evolution, the relative coordinate of two particles changes only in the
$z$ direction. Therefore, what we need to solve is essentially a one-
dimensional two-particle problem, which can be solved exactly in the spirit of
the Bethe ansatz Essler  _et al._ (2005). Here we do not go into details of
analytical solutions, since we can easily diagonalize the Hamiltonian (S25)
numerically for a large system size.

Figure S5: Longitudinal component $K^{\parallel}(\boldsymbol{q},\omega)$ of
the electromagnetic response function for the $\eta$-pairing states in the
Hubbard model in units of $4e^{2}t_{h}C_{M,N}$. (a), (b): Real (a) and
imaginary (b) parts of $K^{\parallel}(\boldsymbol{q},\omega)$ for $U/t_{h}=2$.
(c), (d): Real (c) and imaginary (d) parts of
$K^{\parallel}(\boldsymbol{q},\omega)$ for $U/t_{h}=6$.

In Fig. S5, we plot $K^{\parallel}(\boldsymbol{q},\omega)$ for $U/t_{h}=2$ and
$6$. Compared with the transverse component (Fig. S1), there appear sideband
structures which are shifted by $U$ from the original bands in the
longitudinal component due to the effect of the interaction. Otherwise, both
of them have similar spectral features. In the low-frequency limit, the
longitudinal component vanishes,

$\displaystyle\lim_{\omega\to 0}K^{\parallel}(\boldsymbol{q},\omega)$
$\displaystyle=0,$ (S60)

as required by charge conservation (see Sec. II). In the low-momentum limit,
the longitudinal component agrees with the transverse one,

$\displaystyle\lim_{\boldsymbol{q}\to 0}K^{\parallel}(\boldsymbol{q},\omega)$
$\displaystyle=\lim_{\boldsymbol{q}\to 0}K^{\perp}(\boldsymbol{q},\omega)$
$\displaystyle=8e^{2}t_{h}^{2}C_{M,N}\left(\frac{1}{\omega+U}-\frac{1}{\omega-U}\right),$
(S61)

since the hopping in the $z$ direction is suppressed in the limit of
$\boldsymbol{q}\to 0$ as can be seen from Eq. (S25).

##  II II. Charge conservation

In the Hubbard model [Eq. (1) in the main text], electric charge is conserved
due to the charge $U(1)$ symmetry. This imposes a nontrivial constraint on the
electromagnetic response function Schrieffer (1983). To see this, we introduce
the four-vector form of the electromagnetic response function defined by

$\displaystyle K^{\mu\nu}(\boldsymbol{q},\omega)$
$\displaystyle=-i\theta(t)\langle\psi_{N}|[J^{\mu}(\boldsymbol{R}_{j},t),J^{\nu}(0,0)]|\psi_{N}\rangle$
(S62)

$(\mu,\nu=0,x,y,z)$, where
$J^{\mu}(\boldsymbol{R}_{j})=(c\rho(\boldsymbol{R}_{j}),\boldsymbol{J}(\boldsymbol{R}_{j}))$
is the four-vector current, and

$\displaystyle\rho(\boldsymbol{R}_{j})$
$\displaystyle=e\sum_{\sigma}c_{j\sigma}^{\dagger}c_{j\sigma}$ (S63)

is the local density operator. In the following, we use the metric convention
$\eta^{\mu\nu}={\rm diag}(-,+,+,+)$. The linear response in the four-vector
form reads
$J^{\mu}(\boldsymbol{q},\omega)=-K^{\mu\nu}(\boldsymbol{q},\omega)A_{\nu}(\boldsymbol{q},\omega)$,
where $A_{\nu}=\eta_{\nu\lambda}A^{\lambda}=(-\frac{\phi}{c},\boldsymbol{A})$
and $\phi$ is the scalar potential.

###  II.1 A. Charge response function

The charge response function for the $\eta$-pairing state is given by

$\displaystyle K^{00}(\boldsymbol{R}_{j},t)$
$\displaystyle=-ic^{2}\theta(t)\langle\psi_{N}|[\rho(\boldsymbol{R}_{j},t),\rho(0,0)]|\psi_{N}\rangle.$
(S64)

As before, the density-density correlation function for $N$ particles can be
reduced to the two-particle correlation function. To see this, we define local
$\eta$ operators,

$\displaystyle\eta^{+}(\boldsymbol{R}_{j})$
$\displaystyle:=e^{i\boldsymbol{Q}\cdot\boldsymbol{R}_{j}}c_{j\uparrow}^{\dagger}c_{j\downarrow}^{\dagger},$
(S65) $\displaystyle\eta^{-}(\boldsymbol{R}_{j})$
$\displaystyle:=e^{-i\boldsymbol{Q}\cdot\boldsymbol{R}_{j}}c_{j\downarrow}c_{j\uparrow}.$
(S66)

They satisfy the following commutation relations:

$\displaystyle[\rho(\boldsymbol{R}_{j}),\eta^{\pm}]$ $\displaystyle=\pm
2\eta^{\pm}(\boldsymbol{R}_{j}),$ (S67)
$\displaystyle[\eta^{\pm}(\boldsymbol{R}_{j}),\eta^{\pm}]$ $\displaystyle=0,$
(S68) $\displaystyle[\eta^{\pm}(\boldsymbol{R}_{j}),\eta^{\mp}]$
$\displaystyle=\pm(\rho(\boldsymbol{R}_{j})-1).$ (S69)

Applying the above relations iteratively, we can reduce the $N$-particle
density-density correlation function to

$\displaystyle\langle\psi_{N}|\rho(\boldsymbol{R}_{j},t)\rho(0,0)|\psi_{N}\rangle$
$\displaystyle=4e^{2}C_{M,N}\langle
0|\eta^{-}(\boldsymbol{R}_{j},t)\eta^{+}(0,0)|0\rangle+4e^{2}\frac{\frac{N}{2}(\frac{N}{2}-1)}{M(M-1)}.$
(S70)

From this result, we can evaluate the charge response function as

$\displaystyle K^{00}(\boldsymbol{q},t)$
$\displaystyle=-4i\theta(t)c^{2}e^{2}t_{h}^{2}C_{M,N}$
$\displaystyle\quad\times[\langle\boldsymbol{Q}-\boldsymbol{q},\boldsymbol{r}=0|e^{-iH(\boldsymbol{Q}-\boldsymbol{q})t}|\boldsymbol{Q}-\boldsymbol{q},\boldsymbol{r}=0\rangle$
$\displaystyle\quad-\langle\boldsymbol{Q}+\boldsymbol{q},\boldsymbol{r}=0|e^{iH(\boldsymbol{Q}+\boldsymbol{q})t}|\boldsymbol{Q}+\boldsymbol{q},\boldsymbol{r}=0\rangle],$
(S71)

where

$\displaystyle H(\boldsymbol{Q}\pm\boldsymbol{q})$
$\displaystyle=-t_{h}((1-e^{\pm iq})\Delta_{z}^{+}+(1-e^{\mp
iq})\Delta_{z}^{-})$ $\displaystyle\quad+U(\delta_{\boldsymbol{r},0}-1).$
(S72)

Thus, the problem reduces to solving the two-particle dynamics in the one-
dimensional Hubbard model, which can be diagonalized numerically or
analytically with the Bethe ansatz method.

Figure S6: Charge response function $K^{00}(\boldsymbol{q},\omega)$ for the
$\eta$-pairing states in the Hubbard model in units of
$4c^{2}e^{2}t_{h}C_{M,N}$. (a), (b): Real (a) and imaginary (b) parts of
$K^{00}(\boldsymbol{q},\omega)$ for $U/t_{h}=2$. (c), (d): Real (c) and
imaginary (d) parts of $K^{00}(\boldsymbol{q},\omega)$ for $U/t_{h}=6$.

In Fig. S6, we plot the charge response function for $U/t_{h}=2$ and $6$. We
will see in the next subsection that $K^{00}(\boldsymbol{q},\omega)$ is
related to $K^{\parallel}(\boldsymbol{q},\omega)$ due to the symmetry
constraint. Compared with $K^{\parallel}(\boldsymbol{q},\omega)$, the
contribution of the low-energy sidebands is enhanced in
$K^{00}(\boldsymbol{q},\omega)$, which can also be understood from the
symmetry constraint [Eq. (S79)]. In the low-momentum limit, the charge
response function vanishes,

$\displaystyle\lim_{\boldsymbol{q}\to 0}K^{00}(\boldsymbol{q},\omega)$
$\displaystyle=0,$ (S73)

since $H(\boldsymbol{Q}\pm\boldsymbol{q})=0$ for $q\to 0$ and
$\boldsymbol{r}=0$.

###  II.2 B. Symmetry constraint

Here we see how the symmetry puts a constraint on the electromagnetic response
functions Schrieffer (1983). Our starting point is the continuity equation,

$\displaystyle\frac{d}{dt}\rho(\boldsymbol{R}_{j},t)$
$\displaystyle=-\sum_{\mu=x,y,z}(J^{\mu}(\boldsymbol{R}_{j})-J^{\mu}(\boldsymbol{R}_{j}-\boldsymbol{e}_{\mu})),$
(S74)

which is a direct consequence of the U(1) symmetry of the Hubbard model. After
Fourier transformation, the relation becomes

$\displaystyle-i\omega\rho(\boldsymbol{q},\omega)$
$\displaystyle=-\sum_{\mu=x,y,z}(e^{\frac{i}{2}\boldsymbol{q}\cdot\boldsymbol{e}_{\mu}}-e^{-\frac{i}{2}\boldsymbol{q}\cdot\boldsymbol{e}_{\mu}})J^{\mu}(\boldsymbol{q},\omega).$
(S75)

Taking the expectation value with respect to $|\psi_{N}\rangle$ and
substituting the electromagnetic response function (S62) in Eq. (S75), we
obtain

$\displaystyle-i\frac{\omega}{c}\sum_{\nu=0,x,y,z}K^{0\nu}(\boldsymbol{q},\omega)A_{\nu}(\boldsymbol{q},\omega)$
$\displaystyle=-\sum_{\mu=x,y,z}(e^{\frac{i}{2}\boldsymbol{q}\cdot\boldsymbol{e}_{\mu}}-e^{-\frac{i}{2}\boldsymbol{q}\cdot\boldsymbol{e}_{\mu}})\sum_{\nu=0,x,y,z}K^{\mu\nu}(\boldsymbol{q},\omega)A_{\nu}(\boldsymbol{q},\omega).$
(S76)

Since the above relation must hold for arbitrary
$A_{\nu}(\boldsymbol{q},\omega)$, we conclude that

$\displaystyle
i\frac{\omega}{c}K^{0\nu}(\boldsymbol{q},\omega)-\sum_{\mu=x,y,z}(e^{\frac{i}{2}\boldsymbol{q}\cdot\boldsymbol{e}_{\mu}}-e^{-\frac{i}{2}\boldsymbol{q}\cdot\boldsymbol{e}_{\mu}})K^{\mu\nu}(\boldsymbol{q},\omega)$
$\displaystyle=0.$ (S77)

We assume $\boldsymbol{q}\parallel\boldsymbol{e}_{z}$ without loss of
generality. Then the relation becomes

$\displaystyle K^{0\nu}(\boldsymbol{q},\omega)$
$\displaystyle=\frac{2c\sin\frac{q}{2}}{\omega}K^{z\nu}(\boldsymbol{q},\omega).$
(S78)

Due to Onsager’s reciprocity relation, we also have
$K^{\mu\nu}(\boldsymbol{q},\omega)=K^{\nu\mu}(\boldsymbol{q},\omega)$.
Therefore, $K^{00}(\boldsymbol{q},\omega)$ and
$K^{\parallel}(\boldsymbol{q},\omega)$ are related to each other through

$\displaystyle K^{00}(\boldsymbol{q},\omega)$
$\displaystyle=\frac{4c^{2}\sin^{2}\frac{q}{2}}{\omega^{2}}K^{\parallel}(\boldsymbol{q},\omega).$
(S79)

To summarize, all the components of the electromagnetic response function
$K^{\mu\nu}(\boldsymbol{q},\omega)$ can be expressed in terms of
$K^{\perp}(\boldsymbol{q},\omega)$ and $K^{\parallel}(\boldsymbol{q},\omega)$
as

$\displaystyle K^{\mu\nu}(\boldsymbol{q},\omega)=$
$\displaystyle\begin{pmatrix}\frac{4c^{2}\sin^{2}\frac{q}{2}}{\omega^{2}}K^{\parallel}(\boldsymbol{q},\omega)&0&0&\frac{2c\sin\frac{q}{2}}{\omega}K^{\parallel}(\boldsymbol{q},\omega)\\\
0&K^{\perp}(\boldsymbol{q},\omega)&0&0\\\
0&0&K^{\perp}(\boldsymbol{q},\omega)&0\\\
\frac{2c\sin\frac{q}{2}}{\omega}K^{\parallel}(\boldsymbol{q},\omega)&0&0&K^{\parallel}(\boldsymbol{q},\omega)\end{pmatrix}.$
(S80)

One can check that the symmetry constraint (S78) is consistent with gauge
invariance in the Hubbard model.

##  III III. Dynamical electromagnetic fields

In this section, we give a detailed description of dynamical electromagnetic
fields coupled to the $\eta$-pairing states. We start with the Maxwell
equations,

$\displaystyle\nabla\cdot\boldsymbol{E}=\frac{1}{\varepsilon_{0}}\rho,$ (S81)
$\displaystyle\nabla\times\boldsymbol{E}+\frac{\partial}{\partial
t}\boldsymbol{B}=0,$ (S82) $\displaystyle\nabla\cdot\boldsymbol{B}=0,$ (S83)
$\displaystyle\nabla\times\boldsymbol{B}-\frac{1}{c^{2}}\frac{\partial}{\partial
t}\boldsymbol{E}=\mu_{0}\boldsymbol{j},$ (S84)

where $\boldsymbol{E}$ and $\boldsymbol{B}$ are electric and magnetic fields,
and $\rho$ and $\boldsymbol{j}$ are the charge density and current,
respectively. As usual, we introduce the scalar potential $\phi$ and the
vector potential $\boldsymbol{A}$ through

$\displaystyle\boldsymbol{E}$
$\displaystyle=-\nabla\phi-\frac{\partial}{\partial t}\boldsymbol{A},$ (S85)
$\displaystyle\boldsymbol{B}$ $\displaystyle=\nabla\times\boldsymbol{A}.$
(S86)

In the following, we adopt the Lorenz gauge:

$\displaystyle\frac{1}{c^{2}}\frac{\partial}{\partial
t}\phi+\nabla\cdot\boldsymbol{A}$ $\displaystyle=0.$ (S87)

Then, the equations for $\phi$ and $\boldsymbol{A}$ become

$\displaystyle\frac{1}{c^{2}}\frac{\partial^{2}}{\partial
t^{2}}\phi-\nabla^{2}\phi$ $\displaystyle=\frac{1}{\varepsilon_{0}}\rho,$
(S88) $\displaystyle\frac{1}{c^{2}}\frac{\partial^{2}}{\partial
t^{2}}\boldsymbol{A}-\nabla^{2}\boldsymbol{A}$
$\displaystyle=\mu_{0}\boldsymbol{j}.$ (S89)

To solve these equations, we assume plane-wave solutions,

$\displaystyle\phi(\boldsymbol{r},t)$ $\displaystyle=\phi_{0}e^{-i\omega
t+i\boldsymbol{q}\cdot\boldsymbol{r}},$ (S90)
$\displaystyle\boldsymbol{A}(\boldsymbol{r},t)$
$\displaystyle=\boldsymbol{A}_{0}e^{-i\omega
t+i\boldsymbol{q}\cdot\boldsymbol{r}},$ (S91)

with frequency $\omega$ and momentum $\boldsymbol{q}$. Without loss of
generality, we assume $\boldsymbol{q}\parallel\boldsymbol{e}_{z}$. Combining
the linear-response relation
$j^{\mu}(\boldsymbol{q},\omega)=-K^{\mu\nu}(\boldsymbol{q},\omega)A_{\nu}(\boldsymbol{q},\omega)$
and the gauge condition (S87), we obtain

$\displaystyle-\frac{\omega^{2}}{c^{2}}\phi+q^{2}\phi$
$\displaystyle=\frac{1}{\varepsilon_{0}}K^{00}(q,\omega)\frac{\phi}{c}-\frac{1}{\varepsilon_{0}}K^{0z}(q,\omega)A_{z},$
(S92) $\displaystyle-\frac{\omega^{2}}{c^{2}}A^{x}+q^{2}A^{x}$
$\displaystyle=-\mu_{0}K^{xx}(q,\omega)A_{x},$ (S93)
$\displaystyle-\frac{\omega^{2}}{c^{2}}A^{y}+q^{2}A^{y}$
$\displaystyle=-\mu_{0}K^{yy}(q,\omega)A_{y},$ (S94)
$\displaystyle-\frac{\omega^{2}}{c^{2}}A^{z}+q^{2}A^{z}$
$\displaystyle=-\mu_{0}K^{zz}(q,\omega)A_{z}+\mu_{0}K^{z0}(q,\omega)\frac{\phi}{c},$
(S95) $\displaystyle\frac{\omega}{c^{2}}\phi-qA^{z}$ $\displaystyle=0.$ (S96)

Note that we use the metric $\eta^{\mu\nu}={\rm diag}(-,+,+,+)$ to write down
the above equations. One can see that the transverse ($A^{x},A^{y}$) and
longitudinal ($A^{z}$) components are decoupled.

In the low-energy and long-wavelength limit, the above field equation can be
derived from an effective Lagrangian density

$\displaystyle\mathcal{L}_{\rm eff}$
$\displaystyle=-\frac{1}{4\mu_{0}}F_{\mu\nu}F^{\mu\nu}-V_{\rm eff}(A^{\mu})$
(S97)

with an effective potential

$\displaystyle V_{\rm eff}(A^{\mu})$
$\displaystyle=\frac{1}{2\mu_{0}}m^{2}c^{2}A_{\mu}A^{\mu},$ (S98)

where $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ is the field
strength, and $m$ is the effective mass of the electromagnetic field
corresponding to

$\displaystyle m^{2}$ $\displaystyle=\frac{\mu_{0}}{c^{2}}\lim_{q\to
0}\lim_{\omega\to 0}K^{\perp}(q,\omega)=\frac{\mu_{0}}{\pi c^{2}}D_{s}.$ (S99)

Thus, the squared mass of the electromagnetic field is proportional to the
Meissner weight. In ordinary situations ($D_{s}>0$), the electromagnetic field
acquires a positive squared mass due to the Anderson-Higgs mechanism. In the
opposite case ($D_{s}<0$), the electromagnetic field becomes “tachyonic” with
a negative squared mass.

###  III.1 A. Transverse mode

In order for the transverse modes to exist ($A^{x},A^{y}\neq 0$), the
dispersion must satisfy

$\displaystyle q^{2}-\frac{\omega^{2}}{c^{2}}$
$\displaystyle=-\mu_{0}K^{\perp}(q,\omega),$ (S100)

which corresponds to Eq. (6) in the main text. While the condition (S100)
gives a complicated nonlinear relation between $\omega$ and $q$, the situation
becomes simplified at low momentum.

In the case of $U>0$, Eq. (S100) reduces to

$\displaystyle-\frac{\omega^{2}}{c^{2}}$
$\displaystyle=-8\mu_{0}e^{2}t_{h}^{2}C_{M,N}\left(\frac{1}{\omega+U}-\frac{1}{\omega-U}\right)$
(S101)

in the limit of $q\to 0$ [see Eq. (S61)]. The solution for $\omega^{2}$ is
given by

$\displaystyle\omega^{2}$
$\displaystyle=\frac{1}{2}\left(U^{2}\pm\sqrt{U^{4}-64\mu_{0}c^{2}e^{2}t_{h}^{2}C_{M,N}U}\right).$
(S102)

In order for $\omega$ to take a real value, the interaction strength $U$ must
satisfy

$\displaystyle U^{3}\geq 64\mu_{0}c^{2}e^{2}t_{h}^{2}C_{M,N},$ (S103)

which is exactly the condition derived in the main text. In the thermodynamic
limit ($M,N\to\infty$ with $N/M$ being fixed), $C_{M,N}$ approaches
$\rho(1-\rho)$, where $\rho:=(N/2)/M$ is the doublon density. In this limit,
the condition (S103) becomes

$\displaystyle\frac{U}{t_{h}}\geq
4\left(\frac{\mu_{0}c^{2}e^{2}\rho(1-\rho)}{t_{h}}\right)^{\frac{1}{3}}.$
(S104)

Figure S7: Parameter space (shaded region) where the transverse
electromagnetic field coupled to the $\eta$-pairing state in the Hubbard model
becomes unstable in the long-wavelength limit ($\boldsymbol{q}\to 0$).

In Fig. S7, we plot the parameter space where the frequency $\omega$ has an
imaginary part. One can see that a wide range of the parameter region shows a
dynamical instability of the electromagnetic field coupled to the
$\eta$-pairing state in the long-wavelength limit. The real part of the
frequency at $q\to 0$ is given by

$\displaystyle{\rm Re}\,\omega$
$\displaystyle=\frac{1}{2}\sqrt{\sqrt{64\mu_{0}c^{2}e^{2}t_{h}^{2}C_{M,N}U}+U^{2}},$
(S105)

which is proportional to $U^{1/4}$ at small $U$. The imaginary part of the
frequency at $q\to 0$ is given by

$\displaystyle{\rm Im}\,\omega$
$\displaystyle=\frac{1}{2}\sqrt{\sqrt{64\mu_{0}c^{2}e^{2}t_{h}^{2}C_{M,N}U}-U^{2}}$
(S106)

for $U^{3}<64\mu_{0}c^{2}e^{2}t_{h}^{2}C_{M,N}$.

More generally, if we take into account arbitrary $\boldsymbol{q}$ modes, we
can prove that a dynamical instability exists for all $U$ and $\rho$. First,
we observe that $K^{\perp}(\boldsymbol{q},\omega)$ in Eq. (S31) takes a real
value if and only if $|\omega+U|>4t_{h}\sin\frac{q}{2}$ and
$|\omega-U|>4t_{h}\sin\frac{q}{2}$. Let us first consider the case of
$0<\omega<U-4t_{h}\sin\frac{q}{2}$. Using the result for
$K^{\perp}(\boldsymbol{q},\omega)$ in Eq. (S31), the dispersion relation
(S100) can be written as

$\displaystyle q^{2}-\frac{\omega^{2}}{c^{2}}$
$\displaystyle=-8\mu_{0}e^{2}t_{h}^{2}C_{M,N}\left\\{\frac{1}{\sqrt{(\omega+U)^{2}-16t_{h}^{2}\sin^{2}\frac{q}{2}}}\right.$
$\displaystyle\quad\left.+\frac{1}{\sqrt{(\omega-U)^{2}-16t_{h}^{2}\sin^{2}\frac{q}{2}}}\right\\},$
(S107)

which is negative definite. Therefore, we have $\omega>cq$. For the other case
of $\omega>U+4t_{h}\sin\frac{q}{2}$, the condition (S100) reads

$\displaystyle q^{2}-\frac{\omega^{2}}{c^{2}}$
$\displaystyle=-8\mu_{0}e^{2}t_{h}^{2}C_{M,N}\left\\{\frac{1}{\sqrt{(\omega+U)^{2}-16t_{h}^{2}\sin^{2}\frac{q}{2}}}\right.$
$\displaystyle\quad\left.-\frac{1}{\sqrt{(\omega-U)^{2}-16t_{h}^{2}\sin^{2}\frac{q}{2}}}\right\\},$
(S108)

which is positive definite since
$(\omega+U)^{2}-16t_{h}^{2}\sin^{2}\frac{q}{2}>(\omega-U)^{2}-16t_{h}^{2}\sin^{2}\frac{q}{2}$.
Therefore, we have $\omega<cq$.

Figure S8: Schematic illustration of the range of $(q,\omega)$ in which the
real solution for Eq. (S100) with $U>0$ is allowed as shown by the shaded
region for $qa\ll 1$. The dashed lines indicate that the boundary is not
included.

In Fig. S8, we plot the range of $(q,\omega)$ in which the real solution for
Eq. (S100) is allowed. Since the allowed region is separated into two disjoint
islands, it is clear from a topological point of view that the real band
dispersion $\omega=\omega(q)$ for all $q$ is not possible. This means that
there must always be a region in $(q,\omega)$ where the solution for Eq.
(S100) becomes complex. Thus, there is a dynamical instability for arbitrary
$U(>0)$ and $\rho$.

In the attractive case ($U<0$), the mode equation in the long-wavelength limit
is similarly given by Eq. (S101). The solution for $\omega^{2}$ is the same as
Eq. (S102). When $U<0$, $\omega^{2}$ is always real. However, there exists a
solution with $\omega^{2}<0$ when one chooses the minus sign in Eq. (S102).
Hence the frequency becomes imaginary for $\boldsymbol{q}\to 0$. The inverse
of the imaginary part of $\omega$ is given by

$\displaystyle\frac{1}{{\rm Im}\,\omega}$
$\displaystyle=\frac{1}{\sqrt{\frac{1}{2}\left(\sqrt{U^{4}-64\mu_{0}c^{2}e^{2}t_{h}^{2}C_{M,N}U}-U^{2}\right)}}.$
(S109)

In the attractive case, the electromagnetic field is dynamically unstable
against the $\boldsymbol{q}=0$ mode. This corresponds to the fact that the
electromagnetic field has a negative squared mass. Physically, the
electromagnetic field with a long wavelength penetrates deeply inside the
$\eta$ pairing state, transferring the kinetic energy of doublons to the
electromagnetic field. The time scale of the growth of this instability is
determined by Eq. (S109).

Combining the arguments for the two cases ($U>0$ and $U<0$), we have
established that the electromagnetic field coupled to the $\eta$-pairing state
is always dynamically unstable for all $U$ and $\rho$.

###  III.2 B. Longitudinal mode

The mode equation for the longitudinal components can be derived from Eqs.
(S92), (S95), and (S96). To simplify the situation, we focus on the low-
momentum region ($qa\ll 1$). In this region, the dispersion is determined by

$\displaystyle q^{2}-\frac{\omega^{2}}{c^{2}}$
$\displaystyle=-\mu_{0}K^{\parallel}(q,\omega)\left(1-\frac{c^{2}q^{2}}{\omega^{2}}\right),$
(S110)

where we have used the relation (S80). In the limit of $q\to 0$, the mode
equation becomes

$\displaystyle-\frac{\omega^{2}}{c^{2}}$ $\displaystyle=-\mu_{0}\lim_{q\to
0}K^{\parallel}(q,\omega).$ (S111)

If we recall the relation (S61), the dispersion (S111) is the same as that of
the transverse mode (S100). Therefore, the longitudinal mode has the same
dynamical instability as the transverse one at low momentum.

##  IV IV. Magnetic properties of tachyonic superconductors

In this section, we describe static magnetic properties of tachyonic
superconductors realized as the $\eta$-pairing states in the Hubbard model
with $U<0$. As shown in the main text and in the preceding section, the
tachyonic superconductors are dynamically unstable. Here we focus on the
response of the $\eta$-pairing states against static magnetic fields within
the linear-response regime, and do not consider their decay dynamics.

The static magnetic field $\boldsymbol{B}$ obeys the following Maxwell
equations:

$\displaystyle\nabla\cdot\boldsymbol{B}$ $\displaystyle=0,$ (S112)
$\displaystyle\nabla\times\boldsymbol{B}$
$\displaystyle=\mu_{0}\boldsymbol{j}.$ (S113)

We introduce a static vector potential $\boldsymbol{A}$ as
$\boldsymbol{B}=\nabla\times\boldsymbol{A}$. If we take the Coulomb gauge
($\nabla\cdot\boldsymbol{A}=0$), the equation for $\boldsymbol{A}$ becomes

$\displaystyle\nabla^{2}\boldsymbol{A}$
$\displaystyle=-\mu_{0}\boldsymbol{j}.$ (S114)

To solve the equation, we assume a plane-wave form,

$\displaystyle\boldsymbol{A}(\boldsymbol{r})$
$\displaystyle=\boldsymbol{A}_{0}e^{i\boldsymbol{q}\cdot\boldsymbol{r}},$
(S115)

with amplitude $\boldsymbol{A}_{0}$ and wave number $\boldsymbol{q}$. Without
loss of generality, we choose $\boldsymbol{q}=q\boldsymbol{e}_{z}$. We apply
the linear-response theory to the $\eta$-pairing states to obtain the
equations for $\boldsymbol{A}_{0}$:

$\displaystyle\boldsymbol{q}\cdot\boldsymbol{A}_{0}$ $\displaystyle=0,$ (S116)
$\displaystyle-q^{2}A_{0}^{\mu}$
$\displaystyle=-\mu_{0}j^{\mu}(\boldsymbol{q})$
$\displaystyle=\mu_{0}\sum_{\nu}K^{\mu\nu}(\boldsymbol{q},\omega=0)A_{0\nu}\quad(\mu,\nu=x,y,z).$
(S117)

From Eq. (S116), we find $A_{0}^{z}=0$. In order for the solution
$A_{0}^{\mu}\neq 0$ to exist, $q$ must satisfy a nonlinear equation,

$\displaystyle q^{2}$
$\displaystyle=-\mu_{0}K^{\perp}(\boldsymbol{q},\omega=0)$
$\displaystyle=-16\mu_{0}e^{2}t_{h}^{2}C_{M,N}$
$\displaystyle\quad\times\begin{cases}\frac{{\rm
sgn}(U)}{\sqrt{U^{2}-16t_{h}^{2}\sin^{2}\frac{q}{2}}}&\mbox{for}\;|U|>4t_{h}\sin\frac{q}{2};\\\
0&\mbox{for}\;|U|<4t_{h}\sin\frac{q}{2},\end{cases}$ (S118)

where we have used the result (S31).

One can immediately see that a real solution for Eq. (S118) does not exist for
$U>0$. This is nothing but the Meissner effect; that is, a magnetic field
cannot propagate freely into superconductors. In fact, the magnetic field
decays exponentially in space with the penetration depth

$\displaystyle\lambda$
$\displaystyle=\sqrt{\frac{U}{16\mu_{0}e^{2}t_{h}^{2}C_{M,N}}}.$ (S119)

On the other hand, when $U<0$ a real solution is possible as is clear from the
graphical illustration of Eq. (S118) in Fig. S9. Thus, a magnetic field can
penetrate into tachyonic superconductors without decay.

Figure S9: Graphical illustration of Eq. (S118) for (a) $-4t_{h}<U<0$ and (b) $U<-4t_{h}$. The left-hand side of Eq. (S118) is shown by the red curve, while the right-hand side is shown by the blue curve. | # of real solutions for Eq. (S118)
---|---
$U>0$ | 0
$-4t_{h}<U<0$ | $\infty$
$U<-4t_{h}$ | 2
Table 1: List of the number of real solutions for Eq. (S118) with
$16\mu_{0}e^{2}t_{h}a/\hbar^{2}\ll 1$.

From Fig. S9, we can see that the number of real solutions for Eq. (S118)
changes at the boundary of $U=-4t_{h}$. In Table. 1, we list the number of
real solutions for Eq. (S118). This result suggests that the number of modes
of magnetic fields that can propagate inside tachyonic superconductors for
$-4t_{h}<U<0$ is different from that for $U<-4t_{h}$. Following the main text,
we call the former a type-I tachyonic superconductor, and the latter a type-II
tachyonic superconductor.

In the case of the type-I tachyonic superconductor, we further classify the
solutions into two types according to the number of real solutions in the
range of $-\pi<q<\pi$. When $U$ belongs to the range $(-4t_{h}<)$
$U_{\ast}<U<0$ with a certain boundary $U_{\ast}$, real solutions for Eq.
(S118) do not exist in $-\pi<q<\pi$. Since $16\mu_{0}e^{2}t_{h}a/\hbar^{2}\ll
1$ for ordinary materials ($t_{h}\sim$ 1[eV] and $a\sim$ 1[Å]), the threshold
$U_{\ast}$ is approximately given as

$\displaystyle\frac{U_{\ast}}{t_{h}}$
$\displaystyle\approx-2\sqrt{3}(4\mu_{0}e^{2}t_{h}C_{M,N})^{\frac{1}{3}}.$
(S120)

For $-4t_{h}<U<U_{\ast}$, there exist four real solutions in $-\pi<q<\pi$,
which are denoted by $\pm q_{1}$ and $\pm q_{2}$ with $0<q_{1}<q_{2}$. Using
$16\mu_{0}e^{2}t_{h}a/\hbar^{2}\ll 1$ again, we can approximately evaluate
$q_{1}$ and $q_{2}$ as

$\displaystyle q_{1}$
$\displaystyle\approx\sqrt{\frac{16\mu_{0}e^{2}t_{h}^{2}C_{M,N}}{-U}},$ (S121)
$\displaystyle q_{2}$ $\displaystyle\approx
2\sin^{-1}\left(\frac{-U}{4t_{h}}\right).$ (S122)

Note that $q_{1}$ corresponds to the inverse of the analytically continued
London’s penetration depth ($q_{1}=|\lambda|^{-1}$). If we extend the range of
$q$ to $-\infty<q<\infty$, there are infinitely many real solutions. They are
approximately given as $\pm q_{1}$ and $\pm q_{2}+2n\pi$ ($n\in\mathbb{Z}$).

Figure S10: Log plot of $q_{1}^{-1}$, $q_{2}^{-1}$, and $q_{3}^{-1}$ for
$t_{h}=$ 1[eV], $a=$ 1[Å], and $\rho=0.5$ as a function of $-U/t_{h}$.

In the case of the type-II tachyonic superconductor, there are only two real
solutions $\pm q_{1}$ approximately given by (S121). Instead, there emerge
infinitely many complex solutions approximately given by $\pm iq_{3}+2n\pi$
($n\in\mathbb{Z}$) with

$\displaystyle q_{3}$ $\displaystyle\approx
2\cosh^{-1}\left(\frac{-U}{4t_{h}}\right).$ (S123)

Physically, these solutions correspond to a magnetic field localized near the
surface of a tachyonic superconductor. The localization length $q_{3}^{-1}$
diverges at $U=U_{c}^{\prime}=-4t_{h}$ as
$q_{3}^{-1}\sim|U-U_{c}^{\prime}|^{-1/2}$.

In Fig. S10, we plot $q_{1}^{-1}$, $q_{2}^{-1}$, and $q_{3}^{-1}$ for typical
parameters. The $q_{1}$ mode produces a long-period magnetic structure with
the period of the order of $|\lambda|\gg a$. On the other hand, the $q_{2}$
mode provides a short-period magnetic structure with the period length of the
order of the lattice constant $a$. In type-I tachyonic superconductors, both
the long- and short-period structures are allowed to exist, whereas in type-II
tachyonic superconductors the short-period magnetic structure is screened, and
it can penetrate only near the surface.

In Fig. S11, we show examples of magnetic fields that can be realized in
type-I [Fig. S11(a)] and type-II (b) tachyonic superconductors. Even when
there is no magnetic field outside of tachyonic superconductors, nonzero
magnetic fields can be trapped statically inside tachyonic superconductors. To
support those magnetic fields, constant electric currents are flowing in the
bulk of tachyonic superconductors. In the type-II tachyonic superconductor,
there is also a surface current to satisfy the boundary condition.

Figure S11: Examples of static magnetic fields $B=B(x)$ that can exist inside
(a) type-I and (b) type-II tachyonic superconductors shown by shaded regions.

Figure S12: Examples of (a) stripe-like and (b) vortex-antivortex-like
magnetic structures in tachyonic superconductors shown in a two-dimensional
plane.

Various $q$ modes can be linearly superposed in several different directions.
In Fig. S12, we show two examples of magnetic structures that can be realized
in tachyonic superconductors. If one only takes a single $q_{1}$ mode, it
gives a stripe-like structure as shown in Fig. S12(a). Here we neglect short-
period structures ($q_{2}$ modes). If one superposes two $q_{1}$ modes in $x$
and $y$ directions, one obtains a square lattice with alternating vortex and
antivortex structures as shown in Fig. S12. This is to be contrasted with
Abrikosov’s triangular lattice of vortices in type-II superconductors. One can
also superimpose three $q_{1}$ modes in three different directions, creating a
three-dimensional magnetic structure (not shown). In this way, various
configurations of magnetic fields can be trapped in tachyonic superconductors.
We note, however, that these structures are not dynamically stable as shown in
the main text and in the preceding section.

##  V V. Spontaneous light emission

In this section, we evaluate the rate of spontaneous light emission for
$\eta$-pairing states in the Hubbard model with $U>0$. In the repulsive case,
doublons can decay spontaneously into pairs of single particles by emitting
light with frequency $\omega=U$. After emitting light, the $\eta$-pairing
state $|\psi_{N}\rangle$ is transformed into

$\displaystyle J^{\mu}(\boldsymbol{q}=0)|\psi_{N}\rangle$
$\displaystyle=\frac{e}{\sqrt{\mathcal{N}_{N}}}N(\eta^{+})^{\frac{N}{2}-1}\sum_{\boldsymbol{k}}v^{\mu}(\boldsymbol{k})c_{\boldsymbol{k}\uparrow}^{\dagger}c_{\boldsymbol{Q}-\boldsymbol{k}\downarrow}^{\dagger}|0\rangle,$
(S124)

where
$v^{\mu}(\boldsymbol{k})=\frac{\partial\varepsilon_{\boldsymbol{k}}}{\partial
k^{\mu}}=2t_{h}\sin k^{\mu}$ is the group velocity
($\varepsilon_{\boldsymbol{k}}=-2t_{h}\sum_{\mu=x,y,z}\cos k^{\mu}$ is the
band dispersion).

Let us define a one-doublon-broken state Yang (1989)

$\displaystyle|\zeta_{N,\boldsymbol{a}}\rangle$
$\displaystyle=\frac{1}{\sqrt{\mathcal{N}_{N,\boldsymbol{a}}}}(\eta^{+})^{\frac{N}{2}-1}\eta_{\boldsymbol{a}}^{+}|0\rangle,$
(S125)

where $\mathcal{N}_{N,\boldsymbol{a}}$ is the normalization constant (such
that
$\langle\zeta_{N,\boldsymbol{a}}|\zeta_{N,\boldsymbol{a}^{\prime}}\rangle=\delta_{\boldsymbol{a},\boldsymbol{a}^{\prime}}$),
$\boldsymbol{a}$ represents a lattice-site coordinate, and

$\displaystyle\eta_{\boldsymbol{a}}^{+}$
$\displaystyle=\sum_{\boldsymbol{k}}e^{-i\boldsymbol{k}\cdot\boldsymbol{a}}c_{\boldsymbol{k}\uparrow}^{\dagger}c_{\boldsymbol{Q}-\boldsymbol{k}\downarrow}^{\dagger}.$
(S126)

The state $|\zeta_{N,\boldsymbol{a}}\rangle$ consists of $\frac{N}{2}-1$
doublons with momentum $\boldsymbol{Q}$ and two unpaired particles with the
lattice spacing $\boldsymbol{a}$. One can show that
$|\zeta_{N,\boldsymbol{a}}\rangle$ with $\boldsymbol{a}\neq 0$ is an exact
eigenstate of the Hubbard model [Eq. (1) in the main text] with the
eigenenergy $-U$. At $\boldsymbol{a}=0$, we have
$|\zeta_{N,\boldsymbol{a}}\rangle=|\psi_{N}\rangle$. Using
$|\zeta_{N,\boldsymbol{a}}\rangle$, we can write the one-photon emitted state
(S124) as

$\displaystyle J^{\mu}(\boldsymbol{q}=0)|\psi_{N}\rangle$
$\displaystyle=e\sqrt{\frac{\mathcal{N}_{N,\boldsymbol{a}}}{\mathcal{N}_{N}}}\frac{N}{M}\sum_{\boldsymbol{k}}v^{\mu}(\boldsymbol{k})\sum_{\boldsymbol{a}\neq
0}e^{i\boldsymbol{k}\cdot\boldsymbol{a}}|\zeta_{N,\boldsymbol{a}}\rangle.$
(S127)

Therefore, all the states that are accessible by one-photon emission are
covered by the eigenstates $|\zeta_{N,\boldsymbol{a}}\rangle$
($\boldsymbol{a}\neq 0$).

The rate of spontaneous emission is given by Einstein’s A coefficient Loudon
(2000):

$\displaystyle\Gamma$
$\displaystyle=\frac{\omega^{3}}{3\pi\varepsilon_{0}c^{3}\hbar}\sum_{\mu}|\langle
e|P^{\mu}|\sigma\rangle|^{2},$ (S128)

where $P^{\mu}$ is the polarization operator, $|\sigma\rangle$ is an initial
state, and $|e\rangle$ is a one-phonon emitted state. In the present case, we
take $|\sigma\rangle=|\psi_{N}\rangle$ and
$|e\rangle=|\zeta_{N,\boldsymbol{a}}\rangle$. Since
$\frac{d}{dt}P^{\mu}=J^{\mu}$, the rate $\Gamma$ is rewritten as

$\displaystyle\Gamma$
$\displaystyle=\frac{\omega}{3\pi\varepsilon_{0}c^{3}\hbar}\sum_{\mu}\sum_{\boldsymbol{a}}|\langle\zeta_{N,\boldsymbol{a}}|J^{\mu}|\psi_{N}\rangle|^{2}$
$\displaystyle=\frac{e^{2}\omega}{3\pi\varepsilon_{0}c^{3}\hbar}\sum_{\mu}\sum_{\boldsymbol{a}}\frac{\mathcal{N}_{N,\boldsymbol{a}}}{\mathcal{N}_{N}}\frac{N^{2}}{M^{2}}\left|\sum_{\boldsymbol{k}}v^{\mu}(\boldsymbol{k})e^{i\boldsymbol{k}\cdot\boldsymbol{a}}\right|^{2}.$
(S129)

Using Eq. (S13) and

$\displaystyle\mathcal{N}_{N,\boldsymbol{a}}$
$\displaystyle=\frac{M(M-2)!(\frac{N}{2}-1)!}{(M-\frac{N}{2}-1)!}\quad(\boldsymbol{a}\neq
0),$ (S130)

we obtain

$\displaystyle\Gamma$ $\displaystyle=\frac{8e^{2}\omega
t_{h}^{2}a^{2}dC_{M,N}}{3\pi\varepsilon_{0}c^{3}\hbar^{3}}M,$ (S131)

where $d$ is the dimension of the system. One can see that the rate $\Gamma$
is proportional to the system size, which is natural because the doublon decay
can take place at any lattice site with equal probability.

For ordinary three-dimensional materials, we substitute $t_{h}=$ 1[eV],
$\omega=U=$ 1[eV], $a=$ 1[Å], $\rho=0.5$, and $d=3$ in Eq. (S131), obtaining

$\displaystyle\frac{\Gamma}{M}$ $\displaystyle=2.3\times 10^{7}\;[{\rm
s}^{-1}].$ (S132)