J. Phys. Soc. Jpn.
Downloaded from journals.jps.jp by 157.1.136.5 on 10/10/24

Journal of the Physical Society of Japan 90, 044713 (2021)
https://doi.org/10.7566/JPSJ.90.044713

Mechanism for Synchronization of Charge Oscillations in Dimer Lattices
Kenji Yonemitsu1,2+ and Philipp Werner2
2

1
Department of Physics, Chuo University, Bunkyo, Tokyo 112-8551, Japan
Department of Physics, University of Fribourg, 1700 Fribourg, Switzerland(Received January 11, 2021; accepted March 5, 2021; published online March 30, 2021)

We discuss the mechanism and the conditions for the appearance of synchronized charge oscillations which have been
observed experimentally and theoretically after strong photoexcitation of dimerized systems. In the Hubbard model withon-site repulsion, the Bloch equations for a wave-number-dependent pseudospin — whose components describe thecharge-density diﬀerence, current density, and bond density between the two sublattices — involve an alternatingly tilted
pseudomagnetic ﬁeld, which assists the synchronization of pseudospins with diﬀerent wave numbers, irrespective of theinitial condition. This fact is numerically conﬁrmed by the dynamics in ﬁnite lattices based on the exact diagonalization
method. In the presence of nearest-neighbor repulsion, however, the synchronization can be hindered by excitons.
Therefore, the excitation of a suﬃciently large density of free electron–hole pairs, but low density of excitons, is neededto achieve synchronization.

1.

Introduction

Among various nonequilibrium phenomena, the manyelectron dynamics in far-from-equilibrium lattice systems has
attracted much attention.1–7) Many interesting eﬀects, including photoinduced phase transitions, are nonlinear and caused
by interactions.8–31) With progress in the experimental
techniques, ultrafast phenomena that necessarily involve
high-energy processes become increasingly important. Whichelectronic phase is realized in a correlated electron material isoften determined by a subtle competition between the kinetic
and interaction terms in the Hamiltonian. Thus, low-energy
properties are linked to high-energy processes, which makes
ultrafast control possible in some cases.
During photoexcitation, high harmonic generation isobserved in solids and reﬂects properties of their electronic
states.32–40) Recently, second harmonic generation has been
observed in a centrosymmetric organic superconductor
κ-(bis[ethylenedithio]tetrathiafulvalene)2Cu[N(CN)2]Br [κ(BEDT-TTF)2Cu[N(CN)2]Br] through a nonlinear petahertz
current before substantial scattering processes occur.41) In thiscompound,42–51) which is known to have a dimer lattice,stimulated emission is observed after strong photoexcitation.52) This stimulated emission is also caused by a nonlinear
charge oscillation — an electronic breathing mode — as has
been shown by numerical calculations based on the exact
diagonalization method.53) When the photoexcitation is weak,i.e., when the optical ﬁeld amplitude is small, diﬀerent charge
oscillations appear whose frequencies correspond to the peakenergies in the optical conductivity spectrum. When photoexcitation is strong, however, these charge oscillations
synchronize to produce the electronic breathing mode. The
synchronization can be demonstrated numerically by introducing randomness into the transfer integrals to suppress it
and showing that suﬃciently strong on-site repulsion overcomes the randomness and recovers the eﬀect.54) It has also
been discussed that strong on-site attraction produces a pairanalog of the electronic breathing mode.55)
A widely studied lattice with two atoms in a unit cell is thehoneycomb lattice, on which a synchronization transition andresultant coherent oscillations have been reported to occur ina recent mean-ﬁeld investigation.39) Here, it was pointed outthat the equations of motion are somewhat similar to thosein the Kuramoto model,56–58) where the presence of asynchronization transition is established. This analogy ishighly suggestive, but it implies that whether or not the inphase synchronization is realized depends on the sign ofthe interaction. Hence, the detailed form of the equationsof motion merits further investigations. In the present work,
we will adopt a pseudospin representation and show thatthe synchronization is indeed sensitive to the sign of theinteraction, a conclusion that will be conﬁrmed numericallyusing the exact diagonalization method. We hope that these
insights into the mechanism and conditions for thesynchronization of charge oscillations will be helpful for
future experiments related to ultrafast control of charge
motion.
Within a mean-ﬁeld approximation, charge-oscillation
dynamics is described by the Bloch equations for pseudospins in a similar manner to those for photoinduced
magnetization59) and pairing60) dynamics, or photoinduced
excitonic condensation dynamics,25) which allows us tointuitively grasp the essence of the collective evolution.
The time-dependent BCS pairing problem is known to be
integrable,60) and its similarity to and diﬀerence from thepresent problem needs to be clariﬁed. On short timescales,
the mean-ﬁeld picture basically holds even when electron
correlations are taken into account.31,59)
The numerical calculations in this paper consider onedimensional dimerized lattices. They are bipartite and consistof two sublattices. However, the synchronization mechanism
described by the Bloch equations is independent of thedimensionality of the system, and the synchronization
phenomena have been numerically observed in one- andtwo-dimensional systems.53–55) As representative systems, we
treat the Hubbard and a spinless fermion model. The former
represents a model with repulsive interactions within asublattice, while the latter possesses repulsive interactions
between the sublattices, and their synchronization conditions
are diﬀerent. To show the sensitivity to the initial condition in
the latter model, we prepare initial nonequilibrium states in
diﬀerent ways, by photoexcitation and quenching.
2.

Pseudospin Representation

To describe collective dynamics, the pseudospin representation is often useful. Here, following Ref. 39, we employ

044713-1

Societyof Japan©2021 The Physical©2021
The Author(s)This article is published by the Physical Society of Japan under the terms of the Creative Commons Attribution 4.0 License. Any further distribution of this work must
maintain attribution to the author(s) and the title of the article, journal citation, and DOI.

J. Phys. Soc. Jpn.
Downloaded from journals.jps.jp by 157.1.136.5 on 10/10/24

J. Phys. Soc. Jpn. 90, 044713 (2021)

K. Yonemitsu and P. Wernera mean-ﬁeld approximation to treat charge-oscillation
dynamics in this representation and investigate the mechanism for synchronization. We consider a Hubbard model withon-site repulsion U and a spinless fermion model withnearest-neighbor repulsion V on dimer lattices. Although thediscussion in this section does not depend on the dimensionality of the system, we use one-dimensional lattices withalternating transfer integrals, t1 and t2 , and periodic boundary
conditions when the kinetic term is explicitly shown, to makethe relationship with the results in the next section clear. The
distance between neighboring sites is set to be equal andunity when the system is photoexcited.53) The Hubbard
model can be written as
N1 X
X
HU ¼
½t1 ðeyj oj þ oyj ej Þj¼0þ t2 ðoyj ejþ1 þ eyjþ1 oj Þþ Uðeyj" ej" eyj# ej# þ oyj" oj" oyj# oj# Þ;ð1Þ

