# Multi-Triplon Excitations of Hubbard Ladders with Site-Dependent Potentials

Nobuya Maeshima1,2 and Ken-ichi Hino2,1<EMAIL_ADDRESS>for
Computational Sciences1Center for Computational Sciences University of Tsukuba
University of Tsukuba Tsukuba Tsukuba Ibaraki 305-8577 Ibaraki 305-8577 Japan
2Institute of Pure and Applied Sciences Japan
2Institute of Pure and Applied Sciences University of Tsukuba University of
Tsukuba Tsukuba Tsukuba Ibaraki 305-8573 Ibaraki 305-8573 Japan Japan

###### Abstract

We study low-lying spin-singlet excitations of two-leg Hubbard ladders with
site-dependent potentials. Using general formulas of the charge
disproportionation induced by the site-dependent potentials, we derive the
contributions of spin degrees of freedom to the spectral functions such as the
dynamical charge structure factor $N(\bm{k},\omega)$ and the optical
conductivity $\sigma_{\gamma}(\omega)$ along the $\gamma(=x,y)$-direction of
the two-leg ladders. Numerical results obtained by the Lanczos diagonalization
method have clarified that the multi-triplon singlet states, including two-
and three-triplon excitations, can be detected by observing these quantities.
Furthermore, we have found that ladders with random potentials also have non-
negligible contributions of these excitations to $\sigma_{x}(\omega)$.

##  1 Introduction

Quantum spin systems have long been regarded an ideal arena for exploring
quantum many-body phenomena. In particular, $S=1/2$ two-leg spin ladders are
notable for their tractability and typical low-dimensional features of quantum
systems, such as disordered ground states, finite spin gaps, and non-trivial
excitations [1, 2]. The elementary excitations in $S=1/2$ two-leg ladders are
quasiparticles called ”triplons,” each carrying a spin of $S=1$ [3]. The
triplon state is regarded as a quasiparticle in a wide class of quantum spin
systems with finite spin gaps, and its dynamical properties have been examined
in various materials [4, 5, 6, 7, 8]. Also, in two-leg ladders, experimental
studies using inelastic neutron scattering (INS) have successfully identified
the finite energy gap and the single-triplon dispersion [9, 10].

Recent experimental achievements strongly indicate that two-leg spin ladders
are excellent platforms for studying two-triplon states; high-resolution
analysis of INS experiments has confirmed the existence of a two-triplon
continuum and bound states in the spin-triplet ($S=1$) sector of the excited
states [11, 12]. Besides these two-triplon states with $S=1$, which are
detectable with INS, two-triplon singlet ($S=0$) excitations have also been
observed using Raman spectroscopy [13, 14] and phonon-assisted optical
absorption [15]. In addition, resonant inelastic x-ray scattering (RIXS) has
been employed to explore two-triplon $S=0$ excitations [16, 17]. These
experimental results have been compared with theoretical studies on two-
triplon states [18, 19, 20, 21, 22, 23].

Turning our attention to recent studies of other quasiparticles in quantum
magnets, we can see various experimental challenges to observe excited states
having three or more numbers of quasiparticles, contributing to the spectral
properties of these systems in high-energy regions. The majority of these
studies are devoted to bound states formed by strong inter-quasiparticle
interactions and are investigated using INS [24, 25], THz spectroscopy [26,
27], and infrared magnetospectroscopy [28]. Additionally, multi-quasiparticle
continuum states composed of weakly interacting (or free) quasiparticles have
also been examined in relation to the high-energy spectral components of INS
or RIXS in 1D chains [29, 30], or those of RIXS in a kagome antiferromagnet
[31]. As for the triplon states, in contrast, it has been only proposed
theoretically that three-triplon states of the ladder can appear in the RIXS
intensity [32]. Thus, methods to detect multi-triplon states are still
lacking.

We present an alternative framework to investigate multi-triplon excitations
of the two-leg ladders in the $S=0$ sector. Our approach is based on an
analytical treatment for the one-dimensional (1D) ionic Hubbard model, a
single-band Hubbard chain with an alternating potential [33], where the charge
disproportionation at each site is represented in terms of the spin operators
of the corresponding Heisenberg model. In this work, the analytical method is
extended to Hubbard models on general lattices and applied to the two-leg
Hubbard ladders with site-dependent potentials. We have analytically proved
that the spin-degrees of freedom of the two-leg ladders can have finite
contributions to the dynamical charge structure factor $N(\bm{k},\omega)$ and
the optical conductivity $\sigma_{\gamma}(\omega)$ along the
$\gamma(=x,y)$-direction of the ladders.

Then, we numerically demonstrate that these spectral functions probe the two
and three-triplon excitations of two-leg ladders with periodic potentials and
with strong rung coupling, where spin dimers called “rungs” weakly couple
along one direction. The observed excitations can be selected by tuning the
periodicity of the site-dependent potential or the wave number perpendicular
to the leg. Obtained results of $N(\bm{k},\omega)$ and
$\sigma_{\gamma}(\omega)$ calculated with the Lanczos diagonalization
elucidate that their spectral peaks lie in the two- and three-triplon
continuum and the two-triplon singlet bound state. Also, analytical results
with bond operator method [34, 35, 36] are found to be consistent with our
calculated results. Furthermore, we investigate the optical conductivity
$\sigma_{x}(\omega)$ with random potential and show that observed spectral
signatures result from these multi-triplon excitations.

This paper is organized as follows. In Sec. 2, our theoretical framework for
the single chain is extended to general lattices and then applied to two-leg
ladders. In Sec. 3, our numerical and analytical results are presented, and
Section 4 is devoted to the conclusion.

##  2 Theory

###  2.1 General Hubbard models

Here, we explain the method, which is an extension of our previous framework
for the 1D ionic Hubbard model [33]. We employ the single-band Hubbard model
at half-filling with a site-dependent potential. The Hamiltonian is given by

$\displaystyle{\cal
H}=-\sum_{jl\tau}t_{jl}c^{\dagger}_{j\tau}c_{l\tau}+U\sum_{l}n_{l\uparrow}n_{l\downarrow}+\sum_{l}e_{l}(n_{i}-1),$
(1)

where $c^{\dagger}_{j\tau}$ ($c_{j\tau}$) is the creation (annihilation)
operator of an electron with spin $\tau$ at site $j$,
$n_{j\tau}=c^{\dagger}_{j\tau}c_{j\tau}$, and
$n_{j}=n_{j\uparrow}+n_{j\downarrow}$. The parameter $U$ gives the on-site
Coulomb interaction, and $t_{jl}$ represents the transfer integral between
sites $j$ and $l$. In addition, $e_{j}$ denotes the site-dependent potential
at site $j$. Here we set $t_{jj}=0$ and $e_{j}=\epsilon_{j}+\phi_{j}$ at site
$j$, where $\epsilon_{j}$ is an intrinsic site-dependent potential and
$\phi_{j}$ is a virtually introduced scalar field. For the large-$U$ case of
this model, we derive the Heisenberg Hamiltonian as [37, 33]

${\cal H}^{s}=\sum_{jl}K_{jl}\left(\bm{S}_{j}\cdot{\bf\it
S}_{l}-\frac{1}{4}\right),$ (2)

with the exchange interaction $K_{jl}=2t_{jl}^{2}/(U-e_{j}+e_{l})$. The charge
disproportionation $\delta n_{l}\equiv n_{l}-1$ from the half-filling for
$\phi_{l}=0$ is evaluated as

$\displaystyle\delta n^{s}_{l}=\left.\frac{\partial{\cal
H}^{s}}{\partial\phi_{l}}\right|_{\phi_{l}=0}=\sum_{j}\eta_{lj}\left(\bm{S}_{l}\cdot\bm{S}_{j}-\frac{1}{4}\right),$
(3)

where $\bm{S}_{j}$ is the spin operator at site $j$, and the coefficient
$\eta_{lj}$ associating $\delta n^{s}_{l}$ with the spin operators is given by

$\displaystyle\eta_{lj}=\frac{8Ut^{2}_{lj}(\epsilon_{l}-\epsilon_{j})}{[U^{2}-(\epsilon_{l}-\epsilon_{j})^{2}]^{2}}.$
(4)

Here, the superscript $s$ of quantities means the contribution from the spin
degrees of freedom.

Next, let us consider the dynamical charge structure factor $N(\bm{k},\omega)$
of the Hubbard model (1) defined by

$N(\bm{k},\omega)=\frac{1}{N}\sum_{\alpha\neq
0}|\langle\alpha|n_{\bm{k}}|0\rangle|^{2}\delta(\omega-E_{\alpha}+E_{0}),$ (5)

where $n_{\bm{k}}=\sum_{l}e^{-i{\bm{k}}\cdot\bm{R}_{l}}n_{l}$, ${\bm{k}}$
denotes the wave number, $N$ is the system size, and $\bm{R}_{l}$ represents
the location of site $l$. Its counterpart $N^{s}(\bm{k},\omega)$ of the spin
degrees of freedom is obtained by replacing $n_{\bm{k}}$ in Eq. (5) with

$\delta n^{s}_{\bm{k}}=\sum_{l}e^{-i\bm{k}\cdot\bm{R}_{l}}\delta
n^{s}_{l}=\sum_{l}e^{-i\bm{k}\cdot\bm{R}_{l}}\sum_{j}\eta_{lj}\left(\bm{S}_{l}\cdot\bm{S}_{j}-\frac{1}{4}\right).$
(6)

At the same time, the eigenstate $\ket{\alpha}$ and the eigenenergy
$E_{\alpha}$ of the Hamiltonian (1) are also replaced by their counterparts,
$|\alpha\rangle^{s}$ and $E_{\alpha}^{s}$, of the Heisenberg Hamiltonian (2).

Now, we discuss the optical conductivity of the Hubbard model (1) given by

$\displaystyle\sigma_{\gamma}(\omega)=\frac{\pi}{N\omega}\sum_{\alpha\neq
0}|\langle\alpha|J_{\gamma}|0\rangle|^{2}\delta(\omega-E_{\alpha}+E_{0})$ (7)

and $\sigma^{s}_{\gamma}(\omega)$ in terms of the spin degrees of freedom,
where $J_{\gamma}$ is the $\gamma$ component of the current operator of the
model (1). Following the procedure to derive $N^{s}(\bm{k},\omega)$, we
represent $J_{\gamma}$ with respect to the operators of spin system (2). To
this aim, let us express the electric polarization [37] for $\gamma$ direction
by using $\delta n^{s}_{j}$ of Eq. (3) as