while the Hamiltonian of the spinless fermion model isN1X½t1 ðeyj oj þ oyj ej Þ þ t2 ðoyj ejþ1 þ eyjþ1 oj Þ
HV ¼
j¼0

þ Vðeyj ej oyj oj þ oyj oj eyjþ1 ejþ1 Þ;ð2Þ

where eyj (eyj ) creates an electron with spin σ (a fermion) at
site j of the “even” sublattice and oyj (oyj ) creates an electronwith spin σ (a fermion) at site j of the “odd” sublattice.
In what follows in this section, we mainly discuss the
Hubbard model and its charge-oscillation dynamics. Whennontrivial diﬀerences between the two models exist, they
will be pointed out. We use a mean-ﬁeld approximation
corresponding to the Hartree approximation, which allows
us to capture the essence of pinteraction
ﬃﬃﬃﬃ P ikjeﬀects. Via the
Fouriertransformse¼
ð1=
NÞ j e ej and ok ¼
k
pﬃﬃﬃﬃ P
ð1= NÞ j eikj oj with N denoting the number of dimers,
the mean-ﬁeld Hamiltonian can be written as
!
!
X
ek
hðkÞ
 ðU=2Þn
MF
y
y
HU ¼
; ð3Þ
ek okðU=2Þnokh ðkÞ
kwithn ¼ heyj ej  oyj oj i;for the Hubbard model and asX 2Vn
HVMF ¼eyk oykh ðkÞ
k

hðkÞ

!

2Vn

ek
ok

ð6Þ

hðkÞ ¼ t1 þ t2 eik :

1

ð7Þ
!

As in Ref. 60, we deﬁne the pseudospin components as
r1k ¼ hayk bk þ byk ak i;

ð11Þ

r2k ¼ hiayk bk þ ibyk ak i;

ð12Þ

r3k ¼ hayk ak  byk bk i:

ð13Þ

2
2
2
þ r2k
þ r3k
¼1
r1k

isHere the normalization condition
always satisﬁed for k for which one state is initially2occupied and the other is initially unoccupied, and r1kþ22
r2k þ r3k ¼ 0 is satisﬁed otherwise. In equilibrium withn ¼ 0 and h0 ðkÞ < 0 (i.e., t1 þ t2 < 0), r3k ¼ 1 for k inthe former case. Further, we deﬁne Ω60) as  ¼ ðU=2Þn.
The self-consistency condition leads to
U X
¼
r1k :
ð14Þ
4N k
The pseudospin dynamics is described by the Blochequations,r_1k ¼ 2h0 ðkÞr2k ;ð15Þ
r_2k ¼ 2h0 ðkÞr1k  2r3k ;

ð16Þ

ð17Þ
r_3k ¼ 2r2k ;which are consistent with the normalization condition. For
U > 0, there is no static solution with  ≠ 0.
For P
the spinless fermion model with   2Vn ¼
ðV=NÞ k r1k , the Bloch equations can be obtained byreplacing Ω by  in Eqs. (15)–(17). For V larger than acritical value Vc , there is a static solution with  ≠ 0, and the
dynamical problem becomes equivalent to the time-dependent BCS pairing problem.60)
Equations (15)–(17) describe the Larmor precession of thepseudospin
r_ k ðtÞ ¼ Bk ðtÞ  rk ðtÞ;Bk ðtÞ ¼ t ðB1k ðtÞ; B2k ; B3k Þ ¼ t ð2ðtÞ; 0; 2h0 ðkÞÞ;ð5Þ

;for the spinless fermion model, where constant terms havebeen omitted. On the one-dimensional lattice, the oﬀdiagonal element is given by

1

ð10Þ

ð18Þ

in a pseudo magnetic ﬁeld!

1 yhe ej  oyj oj i;
2 jFollowing Ref. 39, we deﬁne
 y
 
 1
1
ek oyk ¼ ayk byk pﬃﬃﬃ
2 1

with h0 ðkÞ ¼ sgn Reðhð0ÞÞjhðkÞj and n in Eq. (4) as
1 X y
ha bk þ byk ak i:n ¼
2N k k

ð4Þ

with
n ¼

with k ¼ tan1 ½ImðhðkÞÞ=ReðhðkÞÞ and rewrite the meanﬁeld Hamiltonian in Eq. (3) as
!
!
X
ak
 h0 ðkÞ ðU=2Þn
MF
y
y
HU ¼
; ð9Þ
ak bkbk
ðU=2Þnh0 ðkÞkeik =200eik =2

!
;
ð8Þ

ð19Þ

where the time dependence is made explicit here, andB3k ¼ 2h0 ðkÞ > 0 for h0 ðkÞ < 0. The x, y, and z componentsof the pseudospin correspond to the charge-density diﬀerence, current density, and bond density between the “even”
and “odd” sublattices. The rate of rotation around the z axis
depends on the wave number. In what follows, we assume
that U > 0. The rotation axis is tilted by ðtÞ whose sign isthat of nðtÞ, i.e., that of the majority of r1k ðtÞ. When the
majority of r1k ðtÞ take positive (negative) values, ðtÞ ispositive (negative), and the rotation axis is tilted to have apositive (negative) x component, as shown in the left (right)
panel of Fig. 1. By tilting the rotation axis in this way, mostof the pseudospins are located on the same side of therotation axis, those near the bottom t ð0; 0; 1Þ are accelerated, and consequently their dynamics become similar,

044713-2

Societyof Japan©2021 The Physical©2021
The Author(s)J. Phys. Soc. Jpn.
Downloaded from journals.jps.jp by 157.1.136.5 on 10/10/24

J. Phys. Soc. Jpn. 90, 044713 (2021)

K. Yonemitsu and P. WernerFig. 1. (Color online) Precession of a pseudospin when the majority ofthe pseudospins have positive x components (left panel) and negative xcomponents (right panel) for U > 0.