$\displaystyle P^{s}_{\gamma}\equiv\sum_{j}(-e)R_{j}^{\gamma}\delta
n^{s}_{j}=-e\sum_{ij}\eta_{ij}R^{\gamma}_{i}\left({\bf\it S}_{i}\cdot{\bf\it
S}_{j}-\frac{1}{4}\right),$ (8)

where $e$ is the absolute value of the charge of an electron. Then the current
operator is derived as [38]

$\displaystyle J^{s}_{\gamma}$
$\displaystyle\equiv\frac{dP^{s}_{\gamma}}{dt}=-\frac{i}{\hbar}[P^{s}_{\gamma},{\cal
H}^{s}]$
$\displaystyle=-\frac{e}{\hbar}\sum_{ijl}\eta_{ij}\overline{K}_{il}(R^{\gamma}_{i}-R^{\gamma}_{j})\
{\bf\it S}_{i}\cdot({\bf\it S}_{j}\times{\bf\it S}_{l}),$ (9)

where

$\displaystyle\overline{K}_{il}=\frac{4U(t_{il})^{2}}{U^{2}-(\epsilon_{i}-\epsilon_{l})^{2}},$
(10)

and the optical conductivity is defined as

$\displaystyle\sigma^{s}_{\gamma}(\omega)=\frac{\pi}{N\omega}\sum_{\alpha\neq
0}|\,^{s\,}\\!\langle\alpha|J^{s}_{\gamma}|0\rangle^{s}|^{2}\delta(\omega-E^{s}_{\alpha}+E^{s}_{0}),$
(11)

As $n^{s}_{\bm{k}}$ and $J^{s}_{\gamma}$ commute with the total spin ${\bf\it
S}_{\rm tot}=\sum_{l}{\bf\it S}_{l}$, the spectral functions
$N^{s}(\bm{k},\omega)$ and $\sigma^{s}_{\gamma}(\omega)$ can detect the
singlet excitation [37, 38]. Also, it is noted that these formulas can be
applied to the 1D ionic Hubbard model by setting

$\displaystyle t_{jj+1}=t,\ \epsilon_{j}=-(-1)^{j}\frac{\Delta}{2},\ {\rm
and}\ R_{j}=aj,$ (12)

with the lattice constant $a$ and the strength of the potential $\Delta$.
Then, it can be confirmed that the same results of Refs. katsura,yokoi,komaki
are obtained.

We here stress the importance of the quantity $\eta_{lj}$ in our framework.
The point is that $\eta_{lj}$ appears both in $N^{s}(\bm{k},\omega)$ and
$\sigma^{s}_{\gamma}(\omega)$ and that $\eta_{lj}$ is proportional to the
energy difference $(\epsilon_{l}-\epsilon_{j})$ between the neighboring sites
$l$ and $j$. Therefore, the finite site-dependence of the potential
$\epsilon_{l}$ is essential for observing these spectral functions.

###  2.2 Two-leg ladders

Figure 1:  A Hubbard ladder with two legs 0 and 1, and the transfer integrals
$t_{x}$ and $t_{y}$. The arrows are the primitive translational vectors
${\bf\it a}_{x}$ and ${\bf\it a}_{y}$ of this system.

Here, let us focus on two-leg ladders with only nearest-neighbor hoppings and
periodic potentials. Some notations are introduced to discuss the ladder
systems in detail. We label a site in a ladder by its location ${\bf\it
R}=(x,y)$ and use the primitive translational vectors ${\bf\it a}_{x}$ and
${\bf\it a}_{y}$. Here, ${\bf\it a}_{x(y)}$ is parallel to the legs (rungs),
as shown in Fig. 1. A leg in a ladder is labeled by an integer $n=0,1$, and
the location of a site in the leg $n$ is given by ${\bf\it R}=(x,na_{y})$,
where $a_{\gamma}=|{\bf\it a}_{\gamma}|$ is the lattice constant along the
$\gamma$-direction. In addition, a site ${\bf\it R}$ in the leg $n$ is
represented by the short form ${\bf\it R}\in\text{leg}\ n$. Then the transfer
integral $t_{{\bf\it R}{\bf\it R}^{\prime}}$ between ${\bf\it R}$ and ${\bf\it
R}^{\prime}$ is given as

$\displaystyle t_{{\bf\it R}{\bf\it R}^{\prime}}=\begin{cases}t_{x}&{\bf\it
R}^{\prime}={\bf\it R}\pm{\bf\it a}_{x}\\\ t_{y}&\quad({\bf\it
R}\in\text{leg}\ 0\ \text{and}\ {\bf\it R}^{\prime}={\bf\it R}+{\bf\it
a}_{y})\\\ &\text{or}\ ({\bf\it R}\in\text{leg}\ 1\ \text{and}\ {\bf\it
R}^{\prime}={\bf\it R}-{\bf\it a}_{y})\\\ 0&\text{otherwise}\end{cases}$ (13)

and the site-dependent potential is represented by

$\epsilon_{\bm{R}}=-\frac{\Delta}{2}e^{-i\bm{Q}\cdot\bm{R}},$ (14)

where $\Delta>0$ gives the strength of the potential and ${\bf\it
Q}=(Q_{x},Q_{y})$ is a wave vector characterizing the periodicity. Here the
$\gamma$-component $Q_{\gamma}$ is limited to $Q_{\gamma}=0$ or
$\pi/a_{\gamma}$. Then the energy difference $\epsilon_{l}-\epsilon_{j}$ in
Eq. (4) for a neighboring pair $(l,j)$ with $({\bf\it R}_{l},{\bf\it
R}_{j})=({\bf\it R},{\bf\it R}\pm{\bf\it a}_{\gamma})$ along the
$\gamma$-direction is evaluated as

$\displaystyle\epsilon_{{\bf\it R}}-\epsilon_{{\bf\it R}\pm{\bf\it
a}_{\gamma}}$ $\displaystyle=-e^{-i\bm{Q}\cdot\bm{R}}\frac{\Delta(1-e^{\mp
i\bm{Q}\cdot\bm{a}_{\gamma}})}{2}$
$\displaystyle=-e^{-i\bm{Q}\cdot\bm{R}}\Delta_{\gamma}({\bf\it Q}),$ (15)

where

$\displaystyle\Delta_{\gamma}({\bf\it
Q})\equiv\begin{cases}\Delta&Q_{\gamma}=\pi/a_{\gamma}\\\
0&Q_{\gamma}=0\end{cases}$ (16)

and the corresponding $\eta_{lj}$ is simplified as

$\displaystyle\eta_{{\bf\it R}{\bf\it R}\pm{\bf\it a}_{\gamma}}({\bf\it
Q})=-e^{-i\bm{Q}\cdot\bm{R}}\eta_{\gamma}({\bf\it Q}),$ (17)

where

$\displaystyle\eta_{\gamma}({\bf\it
Q})\equiv\frac{8Ut^{2}_{\gamma}\Delta_{\gamma}({\bf\it
Q})}{[U^{2}-\Delta^{2}_{\gamma}({\bf\it Q})]^{2}}.$ (18)

In the same manner, $\overline{K}_{il}$ of Eq. (10) for a pair $(i,l)$ with
$({\bf\it R}_{i},{\bf\it R}_{l})=({\bf\it R},{\bf\it R}\pm{\bf\it
a}_{\gamma})$ is obtained as

$\displaystyle\overline{K}_{{\bf\it R}{\bf\it R}\pm{\bf\it
a}_{\gamma}}({\bf\it
Q})=\frac{4Ut^{2}_{\gamma}}{U^{2}-\Delta^{2}_{\gamma}({\bf\it
Q})}\equiv\overline{K}_{\gamma}({\bf\it Q}).$ (19)

In the following, the argument ${\bf\it Q}$ of $\eta_{\gamma}({\bf\it Q})$ and
$\overline{K}_{\gamma}({\bf\it Q})$ is omitted unless necessary. Here, it is
valuable to note that $\overline{K}_{\gamma}$ is equal to the exchange
interaction of the two-leg spin ladder whose Hamiltonian of Eq. (2) is
represented as

${\cal H}^{s}=\overline{K}_{x}\sum_{{\bf\it R}}\bm{S}_{{\bf\it
R}}\cdot\bm{S}_{{\bf\it R}+{\bf\it a}_{x}}+\overline{K}_{y}\sum_{{\bf\it
R}\in\text{leg}\ 0}\bm{S}_{{\bf\it R}}\cdot\bm{S}_{{\bf\it R}+{\bf\it
a}_{y}}.$ (20)

Using these results, $\delta n^{s}_{\bm{k}}$ of a two-leg ladder is expressed
as

$\displaystyle\delta n^{s}_{{\bf\it
k}}={\mathcal{O}}_{x}({\bm{k}}^{\prime})+{\mathcal{O}}_{y}({\bm{k}}^{\prime}),$
(21)

where

$\displaystyle{\mathcal{O}}_{x}({\bm{k}}^{\prime})$
$\displaystyle\equiv-\eta_{x}\sum_{\bm{R}}e^{-i{\bf\it k}^{\prime}\cdot{\bf\it
R}}\sum_{\mathbf{r}=\pm{\bf\it a}_{x}}{\bf\it S}_{{\bf\it R}}\cdot{\bf\it
S}_{{\bf\it R}+{\bf\it r}},$ (22)
$\displaystyle{\mathcal{O}}_{y}({\bm{k}}^{\prime})$
$\displaystyle\equiv-\eta_{y}\sum_{\bm{R}}e^{-i{\bf\it k}^{\prime}\cdot{\bf\it
R}}{\bf\it S}_{{\bf\it R}}\cdot{\bf\it S}_{{\bf\it R}+{\bf\it a}_{y}},$ (23)

and

$\displaystyle{\bf\it k}^{\prime}\equiv{\bf\it k}+{\bf\it Q}.$ (24)

Here, the constant (-1/4) is omitted because it has no contribution to
$N^{s}(\bm{k},\omega)$. The derivation of Eq. (21) is shown in Appendix A.
These results show that $N^{s}(\bm{k},\omega)$ probes spin excitations with
the wave number $\bm{k}^{\prime}$. Therefore we represent and plot
$N^{s}({\bf\it k},\omega)$ as a function of ${\bf\it k}^{\prime}$ in the
following. This is the generalization of the property of the 1D ionic Hubbard
model, where $N^{s}(k,\omega)$ probes excited states with $k^{\prime}=k+\pi$
[33].

Next, consider the current operator (9). Following the same procedure as that
of $\delta n^{s}_{{\bf\it k}}$, $J^{s}_{x}$ is simplified as