which assists the synchronization of charge oscillations.
These considerations suggest that charge oscillations aresynchronized regardless of the initial distribution of r1k .
On the right hand side of Eq. (16), the sign of h0 ðkÞ andthat of r3k are the same, as is obvious from Eqs. (9) and (13),at least near equilibrium and the sign of the majority of r1k ðtÞand the sign of ðtÞ are the same, as already pointed out.
Thus, a U-driven [see Eq. (14)] but k-independent force
component is applied in the direction of the majority of the
k-dependent current ﬂow, which is equivalent to the above
picture based on the tilted rotation axis.
A solution to the Bloch equations can be obtained,
according to Ref. 60, by the ansatz
r1k ¼ Ak ;_
r2k ¼ Bk ;

r3k ¼ Ck 2  Dk :

ð20Þ

The Bloch equations and the normalization condition lead to
Ak ¼ 2h0 ðkÞBk , Bk ¼ Ck , and
2DkD2  12
4
0
2_
 þ  þ ð2h ðkÞÞ 
2 þ k 2 ¼ 0: ð21ÞCk
Ck
If the above equation is independent of k, and of the form
_ 2 þ 4 þ ð2  2 Þ2  2 2 ¼ 0;


þ
 þ

ð22Þ

it has a solution described by the elliptic cosine function,!
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃþðtÞ ¼ þ cn
2þ þ 2 t; pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ :ð23Þ
2þ þ 2When
we ﬃ set the amplitude þ ¼  and the frequency
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
p0
2þ þ 2 ¼ 2h
ð0Þ for ﬃ t1 ; t2 < 0,53) the elliptic modulus
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
2
becomes þ = þ þ 2 ¼ =ð2h0 ð0ÞÞ. In order for
Eq. (21) to become Eq. (22), the k-dependent factors mustsatisfy1

Ck ¼ 2½ðð2h0 ð0ÞÞ2  ð2h0 ðkÞÞ2 Þ2 þ 4ð2h0 ðkÞÞ2 2  2 ; ð24Þ
Ck
:
ð25Þ
2
The negative sign in Eq. (24) is due to the self-consistency
condition,
P Eq. (14), which can be rewritten as 1=U ¼
1=ð4NÞ 0k 2h0 ðkÞCk , where the summation is over the k
for which the pseudospin has a nonzero magnitude. In thelimit of vanishing oscillation amplitude  ! 0, Dk ! 1, andr3k ! 1, i.e., the solution indeed approaches the equilibrium state. Although the initial condition must satisfy
Eqs. (24) and (25) for the pseudospin dynamics to be exactly
described by the elliptic cosine function and this is not
generally achieved by a photoexcitation, this solution
strongly suggests that charge oscillations are synchronized
by a suﬃciently large U, which has indeed been observed
Dk ¼ ½ð2h0 ð0ÞÞ2  ð2h0 ðkÞÞ2  22 previously53,54) and will also be demonstrated in the next
section. The electronic breathing mode is more easily
described in real space rather than in momentum space, but
it is controlled by a set of equations analogous to Eqs. (15)–(17), as is explained in the Appendix.
For the spinless fermion model, the situation is quite
diﬀerent from the Hubbard model. For V > Vc , a static
charge ordered state is realized in equilibrium. As mentioned
before, its nonequilibrium dynamics is known60) to be
described by another elliptic function, the delta amplitude,
whose small-amplitude limit corresponds to the Higgs
amplitude mode. For 0 < V < Vc , the excitonic eﬀect stems
from the coupling among electron–hole pairs ayk bk with
diﬀerent k’s through the element 2Vn of the mean-ﬁeld
Hamiltonian. An electron and a hole are bound by V to form
an exciton, which is described by a linear combination of
ayk bk . If the initial state contains such excitons, the chargeoscillation dynamics contains a slow component and thesynchronization is diﬃcult to achieve. In the pseudospin
representation, the rotation axis is tilted in the direction
opposite to the Hubbard case. Thus, a V-driven k-independent
force component is applied in the direction opposite to thecurrent ﬂow, which allows the existence of a charge-orderedstate for V > Vc . For 0 < V < Vc , however, this k-independent force component may lead to a synchronization if theinitial state does not have a signiﬁcant excitonic component.
This fact is numerically demonstrated in the next section.
Note that a similar synchronization phenomenon ofpseudospins has been discussed in the context of the Higgs
mode in superconductors.61,62) There, the Bloch equations
were linearized with respect to the time-dependent parts,
from which a resonant precession of pseudospins was
obtained in the long-time asymptotic evolution of the orderparameter. In the present case, on the other hand, theequilibrium state has no order parameter. After the strong
photoexcitation, the system has no external ﬁeld and its initialstate is far from equilibrium. Although Eq. (24) has a factor
that appears associated with resonance, the synchronization
mechanism is diﬀerent from that in the resonant precession
discussed in Ref. 62.
3.

Dynamics from Charge Disproportionation

Employing one-dimensional lattices, Eqs. (1) and (2), withperiodic boundary conditions, we numerically solve the timedependent Schrödinger equation to investigate the chargeoscillation dynamics from far-from-equilibrium states using
the exact diagonalization method.63) For the transfer integrals,
we use t1 ¼ 0:3 and t2 ¼ 0:1, from which the frequencyof the electronic breathing mode is given by53)!osc ¼ 2ðjt1 j þ jt2 jÞ ¼ 0:8:

ð26ÞNote that this frequency is independent of U in the Hubbard
model and independent of V in the one-dimensional spinlessfermion and extended Hubbard models, although it islowered by intersite repulsion Vij in the two-dimensional
extended Hubbard model for κ-(BEDT-TTF)2Cu[N(CN)2]Br.53) In the noninteracting case, the electronic breathingmode is merely one of the optically active modes in dimer
lattices. The other modes are suppressed by U or V afterstrong photoexcitation. The time slice employed in the
numerical solutions of the time-dependent Schrödinger

044713-3Societyof Japan©2021 The Physical©2021
The Author(s)J. Phys. Soc. Jpn.
Downloaded from journals.jps.jp by 157.1.136.5 on 10/10/24

J. Phys. Soc. Jpn. 90, 044713 (2021)

K. Yonemitsu and P. Werner

equation is dt ¼ 0:02 for small U or V, while smaller dt
values are used for large U or V to ensure the conservation of
the total energy and of the norm. The number of sites isdenoted by L (L ¼ 2N). We use system sizes up to L ¼ 16 for
the Hubbard model at quarter ﬁlling and up to L ¼ 24 for the
spinless fermion model at half ﬁlling. When photoexcitations
are used to produce nonequilibrium states, the procedure is
the same as in previous studies,53,54) i.e., we use symmetric
one-cycle electric-ﬁeld pulses, AðtÞ ¼ ðcF=!c Þ½cosð!c tÞ 1ðt þ ð2=!c ÞÞðtÞ, with central frequency !c and ﬁeld
amplitude F. P
In the quench calculations, we add aP
staggeredpotential  j ðeyj ej  oyj oj Þ to Eq. (1) or  j ðeyj ej oyj oj Þ to Eq. (2) to obtain the ground state and use this state as
the initial state at t ¼ 0. The charge-oscillation dynamics isthen computed for t > 0 with  ¼ 0. Fourier spectra of thecharge density denote the absolute values of the Fouriertransforms of the time proﬁles (0 < t < 103 ) of the charge
density immediately after setting  ¼ 0 or immediately after
the photoexcitation.
Before showing numerical results, we consider the
Hubbard model in the limit of inﬁnitely large repulsion,
U ¼ 1. In this limit, it becomes equivalent to a spinlessfermion model at doubled ﬁlling without interaction,
!
!
X
ak
0
 h0 ðkÞyy
HV¼0 ¼
;
ð27Þ
ak bk
bk
0
h0 ðkÞ
kpﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
with h0 ðkÞ ¼  t12 þ t22 þ 2t1 t2 cos k, i.e., to an integrable
system; thus, independent of the frequency, the charge
oscillations do not decay. At quarter ﬁlling, (N=2) spin-upand (N=2) spin-down electrons are present. When an electron
is virtually moved through a distance of L sites, it passes over(N=2) fermions with the opposite spin. If (N=2) is odd,
the outcome is diﬀerent from the similar process in thenoninteracting spinless fermion model or in the fullypolarized Hubbard model by a factor (−1). Thus, the allowed
wave numbers are k ¼ ð2=NÞj with j being an integer foreven (N=2) and k ¼ ð=NÞð2j þ 1Þ for odd (N=2). On theother hand, the Hubbard model with U ¼ 0 is equivalent to
the noninteracting spinless fermion model at the same ﬁlling;
thus, the allowed wave numbers are always k ¼ ð2=NÞj. The
charge-density diﬀerence 2n can be rewritten as
1X y
1X y
2n ¼
hej ej  oyj oj i ¼
ha bk þ byk ak i:
ð28Þ
N j
N k k
Q
Thus, it is maximized by the state k p1ﬃﬃ2 ðayk þ byk Þj0i, which
is set to be the initial state of jðtÞi, so that
Y 1
00
pﬃﬃﬃ ðeih ðkÞt ayk þ eih ðkÞt byk Þj0i:ð29Þ
jðtÞi ¼
2k
In this state, the time evolution of 2n is given by1X
cos 2h0 ðkÞt:2nðtÞ ¼
N k

ð30Þ

Itspﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
Fourier
spectrum
has
peaks
at
!k ¼
ﬃ
2 t12 þ t22 þ 2t1 t2 cos k with the k values described above,i.e., k ¼ 0, π for U ¼ 0 and =2 for U ¼ 1 if L ¼ 4(N ¼ 2), k ¼ 0; =2;  for U ¼ 0 and 1 if L ¼ 8 (N ¼ 4),
etc.
3.1 Hubbard model at quarter ﬁllingFirst we show in Fig. 2 that the Fourier spectrum of theFig. 2. (Color online) Optical conductivity spectrum in the ground stateand Fourier spectra of the charge density from near-equilibrium ( ¼ 102 )and far-from-equilibrium ( ¼ 1) initial states in the Hubbard model with
L ¼ 16 and U ¼ 0:8. The inset illustrates the breathing mode on the
dimerized lattice.

PPcharge density,  heyj ej i or  hoyj oj i, evolving from anear-equilibrium initial state, has peaks at energies where theoptical conductivity spectrum has peaks. Here, the opticalconductivity spectrum is calculated for the ground state as inprevious studies.49) It is also shown in this ﬁgure that theFourier spectrum of the charge density evolving from thenear-equilibrium state with small charge disproportionation( ¼ 102 ) is quite diﬀerent from that evolving from a statewith large charge disproportionation ( ¼ 1). The latter isdominated by the electronic breathing mode at ! ¼ 0:8,
although this mode is signiﬁcantly broadened, reﬂecting its
short lifetime.
The Fourier spectra of the charge density evolving fromthe state with large charge disproportionation ( ¼ 1) are
shown in Figs. 3(a) to 3(c) for diﬀerent system sizes fromL ¼ 4 to 16. For very weak repulsion (U ¼ 102 ), the spectra
have discrete peaks whose energies are given by !k withk ¼ ð2=NÞj. For very strong repulsion (U ¼ 8:0), thesystems become close to the noninteracting spinless fermion
model with periodic (for even N=2) or antiperiodic (for odd
N=2) boundary conditions; thus, the spectra have discrete
peaks whose energies are given by !k with k ¼ ð2=NÞj orð=NÞð2j þ 1Þ, respectively. In both cases, the systems areclose to integrable ones, and all the charge oscillations are
long-lived, i.e., the peaks are narrow and high. In the ﬁguresfor larger systems, the corresponding peak heights are
multiplied by factors smaller than unity, which allows them
to be compared with intermediate repulsion (U ¼ 0:8) cases.
Correlation eﬀects become signiﬁcant when the repulsion
strength is intermediate. In this case, charge oscillations aresynchronized53,54) until they ﬁnally decay due to dephasing.
The resultant electronic breathing mode has a broad peaknear ! ¼ !osc . The peak position is slightly below !osc for
L ¼ 16, but it rapidly approaches !osc with increasing L, and
it is expected to be at !osc in the thermodynamic limit. As thesystem size L increases, the diﬀerence between the very
weak=strong and the intermediate repulsion cases, i.e., thediﬀerence between the nearly integrable and the strongly
correlated cases, becomes apparent in Fig. 3.
For the largest system (L ¼ 16) with intermediate
repulsion (U ¼ 0:8) we have calculated the Fourier spectra

044713-4

Societyof Japan©2021 The Physical©2021
The Author(s)J. Phys. Soc. Jpn.
Downloaded from journals.jps.jp by 157.1.136.5 on 10/10/24

J. Phys. Soc. Jpn. 90, 044713 (2021)

K. Yonemitsu and P. Werner

Δ=1ωc=0.175, F=0.5 (×5)
ωc=0.35, F=0.5 (×3)
ωc=0.7, F=0.5 (×2)

250

(a)

FT of density

200

16-site Hubbard, U=0.8
150
100
50
0
0.2

0.4

0.6

0.8

1

1.2

ω
Fig. 4. (Color online) Fourier spectra of the charge density from far-fromequilibrium initial states prepared by  ¼ 1 and by strong (F ¼ 0:5)
photoexcitations with !c ¼ 0:175, 0.35, and 0.7 in the Hubbard model with
L ¼ 16 and U ¼ 0:8.