$\displaystyle J^{s}_{x}=-\frac{\eta_{x}a_{x}e}{\hbar}\sum_{{\bf\it
R}}e^{-i{\bf\it Q}\dotproduct{\bf\it R}}\left[2\overline{K}_{x}{\bf\it
S}_{{\bf\it R}}\cdot({\bf\it S}_{{\bf\it R}+{\bf\it a}_{x}}\times{\bf\it
S}_{{\bf\it R}-{\bf\it a}_{x}})\right.$ $\displaystyle+\overline{K}_{y}{\bf\it
S}_{{\bf\it R}}\cdot({\bf\it S}_{{\bf\it R}+{\bf\it a}_{x}}\times{\bf\it
S}_{{\bf\it R}+{\bf\it a}_{y}})\left.-\overline{K}_{y}{\bf\it S}_{{\bf\it
R}}\cdot({\bf\it S}_{{\bf\it R}-{\bf\it a}_{x}}\times{\bf\it S}_{{\bf\it
R}+{\bf\it a}_{y}})\right].$ (25)

This tells us that $\sigma^{s}_{x}(\omega)$ probes excited states with
${\bf\it Q}$. For the $y$-direction, we need to be more careful. The current
operator is evaluated as

$\displaystyle
J^{s}_{y}=-\frac{\eta_{y}a_{y}e}{\hbar}\overline{K}_{x}\sum_{{\bf\it
R}}\sum_{{\bf\it r}=\pm{\bf\it a}_{x}}e^{-i\widetilde{{\bf\it
Q}}\dotproduct{\bf\it R}}{\bf\it S}_{{\bf\it R}}\cdot({\bf\it S}_{{\bf\it
R}+{\bf\it a}_{y}}\times{\bf\it S}_{{\bf\it R}+{\bf\it r}}),$ (26)

where

$\displaystyle\widetilde{{\bf\it Q}}={\bf\it Q}+\left(0,\pi/a_{y}\right).$
(27)

These relations suggest that $\sigma^{s}_{y}(\omega)$ probes excited states of
spin systems with $\widetilde{{\bf\it Q}}$, not ${\bf\it Q}$. The additional
wave number $\pi/a_{y}$ results from the boundary conditions in the
$y$-direction; for ${\bf\it R}\in$ leg 1, ${\bf\it R}+{\bf\it a}_{y}={\bf\it
R}-{\bf\it a}_{y}$. It should also be noted that

$\displaystyle{\mathcal{O}}_{x}({\bm{k}}^{\prime})=J^{s}_{x}=0$
$\displaystyle\ \text{for}\ Q_{x}=0\ \text{and}$
$\displaystyle{\mathcal{O}}_{y}({\bm{k}}^{\prime})=J^{s}_{y}=0$
$\displaystyle\ \text{for}\ Q_{y}=0,$ (28)

which are deduced from Eq. (18) and the definitions of the operators. The
derivation of Eqs. (25) and (26) is also explained in Appendix A.

##  3 Results

Figure 2:  Two-leg ladders with the periodic potentials characterized by
${\bf\it Q}=(\pi,\pi)$ and ${\bf\it Q}=(\pi,0)$. The filled and open circles
represent the periodic potentials.

Before discussing the details of our results, we would like to comment on our
numerical calculation. We set $e$, $\hbar$, $a_{x}$, and $a_{y}$ to unity for
simplicity, and we limit $\bm{Q}$ in two cases; $\bm{Q}=(\pi,\pi)$ and
$\bm{Q}=(\pi,0)$ shown in Fig. 2. For the spectral functions, the
$\delta$-function is replaced with the Lorentzian with finite broadening
$\epsilon=0.01t_{x}$, and we impose the periodic boundary condition along the
$x$-direction unless otherwise noted. For numerical calculations of physical
quantities introduced here, we use Lanczos diagonalization. We also note that
quantities with superscript $s$ are those of the two-leg spin ladders.
Quantities without superscripts are from the Hubbard ladders.

First of all, we examine that $N^{s}({\bf\it k}^{\prime},\omega)$ and
$\sigma^{s}_{\gamma}(\omega)$ of the spin ladder well reproduce their
counterparts $N({\bf\it k}^{\prime},\omega)$ and $\sigma_{\gamma}(\omega)$ of
the Hubbard ladder in the low-energy region. Figure 3 shows calculated results
of $N({\bf\it k}^{\prime},\omega)$ and $N^{s}({\bf\it k}^{\prime},\omega)$ for
the ladder with ${\bf\it Q}=(\pi,\pi)$, $t_{y}/t_{x}=1$, $U/t_{x}=20$,
$\Delta/t_{x}=2$, and $N=16$. It can be seen that $N^{s}({\bf\it
k}^{\prime},\omega)$ well reproduce the peak structure of $N({\bf\it
k}^{\prime},\omega)$ for each ${\bf\it k}^{\prime}$. The optical conductivity
spectra $\sigma^{s}_{\gamma}(\omega)$ for $\gamma=x,y$ and their counterpart
$\sigma_{\gamma}(\omega)$, plotted in Fig. 4, also demonstrate that the
contributions from the spin degrees of freedom to these spectral functions
well explain the low-energy behaviors of the original ones.

Figure 3:  (Color online) dynamical charge structure factors, $N({\bf\it
k}^{\prime},\omega)$ (green dashed lines) of the two-leg Hubbard ladder and
$N^{s}({\bf\it k}^{\prime},\omega)$ (black solid lines) of the spin ladder
with ${\bf\it Q}=(\pi,\pi)$, $t_{y}/t_{x}=1$, $U/t_{x}=20$, $\Delta/t_{x}=2$,
and $N=16$.  Figure 4:  (Color online) Optical conductivity
$\sigma_{\gamma}(\omega)$ (green dashed lines) of the two-leg Hubbard ladder
and $\sigma^{s}_{\gamma}(\omega)$ of the spin ladder for $\gamma=x,y$ (black
solid lines) with ${\bf\it Q}=(\pi,\pi)$, $t_{y}/t_{x}=1$, $U/t_{x}=20$,
$\Delta/t_{x}=2$, and $N=16$. The background resulting from the charge degrees
of freedom is subtracted from the result of $\sigma_{\gamma}(\omega)$.

Let us discuss how the multi-triplon excited states appear in these spectral
functions. For this purpose, we consider ladders with large $t_{y}$, where the
series expansion method works well to evaluate the dispersion curves of multi-
triplon excited states [18, 39]. Figure 5 displays $N^{s}({\bf\it
k}^{\prime},\omega)$ of the ladder with ${\bf\it Q}=(\pi,0)$, $t_{y}/t_{x}=2$,
$U/t_{x}=20$, $\Delta/t_{x}=2$, and $N=28$. In addition, spectra of the multi-
triplon states are also plotted. For details on these multi-triplon states,
please see Appendix B. In the upper panel of this figure, we compare
$N^{s}(k^{\prime}_{x},k^{\prime}_{y}=0,\omega)$ with the dispersions of the
two-triplon continuum and the two-triplon singlet bound state. The result of
$\sigma_{x}^{s}(\omega)$ is also shown at the wave number $(\pi,0)$ because
$\sigma_{x}(\omega)$ probes excited states with ${\bf\it Q}$ as discussed
before. These results demonstrate that
$N^{s}(k^{\prime}_{x},k^{\prime}_{y}=0,\omega)$ and $\sigma_{x}^{s}(\omega)$
well observe the above-mentioned two-triplon excitations, including the
continuum and the singlet bound state. In addition, a significant part of the
spectral weight of $N^{s}(k^{\prime}_{x},k^{\prime}_{y}=0,\omega)$ lies close
to the lower boundary of the 2-triplon continuum. In contrast,
$N^{s}(k^{\prime}_{x},k^{\prime}_{y}=\pi,\omega)$, plotted in the lower panel
of Fig. 5, is not so concentrated near the lower boundary of the three-triplon
continuum; its spectral weight lies mainly in two elliptical regions, one
centered at $(k^{\prime}_{x}/\pi,\omega/\overline{K}_{y})=(0.5,2.8)$ and the
other at $(k^{\prime}_{x}/\pi,\omega/\overline{K}_{y})=(1.5,2.8)$. Here, we
make a comment on the relation between our results and those of Ref.
rixs_3tri; the cited paper theoretically shows that the three-triplon bound
states formed by the three-body interaction have significant intensities in
the RIXS spectrum of the two-leg ladder with
$\overline{K}_{x}>\overline{K}_{y}$ [32]. On the other hand,
$N^{s}(k^{\prime}_{x},k^{\prime}_{y}=\pi,\omega)$ of our study does not detect
the three-triplon bound state. This is because the ladder treated here has a
weaker $\overline{K}_{x}$ than $\overline{K}_{y}$, where the three-triplon
bound states are not clearly visible in the spectral functions [32].
Therefore, $N^{s}(k^{\prime}_{x},k^{\prime}_{y}=\pi,\omega)$ would show some
signature of the bound state for $\overline{K}_{x}>\overline{K}_{y}$, which is
beyond the scope of this paper.

These behaviors of $N^{s}({\bf\it k},\omega)$ are qualitatively reproduced
using the bond operator method [34, 35, 36], where the spin operators on a
rung are represented by boson operators, rung-singlet and rung-triplet bosons.
We use the procedure of Ref. gopalan for the two-leg ladder in order to
evaluate $N^{s}({\bf\it k}^{\prime},\omega)$ of the ladder with ${\bf\it
Q}=(\pi,0)$, and the results are shown in Fig. 6. The analytical form of
$N^{s}({\bf\it k}^{\prime},\omega)$ is shown in Appendix C. We can see that
the spectral weight of $N^{s}(k^{\prime}_{x},k^{\prime}_{y}=0,\omega)$ lies in
the region enclosed between the lower and the upper boundaries of the two-
triplon continuum. In particular, the dominant part exists in the vicinity of
the lower boundary. It can also be seen that there is no signal of the bound
state below the continuum, and a non-negligible amount of the weight is found
near the upper boundary of the continuum. This is because the bond operator
theory of Ref. gopalan does not contain the inter-triplon interaction [33].
The result of $N^{s}(k^{\prime}_{x},k^{\prime}_{y}=\pi,\omega)$ also proves
that the two ellipsoidal-shaped structures emerge in the three-triplon
continuum.

We also find that $N^{s}(k^{\prime}_{x},k^{\prime}_{y}=0,\omega)$ vanishes
completely at $k^{\prime}_{x}=\pi$, which is explained as follows. Equation
(28) shows that ${\mathcal{O}}_{y}({\bm{k}}^{\prime})=0$ at ${\bf\it
Q}=(\pi,0)$. Thus $\delta n^{s}_{{\bf\it
k^{\prime}}}={\mathcal{O}}_{x}({\bm{k}}^{\prime})$ in this case. In addition,
${\mathcal{O}}_{x}({\bm{k}}^{\prime})$ is reformulated as

$\displaystyle{\mathcal{O}}_{x}({\bm{k}}^{\prime})=-\eta_{x}(1+e^{-ik_{x}^{\prime}})\sum_{\bm{R}}e^{-i{\bf\it
k}^{\prime}\cdot{\bf\it R}}{\bf\it S}_{{\bf\it R}}\cdot{\bf\it S}_{{\bf\it
R}+{\bf\it a}_{x}},$ (29)

which indicates that ${\mathcal{O}}_{x}({\bm{k}}^{\prime})=0$ at $k_{x}=\pi$.
Hence $N^{s}({\bf\it k}^{\prime},\omega)$ for ${\bf\it Q}=(\pi,0)$ must be
zero at ${\bf\it k}^{\prime}=(\pi,k_{y})$ and around its vicinity
$N^{s}({\bf\it k}^{\prime},\omega)$ becomes small. On the other hand,
$\sigma_{x}^{s}(\omega)$ has a finite spectral weight at ${\bf\it
k}^{\prime}=(\pi,0)$ and detects the two-triplon singlet bound states at this
point. The results of $N^{s}({\bf\it k}^{\prime},\omega)$ and
$\sigma_{x}^{s}(\omega)$ show that they complement each other well for
observing the two-triplon singlet states.

Figure 5:  (Color online) dynamical charge structure factors $N^{s}({\bf\it
k}^{\prime},\omega)$ (black solid lines) and optical conductivity
$\sigma^{s}_{x}(\omega)$ (green solid line) of the two-leg spin ladder with
${\bf\it Q}=(\pi,0)$, $t_{y}/t_{x}=2$, $U/t_{x}=20$, $\Delta/t_{x}=2$, and
$N=28$. Light blue regions show the two-triplon continuum for
$k^{\prime}_{y}=0$ and the three-triplon continuum for $k^{\prime}_{y}=\pi$.
The blue solid line shows the dispersion of the two-triplon singlet bound
state.  Figure 6:  (Color online) $N^{s}({\bf\it k}^{\prime},\omega)$ of the
ladder with ${\bf\it Q}=(\pi,0)$, $t_{y}/t_{x}=2$, $U/t_{x}=20$, and
$\Delta/t_{x}=2$, calculated with the bond operator method. Green solid lines
show upper and lower bounds of the two-triplon continuum for
$k^{\prime}_{y}=0$ and the three-triplon continuum for $k^{\prime}_{y}=\pi$.

Next, we focus on the results of the ladder with ${\bf\it Q}=(\pi,\pi)$.
Figure 7 shows $N^{s}({\bf\it k}^{\prime},\omega)$, $\sigma^{s}_{x}(\omega)$,
$\sigma^{s}_{y}(\omega)$, and dispersions of two-triplon and three-triplon
continuum for ${\bf\it Q}=(\pi,\pi)$. As shown in Eqs. (26) and (27),
$\sigma_{y}(\omega)$ is plotted at $\widetilde{{\bf\it Q}}=(\pi,0)$, while
$\sigma_{x}^{s}(\omega)$ is shown at ${\bf\it Q}=(\pi,\pi)$. Physical
parameters except ${\bf\it Q}$ are the same as those of Fig. 5. The obtained
results demonstrate that the two-triplon and three-triplon excitations are
detected using these spectral quantities. Especially, it is noteworthy that by
changing the wavenumber ${\bf\it Q}$, one can select a multi-triplon state
observable in terms of optical conductivity $\sigma_{x}^{s}(\omega)$.

It should be also found that $N^{s}(k^{\prime}_{x},k^{\prime}_{y}=\pi,\omega)$
for ${\bf\it Q}=(\pi,\pi)$ is quite similar to that of ${\bf\it Q}=(\pi,0)$,
although $N^{s}(k^{\prime}_{x},k^{\prime}_{y}=0,\omega)$ of ${\bf\it
Q}=(\pi,\pi)$ is clearly different from that of ${\bf\it Q}=(\pi,0)$. The
former similarity at $k^{\prime}_{y}=\pi$ results from the fact that
${\mathcal{O}}_{y}({\bm{k}}^{\prime})=0$ for $k^{\prime}_{y}=\pi$ since

$\displaystyle{\mathcal{O}}_{y}({\bm{k}}^{\prime})=-\eta_{y}(1+e^{ik^{\prime}_{y}})\sum_{{\bm{R}}\in\
\text{leg}\ 0}e^{-i{\bf\it k}^{\prime}\cdot{\bf\it R}}{\bf\it S}_{{\bf\it
R}}\cdot{\bf\it S}_{{\bf\it R}+{\bf\it a}_{y}}.$ (30)

Hence, for $k^{\prime}_{y}=\pi$ we obtain

$\displaystyle\delta n^{s}_{{\bf\it k}}={\mathcal{O}}_{x}({\bm{k}}^{\prime})$
(31)

and thus $N^{s}(k^{\prime}_{x},k^{\prime}_{y}=\pi,\omega)$ for ${\bf\it
Q}=(\pi,\pi)$ shows almost the same $({\bf\it k}^{\prime},\omega)$-dependence
as that of $N^{s}(k^{\prime}_{x},k^{\prime}_{y}=\pi,\omega)$ for ${\bf\it
Q}=(\pi,0)$.

The latter difference of $N^{s}(k^{\prime}_{x},k^{\prime}_{y}=0,\omega)$ is a
little bit complicated; the spectral weight of
$N^{s}(k^{\prime}_{x},k^{\prime}_{y}=0,\omega)$ for ${\bf\it Q}=(\pi,\pi)$ is
concentrated around $k^{\prime}_{x}=\pi$ while
$N^{s}(k^{\prime}_{x},k^{\prime}_{y}=0,\omega)$ for ${\bf\it Q}=(\pi,0)$
vanishes at $k^{\prime}_{x}=\pi$. This point originates from the finite
${\mathcal{O}}_{y}$ for ${\bf\it Q}=(\pi,\pi)$. In addition,
$N^{s}(k^{\prime}_{x},k^{\prime}_{y}=0,\omega)$ of ${\bf\it Q}=(\pi,\pi)$ has
a small contribution around $k^{\prime}_{x}=0$, which is understood as
follows. It can be found that $\Delta_{\gamma}({\bf\it Q})=\Delta$ for
${\bf\it Q}=(\pi,\pi)$ and thus

$\displaystyle\eta_{\gamma}({\bf\it
Q})=\frac{8Ut^{2}_{\gamma}\Delta}{[U^{2}-\Delta]^{2}}=2\Delta
U\overline{K}_{\gamma}.$ (32)

Consequently, we obtain

$\displaystyle\delta n^{s}_{{\bf\it k}^{\prime}=(0,0)}=-4\Delta U{\cal
H}^{s},$ (33)

which leads to $\langle\alpha|\delta n^{s}_{{\bf\it
k}^{\prime}=(0,0)}|0\rangle=0$ for $|\alpha\rangle\neq|0\rangle$ and
$N^{s}(k^{\prime}_{x}=k^{\prime}_{y}=0,\omega)=0$.

Figure 7:  (Color online) $N^{s}({\bf\it k}^{\prime},\omega)$ (black solid
lines), $\sigma^{s}_{x}(\omega)$ (green solid line), and
$\sigma^{s}_{y}(\omega)$ (red solid line) of the two-leg spin ladder with
${\bf\it Q}=(\pi,\pi)$, $t_{y}/t_{x}=2$, $U/t_{x}=20$, $\Delta/t_{x}=2$, and
$N=28$. Light blue regions show the two-triplon continuum for
$k^{\prime}_{y}=0$ and the three-triplon continuum for $k^{\prime}_{y}=\pi$.
The blue solid line shows the dispersion of the two-triplon singlet bound
state.

Now, we discuss how to observe these spectral quantities experimentally.
Unfortunately, not so many materials are thought to realize the Hubbard models
with periodic potentials. Besides TTF-CA and its derivatives [40, 41, 42],
which are quasi-one-dimensional organic compounds with alternating potentials,
cold atom systems have been considered as promising realizations [43].
Although spectral quantities of cold atom systems are also observed
experimentally [44, 45], their spectral resolution seems insufficient to
detect the spectral weights associated with the spin degrees of freedom
discussed here.

Here, we explore materials with randomness as alternative systems with site-
dependent potentials. As a typical example, it has been pointed out that an
organic ladder material ($\alpha$-DT-TTF)2[Au(mnt)2] has a weak disorder due
to the donor molecule $\alpha$-DT-TTF [46, 47]. This molecule has two types of
structures, called cis- and trans-forms, which are randomly arranged in the
ladder material [46]. The difference in the structures naturally leads to the
difference in the energy level, which can lead to randomness in the site-
dependent potential in the form that two types of sites with different
potential energies are randomly arranged. Furthermore, X-ray irradiation has
been known to introduce randomness into organic correlated systems [48], which
may allow us to treat various organic ladders with random potentials.

With ($\alpha$-DT-TTF)2[Au(mnt)2] in mind, we treat a two-leg ladder with the
simple random potential expressed by

$\epsilon_{l}=-\frac{\Delta}{2}\xi_{l},$ (34)

where $\xi_{l}$ is the random variable which takes 1 or -1 with the same
probability. Figure 8 (a) displays the calculated results of the optical
conductivity $\sigma_{x}^{s}(\omega)$ of the ladder with the random potential.
The physical parameters are set to $t_{y}/t_{x}=2$, $U/t_{x}=20$,
$\Delta/t_{x}=2$, $N=28$, and we take the average for 200 samples. In
addition, we use the same energy scale
$\overline{K}_{y}=8Ut^{2}_{y}/(U^{2}-\Delta)^{2}$ as that of the case with the
periodic potential of $\mathbf{Q}=(\pi,\pi)$ because of its usefulness
although $\overline{K}_{y}$ is no longer the exact exchange interaction along
the rung. We have confirmed that the ladder with the random potential has the
characteristic spectral structures in $\sigma_{x}^{s}(\omega)$; there are a
sharp peak at $\omega/\overline{K}_{y}\sim 1.8$, a shoulder structure at
around $\omega/\overline{K}_{y}=2.1$, and a small broad peak at around
$\omega/\overline{K}_{y}=2.8$. Judging from the energy scales where these
structures appear, the sharp peak and shoulder are thought to originate from
the two-triplon states, and the broad peak from the three-triplon states. For
further examination, the semi-log plot of the same data is shown in Fig. 8
(b), where numerical results of ladders with periodic potentials are also
plotted. Comparisons among these data suggest that the sharp peak stems from
the two-triplon bound state and the shoulder from the continuum. In addition,
the broad peak corresponds to the three-triplon continuum.