(b)

200
Conductivity/FT of density

180

(c)

160
140

Conductivity (×10)ωc=0.7, F=0.01 (×0.2)
ωc=0.7, F=0.524-site spinless, V=0.2

120
100
80
60
40
20
0
0.2

0.4

0.6

0.8

1

1.2

ω
Fig. 5. (Color online) Optical conductivity spectrum in the ground stateand Fourier spectra of the charge density after weak (F ¼ 102 ) and strong(F ¼ 0:5) photoexcitations with !c ¼ 0:7 in the spinless fermion model with
L ¼ 24 and V ¼ 0:2.

Fig. 3. (Color online) Fourier spectra of the charge density from far-fromequilibrium ( ¼ 1) initial states in the Hubbard model with (a) L ¼ 4,(b) L ¼ 8, (c) L ¼ 16 for U ¼ 102 , 0.8, and 8.0.

with diﬀerent initial states, as shown in Fig. 4. The initialstates with large charge disproportionation are obtained either
as the ground state with large alternating site energies( ¼ 1) or by the application of one-cycle large-amplitude(F ¼ 0:5) electric-ﬁeld pulses with diﬀerent central frequencies !c . Although the detailed structures are diﬀerent, all the
spectra are similar and dominated by the electronic breathingmode at ! ¼ 0:8. Thus, the synchronization of charge
oscillations is independent of the preparation of the initialstate. It is universally observed for states with large charge
disproportionation in the Hubbard model with intermediate
repulsion.
3.2 Spinless fermion model at half ﬁlling
In the previous section, we mentioned that the excitonic

eﬀect should be taken into account in the model with nearestneighbor repulsion. This implies that the charge-oscillation
dynamics is sensitive to how the initial state is prepared.
Hence, we mainly use one-cycle electric-ﬁeld pulses with acentral frequency near the absorption edge at ! ¼ 0:8,
!c ¼ 0:7, which allow excitations of many free electron–holepairs. We show in Fig. 5 that the Fourier spectrum of thecharge density, heyj ej i or hoyj oj i, evolving from a nearequilibrium state, has peaks at energies where the opticalconductivity spectrum has peaks. Both spectra are dominatedby an excitonic peak at ! ’ 0:4 even if the near-equilibrium
initial state is prepared by a weak photoexcitation (F ¼ 102 )with !c ¼ 0:7. This ﬁgure also shows that, if the initial stateis prepared by a strong photoexcitation (F ¼ 0:5) to be far
from equilibrium, the Fourier spectrum is dominated by theelectronic breathing mode at ! ¼ 0:8, which is in contrast to
the Fourier spectrum after a weak photoexcitation.
The Fourier spectra of the charge density after strong
photoexcitations (F ¼ 0:5) are shown in Fig. 6 for diﬀerent
repulsion strengths V below the critical value Vc . For V ¼ 0,
the system is integrable, and all the charge oscillations haveinﬁnitely long lifetimes. The corresponding discrete peaks

044713-5

Societyof Japan©2021 The Physical©2021
The Author(s)J. Phys. Soc. Jpn.
Downloaded from journals.jps.jp by 157.1.136.5 on 10/10/24

J. Phys. Soc. Jpn. 90, 044713 (2021)

K. Yonemitsu and P. Werner200

V=0.0 (×0.1)
V=0.1
24-site spinless, ωc=0.7, F=0.5 V=0.2
200
V=0.3
V=0.4

180
160
FT of density

FT of density

250

150
100
50

140

24-site spinless
V=0.2

Δ=1 (×0.4)
ωc=0.175, F=0.5 (×4)
ωc=0.35, F=0.5 (×4)
ωc=0.7, F=0.5

120
100
80
60
40
20

0

0
0.2

0.4

0.6

0.8

1

1.2

0.2

ω

0.4

0.6

0.8

1

1.2

ω

Fig. 6. (Color online) Fourier spectra of the charge density after strong(F ¼ 0:5) photoexcitations with !c ¼ 0:7 in the spinless fermion model with
L ¼ 24 for V values below Vc .

Fig. 7. (Color online) Fourier spectra of the charge density from far-fromequilibrium initial states prepared by  ¼ 1 and by strong (F ¼ 0:5)
photoexcitations with !c ¼ 0:175, 0.35, and 0.7 in the spinless fermion
model with L ¼ 24 and V ¼ 0:2.

appear at !k with k ¼ ð=6Þj, and their peak heights are
multiplied by 0.1 in order for them to be comparable with thenonintegrable cases. For small V (V ¼ 0:1), some structuresin the Fourier spectrum reﬂect the peaks at V ¼ 0, but they
are signiﬁcantly broadened. For V ¼ 0:1, 0.2, and 0.3, theFourier spectra are dominated by the electronic breathingmode at ! ¼ 0:8, but its peak height decreases withincreasing V. For V ¼ 0:4, the corresponding peak appears
slightly blueshifted, although V is still smaller than Vc . In
general, the peak associated with the electronic breathingmode is considerably broadened when V approaches Vc ’2jt1 j ¼ 0:6, and it disappears in the charge-ordered phase at
V > Vc .53)
For V ¼ 0:2, we have also calculated the Fourier spectrawith diﬀerent initial states with large charge disproportionation, as shown in Fig. 7. They are sensitive to the initialcondition. When the initial state is prepared by a photoexcitation with !c ¼ 0:7, the spectrum is dominated by theelectronic breathing mode. However, when !c is lowered toallow the exciton at ! ’ 0:4 (Fig. 5) to be excited, theFourier spectra have large contributions from the exciton andsmaller contributions from the electronic breathing oscillation. When the initial sate is prepared by  ¼ 1, the Fourier
spectrum becomes more complex, presumably because
various electron–hole pairs are involved. Thus, to synchronize charge oscillations, the initial state must be strongly
photoexcited with a frequency that produces many free andfew bound electron–hole pairs.