Figure 8:  (Color online) (a) $\sigma^{s}_{x}(\omega)$ of the two-leg ladder
with $t_{y}/t_{x}=2$, $U/t_{x}=20$, $\Delta/t_{x}=2$ and $N=28$ under the
random potential. (b) The semi-log plot of $\sigma^{s}_{x}(\omega)$ shown in
(a) (black solid lines). The results for ladders with periodic potentials,
$\sigma^{s}_{x}(\omega)$ for ${\bf\it Q}=(\pi,0)$ (green dot-dashed line),
$\sigma^{s}_{y}(\omega)$ for ${\bf\it Q}=(\pi,\pi)$ (red dashed line), and
$\sigma^{s}_{x}(\omega)$ for ${\bf\it Q}=(\pi,\pi)$ (blue dotted line) are
also shown here. The wave number ${\bf\it Q}$ for each quantity is depicted as
the subscript.

Whether the spectral signature of the three-triplon state is well separated
from that of the two-triplon states depends strongly on the competition
between two energy scales, $\overline{K}_{y}$ and $\overline{K}_{x}$. The
former is roughly equal to the creation energy of one triplon, and the latter
is proportional to the bandwidth of the triplon band. Figure 9, which shows
$\sigma^{s}_{x}(\omega)$ for $U/t_{y}=10$, $N=28$, and $\Delta/t_{y}=1$,
demonstrate how these spectral structures change with increasing
$\overline{K}_{x}/\overline{K}_{y}=(t_{x}/t_{y})^{2}$. For small
$\overline{K}_{x}/\overline{K}_{y}$, the peak of the three-triplon band around
$\omega\sim 3\overline{K}_{y}$ is far from that of the two-triplon band for
$\omega\sim 2\overline{K}_{y}$. As $\overline{K}_{x}/\overline{K}_{y}$
increases, the shoulder structure at $\omega/\overline{K}_{y}\sim 2$
originating from the two-triplon continuum becomes larger and broader. At
$\overline{K}_{x}/\overline{K}_{y}=1/2$, the shoulder merges with the peak of
the three-triplon continuum, and the resulting broad shoulder ranging from
$\omega/\overline{K}_{y}\sim 2$ to $\omega/\overline{K}_{y}\sim 3.5$ is
reminiscent of the multi-triplon state. At
$\overline{K}_{x}/\overline{K}_{y}=1$, it becomes difficult to find out the
signature of the three-triplon state. Since the ladder material ($\alpha$-DT-
TTF)2[Au(mnt)2] is considered to be $(t_{x}/t_{y})^{2}=0.54$, [46] the broad
shoulder structure state might be detected experimentally. Furthermore,
materials with smaller $\overline{K}_{x}/\overline{K}_{y}$, such as
(C5H12N)2CuBr4 ($\overline{K}_{x}/\overline{K}_{y}\sim 1/4$) [12, 49],
irradiated with X-rays could be much better ones to observe these multi-
triplon states.

Figure 9:  (Color online) optical conductivity $\sigma^{s}_{x}(\omega)$ of the
two-leg ladders under the random potential with $U/t_{y}=10$,
$\Delta/t_{y}=1$, $N=28$, with increasing
$\overline{K}_{x}/\overline{K}_{y}=(t_{x}/t_{y})^{2}$.

##  4 Summary

In this work, we have studied the low-energy spectral properties of the two-
leg Hubbard ladders with site-dependent potentials. Our approach combining the
analytical treatment for the Hubbard models with the numerical diagonalization
demonstrated that the spectral functions $N^{s}(\bm{k}^{\prime},\omega)$ and
$\sigma^{s}_{\gamma}(\omega)$, the contributions of the spin degrees of
freedom to the dynamical charge structure factor and the optical conductivity
along the $\gamma=x,y$ direction, detect the two- and three-triplon
excitations of the two-leg spin ladders. The key ingredient to detect these
multi-triplon excitations is the finite site-dependent potential of the
ladder. In our framework, this point is expressed in terms of the quantity
$\eta_{lj}$, which appears in the definitions of
$N^{s}(\bm{k}^{\prime},\omega)$ and $\sigma^{s}_{\gamma}(\omega)$. Numerical
results obtained using the Lanczos diagonalization are consistent with the
results of bond operator theory. We have also examined that, by tuning the
periodicity of the potential, $\sigma^{s}_{\gamma}(\omega)$ can select
detectable excited states; for example, $\sigma^{s}_{x}(\omega)$ probes either
the two-triplon bound state or the three-triplon continuum. Furthermore, we
confirmed that the two-leg ladder with the random potential also has finite
spectral weights in $\sigma^{s}_{x}(\omega)$ stemming from the two- and three-
triplon states, which suggest that ladder materials with intrinsic or acquired
randomness and strong intra-rung couplings are good candidates for
experimentally observing and identifying these multi-triplon states.

## Acknowledgment

This work was supported by a financial support from the University of Tsukuba,
Pre-Strategic Initiatives ’Development Center for High-Function and High-
Performance Organic-Inorganic Spin Electronics’.

##  Appendix A Operators of Two-leg Ladders

Here, we derive some relations between operators used in this work. Now let us
start with the derivation of Eqs. (21-24) for the two-leg ladders with the
nearest-neighbor hoppings and the periodic potentials. Since the summation in
Eq. (6) with respect to $l$ and $j$ is restricted to the nearest-neighbor
pairs, we obtain

$\displaystyle\delta n^{s}_{\bm{k}}=$
$\displaystyle-\eta_{x}\sum_{\bm{R}}e^{-i{\bf\it k}^{\prime}\cdot{\bf\it
R}}\sum_{\mathbf{r}=\pm{\bf\it a}_{x}}{\bf\it S}_{{\bf\it R}}\cdot{\bf\it
S}_{{\bf\it R}+{\bf\it r}}$
$\displaystyle-\eta_{y}\sum_{{\bm{R}}\in\text{leg}\ 0}e^{-i{\bf\it
k}^{\prime}\cdot{\bf\it R}}{\bf\it S}_{{\bf\it R}}\cdot{\bf\it S}_{{\bf\it
R}+{\bf\it a}_{y}}$ $\displaystyle-\eta_{y}\sum_{{\bm{R}}\in\text{leg}\
1}e^{-i{\bf\it k}^{\prime}\cdot{\bf\it R}}{\bf\it S}_{{\bf\it R}}\cdot{\bf\it
S}_{{\bf\it R}-{\bf\it a}_{y}},$ (35)

where the coefficient $\eta_{lj}$ for the pair $(lj)$ is replaced by
$\eta_{{\bf\it R}{\bf\it R}\pm{\bf\it a}_{\gamma}}$ of Eq. (17), and the
constant (-1/4) in Eq. (3) is omitted as noted above. Then using the periodic
boundary condition for the $y$-direction, ${\bf\it R}+{\bf\it a}_{y}={\bf\it
R}-{\bf\it a}_{y}$ for ${\bf\it R}\in$ leg 1, the second and third terms are
summarized as

$\displaystyle\delta n^{s}_{{\bf\it k}}$
$\displaystyle=-\eta_{x}\sum_{\bm{R}}e^{-i{\bf\it k}^{\prime}\cdot{\bf\it
R}}\sum_{\mathbf{r}=\pm{\bf\it a}_{x}}{\bf\it S}_{{\bf\it R}}\cdot{\bf\it
S}_{{\bf\it R}+{\bf\it r}}$ $\displaystyle-\eta_{y}\sum_{\bm{R}}e^{-i{\bf\it
k}^{\prime}\cdot{\bf\it R}}{\bf\it S}_{{\bf\it R}}\cdot{\bf\it S}_{{\bf\it
R}+{\bf\it a}_{y}}.$ (36)

As for the current operator $J^{s}_{x}$, Eq. (9) for $\gamma=x$ is evaluated
as

$\displaystyle J^{s}_{x}$
$\displaystyle=\frac{e}{\hbar}\eta_{x}\overline{K}_{x}\left[\sum_{{\bm{R}}}e^{-i{\bf\it
Q}\cdot{\bf\it R}}(-a_{x}){\bf\it S}_{{\bf\it R}}\cdot({\bf\it S}_{{\bf\it
R}+{\bf\it a}_{x}}\times{\bf\it S}_{{\bf\it R}-{\bf\it a}_{x}})\right.$
$\displaystyle+\left.\sum_{{\bm{R}}}e^{-i{\bf\it Q}\cdot{\bf\it
R}}(a_{x}){\bf\it S}_{{\bf\it R}}\cdot({\bf\it S}_{{\bf\it R}-{\bf\it
a}_{x}}\times{\bf\it S}_{{\bf\it R}+{\bf\it a}_{x}})\right]$
$\displaystyle+\frac{e}{\hbar}\eta_{x}\overline{K}_{y}\left[\sum_{{\bm{R}}\in\text{leg}\
0}e^{-i{\bf\it Q}\cdot{\bf\it R}}(-a_{x}){\bf\it S}_{{\bf\it R}}\cdot({\bf\it
S}_{{\bf\it R}+{\bf\it a}_{x}}\times{\bf\it S}_{{\bf\it R}+{\bf\it
a}_{y}})\right.$ $\displaystyle+\sum_{{\bm{R}}\in\text{leg}\ 0}e^{-i{\bf\it
Q}\cdot{\bf\it R}}(a_{x}){\bf\it S}_{{\bf\it R}}\cdot({\bf\it S}_{{\bf\it
R}-{\bf\it a}_{x}}\times{\bf\it S}_{{\bf\it R}+{\bf\it a}_{y}})$
$\displaystyle+\sum_{{\bm{R}}\in\text{leg}\ 1}e^{-i{\bf\it Q}\cdot{\bf\it
R}}(-a_{x}){\bf\it S}_{{\bf\it R}}\cdot({\bf\it S}_{{\bf\it R}+{\bf\it
a}_{x}}\times{\bf\it S}_{{\bf\it R}-{\bf\it a}_{y}})$
$\displaystyle+\left.\sum_{{\bm{R}}\in\text{leg}\ 1}e^{-i{\bf\it
Q}\cdot{\bf\it R}}(a_{x}){\bf\it S}_{{\bf\it R}}\cdot({\bf\it S}_{{\bf\it
R}-{\bf\it a}_{x}}\times{\bf\it S}_{{\bf\it R}-{\bf\it a}_{y}})\right].$ (37)

Now, using the same boundary condition for $y$-direction, we obtain