number, while the term derived from the interaction is
proportional to a wave number independent and thus
universal factor. The latter facilitates the synchronization of
charge oscillations. In the Hubbard model with U > 0, more
precisely speaking, when the interaction within a sublattice
is repulsive, the kinetic term and the interaction term havethe same sign and constructively work together. From theviewpoint of pseudospins, the rotation axis describing theLarmor precession is alternatingly tilted in such a way thatit leads pseudospins of diﬀerent wave numbers to besynchronized. In the spinless fermion model with V > 0,
more precisely speaking, when the interaction between
the two sublattices is repulsive, the kinetic term and theinteraction term have opposite signs. This allows theexistence of static charge order for V > Vc . For V < Vc , thesynchronization of charge oscillations is possible but requires
suitable excitation protocols.
The above insights based on the analytic form of the Blochequations are all consistent with the Fourier analyses of thenumerically obtained charge-oscillation dynamics. To prepare initial states with charge disproportionation, we either
quenched a staggered potential to zero or applied one-cycle
electric-ﬁeld pulses. In the Hubbard model with intermediate
repulsion, i.e., suﬃciently away from the integrable limits of
U ¼ 0 and 1, the synchronization is achieved irrespective of
how the initial state is prepared. In the spinless fermion
model with V < Vc , a suﬃciently large number of freeelectron–hole pairs must be excited to achieve the synchronization, which is hindered by excitons and is not achievedby the quenching. These results clarify the conditions for thesynchronization of charge oscillations on dimer lattices,which may guide the study of synchronization phenomena indiﬀerent classes of lattices.

4.

Conclusions

To clarify the mechanism for the previously observed
synchronization of charge oscillations52–54) and investigate
the conditions for the appearance of this phenomenon, we
analytically studied the Bloch equations in the mean-ﬁeld
approximation and numerically studied the charge-oscillation
dynamics using the exact diagonalization method for the
Hubbard and spinless fermion models on one-dimensional
dimerized lattices.
In one of the Bloch equations, the time derivatives ofcurrent densities between the sublattices are determined by
two terms, one of which is of kinetic origin and the other isof interaction origin. The kinetic term depends on the wave

Acknowledgments This work was supported by JSPS KAKENHI Grant
No. JP16K05459, MEXT Q-LEAP Grant No. JPMXS0118067426, and JST
CREST Grant No. JPMJCR1901. PW acknowledges support from ERC
Consolidator Grant No. 724103.

Appendix:

Breathing Mode in the Pseudospin Picture

In this appendix we show that the equations governing thecharge motion in a real-space representation are equivalent
to the pseudo-spin equations (15)–(17) and explain the044713-6

Societyof Japan©2021 The Physical©2021
The Author(s)J. Phys. Soc. Jpn.
Downloaded from journals.jps.jp by 157.1.136.5 on 10/10/24

J. Phys. Soc. Jpn. 90, 044713 (2021)

K. Yonemitsu and P. Werner

two sublattices weighted by transfer integrals. In the last two
terms, we use the mean-ﬁeld approximation to replace theone-body operators with spin  ¼  by their expectation
values. Because the sites ði; jÞ and ðl; mÞ [ða; bÞ and ðc; dÞ]belong to the same sublattice, the last two terms vanish owing
to the translational symmetry. Thus, Eq. (A·3) is equivalent
to Eq. (15) for the charge-density diﬀerence between thesublattices. Using Eq. (A·2), we also obtaindðicyij clm þ icylm cij ÞdtXtll0 ;m ðcyij cl0 m þ cyl0 m cij Þ¼
l0

Fig. A·1. (Color online) Notation used for the transfer integrals betweensites on the square lattice.



X
m0

þ
synchronization leading to the observed breathing mode.
Although the following discussion does not depend on thedimensionality of the system, we consider the electronicbreathing mode in the square lattice shown in Fig. A·1. The
Hubbard model can be written as
X
X
H2D ¼
til; j cyij clj þti;jm cyij cim
ilj

ijm

X y
þU
cij" cij" cyij# cij# ;

ðA:1Þ

ij

i0

þ

ðA:2Þ

Using Eq. (A·2), we easily obtaind yðc clm  cyab ccd Þ
dt ij
!
XX
yytll0 ;m cij cl0 m þtl;mm0 cij clm0 
¼ i
m0

l0

i

X
a0

þi

X
i0

þi

X
c0

taa0 ;b cya0 b ccd þ
tii0 ; j cyi0 j clm þ

X
b0

X

tcc0 ;d cyab cc0 d þj0!

ti;jj 0 cyij 0  clmXd0!
ta;bb0 cyab0  ccd!
tc;dd0 cyab ccd0  iUðcylm clm  cyij cij Þcyij clmþ iUðcycd ccd  cyab cab Þcyab ccd :ðA:3Þ

When we assume that the sites ði; jÞ and ðl; mÞ belong to theeven sublattice and the sites ða; bÞ and ðc; dÞ belong to the
odd sublattice, the sites ða0 ; bÞ, ða; b0 Þ, ðc0 ; dÞ, and ðc; d 0 Þbelong to the even sublattice, and the sites ði0 ; jÞ, ði; j 0 Þ,
ðl 0 ; mÞ, and ðl; m0 Þ belong to the odd sublattice. Thus, the ﬁrst
four terms can be regarded as current densities between thetii0 ; j ðcyi0 j clm þ cylm ci0 j Þti;jj 0 ðcyij 0  clm þ cylm cij 0  Þþ Uðcyij cij  cylm clm Þðcyij clm þ cylm cij Þ;

ðA:4Þ

and
d y
ðc clm þ cylm cij Þ
dt ij
X
tll0 ;m ðicyij cl0 m þ icyl0 m cij Þ
¼
l0

þ

X
m0



X
i0



X
j0

j0

 iUðcylm clm  cyij cij Þcyij clm :X
j0

where cyij creates an electron with spin σ at site ði; jÞ. In the
Heisenberg picture, the time evolution of the one-bodyoperator cyij clm is given byXX
d ycij clm ¼ i
tll0 ;m cyij cl0 m  itl;mm0 cyij clm0 
dt
m0
l0
X
X
y
þi
tii0 ;j ci0 j clm þ iti;jj 0 cyij 0  clmi0Xtl;mm0 ðcyij clm0  þ cylm0  cij Þtl;mm0 ðicyij clm0  þ icylm0  cij Þ
tii0 ;j ðicyi0 j clm þ icylm ci0 j Þti;jj 0 ðicyij 0  clm þ icylm cij 0  Þ Uðcyij cij  cylm clm Þðicyij clm þ icylm cij Þ: ðA:5Þ
When we assume that the site ði; jÞ belongs to the evensublattice and the site ðl; mÞ belongs to the odd sublattice, thesites ðl 0 ; mÞ and ðl; m0 Þ belong to the even sublattice, and thesites ði0 ; jÞ and ði; j 0 Þ belong to the odd sublattice. Then, theﬁrst and second terms in Eqs. (A·4) and (A·5) are one-bodyoperators acting within the even sublattice, and the third andfourth terms are ones acting within the odd sublattice. The
ﬁrst four terms in Eq. (A·4) can be regarded as chargedensity diﬀerences between the two sublattices weighted
by transfer integrals. In the mean-ﬁeld approximation,Uðcyij cij  cylm clm Þ
isreplacedbyUhcyij cij yclm clm i ¼ Un ¼ 2, which is independent of ði; jÞ, ðl; mÞ, The last term in Eq. (A·4) is regarded as (2) timesor .
the bond density between sites ði; jÞ and ðl; mÞ. Thus,
Eq. (A·4) is equivalent to Eq. (16) for the current density
between the sublattices. The ﬁrst and second (third andfourth) terms in Eq. (A·5) can be regarded as current
densities within the even (odd) sublattice weighted bytransfer integrals, which give a small contribution and areunimportant. [In momentum space, they can be eliminated by
the unitary transformation within each sublattice using eik =2
in Eq. (8).] The last term in Eq. (A·5) can be regarded as(2) times the current density between sites ði; jÞ and ðl; mÞ.
Thus, Eq. (A·5) is equivalent to Eq. (17) [though both sides
are multiplied by (−1)] for the bond density between thesublattices.