$\displaystyle J^{s}_{x}$
$\displaystyle=-\frac{e}{\hbar}\eta_{x}a_{x}\sum_{{\bm{R}}}e^{-i{\bf\it
Q}\cdot{\bf\it R}}\left[2\overline{K}_{x}{\bf\it S}_{{\bf\it R}}\cdot({\bf\it
S}_{{\bf\it R}+{\bf\it a}_{x}}\times{\bf\it S}_{{\bf\it R}-{\bf\it
a}_{x}})\right.$ $\displaystyle+\overline{K}_{y}{\bf\it S}_{{\bf\it
R}}\cdot({\bf\it S}_{{\bf\it R}+{\bf\it a}_{x}}\times{\bf\it S}_{{\bf\it
R}+{\bf\it a}_{y}})$ $\displaystyle-\left.\overline{K}_{y}{\bf\it S}_{{\bf\it
R}}\cdot({\bf\it S}_{{\bf\it R}-{\bf\it a}_{x}}\times{\bf\it S}_{{\bf\it
R}+{\bf\it a}_{y}})\right].$ (38)

For $J^{s}_{y}$, we need a slightly different treatment. The definition (9)
for $\gamma=y$ is evaluated as

$\displaystyle J^{s}_{y}=$
$\displaystyle\frac{e\eta_{y}\overline{K}_{x}}{\hbar}\left[\sum_{{\bf\it
R}\in\text{leg}\ 0}e^{-i{\bf\it Q}\dotproduct{\bf\it R}}\sum_{{\bf\it
r}=\pm{\bf\it a}_{x}}(-a_{y}){\bf\it S}_{{\bf\it R}}\cdot({\bf\it S}_{{\bf\it
R}+{\bf\it a}_{y}}\times{\bf\it S}_{{\bf\it R}+{\bf\it r}})\right.$
$\displaystyle\left.+\sum_{{\bf\it R}\in\text{leg}\ 1}e^{-i{\bf\it
Q}\dotproduct{\bf\it R}}\sum_{{\bf\it r}=\pm{\bf\it a}_{x}}(a_{y}){\bf\it
S}_{{\bf\it R}}\cdot({\bf\it S}_{{\bf\it R}+{\bf\it a}_{x}}\times{\bf\it
S}_{{\bf\it R}+{\bf\it r}})\right].$ (39)

Here we introduce the modified wave number $\widetilde{{\bf\it Q}}={\bf\it
Q}+\left(0,\pi/a_{y}\right)$ of Eq. (27), and then obtain

$\displaystyle e^{-i\widetilde{{\bf\it Q}}\cdot{\bf\it
R}}=\begin{cases}e^{-i{\bf\it Q}\cdot{\bf\it R}}&\text{for}\ {\bf\it R}\ \in\
\text{leg}\ 0\\\ -e^{-i{\bf\it Q}\cdot{\bf\it R}}&\text{for}\ {\bf\it R}\ \in\
\text{leg}\ 1\end{cases}.$ (40)

Therefore, Equation (39) is summarized as

$\displaystyle
J^{s}_{y}=-\frac{\eta_{y}a_{y}e}{\hbar}\overline{K}_{x}\sum_{{\bf\it
R}}\sum_{{\bf\it r}=\pm{\bf\it a}_{x}}e^{-i\widetilde{{\bf\it
Q}}\dotproduct{\bf\it R}}{\bf\it S}_{{\bf\it R}}\cdot({\bf\it S}_{{\bf\it
R}+{\bf\it a}_{y}}\times{\bf\it S}_{{\bf\it R}+{\bf\it r}}).$ (41)

##  Appendix B Dispersion of Multi-Triplon States

In this part, we make some comments on the dispersion relations of multi-
triplon states. As for the one-triplon dispersion, we use the result of Ref.
oitmaa. The dispersion $\omega_{1}(k_{x})$ of the one-triplon with the wave
number $k_{x}$ along $x$-direction is given by

$\displaystyle\omega_{1}(k_{x})=\overline{K}_{y}\sum_{n=0}^{8}\sum_{m=0}^{8}a_{n,m}\lambda^{n}\cos(mk_{x}),$
(42)

where $\lambda=\overline{K}_{x}/\overline{K}_{y}$ is the expansion parameter
and the coefficient $a_{n,m}$ is shown in Table III of Ref. oitmaa.

An eigenvalue of an $n$-triplon unbound state with the wave number $p_{\rm
tot}$ is given by the following simple summation

$\displaystyle\omega^{C}_{nt}(p_{\rm
tot})=\sum_{i=1,\cdots,n}\omega_{1}(p_{i})$ (43)

under the limitation of $\sum_{i=1,\cdots,n}p_{i}=p_{\rm tot}$. By taking all
allowed combinations of $(p_{1},\cdots,p_{n})$, we obtain the $n$-triplon
continuum at $k_{x}=p_{\rm tot}$. Here, we note that one triplon carries the
wave number $k_{y}=\pi$ along $y$-direction. Hence $n$-triplon states have the
wave number $k_{y}=n\pi$ (mod $2\pi$).

As for the two-triplon singlet bound state, we follow the procedure of Ref.
kumar; we employ the results of Ref. zheng for the dispersion relation
$\omega^{S}_{2t}(k_{x})$ of the two-triplon singlet bound state, which is
given by

$\displaystyle\omega^{S}_{2t}(k_{x})=$
$\displaystyle\overline{K}_{y}\left[2-\frac{3}{2}\lambda+\frac{19}{16}\lambda^{2}-\frac{9}{32}\lambda^{3}\right.$
$\displaystyle-\left(\frac{1}{2}\lambda-\frac{1}{8}\lambda^{2}+\frac{51}{128}\lambda^{3}\right)\cos(k_{x})$
$\displaystyle-\left(\frac{5}{16}\lambda^{2}+\frac{21}{32}\lambda^{3}\right)\cos(2k_{x})$
$\displaystyle\left.-\frac{37}{128}\lambda^{3}\cos(3k_{x})\right].$ (44)

##  Appendix C Bond Operator Theory

We here review the bond operator theory [34, 35, 36], which was developed to
investigate the ground state phase diagram of $S=1/2$ two-dimensional
frustrated Heisenberg model [34]. The application to the two-leg ladder was
made by Gopalan et al. [35], and more excellent treatment of the ladder is
discussed later [36]. Here, we follow the procedure of Ref. gopalan and pick
up only essential points to reproduce our calculations for the convenience of
the reader.

Let us introduce the bond operators. The $\alpha$-component of the spin
operator of the site ${\bf\it R}=(l,n)$ on the leg $n=0$ or $1$ is represented
as

$\displaystyle S^{\alpha}_{0l}$ $\displaystyle\equiv S^{\alpha}_{{\bf\it
R}=(l,0)}=\frac{1}{2}(s^{\dagger}_{l}t_{l\alpha}+t^{\dagger}_{l\alpha}s_{l}-\varepsilon_{\alpha\beta\gamma}t^{\dagger}_{l\beta}t_{l\gamma})$
(45)

and

$\displaystyle S^{\alpha}_{1l}$ $\displaystyle\equiv S^{\alpha}_{{\bf\it
R}=(l,1)}=\frac{1}{2}(-s^{\dagger}_{l}t_{l\alpha}-t^{\dagger}_{l\alpha}s_{l}-\varepsilon_{\alpha\beta\gamma}t^{\dagger}_{l\beta}t_{l\gamma}),$
(46)

where $s_{l}$ and $t_{l\alpha}(\alpha=x,y,z)$ are bond operators for the
singlet and triplet states of a rung labeled by $l$. Here, we use the
summation convention for repeated indices and the totally antisymmetric tensor
$\varepsilon_{\alpha\beta\gamma}$. These operators have bosonic statistics to
satisfy the algebra of the spin operator $S^{\alpha}_{{\bf\it R}}$. In
addition, they follow the hard-core condition

$\displaystyle s^{\dagger}_{l}s_{l}+t^{\dagger}_{l\alpha}t_{l\alpha}=1$ (47)

for each rung.

Using the bond operators, we rewrite the Hamiltonian of the two-leg ladder
(20) as

$\displaystyle{\cal H}^{s}=H_{0}+\lambda H_{1},$ (48)

where

$\displaystyle H_{0}$
$\displaystyle=\sum_{l}\overline{K}_{y}\left(-\frac{3}{4}s^{\dagger}_{l}s_{l}+\frac{1}{4}t^{\dagger}_{l\alpha}t_{l\alpha}\right)$
$\displaystyle-\sum_{l}\mu_{l}(s^{\dagger}_{l}s_{l}+t^{\dagger}_{l\alpha}t_{l\alpha}-1)$
(49)

and

$\displaystyle
H_{1}=\frac{\overline{K}_{y}}{2}\sum_{l}\left(t^{\dagger}_{l\alpha}t_{l+1\alpha}s^{\dagger}_{l+1}s_{l}+t^{\dagger}_{l\alpha}t^{\dagger}_{l+1\alpha}s_{l}s_{l+1}+\text{H.c.}\right).$
(50)

Here we note that the chemical potential $\mu_{l}$ is introduced to impose the
hard-core condition (47) and the quartic term with respect to
$t^{(\dagger)}_{l}$ ($H_{2}$ in Ref. gopalan ) is omitted here following the
theoretical treatment [35].

Now, the replacement $s^{(\dagger)}_{l}\to\langle
s^{(\dagger)}_{l}\rangle=\bar{s}$ using the mean field $\bar{s}$ is introduced
to diagonalize the Hamiltonian (48). In addition, the chemical potential
$\mu_{l}$ is replaced by the site-independent one $\mu$ because of the
translational invariance of the ladder system. Applying the Fourier
transformation of the operations
$t_{l\alpha}=N_{x}^{-1}\sum_{k}t_{k\alpha}e^{ikl}$, we obtain the mean-field
Hamiltonian

$\displaystyle H_{m}(\mu,\bar{s})=$ $\displaystyle
N_{x}\left(-\frac{3}{4}\overline{K}_{y}\bar{s}^{2}-\mu\bar{s}^{2}+\mu\right)$
$\displaystyle+\sum_{k}[\Lambda_{k}t^{\dagger}_{k\alpha}t_{k\alpha}+\Delta_{k}(t^{\dagger}_{k\alpha}t^{\dagger}_{-k\alpha}+t_{k\alpha}t_{-k\alpha})]$
(51)

where

$\displaystyle\Lambda_{k}=\frac{1}{4}-\mu+\lambda\overline{K}_{y}\bar{s}^{2}\cos
k,\quad\Delta_{k}=\frac{1}{2}\overline{K}_{y}\bar{s}^{2}\cos k,$ (52)

and $N_{x}$ is the number of rungs. Using the Bogoliubov transformation

$\displaystyle\gamma_{k\alpha}=\cosh\theta_{k}t_{k\alpha}+\sinh\theta_{k}t^{\dagger}_{-k\alpha}$
(53)

with