044713-7

Societyof Japan©2021 The Physical©2021
The Author(s)J. Phys. Soc. Jpn.
Downloaded from journals.jps.jp by 157.1.136.5 on 10/10/24

J. Phys. Soc. Jpn. 90, 044713 (2021)

K. Yonemitsu and P. Werner

It is now clear that a common factor Un (¼ 2) appears
independently of the relative position of sites ði; jÞ and ðl; mÞ.
When U is positive, the force in the last term of Eq. (A·4)
is applied in the direction of the current ﬂow. [The current
direction is easily checked by setting ðl; mÞ ¼ ði; jÞ in
Eq. (A·2).] Thus, the charge-density diﬀerence between thesublattices synchronizes the charge motion on diﬀerent bonds
to achieve the electronic breathing mode.

+

kxy@phys.chuo-u.ac.jp
1) D. N. Basov, R. D. Averitt, D. van der Marel, M. Dressel, and K.
Haule, Rev. Mod. Phys. 83, 471 (2011).
2) D. Nicoletti and A. Cavalleri, Adv. Opt. Photonics 8, 401 (2016).
3) C. Giannetti, M. Capone, D. Fausti, M. Fabrizio, F. Parmigiani, and D.
Mihailovic, Adv. Phys. 65, 58 (2016).
4) S. Kaiser, Phys. Scr. 92, 103001 (2017).
5) M. Chergui and E. Collet, Chem. Rev. 117, 11025 (2017).
6) Y. Kawakami, H. Itoh, K. Yonemitsu, and S. Iwai, J. Phys. B 51,
174005 (2018).
7) S. Ishihara, J. Phys. Soc. Jpn. 88, 072001 (2019).
8) S. Iwai, M. Ono, A. Maeda, H. Matsuzaki, H. Kishida, H. Okamoto,and Y. Tokura, Phys. Rev. Lett. 91, 057401 (2003).
9) D. Fausti, R. I. Tobey, N. Dean, S. Kaiser, A. Dienst, M. C. Hoﬀmann,
S. Pyon, T. Takayama, H. Takagi, and A. Cavalleri, Science 331, 189(2011).
10) S. Wall, D. Brida, S. R. Clark, H. P. Ehrke, D. Jaksch, A. Ardavan, S.
Bonora, H. Uemura, Y. Takahashi, T. Hasegawa, H. Okamoto, G.
Cerullo, and A. Cavalleri, Nat. Phys. 7, 114 (2011).
11) N. Tsuji, T. Oka, P. Werner, and H. Aoki, Phys. Rev. Lett. 106, 236401(2011).
12) N. Tsuji, T. Oka, H. Aoki, and P. Werner, Phys. Rev. B 85, 155124(2012).
13) M. Gao, C. Lu, H. Jean-Ruel, L. C. Liu, A. Marx, K. Onda, S.
Koshihara, Y. Nakano, X. Shao, T. Hiramatsu, G. Saito, H. Yamochi,
R. R. Cooney, G. Moriena, G. Sciaini, and R. J. D. Miller, Nature 496,
343 (2013).
14) J. Ohara, Y. Kanamori, and S. Ishihara, Phys. Rev. B 88, 085107(2013).
15) T. Ishikawa, Y. Sagae, Y. Naitoh, Y. Kawakami, H. Itoh, K.
Yamamoto, K. Yakushi, H. Kishida, T. Sasaki, S. Ishihara, Y.
Tanaka, K. Yonemitsu, and S. Iwai, Nat. Commun. 5, 5528 (2014).
16) S. Kaiser, C. R. Hunt, D. Nicoletti, W. Hu, I. Gierz, H. Y. Liu, M. L.
Tacon, T. Loew, D. Haug, B. Keimer, and A. Cavalleri, Phys. Rev. B
89, 184516 (2014).
17) L. Stojchevska, I. Vaskivskyi, T. Mertelj, P. Kusar, D. Svetin, S.
Brazovskii, and D. Mihailovic, Science 344, 177 (2014).
18) R. Fukaya, Y. Okimoto, M. Kunitomo, K. Onda, T. Ishikawa, S.
Koshihara, H. Hashimoto, S. Ishihara, A. Isayama, H. Yui, and T.
Sasagawa, Nat. Commun. 6, 8519 (2015).
19) T.-R. T. Han, F. Zhou, C. D. Malliakas, P. M. Duxbury, S. D. Mahanti,
M. G. Kanatzidis, and C.-Y. Ruan, Sci. Adv. 1, e1400173 (2015).
20) M. Mitrano, A. Cantaluppi, D. Nicoletti, S. Kaiser, A. Perucchi, S.
Lupi, P. D. Pietro, D. Pontiroli, M. Riccò, S. R. Clark, D. Jaksch, and
A. Cavalleri, Nature 530, 461 (2016).
21) L. Rettig, R. Cortés, J.-H. Chu, I. Fisher, F. Schmitt, R. Moore, Z.-X.
Shen, P. Kirchmann, M. Wolf, and U. Bovensiepen, Nat. Commun. 7,10459 (2016).
22) A. Singer, S. K. K. Patel, R. Kukreja, V. Uhlíř, J. Wingert, S. Festersen,
D. Zhu, J. M. Glownia, H. T. Lemke, S. Nelson, M. Kozina, K.
Rossnagel, M. Bauer, B. M. Murphy, O. M. Magnussen, E. E.
Fullerton, and O. G. Shpyrko, Phys. Rev. Lett. 117, 056401 (2016).
23) S. Mor, M. Herzog, D. Golež, P. Werner, M. Eckstein, N. Katayama,
M. Nohara, H. Takagi, T. Mizokawa, C. Monney, and J. Stähler, Phys.
Rev. Lett. 119, 086401 (2017).
24) A. Ono and S. Ishihara, Phys. Rev. Lett. 119, 207202 (2017).
25) Y. Murakami, D. Golež, M. Eckstein, and P. Werner, Phys. Rev. Lett.
119, 247601 (2017).