$\displaystyle\cosh\theta_{k}$
$\displaystyle=\frac{1}{\sqrt{2}}\sqrt{\frac{\Lambda_{k}}{\omega_{k}}+1}$ (54)

and

$\displaystyle\sinh\theta_{k}$
$\displaystyle=\frac{1}{\sqrt{2}}\text{sgn}(\Delta_{k})\sqrt{\frac{\Lambda_{k}}{\omega_{k}}-1},$
(55)

we diagonalize the mean-field Hamiltonian as

$\displaystyle H_{m}(\mu,\bar{s})=$ $\displaystyle
N_{x}\left(-\frac{3}{4}\overline{K}_{y}\bar{s}^{2}-\mu\bar{s}^{2}+\mu\right)-\frac{N_{x}}{2}\left(\frac{\overline{K}_{y}}{4}-\mu\right)$
$\displaystyle+\sum_{k}\omega_{k}\left(\gamma^{\dagger}_{k\alpha}\gamma_{k\alpha}+\frac{1}{2}\right).$
(56)

Here, the quantity

$\displaystyle\omega_{k}=\sqrt{\Lambda^{2}_{k}-(2\Delta_{k})^{2}}$ (57)

is the dispersion of the Bogoliubov boson, which corresponds to the triplon of
the two-leg ladder.

The next step is to obtain physical quantities within this mean-field theory.
Following the method in the original paper [35], we calculate the mean fields
$\bar{s}$ and $\mu$ at the ground state by numerically solving the saddle-
point equations

$\displaystyle\left\langle\frac{\partial H_{m}}{\partial\mu}\right\rangle=0\
\text{and}\ \left\langle\frac{\partial
H_{m}}{\partial\bar{s}}\right\rangle=0.$ (58)

Then, we evaluate the dynamical charge structure factor
$N^{s}(\bm{k}^{\prime},\omega)$ for ${\bf\it Q}=(\pi,0)$, where the charge
disproportionation $\delta n^{s}_{{\bm{k}}^{\prime}}$ of Eq. (21) is given by

$\displaystyle\delta n^{s}_{{\bf\it k}}$
$\displaystyle={\mathcal{O}}_{x}({\bm{k}}^{\prime})$
$\displaystyle=-\eta_{x}\sum_{\bm{R}}e^{-i{\bf\it k}^{\prime}\cdot{\bf\it
R}}\sum_{\mathbf{r}=\pm{\bf\it a}_{x}}{\bf\it S}_{{\bf\it R}}\cdot{\bf\it
S}_{{\bf\it R}+{\bf\it r}}.$ (59)

For this purpose, we represent

$\displaystyle\delta n^{s}_{{{\bf\it k}}^{\prime}}$ $\displaystyle=\delta
n^{0}_{k_{x}^{\prime}}\quad\text{for}\quad\bm{k}^{\prime}=(k_{x}^{\prime},0)$
(60) $\displaystyle=\delta
n^{\pi}_{k_{x}^{\prime}}\quad\text{for}\quad\bm{k}^{\prime}=(k_{x}^{\prime},\pi),$
(61)

and then express these quantities with the boson operators.

First, we deal with

$\displaystyle\delta
n^{0}_{k_{x}^{\prime}}=-\eta_{x}\sum_{\bm{l}}e^{-ik_{x}^{\prime}l}\sum_{r=\pm
1}({\bf\it S}_{0l}\cdot{\bf\it S}_{0l+r}+{\bf\it S}_{1l}\cdot{\bf\it
S}_{1l+r}),$ (62)

which is necessary to calculate
$N^{s}(k^{\prime}_{x},k^{\prime}_{y}=0,\omega)$. Here, we use the boson
representation

$\displaystyle\bm{S}_{0l}\cdot\bm{S}_{0j}+\bm{S}_{1l}\cdot\bm{S}_{1j}$
$\displaystyle=\frac{\bar{s}^{2}}{2}\left(t^{\dagger}_{l\alpha}t_{j\alpha}+t^{\dagger}_{l\alpha}t^{\dagger}_{j\alpha}+t^{\dagger}_{j\alpha}t_{l\alpha}+t_{j\alpha}t_{l\alpha}\right),$
(63)

where the replacement $s^{(\dagger)}_{i}\to\bar{s}$ is introduced, and only
the two-body terms of the boson operators are left. Using the Fourier
transformation and the Bogoliubov boson $\gamma^{(\dagger)}_{k\alpha}$,
$\delta n^{0}_{k_{x}^{\prime}}$ is expressed as

$\displaystyle\delta n^{0}_{k_{x}^{\prime}}$
$\displaystyle=\eta_{x}\bar{s}^{2}\sum_{p}\left\\{-[\cos
p+\cos(p-k_{x}^{\prime})]\cosh\theta_{p-k_{x}^{\prime}}\sinh\theta_{p}\right.$
$\displaystyle\times\gamma^{\dagger}_{p-k_{x}^{\prime}\alpha}\gamma^{\dagger}_{-p\alpha}+\cos
p\cosh\theta_{p+k_{x}^{\prime}}\cosh\theta_{p}\gamma^{\dagger}_{-p-k_{x}^{\prime}\alpha}\gamma^{\dagger}_{p\alpha}$
$\displaystyle+\cos(p-k_{x}^{\prime})\sinh\theta_{p-k_{x}^{\prime}}\sinh\theta_{p}\gamma^{\dagger}_{p-k_{x}^{\prime}\alpha}\gamma^{\dagger}_{-p\alpha}$
$\displaystyle-[\cos
p+\cos(p-k_{x}^{\prime})]\sinh\theta_{p-k_{x}^{\prime}}\cosh\theta_{p}\gamma_{-p+k_{x}^{\prime}\alpha}\gamma_{p\alpha}$
$\displaystyle+\cos
p\sinh\theta_{p+k_{x}^{\prime}}\sinh\theta_{p}\gamma_{p+k_{x}^{\prime}\alpha}\gamma_{-p\alpha}$
$\displaystyle+\left.\cos(p-k_{x}^{\prime})\cosh\theta_{p-k_{x}^{\prime}}\cosh\theta_{p}\gamma_{-p+k_{x}^{\prime}\alpha}\gamma_{p\alpha}\right\\},$
(64)

where the two-boson creation or annihilation terms are left to treat the two-
triplon continuum contributing to
$N^{s}(k^{\prime}_{x},k^{\prime}_{y}=0,\omega)$. Then
$N^{s}(k^{\prime}_{x},k^{\prime}_{y}=0,\omega)$ is represented as

$\displaystyle N^{s}(k^{\prime}_{x},k^{\prime}_{y}=0,\omega)=\sum_{\alpha\neq
0}|\langle\alpha|\delta n^{0}_{k^{\prime}_{x}}|0\rangle|^{2}\delta(\omega-
E_{\alpha}+E_{0})$
$\displaystyle=6\eta_{x}^{2}\bar{s}^{4}\sum_{p}\left[\cos\left(p+\frac{k^{\prime}_{x}}{2}\right)\cos\left(\frac{k^{\prime}_{x}}{2}\right)\sinh(\theta_{p}+\theta_{p+k^{\prime}_{x}})\right.$
$\displaystyle-\left.\cos p\cosh(\theta_{p}+\theta_{p+k^{\prime}_{x}})\right]$
$\displaystyle\times\left[\cos\left(p+\frac{k^{\prime}_{x}}{2}\right)\cos\left(\frac{k^{\prime}_{x}}{2}\right)\sinh(\theta_{p}+\theta_{p+k^{\prime}_{x}})\right.$
$\displaystyle-\left.\cos\left(p+\frac{k^{\prime}_{x}}{2}\right)\cos\left(\frac{k^{\prime}_{x}}{2}\right)\cosh(\theta_{p}+\theta_{p+k^{\prime}_{x}})\right]\delta(\omega-\omega_{p}-\omega_{p+k^{\prime}_{x}}).$
(65)

As for

$\displaystyle\delta
n^{\pi}_{k_{x}^{\prime}}=-\eta_{x}\sum_{\bm{l}}e^{-ik_{x}^{\prime}l}\sum_{r=\pm
1}({\bf\it S}_{0l}\cdot{\bf\it S}_{0l+r}-{\bf\it S}_{1l}\cdot{\bf\it
S}_{1l+r}),$ (66)

we use the same notation as

$\displaystyle\bm{S}_{0l}\cdot\bm{S}_{0j}-\bm{S}_{1l}\cdot\bm{S}_{1j}$
$\displaystyle=-\frac{i\bar{s}}{2}\varepsilon_{\alpha\beta\delta}\left(t^{\dagger}_{l\alpha}t^{\dagger}_{j\beta}t_{j\delta}+t_{l\alpha}t^{\dagger}_{j\beta}t_{j\delta}+t^{\dagger}_{j\alpha}t^{\dagger}_{l\beta}t_{l\delta}+t_{j\alpha}t^{\dagger}_{l\beta}t_{l\delta}\right),$
(67)

and then obtain

$\displaystyle\delta n^{\pi}_{k_{x}^{\prime}}$
$\displaystyle=-\frac{i\bar{s}}{2}\eta_{x}\varepsilon_{\alpha\beta\delta}(1+e^{-ik_{x}^{\prime}})\frac{1}{\sqrt{N_{x}}}\sum_{pq}\left[e^{i(k_{x}^{\prime}+p)}+e^{-p}\right]$
$\displaystyle\times(\cosh\theta_{p}-\sin\theta_{p})(\cosh\theta_{q}\sinh\theta_{k_{x}^{\prime}+p+q}\gamma^{\dagger}_{p\alpha}\gamma^{\dagger}_{q\beta}\gamma^{\dagger}_{-k_{x}^{\prime}-p-q\delta}$
$\displaystyle+\sinh\theta_{q}\cosh\theta_{k_{x}^{\prime}+p+q}\gamma_{-p\alpha}\gamma_{-q\beta}\gamma_{k_{x}^{\prime}+p+q\delta}).$
(68)

It should be noted that the three-boson creation or annihilation terms are
left in Eq. (68). Finally, $N^{s}(k^{\prime}_{x},k^{\prime}_{y}=\pi,\omega)$
is expressed as

$\displaystyle N^{s}(k^{\prime}_{x},k^{\prime}_{y}=\pi,\omega)$
$\displaystyle=\sum_{\alpha\neq 0}|\langle\alpha|\delta
n^{\pi}_{k^{\prime}_{x}}|0\rangle|^{2}\delta(\omega-E_{n}+E_{n})$
$\displaystyle=\frac{12\eta_{x}^{2}\bar{s}^{2}}{N_{x}}(1+\cos
k^{\prime}_{x})\sum_{pq}\cos\left(p+\frac{k^{\prime}_{x}}{2}\right)$
$\displaystyle\times\cosh\theta_{q}\sinh\theta_{p+q+k^{\prime}_{x}}(\cosh\theta_{p}-\sinh\theta_{p})$
$\displaystyle\times\left[\cos\left(p+\frac{k^{\prime}_{x}}{2}\right)(\cosh\theta_{p}-\sinh\theta_{p})\sinh(\theta_{p+q+k^{\prime}_{x}+}-\theta_{q})\right.$
$\displaystyle+\cos\left(p+\frac{k^{\prime}_{x}}{2}\right)(\cosh\theta_{q}-\sinh\theta_{q})\sinh(\theta_{q}-\theta_{p+q+k^{\prime}_{x}+})$
$\displaystyle+\left.\cos\left(p+q+\frac{k^{\prime}_{x}}{2}\right)(\cosh\theta_{p+q+k^{\prime}_{x}}-\sinh\theta_{p+q+k^{\prime}_{x}})\sinh(\theta_{q}-\theta_{p})\right]$
$\displaystyle\times\delta(\omega-\omega_{p}-\omega_{q}-\omega_{p+q+k^{\prime}_{x}}).$
(69)

## References

  * [1] E. Dagotto and T. M. Rice, Science 271, 618 (1996).
  * [2] E. Dagotto, Rep. Prog. Phys. 62, 1525 (1999).
  * [3] K. P. Schmidt and G. S. Uhrig, Phys. Rev. Lett. 90, 227204 (2003).
  * [4] J. P. Renard, M. Verdaguer, L. P. Regnault, W. A. C. Erkelens, J. Rossat-Mignod, and W. G. Stirling, Europhys. Lett. 3, 945 (1987).
  * [5] M. Nishi, O. Fujita, and J. Akimitsu, Phys. Rev. B 50, 6508(R) (1994).
  * [6] G. Xu, C. Broholm, D. H. Reich, and M. A. Adams, Phys. Rev. Lett. 84, 4465 (2000).
  * [7] H. Kageyama, M. Nishi, N. Aso, K. Onizuka, T. Yosihama, K. Nukui, K. Kodama, K. Kakurai, and Y. Ueda, Phys. Rev. Lett. 84, 5876 (2000).
  * [8] A. Oosawa, T. Kato, H. Tanaka, K. Kakurai, M. Müller, and H.-J. Mikeska, Phys. Rev. B 65, 094426 (2002).
  * [9] R. S. Eccleston, M. Uehara, J. Akimitsu, H. Eisaki, N. Motoyama, and S.-i. Uchida, Phys. Rev. Lett. 81, 1702 (1998).
  * [10] M. Matsuda, K. Katsumata, R. S. Eccleston, S. Brehmer, and H.-J. Mikeska, Phys. Rev. B 62, 8903 (2000).
  * [11] S. Notbohm, P. Ribeiro, B. Lake, D. A. Tennant, K. P. Schmidt, G. S. Uhrig, C. Hess, R. Klingeler, G. Behr, B. Büchner, M. Reehuis, R. I. Bewley, C. D. Frost, P. Manuel, and R. S. Eccleston, Phys. Rev. Lett. 98, 027403 (2007).
  * [12] S. Ward, M. Mena, P. Bouillot, C. Kollath, T. Giamarchi, K. P. Schmidt, B. Normand, K. W. Krämer, D. Biner, R. Bewley, T. Guidi, M. Boehm, D. F. McMorrow, and Ch. Rüegg, Phys. Rev. Lett. 118, 177202 (2017).
  * [13] A. Gozar, G. Blumberg, B. S. Dennis, B. S. Shastry, N. Motoyama, H. Eisaki, and S. Uchida, Phys. Rev. Lett. 87, 197202 (2001).
  * [14] K. P. Schmidt, A. Gössling, U. Kuhlmann, C. Thomsen, A. Löffert, C. Gross, and W. Assmus, Phys. Rev. B 72, 094419 (2005).
  * [15] M. Windt, M. Grüninger, T. Nunner, C. Knetter, K. P. Schmidt, G. S. Uhrig, T. Kopp, A. Freimuth, U. Ammerahl, B. Büchner, and A. Revcolevschi, Phys. Rev. Lett. 87, 127002 (2001).
  * [16] J. Schlappa, T. Schmitt, F. Vernay, V. N. Strocov, V. Ilakovac, B. Thielemann, H. M. Rønnow, S. Vanishri, A. Piazzalunga, X. Wang, L. Braicovich, G. Ghiringhelli, C. Marin, J. Mesot, B. Delley, and L. Patthey, Phys. Rev. Lett. 103, 047401 (2009).
  * [17] Y. Tseng, E. Paris, K. P. Schmidt, W. Zhang, T. C. Asmara, R. Bag, V. N. Strocov, S. Singh, J. Schlappa, H. M. Rønnow, and T. Schmitt, Commun. Phys. 6, 138 (2023).
  * [18] W. Zheng, C. J. Hamer, R. R. P. Singh, S. Trebst, and H. Monien, Phys. Rev. B 63, 144410 (2001).
  * [19] K. P. Schmidt and G. S. Uhrig, Mod. Phys. Lett. B 19, 1179 (2005).
  * [20] T. S. Nunner, P. Brune, T. Kopp, M. Windt, and M. Grüninger, Phys. Rev. B 66, 180404(R) (2002).
  * [21] A. Nocera, N. D. Patel, J. Fernandez-Baca, E. Dagotto, and G. Alvarez, Phys. Rev. B 94, 205145 (2016).
  * [22] T. Nagao and J.-i. Igarashi, Phys. Rev. B 85, 224436 (2012).
  * [23] U. Kumar, A. Nocera, E. Dagotto, and S. Johnston, Phys. Rev. B 99, 205130 (2019).
  * [24] R. L. Dally, A. J. R. Heng, A. Keselman, M. M. Bordelon, M. B. Stone, L. Balents, and S. D. Wilson, Phys. Rev. Lett. 124, 197203 (2020).
  * [25] X. Bai, S.-S. Zhang, H. Zhang, Z. Dun, W. A. Phelan, V. O. Garlea, M. Mourigal, C. D. Batista, Nat. Commun. 14, 4199 (2023).
  * [26] Z. Wang, J. Wu, W. Yang, A. K. Bera, D. Kamenskyi, A. T. M. N. Islam, S. Xu, J. M. Law, B. Lake, C. Wu, and A. Loid, Nature 554, 219 (2018).
  * [27] A. Legros, S.-S. Zhang, X. Bai, H. Zhang, Z. Dun, W. A. Phelan, C. D. Batista, M. Mourigal, and N. P. Armitage, Phys. Rev. Lett. 127, 267201 (2021).
  * [28] F. Le Mardelé, A. El Mendili, M. E. Zhitomirsky, I. Mohelsky, D. Jana, I. Plutnarova, Z. Sofer, C. Faugeras, M. Potemski, and M. Orlita, arXiv:2401.05557.
  * [29] M. Enderle, B. Fåk, H.-J. Mikeska, R. K. Kremer, A. Prokofiev, and W. Assmus, Phys. Rev. Lett. 104, 237207 (2010).
  * [30] J. Schlappa, U. Kumar, K.J. Zhou, S. Singh, M. Mourigal, V.N. Strocov, A. Revcolevschi, L. Patthey, H.M. Rønnow, S. Johnston, and T. Schmitt, Nat. Commun. 9, 5394 (2018).
  * [31] R. W. Smaha, J. Pelliciari, I. Jarrige, V. Bisogni, A. T. Breidenbach, J. M. Jiang, J. Wen, H.-C. Jiang, and Y. S. Lee, Phys. Rev. B 107, L060402 (2023).
  * [32] G. Schmiedinghoff, L. Müller, U. Kumar, G. S. Uhrig, and B. Fauseweh, Commun. Phys. 5, 218 (2022).
  * [33] Y. Komaki, Y. Iwase, S. Yanagimatsu, Y. Muta, N. Maeshima, and K. Hino, J. Phys. Soc. Jpn. 88, 054709 (2019).
  * [34] S. Sachdev and R. N. Bhatt, Phys. Rev. B 41, 9323 (1990).
  * [35] S. Gopalan, T. M. Rice, and M. Sigrist, Phys. Rev. B 49, 8901 (1994).
  * [36] B. Normand and C. Rüegg, Phys. Rev. B 83, 054415 (2011).
  * [37] H. Katsura, M. Sato, T. Furuta, and N. Nagaosa, Phys. Rev. Lett. 103, 177402 (2009).
  * [38] K. Yokoi. N. Maeshima, and K. Hino, J. Phys. Soc. Jpn. 86, 104708 (2017).
  * [39] J. Oitmaa, R. R. P. Singh, and Z. Weihong, Phys. Rev. B 54, 1009 (1996).
  * [40] J. B. Torrance, J. E. Vazquez, J. J. Mayerle, and V. Y. Lee, Phys. Rev. Lett. 46, 253 (1981).
  * [41] Y. Tokura, S. Koshihara, Y. Iwasa, H. Okamoto, T. Komatsu, T. Koda, N. Iwasawa, and G. Saito, Phys. Rev. Lett. 63, 2405 (1989).
  * [42] S. Horiuchi, Y. Okimoto, R. Kumai, and Y. Tokura, J. Am. Chem. Soc. 123, 665 (2001).
  * [43] M. Messer, R. Desbuquois, T. Uehlinger, G. Jotzu, S. Huber, D. Greif, and T. Esslinger, Phys. Rev. Lett. 115, 115303 (2015).
  * [44] T. L. Yang, P. Griśins, Y. T. Chang, Z. H. Zhao, C. Y. Shih, T. Giamarchi, and R. G. Hulet, Phys. Rev. Lett. 121, 103001 (2018).
  * [45] R. Anderson, F. Wang, P. Xu, V. Venu, S. Trotzky, F. Chevy, and J. H. Thywissen, Phys. Rev. Lett. 122, 153602 (2019).
  * [46] R. A. L. Silva, A. I. S. Neves, E. B. Lopes, I. C. Santos, J. T. Coutinho, L. C. J. Pereira, C. Rovira, M. Almeida, and D. Belo, Inorg. Chem. 52, 5300 (2013).
  * [47] R. A. L. Silva and M. Almeida, J. Mater. Chem. C 9, 10573 (2021).
  * [48] T. Sakaki, Crystals 2, 374 (2012).
  * [49] P. Bouillot, C. Kollath, A. M. Läuchli, M. Zvonarev, B. Thielemann, C. Rüegg, E. Orignac, R. Citro, M. Klanjšek, C. Berthier, M. Horvatić, and T. Giamarchi, Phys. Rev. B 83, 054407 (2011).