26) Y. Tanaka, M. Daira, and K. Yonemitsu, Phys. Rev. B 97, 115105(2018).
27) Y. Kinoshita, N. Kida, Y. Magasaki, T. Morimoto, T. Terashige, T.
Miyamoto, and H. Okamoto, Phys. Rev. Lett. 124, 057402 (2020).
28) K. Gillmeister, D. Golež, C.-T. Chiang, N. Bittner, Y. Pavlyukh, J.
Berakdar, P. Werner, and W. Widdra, Nat. Commun. 11, 4095 (2020).
29) M. Schüler, U. D. Giovannini, H. Hübener, A. Rubio, M. A. Sentef,
T. P. Devereaux, and P. Werner, Phys. Rev. X 10, 041013 (2020).
30) N. Bittner, D. Golež, M. Eckstein, and P. Werner, Phys. Rev. B 102,
235169 (2020).
31) P. Werner and Y. Murakami, Phys. Rev. B 102, 241103(R) (2020).
32) S. Ghimire, A. D. DiChiara, E. Sistrunk, P. Agostini, L. F. DiMauro,and D. A. Reis, Nat. Phys. 7, 138 (2011).
33) O. Schubert, M. Hohenleutner, F. Langer, B. Urbanek, C. Lange, U.
Huttner, D. Golde, T. Meier, M. Kira, S. W. Koch, and R. Huber, Nat.
Photonics 8, 119 (2014).
34) T. T. Luu, M. Garg, S. Y. Kruchinin, A. Moulet, M. T. Hassan, and E.
Goulielmakis, Nature 521, 498 (2015).
35) N. Yoshikawa, T. Tamaya, and K. Tanaka, Science 356, 736 (2017).
36) R. E. F. Silva, I. V. Blinov, A. N. Rubtsov, O. Smirnova, and M.
Ivanov, Nat. Photonics 12, 266 (2018).
37) Y. Murakami, M. Eckstein, and P. Werner, Phys. Rev. Lett. 121,057405 (2018).
38) T. N. Ikeda, K. Chinzei, and H. Tsunetsugu, Phys. Rev. A 98, 063426(2018).
39) T. Nag, R.-J. Slager, T. Higuchi, and T. Oka, Phys. Rev. B 100, 134301(2019).
40) M. Lysne, Y. Murakami, M. Schüler, and P. Werner, Phys. Rev. B 102,
081121(R) (2020).
41) Y. Kawakami, T. Amano, H. Ohashi, H. Itoh, Y. Nakamura, H.
Kishida, T. Sasaki, G. Kawaguchi, H. M. Yamamoto, K. Yamamoto, S.
Ishihara, K. Yonemitsu, and S. Iwai, Nat. Commun. 11, 4138 (2020).
42) F. Kagawa, K. Miyagawa, and K. Kanoda, Nature 436, 534 (2005).
43) Y. Kawakami, S. Iwai, T. Fukatsu, M. Miura, N. Yoneyama, T. Sasaki,
and N. Kobayashi, Phys. Rev. Lett. 103, 066403 (2009).
44) M. Abdel-Jawad, I. Terasaki, T. Sasaki, N. Yoneyama, N. Kobayashi,
Y. Uesu, and C. Hotta, Phys. Rev. B 82, 125119 (2010).
45) M. Naka and S. Ishihara, J. Phys. Soc. Jpn. 79, 063707 (2010).
46) H. Gomi, T. Imai, A. Takahashi, and M. Aihara, Phys. Rev. B 82,
035101 (2010).
47) C. Hotta, Phys. Rev. B 82, 241104(R) (2010).
48) S. Dayal, R. T. Clay, H. Li, and S. Mazumdar, Phys. Rev. B 83,
245106 (2011).
49) K. Yonemitsu, S. Miyashita, and N. Maeshima, J. Phys. Soc. Jpn. 80,
084710 (2011).
50) K. Itoh, H. Itoh, M. Naka, S. Saito, I. Hosako, N. Yoneyama, S.
Ishihara, T. Sasaki, and S. Iwai, Phys. Rev. Lett. 110, 106401 (2013).
51) H. Gomi, T. Kawatani, T. J. Inagaki, and A. Takahashi, J. Phys. Soc.
Jpn. 83, 094714 (2014).
52) Y. Kawakami, T. Amano, Y. Yoneyama, Y. Akamine, H. Itoh, G.
Kawaguchi, H. M. Yamamoto, H. Kishida, K. Itoh, T. Sasaki, S.
Ishihara, Y. Tanaka, K. Yonemitsu, and S. Iwai, Nat. Photonics 12, 474(2018).
53) K. Yonemitsu, J. Phys. Soc. Jpn. 87, 044708 (2018).
54) T. Shimada and K. Yonemitsu, J. Phys. Soc. Jpn. 89, 084701 (2020).
55) K. Yonemitsu, J. Phys. Soc. Jpn. 87, 124703 (2018).
56) Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence(Springer, Berlin, 1984).
57) J. A. Acebrón, L. L. Bonilla, C. J. P. Vicente, F. Ritort, and R. Spigler,Rev. Mod. Phys. 77, 137 (2005).
58) F. A. Rodrigues, T. K. D. Peron, P. Ji, and J. Kurths, Phys. Rep. 610, 1(2016).
59) N. Tsuji, M. Eckstein, and P. Werner, Phys. Rev. Lett. 110, 136404(2013).
60) R. A. Barankov, L. S. Levitov, and B. Z. Spivak, Phys. Rev. Lett. 93,
160401 (2004).
61) R. Matsunaga, N. Tsuji, H. Fujita, A. Sugioka, K. Makise, Y. Uzawa,
H. Terai, Z. Wang, H. Aoki, and R. Shimano, Science 345, 1145
(2014).
62) N. Tsuji and H. Aoki, Phys. Rev. B 92, 064508 (2015).
63) K. Yonemitsu and N. Maeshima, Phys. Rev. B 79, 125118 (2009).

044713-8

Societyof Japan©2021 The Physical©2021
The Author(s)