#  Polarization-dependent magnetic properties of periodically driven
$\alpha$-RuCl3

Naoya Arakawa<EMAIL_ADDRESS>The Institute of Science and
Engineering, Chuo University, Bunkyo, Tokyo 112-8551, Japan  Kenji Yonemitsu
The Institute of Science and Engineering, Chuo University, Bunkyo, Tokyo
112-8551, Japan  Department of Physics, Chuo University, Bunkyo, Tokyo
112-8551, Japan

###### Abstract

We study magnetic properties of a periodically driven Mott insulator with
strong spin-orbit coupling and show some properties characteristic of linearly
polarized light. We consider a $t_{2g}$-orbital Hubbard model driven by
circularly or linearly polarized light with strong spin-orbit coupling and
derive its effective Hamiltonian in the strong-interaction limit for a high-
frequency case. We show that linearly polarized light can change not only the
magnitudes and signs of the exchange interactions, but also their bond
anisotropy even without the bond-anisotropic hopping integrals. Because of
this property, the honeycomb-network spin system could be transformed into
weakly coupled zigzag or step spin chains for the light field polarized along
the $b$\- or $a$-axis, respectively. Then, analyzing how the light fields
affect several magnetic states in a mean-field approximation, we show that
linearly polarized light can change the relative stability of the competing
magnetic states, whereas such a change is absent for circularly polarized
light. We also analyze the effects of both the bond anisotropy of nearest-
neighbor hopping integrals and a third-neighbor hopping integral on the
magnetic states and show that the results obtained in a simple model, in which
the bond-averaged nearest-neighbor hopping integrals are considered, remain
qualitatively unchanged except for the stability of zigzag states in the non-
driven case and the degeneracy lifting of the zigzag or stripy states.

##  I Introduction

In periodically driven systems, magnetic properties can be controlled via a
time-periodic field. In the presence of a time-periodic field the solution to
the Schrödinger equation satisfies the Floquet theorem and the Floquet
Hamiltonian, a time-independent Hamiltonian, can describe the time evolution
in steps of the driving period $T$ Floq1 ; Floq2 ; review1 . Such a
description may be appropriate if the effects of heating due to the driving
field are negligible; such a situation could be realized before the system
approaches an infinite-temperature state Heat-Floq1 ; Heat-Floq2 . Since the
Floquet Hamiltonian usually depends on parameters of the driving field, it is
possible to control magnetic properties of a periodically driven system. For
example, in a single-orbital Mott insulator driven by $E(t)=E_{0}\cos\omega
t$, we can change the magnitude and sign of the antiferromagnetic Heisenberg
interaction by tuning $E_{0}$ and $\omega$ Floquet-1orbMott . This or an
extended method could be used to control magnetic properties of Mott
insulators.

If a periodically driven Mott insulator has strong spin-orbit coupling (SOC),
it may offer possibility for controlling various exchange interactions and
magnetic states. The low-energy excitations of a Mott insulator with strong
SOC can be described by the spin and orbital entangled degrees of freedom
StrongLS-Mott ; StrongLS-Mott-review ; Takagi-review  and its effective
Hamiltonian has not only the Heisenberg interaction, but also the anisotropic
exchange interactions StrongLS-Ex1 ; StrongLS-Ex2 ; Rau-PRL ; Valenti-PRB ;
NA-Jeff . Then, various magnetic states appear, depending on the values of the
Heisenberg interaction and the anisotropic exchange interactions Rau-PRL ;
Valenti-PRB ; PD1 ; PD2 ; PD3 ; PD4 ; PD5 . In our previous paper Floquet-NA
we showed that by applying a circularly polarized light field
$\bm{E}(t)={}^{t}(E_{0}\cos\omega t\ -E_{0}\sin\omega t)$ to a multiorbital
Mott insulator with strong SOC and tuning $\omega$ and $E_{0}$, the magnitudes
and signs of three exchange interactions can be changed simultaneously; these
interactions are the Heisenberg interaction $J$, the Kitaev interaction $K$,
and the off-diagonal symmetric exchange interaction $\Gamma$.

The aim of this paper is to clarify the polarization dependences of magnetic
properties for a periodically driven Mott insulator with strong SOC. First, it
is essential to understand the similarities and differences between the
effects of circularly and linearly polarized light. In general, magnetic
properties of solids depend on the polarization of light review-opt . It is
also necessary to clarify how the changes in exchange interactions due to a
light field affect energies of several magnetic states. These topics are not
discussed in our previous paper Floquet-NA .

In this paper we study the exchange interactions and the energies of several
magnetic states for a periodically driven Mott insulator with strong SOC. We
consider a $t_{2g}$-orbital Hubbard model on the honeycomb lattice with strong
SOC and a field of linearly or circularly polarized light as a model of
periodically driven $\alpha$-RuCl3. To analyze the effects of both the bond
anisotropy of nearest-neighbor hopping integrals and a third-neighbor hopping
integral, we consider five cases of our model, including a simple case
considered in our previous paper Floquet-NA ; the differences among them are
about the hopping integrals. Treating the effects of one of the light fields
as Peierls phase factors and using the Floquet theory Floquet-1orbMott ;
Floquet-MultiMott ; Floquet-NA , we derive the effective Hamiltonian of
periodically driven $\alpha$-RuCl3 in the strong-interaction limit for a high-
frequency case. Evaluating the exchange interactions numerically in the first
two cases of our model, we show that linearly polarized light can be used to
change not only the magnitudes and signs of $J$, $K$, and $\Gamma$, but also
their bond anisotropy even without the bond anisotropy of the hopping
integrals. This property, which is distinct from the effects of circularly
polarized light Floquet-NA , could be used to transform the honeycomb-network
spin system into weakly coupled zigzag or step spin chains. We also show that
the effects of the bond anisotropy of the nearest-neighbor hopping integrals
are just quantitative in the sense that it only induces weak bond anisotropy
of the exchange interactions. Then, by using a mean-field approximation (MFA),
we evaluate the expectation value of our effective Hamiltonian. Analyzing the
effects of the light fields on the energies of several magnetic states in the
five cases of our model, we show that linearly polarized light can change the
relative stability of the competing magnetic states, whereas circularly
polarized light does not. Furthermore, we show that the results obtained in
the simple case remain qualitatively unchanged except that the stability of
the zigzag states or the ferromagnetic state depends on the strength of the
bond anisotropy of the hopping integrals and the value of the third-neighbor
hopping integral and that the degeneracy of the zigzag or stripy states is
lifted not only by the fields of linearly polarized light, but also by the
bond anisotropy of the hopping integrals.

The rest of the paper is organized as follows. In Sec. II we introduce the
Hamiltonian of our model. It consists of the hopping integrals of the
$t_{2g}$-orbital electrons on the honeycomb lattice with a light field
$\bm{E}(t)$, their $LS$ coupling, and the $t_{2g}$-orbital Hubbard
interactions. There are three $\bm{E}(t)$’s considered:
$\bm{E}_{\textrm{circ}}(t)$, $\bm{E}_{\textrm{linear-}b}(t)$, and
$\bm{E}_{\textrm{linear-}a}(t)$ [Eqs. (2)–(4)]. We consider five cases of our
model, one of which is used in our previous paper Floquet-NA ; the others are
used to study the effects of the bond-anisotropic nearest-neighbor hopping
integrals and the third-neighbor hopping integral. In Sec. III A we express
the effective Hamiltonian of our periodically driven Mott insulator in terms
of the isotropic and the anisotropic exchange interactions. In Sec. III B we
present the $|u_{ij}|(=eE_{0}/\omega)$ dependences of the exchange
interactions estimated numerically in the first two cases of our model. We
compare the results for the three $\bm{E}(t)$’s and discuss the properties
characteristic of the linearly polarized light. In Sec. IV A, by applying the
MFA to our effective Hamiltonian, we derive an expression of the energy of a
magnetic state characterized by the ordering vector. We also explain the
characteristics of the magnetic states considered. In Sec. IV B we present the
$|u_{ij}|$ dependences of the numerically evaluated energies of several
magnetic states in the five cases of our model. Comparing the results for the
three $\bm{E}(t)$’s, we discuss the similarities and differences between the
effects of the fields of circularly and linearly polarized light. We also
analyze how the bond anisotropy of the nearest-neighbor hopping integrals and
the third-neighbor hopping integral affect the magnetic states. In Sec. V we
discuss the validity of our model, the effects of heating, and an experimental
observation of our results. Furthermore, we remark on a property induced by a
field of circularly polarized light and several directions for further
relevant research. In Sec. VI we summarize the main results and their
implications.

Throughout this paper, we set $\hbar=1$ and $a_{\textrm{NN}}=1$, where
$a_{\textrm{NN}}$ denotes the distance between nearest-neighbor sites on the
honeycomb lattice (Fig. 1). This choice of $a_{\textrm{NN}}$ leads to
$a_{\textrm{2nd}}=\sqrt{3}$ and $a_{\textrm{3rd}}=2$, where $a_{\textrm{2nd}}$
and $a_{\textrm{3rd}}$ denote the distances between second-neighbor and third-
neighbor, respectively, sites. For simplicity, we neglect a small difference
between the lengths of a $Z$ bond and of an $X$ or $Y$ bond.

##  II Model

Figure 1:  Structure of the honeycomb lattice. Blue, green, and red lines
represent three nearest-neighbor bonds, i.e., $X$, $Y$, and $Z$ bonds,
respectively; blue, green, and red dashed lines represent three third-neighbor
bonds, i.e., $X_{3}$, $Y_{3}$, and $Z_{3}$ bonds, respectively. Black and
white circles denote the $A$ and $B$ sublattices, respectively. This
sublattice structure is necessary because the honeycomb lattice, which is not
a Bravais lattice, can be represented as a triangular Bravais lattice with a
two-point basis Ashcro-Merm . The $a$ and $b$ axes of the crystal are also
shown.

Our model Hamiltonian consists of three parts:

$\displaystyle H=H_{\textrm{KE}}+H_{\textrm{SOC}}+H_{\textrm{int}},$ (1)

where $H_{\textrm{KE}}$ represents the kinetic energy, $H_{\textrm{SOC}}$
represents the atomic SOC StrongLS-Mott-review , and $H_{\textrm{int}}$
represents the Coulomb interactions Kanamori .

Table 1: Five cases of our model. Bond anisotropy of the nearest-neighbor hopping integrals is absent in the first and third cases and present in the other three cases. The third nearest-neighbor hopping integral is neglected in the first two cases and considered in the last three cases.  Parameter |  $1$st |  $2$nd |  $3$rd |  $4$th |  $5$th
---|---|---|---|---|---
Bond anisotropy | Absent | Present | Absent | Present | Present
$t_{\textrm{3rd}}$ (meV) | $0$ | $0$ | $-40$ | $-40$ | $-60$

The kinetic energy is given by the hopping integrals of the $t_{2g}$-orbital
electrons on the honeycomb lattice (Fig. 1) in the presence of a light field
$\bm{E}(t)={}^{t}(E_{\bar{b}}(t)\ E_{a}(t))$. Here $E_{\bar{b}}(t)$ is the
component antiparallel to the $b$ axis and $E_{a}(t)$ is that parallel to the
$a$ axis. The $\bm{E}(t)$ for circularly polarized light is given by

$\displaystyle\bm{E}_{\textrm{circ}}(t)={}^{t}(E_{0}\cos\omega t\
-E_{0}\sin\omega t),$ (2)

and that for linearly polarized light is given by

$\displaystyle\bm{E}_{\textrm{linear-}b}(t)={}^{t}(E_{0}\sin\omega t\ 0),$ (3)

or

$\displaystyle\bm{E}_{\textrm{linear-}a}(t)={}^{t}(0\ E_{0}\sin\omega t).$ (4)

We do not consider the helicity of circularly polarized light because the
magnetic properties shown in this paper for
$\bm{E}(t)=\bm{E}_{\textrm{circ}}(t)$ remain unchanged even for
$\bm{E}(t)={}^{t}(E_{0}\cos\omega t\ E_{0}\sin\omega t)$. Then, we treat the
effects of $\bm{E}(t)$ as Peierls phase factors:

$\displaystyle
H_{\textrm{KE}}=\sum_{i,j}\sum_{a,b}\sum_{\sigma}t_{iajb}e^{-ie(\bm{R}_{i}-\bm{R}_{j})\cdot\bm{A}(t)}c_{ia\sigma}^{\dagger}c_{jb\sigma},$
(5)

where $\bm{A}(t)={}^{t}(A_{\bar{b}}(t)\ A_{a}(t))$ for
$\bm{E}(t)=\bm{E}_{\textrm{circ}}(t)$, $\bm{E}_{\textrm{linear-}b}(t)$, or
$\bm{E}_{\textrm{linear-}a}(t)$ is given by

$\displaystyle\bm{A}_{\textrm{circ}}(t)={}^{t}\Bigl{(}-\frac{E_{0}}{\omega}\sin\omega
t\ -\frac{E_{0}}{\omega}\cos\omega t\Bigr{)},$ (6)
$\displaystyle\bm{A}_{\textrm{linear-}b}(t)={}^{t}\Bigl{(}\frac{E_{0}}{\omega}\cos\omega
t\ 0\Bigr{)},$ (7)

or

$\displaystyle\bm{A}_{\textrm{linear-}a}(t)={}^{t}\Bigl{(}0\
\frac{E_{0}}{\omega}\cos\omega t\Bigr{)},$ (8)

respectively. As for $t_{iajb}$’s, we consider five cases (TABLE 1). In the
first case $t_{iajb}$’s are parametrized by three nearest-neighbor hopping
integrals Rau-PRL ; Floquet-NA : the finite $t_{iajb}$’s for the $Z$ bonds
(Fig. 1) are

$\displaystyle t_{id_{yz}jd_{yz}}=t_{id_{zx}jd_{zx}}=t_{1},$ (9)
$\displaystyle t_{id_{yz}jd_{zx}}=t_{id_{zx}jd_{yz}}=t_{2},$ (10)
$\displaystyle t_{id_{xy}jd_{xy}}=t_{3},$ (11)

and those for the $X$ and the $Y$ bonds are obtained by replacing $(x,y,z)$ in
Eqs. (9)–(11) by $(y,z,x)$ and $(z,x,y)$, respectively. This case, which is
used also in our previous paper Floquet-NA , corresponds to a minimal model of
$\alpha$-RuCl3. In the second case the bond anisotropy of $t_{1}$, $t_{2}$,
and $t_{3}$ is considered: the finite $t_{iajb}$’s for the $Z$ bonds are the
same as Eqs. (9)–(11), whereas those for the $X$ or the $Y$ bonds are

$\displaystyle t_{id_{zx}jd_{zx}}=t_{id_{xy}jd_{xy}}=t_{1}^{\prime},$ (12)
$\displaystyle t_{id_{zx}jd_{xy}}=t_{id_{xy}jd_{zx}}=t_{2}^{\prime},$ (13)
$\displaystyle t_{id_{yz}jd_{yz}}=t_{3}^{\prime},$ (14)

or

$\displaystyle t_{id_{xy}jd_{xy}}=t_{id_{yz}jd_{yz}}=t_{1}^{\prime},$ (15)
$\displaystyle t_{id_{xy}jd_{yz}}=t_{id_{yz}jd_{xy}}=t_{2}^{\prime},$ (16)
$\displaystyle t_{id_{zx}jd_{zx}}=t_{3}^{\prime},$ (17)

respectively. This case is used to study the effects of the weak bond
anisotropy Valenti-PRB  of $\alpha$-RuCl3 on the magnetic properties. In the
third or the fourth case we consider not only the nearest-neighbor hopping
integrals considered in the first or the second case, respectively, but also a
third-neighbor hopping integral $t_{\textrm{3rd}}$. The $t_{\textrm{3rd}}$ for
the $Z_{3}$, $X_{3}$, or $Y_{3}$ bonds (Fig. 1) is the intraorbital hopping
integral of the $d_{xy}$, $d_{yz}$, or $d_{zx}$ orbital, respectively. (We
consider only this among the third neighbor hopping integrals because it is
the largest Valenti-PRB .) The fifth case is almost the same as the fourth
case except for the value of $t_{\textrm{3rd}}$ (see the first paragraph of
Sec. IV B). The last three cases are used to analyze the effects of
$J_{\textrm{3rd}}$, the Heisenberg interaction between third neighbors, on
several magnetic states. In Sec. V we will compare our choices of $t_{iajb}$’s
with the result obtained in first-principles calculations Valenti-PRB .

The atomic SOC is given by the $LS$ coupling of the $t_{2g}$-orbital electrons
StrongLS-Mott-review ; Takagi-review . Because of its nonperturbative effect,
the $t_{2g}$-orbital states with spin degrees of freedom are converted into
the $j_{\textrm{eff}}=1/2$ doublet and the $j_{\textrm{eff}}=3/2$ quartet. (In
this argument we have omitted the components of the $LS$ coupling between the
$t_{2g}$ and the $e_{g}$ orbitals because the crystal field energy between
them is supposed to be sufficiently large; this is the reason why the total
angular momentum is not $j$, but $j_{\textrm{eff}}$.) For $\alpha$-RuCl3 the
low-energy properties can be described by the $j_{\textrm{eff}}=1/2$ doublet
Jeff-RuCl3-1 ; Jeff-RuCl3-2 , which is occupied by an electron (or a hole) per
site,

$\displaystyle|+\rangle_{i}=\frac{1}{\sqrt{3}}(c_{id_{yz}\downarrow}^{\dagger}+ic_{id_{zx}\downarrow}^{\dagger}+c_{id_{xy}\uparrow}^{\dagger})|0\rangle,$
(18)
$\displaystyle|-\rangle_{i}=\frac{1}{\sqrt{3}}(c_{id_{yz}\uparrow}^{\dagger}-ic_{id_{zx}\uparrow}^{\dagger}-c_{id_{xy}\downarrow}^{\dagger})|0\rangle.$
(19)

In these states the spin and the orbital are entangled.

The Coulomb interactions are given by the multiorbital Hubbard interactions
Kanamori  of the $t_{2g}$-orbital electrons:

$\displaystyle
H_{\textrm{int}}=\sum_{i}\Bigl{\\{}\sum_{a,b}c_{ia\uparrow}^{\dagger}c_{ia\downarrow}^{\dagger}[U\delta_{a,b}+J^{\prime}(1-\delta_{a,b})]c_{ib\downarrow}c_{ib\uparrow}$
$\displaystyle+\sum_{\begin{subarray}{c}a,b\\\
a>b\end{subarray}}\sum_{\sigma,\sigma^{\prime}}c_{ia\sigma}^{\dagger}c_{ib\sigma^{\prime}}^{\dagger}(U^{\prime}c_{ib\sigma^{\prime}}c_{ia\sigma}-J_{\textrm{H}}c_{ib\sigma}c_{ia\sigma^{\prime}})\Bigr{\\}},$
(20)

where $U$, $J^{\prime}$, $U^{\prime}$, and $J_{\textrm{H}}$ are the
intraorbital Hubbard interaction, the pair hopping, the interorbital Hubbard
interaction, and the Hund’s coupling, respectively. In deriving the exchange
interactions of periodically driven $\alpha$-RuCl3 (Sec. III A) we use
$H_{\textrm{int}}$ expressed in terms of the irreducible representations of
doubly occupied states NA-Jeff ; Floquet-NA ; Ishihara :

$\displaystyle
H_{\textrm{int}}=\sum_{i}\sum_{\Gamma,g_{\Gamma}}U_{\Gamma}|i;\Gamma,g_{\Gamma}\rangle\langle
i;\Gamma,g_{\Gamma}|,$ (21)

where $U_{\Gamma}$’s are given by

$\displaystyle U_{A_{1}}=U+2J^{\prime},$ (22) $\displaystyle
U_{E}=U-J^{\prime},$ (23) $\displaystyle U_{T_{1}}=U^{\prime}-J_{\textrm{H}},$
(24) $\displaystyle U_{T_{2}}=U^{\prime}+J_{\textrm{H}},$ (25)

and $|i;\Gamma,g_{\Gamma}\rangle$’s are given by

$\displaystyle|i;A_{1}\rangle=\frac{1}{\sqrt{3}}(c_{id_{yz}\uparrow}^{\dagger}c_{id_{yz}\downarrow}^{\dagger}+c_{id_{zx}\uparrow}^{\dagger}c_{id_{zx}\downarrow}^{\dagger}+c_{id_{xy}\uparrow}^{\dagger}c_{id_{xy}\downarrow}^{\dagger})|0\rangle,$
(26)
$\displaystyle|i;E,u\rangle=\frac{1}{\sqrt{6}}(c_{id_{yz}\uparrow}^{\dagger}c_{id_{yz}\downarrow}^{\dagger}+c_{id_{zx}\uparrow}^{\dagger}c_{id_{zx}\downarrow}^{\dagger}$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
-2c_{id_{xy}\uparrow}^{\dagger}c_{id_{xy}\downarrow}^{\dagger})|0\rangle,$
(27)
$\displaystyle|i;E,v\rangle=\frac{1}{\sqrt{2}}(c_{id_{yz}\uparrow}^{\dagger}c_{id_{yz}\downarrow}^{\dagger}-c_{id_{zx}\uparrow}^{\dagger}c_{id_{zx}\downarrow}^{\dagger})|0\rangle,$
(28)
$\displaystyle|i;T_{1},\alpha_{+}\rangle=c_{id_{yz}\uparrow}^{\dagger}c_{id_{zx}\uparrow}^{\dagger}|0\rangle,$
(29)
$\displaystyle|i;T_{1},\alpha_{-}\rangle=c_{id_{yz}\downarrow}^{\dagger}c_{id_{zx}\downarrow}^{\dagger}|0\rangle,$
(30)
$\displaystyle|i;T_{1},\alpha\rangle=\frac{1}{\sqrt{2}}(c_{id_{yz}\uparrow}^{\dagger}c_{id_{zx}\downarrow}^{\dagger}+c_{id_{yz}\downarrow}^{\dagger}c_{id_{zx}\uparrow}^{\dagger})|0\rangle,$
(31)
$\displaystyle|i;T_{1},\beta_{+}\rangle=c_{id_{zx}\uparrow}^{\dagger}c_{id_{xy}\uparrow}^{\dagger}|0\rangle,$
(32)
$\displaystyle|i;T_{1},\beta_{-}\rangle=c_{id_{zx}\downarrow}^{\dagger}c_{id_{xy}\downarrow}^{\dagger}|0\rangle,$
(33)
$\displaystyle|i;T_{1},\beta\rangle=\frac{1}{\sqrt{2}}(c_{id_{zx}\uparrow}^{\dagger}c_{id_{xy}\downarrow}^{\dagger}+c_{id_{zx}\downarrow}^{\dagger}c_{id_{xy}\uparrow}^{\dagger})|0\rangle,$
(34)
$\displaystyle|i;T_{1},\gamma_{+}\rangle=c_{id_{xy}\uparrow}^{\dagger}c_{id_{yz}\uparrow}^{\dagger}|0\rangle,$
(35)
$\displaystyle|i;T_{1},\gamma_{-}\rangle=c_{id_{xy}\downarrow}^{\dagger}c_{id_{yz}\downarrow}^{\dagger}|0\rangle,$
(36)
$\displaystyle|i;T_{1},\gamma\rangle=\frac{1}{\sqrt{2}}(c_{id_{xy}\uparrow}^{\dagger}c_{id_{yz}\downarrow}^{\dagger}+c_{id_{xy}\downarrow}^{\dagger}c_{id_{yz}\uparrow}^{\dagger})|0\rangle,$
(37)
$\displaystyle|i;T_{2},\alpha\rangle=\frac{1}{\sqrt{2}}(c_{id_{yz}\uparrow}^{\dagger}c_{id_{zx}\downarrow}^{\dagger}-c_{id_{yz}\downarrow}^{\dagger}c_{id_{zx}\uparrow}^{\dagger})|0\rangle,$
(38)
$\displaystyle|i;T_{2},\beta\rangle=\frac{1}{\sqrt{2}}(c_{id_{zx}\uparrow}^{\dagger}c_{id_{xy}\downarrow}^{\dagger}-c_{id_{zx}\downarrow}^{\dagger}c_{id_{xy}\uparrow}^{\dagger})|0\rangle,$
(39)
$\displaystyle|i;T_{2},\gamma\rangle=\frac{1}{\sqrt{2}}(c_{id_{xy}\uparrow}^{\dagger}c_{id_{yz}\downarrow}^{\dagger}-c_{id_{xy}\downarrow}^{\dagger}c_{id_{yz}\uparrow}^{\dagger})|0\rangle.$
(40)

##  III Exchange interactions

In this section we study the exchange interactions of periodically driven
$\alpha$-RuCl3. In Sec. III A we derive an effective Hamiltonian expressed in
terms of the exchange interactions in the strong-interaction limit for a high-
frequency case of the driving field. This derivation is performed in the first
case of our model, and the changes in the last four cases are also remarked
on. In Sec. III B we present the dependences of the exchange interactions on
the dimensionless parameter $|u_{ij}|$ for some non-resonant values of
$\omega$. The reason why we use non-resonant $\omega$’s is that our theory is
valid only for such $\omega$’s [see the derivation of Eq. (57)]. We also
discuss the effects of circularly or linearly polarized light and the
similarities and differences between them.

###  III.1 Theory

We derive the exchange interactions of our periodically driven Mott insulator
with strong SOC. Since the derivation for circularly polarized light has been
described in our previous paper Floquet-NA , we explain the main points and
some changes for linearly polarized light. Here we derive the expression in
the first case of our model; the expressions in the other four cases can be
obtained from symmetry arguments (see the fifth paragraph of this section).
The derivation consists of three steps.

First, we derive an effective Hamiltonian for the periodically driven Mott
insulator with strong SOC. To do this, we consider the strong-interaction
limit in which $t_{iajb}$’s in Eq. (5) are much smaller than $U_{\Gamma}$’s in
Eq. (21). In this limit the solution to Schrödinger’s equation can be
approximately expressed as
$|\Psi\rangle_{t}\approx|\Psi_{0}\rangle_{t}+|\Psi_{1}\rangle_{t}$ with
$|\Psi_{0}\rangle_{t}$ and $|\Psi_{1}\rangle_{t}$, the states without and
with, respectively, a doubly occupied site. Thus, the solution can be obtained
by solving a set of simultaneous equations,

$\displaystyle
i\partial_{t}|\Psi_{0}\rangle_{t}=\mathcal{P}_{0}H_{\textrm{KE}}|\Psi_{1}\rangle_{t}+H_{\textrm{SOC}}|\Psi_{0}\rangle_{t},$
(41) $\displaystyle
i\partial_{t}|\Psi_{1}\rangle_{t}=H_{\textrm{KE}}|\Psi_{0}\rangle_{t}+(\mathcal{P}_{1}H_{\textrm{KE}}\mathcal{P}_{1}+\tilde{H}_{\textrm{int}})|\Psi_{1}\rangle_{t},$
(42)

where $\mathcal{P}_{0}$ and $\mathcal{P}_{1}$ denote the projections onto the
subspaces without and with, respectively, a doubly occupied site, and
$\tilde{H}_{\textrm{int}}$ is defined as

$\displaystyle\tilde{H}_{\textrm{int}}=H_{\textrm{int}}+H_{\textrm{SOC}}.$
(43)

Then, we suppose that $\omega$ is much larger than $t_{iajb}$’s. In this
situation $\mathcal{P}_{1}H_{\textrm{KE}}\mathcal{P}_{1}$ in Eq. (42) could be
replaced by its time-averaged one $\bar{H}_{\textrm{KE}}$. As derived in
Appendix A, $\bar{H}_{\textrm{KE}}$ is given by

$\displaystyle\bar{H}_{\textrm{KE}}=\mathcal{P}_{1}\sum_{i,j}\sum_{a,b}\sum_{\sigma}t_{iajb}\mathcal{J}_{0}[\tilde{u}_{ij}^{(\textrm{p})}]c_{ia\sigma}^{\dagger}c_{jb\sigma}\mathcal{P}_{1},$
(44)

where $\mathcal{J}_{n}(x)$ is the $n$th Bessel function of the first kind,
$\tilde{u}_{ij}^{(\textrm{p})}$’s for $\bm{A}(t)=\bm{A}_{\textrm{circ}}(t)$,
$\bm{A}_{\textrm{linear-}b}(t)$, and $\bm{A}_{\textrm{linear-}a}(t)$ are given
by

$\displaystyle\tilde{u}_{ij}^{(\textrm{c})}=u_{ij}=\frac{eE_{0}}{\omega}\textrm{sgn}(i-j),$
(45)
$\displaystyle\tilde{u}_{ij}^{(\textrm{lb})}=\begin{cases}\frac{1}{2}u_{ij}\ \
(X\ \textrm{or}\ Y\ \textrm{bonds}),\\\ u_{ij}\ \ \ \ (Z\
\textrm{bonds}),\end{cases}$ (46)

and

$\displaystyle\tilde{u}_{ij}^{(\textrm{la})}=\begin{cases}\frac{\sqrt{3}}{2}u_{ij}\
\ (X\ \textrm{or}\ Y\ \textrm{bonds}),\\\ \ 0\ \ \ \ \ \ \ (Z\
\textrm{bonds}),\end{cases}$ (47)

respectively, and $\textrm{sgn}(i-j)$ is $1$ for $i\in A$ or $-1$ for $i\in
B$. (Note that $i\in A$ or $i\in B$ means that $i$ is on the $A$ sublattice or
on the $B$ sublattice, respectively.) By using the replacement
$\mathcal{P}_{1}H_{\textrm{KE}}\mathcal{P}_{1}\rightarrow\bar{H}_{\textrm{KE}}$,
we can write Eq. (42) as

$\displaystyle(i\partial_{t}-\bar{H}_{\textrm{KE}}-\tilde{H}_{\textrm{int}})|\Psi_{1}\rangle_{t}=H_{\textrm{KE}}|\Psi_{0}\rangle_{t}.$
(48)

This is equivalent to the following equation:

$\displaystyle
i\partial_{t}[e^{i(\bar{H}_{\textrm{KE}}+\tilde{H}_{\textrm{int}})t}|\Psi_{1}\rangle_{t}]=e^{i(\bar{H}_{\textrm{KE}}+\tilde{H}_{\textrm{int}})t}H_{\textrm{KE}}|\Psi_{0}\rangle_{t}.$
(49)

As derived in Appendix B, the solution to Eq. (49) can be expressed as
follows:

$\displaystyle|\Psi_{1}\rangle_{t}=\sum_{i,j,a,b,\sigma}\sum_{n=-\infty}^{\infty}\frac{t_{iajb}\tilde{\mathcal{J}}_{-n\nu_{ij}}^{(\textrm{p})}(u_{ij})e^{-in\omega
t}}{n\omega-\bar{H}_{\textrm{KE}}-\tilde{H}_{\textrm{int}}}c_{ia\sigma}^{\dagger}c_{jb\sigma}|\Psi_{0}\rangle_{t},$
(50)

where $\tilde{\mathcal{J}}_{-n\nu_{ij}}^{(\textrm{p})}(u_{ij})$’s for
$\bm{A}(t)=\bm{A}_{\textrm{circ}}(t)$, $\bm{A}_{\textrm{linear-}b}(t)$, and
$\bm{A}_{\textrm{linear-}a}(t)$ are given by

$\displaystyle\tilde{\mathcal{J}}_{-n\nu_{ij}}^{(\textrm{c})}(u_{ij})=\begin{cases}\
\mathcal{J}_{-n}(u_{ij})e^{-in\frac{5\pi}{3}}\ \ (\nu_{ij}=X),\\\ \
\mathcal{J}_{-n}(u_{ij})e^{-in\frac{\pi}{3}}\ \ \ (\nu_{ij}=Y),\\\ \
\mathcal{J}_{-n}(u_{ij})e^{-in\pi}\ \ \ \ (\nu_{ij}=Z),\end{cases}$ (51)
$\displaystyle\tilde{\mathcal{J}}_{-n\nu_{ij}}^{(\textrm{lb})}(u_{ij})=\begin{cases}\
\mathcal{J}_{-n}(\frac{u_{ij}}{2})e^{+in\frac{\pi}{2}}\ \ \ (\nu_{ij}=X),\\\ \
\mathcal{J}_{-n}(\frac{u_{ij}}{2})e^{+in\frac{\pi}{2}}\ \ \ (\nu_{ij}=Y),\\\ \
\mathcal{J}_{-n}(u_{ij})e^{-in\frac{\pi}{2}}\ \ \ (\nu_{ij}=Z),\end{cases}$
(52)

and

$\displaystyle\tilde{\mathcal{J}}_{-n\nu_{ij}}^{(\textrm{la})}(u_{ij})=\begin{cases}\
\mathcal{J}_{-n}(\frac{\sqrt{3}}{2}u_{ij})e^{-in\frac{\pi}{2}}\ \
(\nu_{ij}=X),\\\ \
\mathcal{J}_{-n}(\frac{\sqrt{3}}{2}u_{ij})e^{+in\frac{\pi}{2}}\ \
(\nu_{ij}=Y),\\\ \ \delta_{n,0}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
(\nu_{ij}=Z),\end{cases}$ (53)

respectively. By substituting Eq. (50) into Eq. (41) and omitting the constant
term $H_{\textrm{SOC}}|\Psi_{0}\rangle_{t}$, we obtain

$\displaystyle
i\partial_{t}|\Psi_{0}\rangle_{t}=H_{\textrm{eff}}(t)|\Psi_{0}\rangle_{t},$
(54)

where

$\displaystyle H_{\textrm{eff}}(t)=$
$\displaystyle\sum_{i,j,i^{\prime},j^{\prime}}\sum_{n,m=-\infty}^{\infty}\mathcal{P}_{0}T_{i^{\prime}j^{\prime}}\tilde{\mathcal{J}}_{m\nu_{i^{\prime}j^{\prime}}}^{(\textrm{p})}(u_{i^{\prime}j^{\prime}})$
$\displaystyle\times\frac{e^{i(m-n)\omega
t}}{n\omega-\bar{H}_{\textrm{KE}}-\tilde{H}_{\textrm{int}}}\tilde{\mathcal{J}}_{-n\nu_{ij}}^{(\textrm{p})}(u_{ij})T_{ij}\mathcal{P}_{0},$
(55)

and

$\displaystyle
T_{ij}=\sum_{a,b}\sum_{\sigma}t_{iajb}c_{ia\sigma}^{\dagger}c_{jb\sigma}.$
(56)

Furthermore, since in the denominator of Eq. (55) $H_{\textrm{int}}$ gives the
largest contribution of $\bar{H}_{\textrm{KE}}$ and
$\tilde{H}_{\textrm{int}}(=H_{\textrm{int}}+H_{\textrm{SOC}})$, we replace
$n\omega-\bar{H}_{\textrm{KE}}-\tilde{H}_{\textrm{int}}$ in Eq. (55) by
$n\omega-H_{\textrm{int}}$; this replacement may be sufficient for non-
resonant $\omega$’s (i.e., the $\omega$’s at which the denominator does not
diverge). As a result, we have

$\displaystyle H_{\textrm{eff}}(t)=$
$\displaystyle\sum_{i,j}\sum_{n,m=-\infty}^{\infty}\mathcal{P}_{0}T_{ji}\tilde{\mathcal{J}}_{m\nu_{ji}}^{(\textrm{p})}(u_{ji})$
$\displaystyle\times\frac{e^{i(m-n)\omega t}}{n\omega-
H_{\textrm{int}}}\tilde{\mathcal{J}}_{-n\nu_{ij}}^{(\textrm{p})}(u_{ij})T_{ij}\mathcal{P}_{0}.$
(57)

Second, we derive the leading term of $H_{\textrm{eff}}(t)$. By expressing
$H_{\textrm{eff}}(t)$ as the Fourier series
$H_{\textrm{eff}}(t)=\sum_{l}e^{il\omega t}H_{l}$ and using a high-frequency
expansion of the Floquet theory review1 , we can write $H_{\textrm{eff}}(t)$
in the form

$\displaystyle
H_{\textrm{eff}}(t)=H_{0}+O\Bigl{(}\frac{J_{\textrm{ex}}^{2}}{\omega}\Bigr{)},$
(58)

where $J_{\textrm{ex}}$ is of the order of the exchange interactions. If
$\omega$ is high enough to satisfy $\omega\gg|J_{\textrm{ex}}|$, the leading
term of $H_{\textrm{eff}}(t)$ comes from the Floquet Hamiltonian
$\bar{H}_{\textrm{eff}}(=H_{0})$:

$\displaystyle\bar{H}_{\textrm{eff}}=$
$\displaystyle\frac{\omega}{2\pi}\int_{0}^{2\pi/\omega}dtH_{\textrm{eff}}(t)$
$\displaystyle=$
$\displaystyle\sum_{i,j}\sum_{n=-\infty}^{\infty}\mathcal{P}_{0}T_{ji}\frac{\mathcal{J}_{n}[\tilde{u}_{ij}^{(\textrm{p})}]^{2}}{n\omega-
H_{\textrm{int}}}T_{ij}\mathcal{P}_{0}.$ (59)

Furthermore, by using Eq. (21) and expressing the projection operators
$\mathcal{P}_{0}$’s in terms of the possible states, we can rewrite Eq. (59)
as follows:

$\displaystyle\bar{H}_{\textrm{eff}}=$
$\displaystyle\sum_{i,j}\sum_{n=-\infty}^{\infty}\sum_{\Gamma,g_{\Gamma}}\mathcal{P}_{0}T_{ji}|i;\Gamma,g_{\Gamma}\rangle\frac{\mathcal{J}_{n}[\tilde{u}_{ij}^{(\textrm{p})}]^{2}}{n\omega-
U_{\Gamma}}$ $\displaystyle\times\langle
i;\Gamma,g_{\Gamma}|T_{ij}\mathcal{P}_{0}$ $\displaystyle=$
$\displaystyle\sum_{i,j}\sum_{n=-\infty}^{\infty}\sum_{\Gamma,g_{\Gamma}}\sum_{\textrm{i},\textrm{f}}\langle\textrm{f}|T_{ji}|i;\Gamma,g_{\Gamma}\rangle\frac{\mathcal{J}_{n}[\tilde{u}_{ij}^{(\textrm{p})}]^{2}}{n\omega-
U_{\Gamma}}$ $\displaystyle\times\langle
i;\Gamma,g_{\Gamma}|T_{ij}|\textrm{i}\rangle|\textrm{f}\rangle\langle\textrm{i}|,$
(60)

where $|\textrm{i}\rangle$ and $|\textrm{f}\rangle$ are restricted to the
$j_{\textrm{eff}}=1/2$ subspace. Equation (60) shows that a light field
affects the effective Hamiltonian through the factor
$\mathcal{J}_{n}[\tilde{u}_{ij}^{(\textrm{p})}]^{2}/(n\omega-U_{\Gamma})$
(i.e., the changes due to the Bessel functions and the energy shifts of
intermediate states). Note that if we replace
$\sum_{n=-\infty}^{\infty}\mathcal{J}_{n}[\tilde{u}_{ij}^{(\textrm{p})}]^{2}/(n\omega-
U_{\Gamma})$ in Eq. (60) by $1/(-U_{\Gamma})$, the resultant equation gives
the effective Hamiltonian in the absence of $\bm{E}(t)$.

Third, we express Eq. (60) in terms of exchange interactions. To do this, we
calculate the possible terms for the $Z$ bonds on the honeycomb lattice; the
other terms can be obtained from symmetry arguments. $T_{ij}$ for the $Z$
bonds, $T_{ij}^{Z}$, is given by

$\displaystyle T_{ij}^{Z}=$ $\displaystyle
t_{1}\sum_{\sigma}(c_{id_{yz}\sigma}^{\dagger}c_{jd_{yz}\sigma}+c_{id_{zx}\sigma}^{\dagger}c_{jd_{zx}\sigma})$
$\displaystyle+t_{2}\sum_{\sigma}(c_{id_{yz}\sigma}^{\dagger}c_{jd_{zx}\sigma}+c_{id_{zx}\sigma}^{\dagger}c_{jd_{yz}\sigma})$
$\displaystyle+t_{3}\sum_{\sigma}c_{id_{xy}\sigma}^{\dagger}c_{jd_{xy}\sigma}.$
(61)

Then, since Eq. (60) is the sum of two-site terms, we express
$|\textrm{i}\rangle$ and $|\textrm{f}\rangle$ as the products of the
$j_{\textrm{eff}}=1/2$ states at two sites:

$\displaystyle\\{|\textrm{i}\rangle,|\textrm{f}\rangle\\}=\\{|+,+\rangle,|+,-\rangle,|-,+\rangle,|-,-\rangle\\},$
(62)

where $|+,+\rangle=|+\rangle_{1}|+\rangle_{2}$,
$|+,-\rangle=|+\rangle_{1}|-\rangle_{2}$,
$|-,+\rangle=|-\rangle_{1}|+\rangle_{2}$, and
$|-,-\rangle=|-\rangle_{1}|-\rangle_{2}$; $|+\rangle_{i}$ and $|-\rangle_{i}$
have been defined in Eqs. (18) and (19). Since $|i;\Gamma,g_{\Gamma}\rangle$’s
and $U_{\Gamma}$’s are given by Eqs. (26)–(40) and Eqs. (22)–(25), we can
write the finite terms of Eq. (60) for the $Z$ bonds in the form

$\displaystyle\sum_{\langle
i,j\rangle_{Z}}\Bigl{[}J_{Z}\bm{S}_{i}\cdot\bm{S}_{j}+K_{Z}S_{i}^{z}S_{j}^{z}+\Gamma_{Z}(S_{i}^{x}S_{j}^{y}+S_{i}^{y}S_{j}^{x})\Bigr{]},$
(63)

where the summation $\sum_{\langle i,j\rangle_{Z}}$ is over all the $Z$ bonds,

$\displaystyle J_{Z}=$
$\displaystyle\sum_{n=-\infty}^{\infty}\frac{4\mathcal{J}_{n}(u_{ij}^{Z})^{2}}{27}\Bigl{[}\frac{(2t_{1}+t_{3})^{2}}{U+2J^{\prime}-n\omega}+\frac{6t_{1}(t_{1}+2t_{3})}{U^{\prime}-J_{\textrm{H}}-n\omega}$
$\displaystyle\ \ \ \
+\frac{2[(t_{1}-t_{3})^{2}-3t_{2}^{2}]}{U-J^{\prime}-n\omega}+\frac{6t_{2}^{2}}{U^{\prime}+J_{\textrm{H}}-n\omega}\Bigr{]},$
(64) $\displaystyle K_{Z}=$
$\displaystyle\sum_{n=-\infty}^{\infty}\frac{4\mathcal{J}_{n}(u_{ij}^{Z})^{2}}{9}\Bigl{[}\frac{4t_{2}^{2}}{U-J^{\prime}-n\omega}-\frac{(t_{1}-t_{3})^{2}+t_{2}^{2}}{U^{\prime}+J_{\textrm{H}}-n\omega}$
$\displaystyle\ \ \ \
-\frac{3t_{2}^{2}-(t_{1}-t_{3})^{2}}{U^{\prime}-J_{\textrm{H}}-n\omega}\Bigr{]},$
(65) $\displaystyle\Gamma_{Z}=$
$\displaystyle\sum_{n=-\infty}^{\infty}\frac{8\mathcal{J}_{n}(u_{ij}^{Z})^{2}}{9}\Bigl{[}\frac{t_{2}(t_{1}-t_{3})}{U^{\prime}-J_{\textrm{H}}-n\omega}-\frac{t_{2}(t_{1}-t_{3})}{U^{\prime}+J_{\textrm{H}}-n\omega}\Bigr{]},$
(66)

and

$\displaystyle u_{ij}^{Z}=\begin{cases}\ u_{ij}\ \ \textrm{for}\
\bm{A}(t)=\bm{A}_{\textrm{circ}}(t),\\\ \ u_{ij}\ \ \textrm{for}\
\bm{A}(t)=\bm{A}_{\textrm{linear-}b}(t),\\\ \ 0\ \ \ \ \textrm{for}\
\bm{A}(t)=\bm{A}_{\textrm{linear-}a}(t).\end{cases}$ (67)

The derivation of Eq. (63) is described in Appendix C. Equations (64)–(67)
show that for $\bm{A}(t)=\bm{A}_{\textrm{circ}}(t)$ and
$\bm{A}_{\textrm{linear-}b}(t)$ $J_{Z}$, $K_{Z}$, and $\Gamma_{Z}$ can be
changed by varying $\omega$, $E_{0}$, or both, whereas for
$\bm{A}(t)=\bm{A}_{\textrm{linear-}a}(t)$ those remain unchanged due to
$\mathcal{J}_{n}(0)=\delta_{n,0}$. This means that the exchange interactions
for the $Z$ bonds are not affected if the light is polarized along the
$a$-axis. Then, replacing $(x,y,z,Z)$ in Eq. (63) by $(y,z,x,X)$ or
$(z,x,y,Y)$, we obtain the possible terms for the $X$ or the $Y$ bonds,
respectively. As a result, $\bar{H}_{\textrm{eff}}$’s for the $X$ and the $Y$
bonds are given by

$\displaystyle\sum_{\langle
i,j\rangle_{X}}\Bigl{[}J_{X}\bm{S}_{i}\cdot\bm{S}_{j}+K_{X}S_{i}^{x}S_{j}^{x}+\Gamma_{X}(S_{i}^{y}S_{j}^{z}+S_{i}^{z}S_{j}^{y})\Bigr{]},$
(68)

and

$\displaystyle\sum_{\langle
i,j\rangle_{Y}}\Bigl{[}J_{Y}\bm{S}_{i}\cdot\bm{S}_{j}+K_{Y}S_{i}^{y}S_{j}^{y}+\Gamma_{Y}(S_{i}^{z}S_{j}^{x}+S_{i}^{x}S_{j}^{z})\Bigr{]},$
(69)

respectively. Here the summations $\sum_{\langle i,j\rangle_{X}}$ and
$\sum_{\langle i,j\rangle_{Y}}$ are over all the $X$ and the $Y$ bonds,
respectively; $J_{X}$, $J_{Y}$, $K_{X}$, $K_{Y}$, $\Gamma_{X}$, and
$\Gamma_{Y}$ are given by

$\displaystyle J_{X}$ $\displaystyle=J_{Y}$
$\displaystyle=\sum_{n=-\infty}^{\infty}\frac{4\mathcal{J}_{n}(u_{ij}^{X})^{2}}{27}\Bigl{[}\frac{(2t_{1}+t_{3})^{2}}{U+2J^{\prime}-n\omega}+\frac{6t_{1}(t_{1}+2t_{3})}{U^{\prime}-J_{\textrm{H}}-n\omega}$
$\displaystyle\ \ \ \
+\frac{2[(t_{1}-t_{3})^{2}-3t_{2}^{2}]}{U-J^{\prime}-n\omega}+\frac{6t_{2}^{2}}{U^{\prime}+J_{\textrm{H}}-n\omega}\Bigr{]},$
(70) $\displaystyle K_{X}$ $\displaystyle=K_{Y}$
$\displaystyle=\sum_{n=-\infty}^{\infty}\frac{4\mathcal{J}_{n}(u_{ij}^{X})^{2}}{9}\Bigl{[}\frac{4t_{2}^{2}}{U-J^{\prime}-n\omega}-\frac{(t_{1}-t_{3})^{2}+t_{2}^{2}}{U^{\prime}+J_{\textrm{H}}-n\omega}$
$\displaystyle\ \ \ \
-\frac{3t_{2}^{2}-(t_{1}-t_{3})^{2}}{U^{\prime}-J_{\textrm{H}}-n\omega}\Bigr{]},$
(71) $\displaystyle\Gamma_{X}$ $\displaystyle=\Gamma_{Y}$
$\displaystyle=\sum_{n=-\infty}^{\infty}\frac{8\mathcal{J}_{n}(u_{ij}^{X})^{2}}{9}\Bigl{[}\frac{t_{2}(t_{1}-t_{3})}{U^{\prime}-J_{\textrm{H}}-n\omega}-\frac{t_{2}(t_{1}-t_{3})}{U^{\prime}+J_{\textrm{H}}-n\omega}\Bigr{]},$
(72)

where

$\displaystyle u_{ij}^{X}=\begin{cases}\ u_{ij}\ \ \ \ \ \ \textrm{for}\
\bm{A}(t)=\bm{A}_{\textrm{circ}}(t),\\\ \ \frac{1}{2}u_{ij}\ \ \ \
\textrm{for}\ \bm{A}(t)=\bm{A}_{\textrm{linear-}b}(t),\\\ \
\frac{\sqrt{3}}{2}u_{ij}\ \ \textrm{for}\
\bm{A}(t)=\bm{A}_{\textrm{linear-}a}(t).\end{cases}$ (73)

Comparing Eqs. (70)–(73) with Eqs. (64)–(67), we see linearly polarized light
can induce the bond anisotropy of the exchange interactions, i.e., the
differences between the exchange interactions for the $Z$ bonds and those for
the $X$ or $Y$ bonds, even without the bond anisotropy of the hopping
integrals. Such light-induced bond anisotropy does not appear for circularly
polarized light. The origin of this light-induced bond anisotropy is the
difference in the argument of the Bessel function [Eqs. (67) and (73)]. In
Sec. III B we will analyze the light-induced bond anisotropy quantitatively.
Then, by combining Eqs. (68) and (69) with Eq. (63), we can express
$\bar{H}_{\textrm{eff}}$ as follows:

$\displaystyle\bar{H}_{\textrm{eff}}=\sum_{\langle
i,j\rangle}\Bigl{[}J_{\delta}\bm{S}_{i}\cdot\bm{S}_{j}+K_{\delta}S_{i}^{\gamma}S_{j}^{\gamma}+\Gamma_{\delta}(S_{i}^{\alpha}S_{j}^{\beta}+S_{i}^{\beta}S_{j}^{\alpha})\Bigr{]},$
(74)

where

$\displaystyle(\alpha,\beta,\gamma,\delta)=\begin{cases}\ (y,z,x,X)\ \ \ \ (X\
\textrm{bonds}),\\\ \ (z,x,y,Y)\ \ \ \ (Y\ \textrm{bonds}),\\\ \ (x,y,z,Z)\ \
\ \ (Z\ \textrm{bonds}).\end{cases}$ (75)

Before showing the results obtained in numerical calculations, we comment on
the exchange interactions in the other four cases of our model. In the second
case, in which the bond anisotropy of $t_{1}$, $t_{2}$, and $t_{3}$ is
considered [Eqs. (12)–(17)], the exchange interactions for the $X$ and $Y$
bonds are given by

$\displaystyle J_{X}$ $\displaystyle=J_{Y}$
$\displaystyle=\sum_{n=-\infty}^{\infty}\frac{4\mathcal{J}_{n}(u_{ij}^{X})^{2}}{27}\Bigl{[}\frac{(2t_{1}^{\prime}+t_{3}^{\prime})^{2}}{U+2J^{\prime}-n\omega}+\frac{6t_{1}^{\prime}(t_{1}^{\prime}+2t_{3}^{\prime})}{U^{\prime}-J_{\textrm{H}}-n\omega}$
$\displaystyle\ \ \ \
+\frac{2[(t_{1}^{\prime}-t_{3}^{\prime})^{2}-3(t_{2}^{\prime})^{2}]}{U-J^{\prime}-n\omega}+\frac{6(t_{2}^{\prime})^{2}}{U^{\prime}+J_{\textrm{H}}-n\omega}\Bigr{]},$
(76) $\displaystyle K_{X}$ $\displaystyle=K_{Y}$
$\displaystyle=\sum_{n=-\infty}^{\infty}\frac{4\mathcal{J}_{n}(u_{ij}^{X})^{2}}{9}\Bigl{[}\frac{4(t_{2}^{\prime})^{2}}{U-J^{\prime}-n\omega}-\frac{(t_{1}^{\prime}-t_{3}^{\prime})^{2}+(t_{2}^{\prime})^{2}}{U^{\prime}+J_{\textrm{H}}-n\omega}$
$\displaystyle\ \ \ \
-\frac{3(t_{2}^{\prime})^{2}-(t_{1}^{\prime}-t_{3}^{\prime})^{2}}{U^{\prime}-J_{\textrm{H}}-n\omega}\Bigr{]},$
(77) $\displaystyle\Gamma_{X}$ $\displaystyle=\Gamma_{Y}$
$\displaystyle=\sum_{n=-\infty}^{\infty}\frac{8\mathcal{J}_{n}(u_{ij}^{X})^{2}}{9}\Bigl{[}\frac{t_{2}^{\prime}(t_{1}^{\prime}-t_{3}^{\prime})}{U^{\prime}-J_{\textrm{H}}-n\omega}-\frac{t_{2}^{\prime}(t_{1}^{\prime}-t_{3}^{\prime})}{U^{\prime}+J_{\textrm{H}}-n\omega}\Bigr{]},$
(78)

rather than by Eqs. (70)–(72); those for the $Z$ bonds are given by Eqs.
(64)–(66). Thus, by comparing the results in this case and the first case, we
can understand how the light-induced bond anisotropy of the exchange
interactions is affected by the bond anisotropy of the nearest-neighbor
hopping integrals. Those results are shown in Sec. III B. In Sec. IV B we will
analyze the energies of several magnetic states not only in those cases, but
also in the additional three cases. In the latter three cases we consider the
third-neighbor hopping integral $t_{\textrm{3rd}}$, as well as the nearest-
neighbor hopping integrals; it leads to the third-neighbor Heisenberg
interaction

$\displaystyle\sum_{\langle\langle\langle
i,j\rangle\rangle\rangle}J_{\delta}^{\textrm{3rd}}\bm{S}_{i}\cdot\bm{S}_{j},$
(79)

where the summation $\sum_{\langle\langle\langle i,j\rangle\rangle\rangle}$ is
over all the $Z_{3}$, $X_{3}$, and $Y_{3}$ bonds (Fig. 1), $\delta$ is $Z$,
$X$, or $Y$ for the $Z_{3}$, $X_{3}$, or $Y_{3}$ bonds, respectively, and
$J^{\textrm{3rd}}_{\delta}$ is given (in the second-order perturbation theory)
by

$\displaystyle J_{\delta}^{\textrm{3rd}}=$
$\displaystyle\sum_{n=-\infty}^{\infty}\frac{4\mathcal{J}_{n}(2u_{ij}^{\delta})^{2}}{27}\Bigl{(}\frac{t_{\textrm{3rd}}^{2}}{U+2J^{\prime}-n\omega}+\frac{2t_{\textrm{3rd}}^{2}}{U-J^{\prime}-n\omega}\Bigr{)}.$
(80)

[For example, $J_{Z}^{\textrm{3rd}}$ is obtained by replacing $u_{ij}^{Z}$,
$t_{3}$, $t_{1}$, and $t_{2}$ in Eq. (64) by $2u_{ij}^{Z}$,
$t_{\textrm{3rd}}$, $0$, and $0$, respectively.] In the third case of our
model the effective Hamiltonian consists of the sum of Eqs. (74) and (79) with
Eqs. (64)–(67), (70)–(73), and (80). In the fourth or the fifth case the
effective Hamiltonian is the same as that in the third case except that the
exchange interactions for the $X$ and the $Y$ bonds are given by Eqs.
(76)–(78). Analyses of several magnetic states in those three cases may be
useful to clarify the role of $J_{\delta}^{\textrm{3rd}}$, which is shown to
be important without $\bm{E}(t)$ Valenti-PRB , in determining the magnetic
states of periodically driven $\alpha$-RuCl3.

Figure 2:  The $|u_{ij}|(=|\frac{eE_{0}}{\omega}|)$ dependences of $J$, $K$,
and $\Gamma$ with $\bm{E}_{\textrm{circ}}(t)$ at (a) $\omega=3$, (b) $2.2$,
and (c) $1.8$ eV in the first case of our model. In this case there is no bond
anisotropy of the exchange interactions: $J_{Z}=J_{X}=J_{Y}=J$,
$K_{Z}=K_{X}=K_{Y}=K$, and $\Gamma_{Z}=\Gamma_{X}=\Gamma_{Y}=\Gamma$.

In our analyses we do not consider the third-neighbor Kitaev interaction
$K_{\delta}^{\textrm{3rd}}$, although it is also induced by
$t_{\textrm{3rd}}$. This is because the value of $K_{\delta}^{\textrm{3rd}}$
is overestimated if the third-neighbor hopping integrals other than
$t_{\textrm{3rd}}$ are omitted (the value of $K_{\delta}^{\textrm{3rd}}$ is
very small in a more realistic situation Valenti-PRB ). In contrast to
$K_{\delta}^{\textrm{3rd}}$, the value of $J_{\delta}^{\textrm{3rd}}$ is
underestimated. For more details about the values of
$J_{\delta}^{\textrm{3rd}}$ and $K_{\delta}^{\textrm{3rd}}$, see Appendix D.

###  III.2 Results

Figure 3:  The $|u_{ij}|(=|\frac{eE_{0}}{\omega}|)$ dependences of
$J_{X}(=J_{Y})$, $J_{Z}$, $K_{X}(=K_{Y})$, $K_{Z}$, $\Gamma_{X}(=\Gamma_{Y})$,
and $\Gamma_{Z}$ with $\bm{E}_{\textrm{circ}}(t)$ in the second case of our
model. The value of $\omega$ is $3$ eV in (a)–(c) and $2.2$ eV in (d)–(f). In
contrast to the first case with $\bm{E}_{\textrm{circ}}(t)$ (Fig. 2), the
exchange interactions become bond-anisotropic due to the bond anisotropy of
the nearest-neighbor hopping integrals.

We numerically evaluate the exchange interactions for some non-resonant
$\omega$’s in the first two cases of our model. To do this, we replace
$\sum_{n=-\infty}^{\infty}$’s in the exchange interactions by
$\sum_{n=-n_{\textrm{max}}}^{n_{\textrm{max}}}$’s and set
$n_{\textrm{max}}=500$; in the first case the exchange interactions are given
by Eqs. (64)–(66) and Eqs. (70)–(72), whereas in the second case those are
given by Eqs. (64)–(66) and Eqs. (76)–(78). Furthermore, we set
$J^{\prime}=J_{\textrm{H}}$, $U^{\prime}=U-2J_{\textrm{H}}$, $U=3$ eV, and
$J_{\textrm{H}}=0.5$ eV. We choose the values of the nearest-neighbor hopping
integrals as follows: in the first case $t_{1}=47$ meV, $t_{2}=160$ meV, and
$t_{3}=-129$ meV; in the second case $t_{1}=51$ meV, $t_{2}=158$ meV,
$t_{3}=-154$ meV, $t_{1}^{\prime}=45$ meV, $t_{2}^{\prime}=162$ meV, and
$t_{3}^{\prime}=-103$ meV. The values in the first case correspond to the
averages Floquet-NA  of the values obtained in the first-principles
calculations (e.g., $t_{1}$ is the average of $t_{1}$, $t_{1a}^{\prime}$, and
$t_{1b}^{\prime}$ of Ref. Valenti-PRB, ); and the values in the second case
are consistent with those obtained in the first-principles calculations
Valenti-PRB  except that $t_{1}^{\prime}$ is the average of $t_{1a}^{\prime}$
and $t_{1b}^{\prime}$ (i.e., the tiny difference between them is neglected in
our analyses).

First, we present the $|u_{ij}|$ dependences of the exchange interactions with
circularly polarized light. Those in the first case at $\omega=3$, $2.2$, and
$1.8$ eV are shown in Figs. 2(a)–(c). We see that only the magnitudes of the
exchange interactions are changed at $\omega=3$ eV, whereas their magnitudes
and signs can be changed at $\omega=2.2$ and $1.8$ eV. This property remains
unchanged in the second case [Figs. 3(a)–(f)]. The main difference between the
results in these two cases is that there is no bond anisotropy of the exchange
interactions in the first case (i.e., $J_{Z}=J_{X}=J_{Y}=J$,
$K_{Z}=K_{X}=K_{Y}=K$, and $\Gamma_{Z}=\Gamma_{X}=\Gamma_{Y}=\Gamma$), whereas
it is induced by the bond anisotropy of the nearest-neighbor hopping integrals
in the second case. Note that at $|u_{ij}|=0$ in the latter case we have
$J_{X}/J_{Z}\sim 0.63$, $K_{X}/K_{Z}\sim 1.7$, and $\Gamma_{X}/\Gamma_{Z}\sim
0.74$, which are consistent with the values obtained in the first-principles
calculations Valenti-PRB  (i.e., $J_{X}/J_{Z}\sim 0.64$, $K_{X}/K_{Z}\sim
1.5$, and $\Gamma_{X}/\Gamma_{Z}\sim 0.74$).

Figure 4:  The $|u_{ij}|(=|\frac{eE_{0}}{\omega}|)$ dependences of
$J_{X}(=J_{Y})$, $J_{Z}$, $K_{X}(=K_{Y})$, $K_{Z}$, $\Gamma_{X}(=\Gamma_{Y})$,
and $\Gamma_{Z}$ with $\bm{E}_{\textrm{linear-}b}(t)$ in the first case of our
model. The value of $\omega$ is $3$ eV in (a)–(c) and $2.2$ eV in (d)–(f). In
contrast to the first case with $\bm{E}_{\textrm{circ}}(t)$ (Fig. 2), the
light field induces the bond anisotropy of the exchange interactions even
without the bond anisotropy of the hopping integrals.  Figure 5:  The
$|u_{ij}|(=|\frac{eE_{0}}{\omega}|)$ dependences of $J_{X}(=J_{Y})$, $J_{Z}$,
$K_{X}(=K_{Y})$, $K_{Z}$, $\Gamma_{X}(=\Gamma_{Y})$, and $\Gamma_{Z}$ with
$\bm{E}_{\textrm{linear-}a}(t)$ in the first case of our model. The value of
$\omega$ is $3$ eV in (a)–(c) and $2.2$ eV in (d)–(f). As well as the first
case with $\bm{E}_{\textrm{linear-}b}(t)$ (Fig. 4), the bond anisotropy of the
exchange interactions is induced by linearly polarized light; this contrasts
with the case with circularly polarized light (Fig. 2). The exchange
interactions for the $Z$ bonds are independent of $|u_{ij}|$ because of Eq.
(67).  Figure 6:  The $|u_{ij}|(=|\frac{eE_{0}}{\omega}|)$ dependences of
$J_{X}(=J_{Y})$, $J_{Z}$, $K_{X}(=K_{Y})$, $K_{Z}$, $\Gamma_{X}(=\Gamma_{Y})$,
and $\Gamma_{Z}$ with $\bm{E}_{\textrm{linear-}b}(t)$ in the second case of
our model. The value of $\omega$ is $3$ eV in (a)–(c) and $2.2$ eV in (d)–(f).
The bond-anisotropic $|u_{ij}|$ dependences in this case are similar to those
obtained in the first case with $\bm{E}_{\textrm{linear-}b}(t)$ (Fig. 4).
Figure 7:  The $|u_{ij}|(=|\frac{eE_{0}}{\omega}|)$ dependences of
$J_{X}(=J_{Y})$, $J_{Z}$, $K_{X}(=K_{Y})$, $K_{Z}$, $\Gamma_{X}(=\Gamma_{Y})$,
and $\Gamma_{Z}$ with $\bm{E}_{\textrm{linear-}a}(t)$ in the second case of
our model. The value of $\omega$ is $3$ eV in (a)–(c) and $2.2$ eV in (d)–(f).
As well as the case with $\bm{E}_{\textrm{linear-}b}(t)$, the bond-anisotropic
$|u_{ij}|$ dependences in this case are similar to those obtained in the first
case with $\bm{E}_{\textrm{linear-}a}(t)$ (Fig. 5).

Before showing the results with linearly polarized light, we explain the
mechanism of the magnitude or sign changes in the exchange interactions. Here
we focus on the $|u_{ij}|$ dependences of $J$ in the first case at $\omega=3$
and $2.2$ eV. This is enough in understanding the magnitude or sign changes
shown above because of the following three facts: $J$, $K$ and $\Gamma$ have
the similar $|u_{ij}|$ dependences; the results at $\omega=1.8$ eV are
essentially the same as those at $\omega=2.2$ eV; and the $|u_{ij}|$
dependences of the exchange interactions in the second case are similar to
those in the first case. As we explain below, the magnitude changes and the
difference in the sign changes can be understood by considering the leading
terms of Eq. (64), the $n=0$ and the $n=1$ terms. At $\omega=3$ eV in the
first case we can express the leading terms of $J(=J_{Z})$ as follows:

$\displaystyle J$
$\displaystyle\approx(J_{1}+J_{2}+J_{3})\mathcal{J}_{0}(u_{ij})^{2}+(c_{1}^{\prime}J_{1}-c_{2}^{\prime}J_{2})\mathcal{J}_{1}(u_{ij})^{2}$
$\displaystyle\approx(J_{2}+J_{3})\mathcal{J}_{0}(u_{ij})^{2}-c_{2}^{\prime}J_{2}\mathcal{J}_{1}(u_{ij})^{2},$
(81)

where $J_{1}=\frac{4(2t_{1}+t_{3})^{2}}{27(U+2J_{\textrm{H}})}$,
$J_{2}=\frac{8(t_{1}-t_{3})^{2}}{27(U-J_{\textrm{H}})}$,
$J_{3}=\frac{8t_{1}(t_{1}+2t_{3})}{9(U-3J_{\textrm{H}})}$,
$c_{1}^{\prime}=\frac{U+2J_{\textrm{H}}}{\delta\omega_{1}^{\prime}}$,
$c_{2}^{\prime}=\frac{U-J_{\textrm{H}}}{\delta\omega_{2}^{\prime}}$, and
$\omega=U-J_{\textrm{H}}+\delta\omega_{2}^{\prime}=U+2J_{\textrm{H}}-\delta\omega_{1}^{\prime}$
(i.e., $\delta\omega_{2}^{\prime}=0.5$ eV and $\delta\omega_{1}^{\prime}=1$ eV
at $\omega=3$ eV). In deriving the second line of Eq. (81) we have used
$J_{1}\ll J_{2},|J_{3}|$, which is satisfied in $\alpha$-RuCl3. Because of
this property, the $|u_{ij}|$ dependence of $J$ is similar to those of $K$ and
$\Gamma$, as described in Ref. Floquet-NA, . Since $J_{2}>0$, $J_{3}<0$, and
$J_{2}+J_{3}<0$ are also satisfied, Eq. (81) shows that $J$ is always
negative, i.e., its sign remains unchanged, although its magnitude is changed
due to the Bessel functions. Then, at $\omega=2.2$ eV in the first case the
leading terms of Eq. (64) are given by

$\displaystyle J$
$\displaystyle\approx(J_{2}+J_{3})\mathcal{J}_{0}(u_{ij})^{2}+(c_{2}J_{2}-c_{3}J_{3})\mathcal{J}_{1}(u_{ij})^{2},$
(82)

where $c_{2}=\frac{U-J_{\textrm{H}}}{\delta\omega_{2}}$,
$c_{3}=\frac{U-3J_{\textrm{H}}}{\delta\omega_{3}}$, and
$\omega=U-3J_{\textrm{H}}+\delta\omega_{3}=U-J_{\textrm{H}}-\delta\omega_{2}$
(i.e., $\delta\omega_{3}=0.7$ eV and $\delta\omega_{2}=0.3$ eV at $\omega=2.2$
eV). [We have used $J_{1}\ll J_{2},|J_{3}|$ again in the derivation of Eq.
(82).] In contrast to Eq. (81), the term including
$\mathcal{J}_{1}(u_{ij})^{2}$ in Eq. (82) gives the positive-sign
contribution. Thus, we can see from Eq. (82) that the competition between the
negative-sign term including $\mathcal{J}_{0}(u_{ij})^{2}$ and the positive-
sign term including $\mathcal{J}_{1}(u_{ij})^{2}$ is the origin of the sign
changes at $|u_{ij}|\sim 0.4$, $3.5$ in Fig. 2(b); and that the magnitude
changes come from the Bessel functions. Since the leading terms at
$\omega=1.8$ eV is also written in the form of Eq. (82), we can similarly
understand the magnitude and sign changes at $\omega=1.8$ eV. Note that a
similar argument is applicable to Eqs. (65), (66), (70)–(72), and (76)–(78).

Figure 8:  Weakly coupled (a) zigzag and (b) step spin chains with
$\bm{E}_{\textrm{linear-}b}$ and $\bm{E}_{\textrm{linear-}a}$, respectively.
The definitions of blue, green, and red lines and black and white circles are
the same as those in Fig. 1. Thicker bonds represent the bonds with the
exchange interactions which are larger in magnitude.

We turn to the $|u_{ij}|$ dependences of the exchange interactions with
linearly polarized light. The $|u_{ij}|$ dependences with
$\bm{E}_{\textrm{linear-}b}(t)$ are shown in Figs. 4 and 6, and those with
$\bm{E}_{\textrm{linear-}a}(t)$ are shown in Figs. 5 and 7. (Note that the
results at $\omega=1.8$ eV are not shown because they are qualitatively the
same as those at $\omega=2.2$ eV.) Some properties are similar to those with
circularly polarized light (Figs. 2 and 3): the sign changes in the exchange
interactions are absent at $\omega=3$ eV and present at $\omega=2.2$ eV; and
the results obtained in the first case remain qualitatively unchanged even in
the second case. The magnitude and sign changes in the exchange interactions
can be understood in a way similar to those with circularly polarized light.
We also see from Figs. 4–7 that linearly polarized light can change the ratios
$J_{X}/J_{Z}$, $K_{X}/K_{Z}$, and $\Gamma_{X}/\Gamma_{Z}$ even without the
bond anisotropy of the hopping integrals. [As explained below Eq. (73), this
property results from the difference between $u_{ij}^{Z}$ and $u_{ij}^{X}$.]
Because of this property, those ratios can have values which cannot be
realized in non-driven $\alpha$-RuCl3: although $|J_{X}|<|J_{Z}|$,
$|K_{X}|>|K_{Z}|$, and $\Gamma_{X}<\Gamma_{Z}$ hold in non-driven
$\alpha$-RuCl3, $|J_{X}|>|J_{Z}|$, $|K_{X}|<|K_{Z}|$, and
$\Gamma_{X}>\Gamma_{Z}$ are possible in $\alpha$-RuCl3 driven by a field of
linearly polarized light. In addition, it is possible to change the signs of
the exchange interactions only for the $Z$ bonds or only for the $X$ and $Y$
bonds; for example, at $|u_{ij}|\sim 0.5$ in Figs. 4(d)–4(f) we can change the
signs for the $Z$ bonds without changing those for the $X$ and $Y$ bonds.
Then, the honeycomb-network spin system could be transformed either into
weakly coupled zigzag spin chains in the case of
$\bm{E}_{\textrm{linear-}b}(t)$ for $|u_{ij}|\approx 0.4{-}0.42$ at
$\omega=2.2$ eV [e.g., see Figs. 4(d)–4(f) and 6(d)–6(f)] or into weakly
coupled step spin chains in the case of $\bm{E}_{\textrm{linear-}a}(t)$ for
$|u_{ij}|\approx 0.48{-}0.5$ at $\omega=2.2$ eV [e.g., see Figs. 5(d)–5(f) and
7(d)–7(f)]. In the weakly coupled zigzag spin chains [Fig. 8(a)]
$J_{X}=J_{Y}$, $K_{X}=K_{Y}$, and $\Gamma_{X}=\Gamma_{Y}$ are dominant and
$J_{Z}$, $K_{Z}$, and $\Gamma_{Z}$ give the weak coupling between zigzag
chains; in the weakly coupled step spin chains [Fig. 8(b)] $J_{Z}$, $K_{Z}$,
and $\Gamma_{Z}$ are dominant and $J_{X}=J_{Y}$, $K_{X}=K_{Y}$, and
$\Gamma_{X}=\Gamma_{Y}$ give the weak coupling between step chains. Note that
in the first case of our model with $\bm{E}_{\textrm{linear-}b}(t)$ for
$|u_{ij}|\approx 0.42$ at $\omega=2.2$ eV $J_{X}/J_{Z}\sim 9.8$,
$K_{X}/K_{Z}\sim 4.6$, and $\Gamma_{X}/\Gamma_{Z}\sim 4.6$; and that in the
first case with $\bm{E}_{\textrm{linear-}a}(t)$ for $|u_{ij}|\approx 0.48$ at
$\omega=2.2$ eV $J_{X}/J_{Z}\sim 0.1$, $K_{X}/K_{Z}\sim 0.2$, and
$\Gamma_{X}/\Gamma_{Z}\sim 0.2$.

##  IV Magnetic states

In this section we study several magnetic states in periodically driven
$\alpha$-RuCl3. In Sec. IV A we evaluate the expectation value of our Floquet
Hamiltonian within the MFA and express it in a quadratic form. Then, we
explain how to obtain the energies and spin configurations of magnetic states.
We also remark on the magnetic states considered in our analyses. In Sec. IV B
we present the $|u_{ij}|$ dependences of the energies of the magnetic states
for some non-resonant $\omega$’s and discuss the effects of one of the light
fields and the differences due to the polarization of light. (The reason why
we use non-resonant $\omega$’s has been explained in Sec. III.)

###  IV.1 Theory

Applying the MFA to our Floquet Hamiltonian, we derive an expression of its
expectation value in a quadratic form. Since the MFA for Mott insulators with
strong SOC has been explained, for example, in Ref. NA-Jeff, , we explain the
main points here. By using the MFA, we can write the expectation value of Eq.
(74) as

$\displaystyle\langle\bar{H}_{\textrm{eff}}\rangle=\sum_{\langle
i,j\rangle}\sum_{\mu,\nu=x,y,z}M_{\mu\nu}^{ij}\langle
S_{i}^{\mu}\rangle\langle S_{j}^{\nu}\rangle,$ (83)

where $M_{\mu\nu}^{ij}$ is given in the first or second case of our model by

$\displaystyle M_{\mu\nu}^{ij}=\begin{cases}\ J_{\delta}+K_{\delta}\ \
(\mu=\nu=\gamma),\\\ \ J_{\delta}\ \ \ \ \ \ \ \ \ (\mu=\nu=\alpha\
\textrm{or}\ \beta),\\\ \ \Gamma_{\delta}\ \ \ \ \ \ \ \ \
(\mu=\alpha,\nu=\beta),\\\ \ \Gamma_{\delta}\ \ \ \ \ \ \ \ \
(\mu=\beta,\nu=\alpha),\\\ \ 0\ \ \ \ \ \ \ \ \ \ \
(\textrm{otherwise}).\end{cases}$ (84)

[Note that $\delta$, $\gamma$, $\alpha$, and $\beta$ have been defined in Eq.
(75).] In the third, fourth, or fifth case of our model the contribution from
Eq. (79) is added to Eq. (83). In the MFA the expectation value of spin
operators should satisfy the hard-spin constraints:

$\displaystyle|\langle\bm{S}_{i}\rangle|^{2}=S^{2},\
|\langle\bm{S}_{j}\rangle|^{2}=S^{2},$ (85)

where $i$ and $j$ belong to the $A$ and the $B$ sublattices (Fig. 1),
respectively, and $S$ is $1/2$ in the case of $\alpha$-RuCl3. Since $\langle
S_{i}^{\mu}\rangle$ and $\langle S_{j}^{\nu}\rangle$ are expressed as

$\displaystyle\langle S_{i}^{\mu}\rangle$
$\displaystyle=\sqrt{\frac{2}{N}}\sum_{\bm{q}}\langle S_{\bm{q}A}^{\mu}\rangle
e^{i\bm{q}\cdot\bm{R}_{i}},$ (86) $\displaystyle\langle S_{j}^{\nu}\rangle$
$\displaystyle=\sqrt{\frac{2}{N}}\sum_{\bm{q}}\langle S_{\bm{q}B}^{\nu}\rangle
e^{i\bm{q}\cdot\bm{R}_{j}},$ (87)

where $N$ is the total number of sites, we can express Eq. (83) in the
following quadratic form:

$\displaystyle\langle\bar{H}_{\textrm{eff}}\rangle=$
$\displaystyle\sum_{\bm{q}}\sum_{\mu,\nu=x,y,z}\sum_{l,l^{\prime}=A,B}\langle
S_{\bm{q}l}^{\mu}\rangle^{\ast}[M_{ll^{\prime}}(\bm{q})]_{\mu\nu}\langle
S_{\bm{q}l^{\prime}}^{\nu}\rangle$ $\displaystyle=$
$\displaystyle\sum_{\bm{q}}\sum_{\mu,\nu=x,y,z}\langle
S_{-\bm{q}A}^{\mu}\rangle[M(\bm{q})]_{\mu\nu}\langle S_{\bm{q}B}^{\nu}\rangle$
$\displaystyle+\sum_{\bm{q}}\sum_{\mu,\nu=x,y,z}\langle
S_{-\bm{q}B}^{\mu}\rangle[M(\bm{q})^{\dagger}]_{\mu\nu}\langle
S_{\bm{q}A}^{\nu}\rangle,$ (88)

where the $\mu\times\nu$ matrix $M(\bm{q})$ is given by the matrix

$\displaystyle\left(\begin{array}[]{@{\,}ccc@{\,}}J(\bm{q})+K_{x}(\bm{q})&\Gamma_{z}(\bm{q})&\Gamma_{y}(\bm{q})\\\\[3.0pt]
\Gamma_{z}(\bm{q})&J(\bm{q})+K_{y}(\bm{q})&\Gamma_{x}(\bm{q})\\\\[3.0pt]
\Gamma_{y}(\bm{q})&\Gamma_{x}(\bm{q})&J(\bm{q})+K_{z}(\bm{q})\end{array}\right),$
(92)

and $J(\bm{q})$, $K_{\mu}(\bm{q})$’s, and $\Gamma_{\mu}(\bm{q})$’s are defined
as

$\displaystyle
J(\bm{q})=\frac{J_{X}}{2}e^{-i\frac{q_{x}}{2}+i\frac{\sqrt{3}}{2}q_{y}}+\frac{J_{Y}}{2}e^{-i\frac{q_{x}}{2}-i\frac{\sqrt{3}}{2}q_{y}}+\frac{J_{Z}}{2}e^{iq_{x}}$
$\displaystyle\ \
+\frac{J_{X}^{\textrm{3rd}}}{2}e^{iq_{x}-i\sqrt{3}q_{y}}+\frac{J_{Y}^{\textrm{3rd}}}{2}e^{iq_{x}+i\sqrt{3}q_{y}}+\frac{J_{Z}^{\textrm{3rd}}}{2}e^{-2iq_{x}},$
(93) $\displaystyle
K_{x}(\bm{q})=\frac{K_{X}}{2}e^{-i\frac{q_{x}}{2}+i\frac{\sqrt{3}}{2}q_{y}},$
(94) $\displaystyle
K_{y}(\bm{q})=\frac{K_{Y}}{2}e^{-i\frac{q_{x}}{2}-i\frac{\sqrt{3}}{2}q_{y}},$
(95) $\displaystyle K_{z}(\bm{q})=\frac{K_{Z}}{2}e^{iq_{x}},$ (96)
$\displaystyle\Gamma_{x}(\bm{q})=\frac{\Gamma_{X}}{2}e^{-i\frac{q_{x}}{2}+i\frac{\sqrt{3}}{2}q_{y}},$
(97)
$\displaystyle\Gamma_{y}(\bm{q})=\frac{\Gamma_{Y}}{2}e^{-i\frac{q_{x}}{2}-i\frac{\sqrt{3}}{2}q_{y}},$
(98) $\displaystyle\Gamma_{z}(\bm{q})=\frac{\Gamma_{Z}}{2}e^{iq_{x}}.$ (99)

The details of the derivation of Eq. (88) are described in Appendix E.
Furthermore, we make some remarks about momentum in the case of the honeycomb
lattice in Appendix F. Note that the MFA can reproduce the phase diagram
obtained in the Luttinger-Tizsa method Luttinger-Tizsa1 ; Luttinger-Tizsa2
[i.e., Fig. 2(a) of Ref. Rau-PRL,  except the region surrounded by the dashed
white line remark-Rau ].

Figure 9:  Spin configurations of the magnetic states considered in our
analyses: (a) the magnetic states with $\bm{Q}=\bm{0}$, the ferromagnetic
state and the antiferromagnetic state, (b) the magnetic states with
$\bm{Q}=\bm{Q}_{\textrm{ZZ-}X}(=\bm{K}_{01}/2)$ or
$\bm{Q}=\bm{Q}_{\textrm{ZZ-}Z}(=\bm{K}_{10}/2)$, the zigzag states and the
stripy states, and (c) the magnetic state with $\bm{Q}=\bm{Q}_{120}$, the
$120^{\circ}$ order state. Specific directions of spins, which are represented
by the arrows, are chosen because the relative angles between the neighboring
spins are essential for understanding the differences in the spin
configurations.

By using Eqs. (85)–(99), we can obtain the energies and spin configurations of
magnetic states. Since Eq. (88) is quadratic in spin variables, we obtain six
eigenvalues and the corresponding eigenvectors for each $\bm{q}$ by
diagonalizing the matrix

$\displaystyle\left(\begin{array}[]{@{\,}cc@{\,}}0&M(\bm{q})\\\\[3.0pt]
M(\bm{q})^{\dagger}&0\end{array}\right),$ (102)

where $M(\bm{q})$ has been defined in Eq. (92) and $0$ represents the zero
matrix (i.e., the $\mu\times\nu$ matrix of which the components are all zero).
At a certain $\bm{q}$, for example $\bm{q}=\bm{Q}$, the eigenvalue which is
the smallest of the six ones gives the energy of a magnetic state
characterized by the ordering vector $\bm{Q}$. Then, its spin configuration
can be obtained by combining the corresponding eigenvectors and Eqs.
(85)–(87); the obtained spin configuration is valid only if it is consistent
with Eq. (85).

In the analyses of Sec. IV B we consider several magnetic states, which can be
classified into three groups. The first group consists of the magnetic states
with $\bm{Q}=\bm{0}$, which include a ferromagnetic state and an
antiferromagnetic state [Fig. 9(a)]. Note that the $\bm{Q}$ of the
antiferromagnetic state becomes $\bm{0}$ in the presence of a sublattice
structure (because of it, the spins on a sublattice are all ferromagnetic,
i.e., parallel). The second group consists of the magnetic states with
$\bm{Q}=\bm{K}/2$, where $\bm{K}$ is the reciprocal lattice vector [Eq.
(162)]; they include the zigzag states and the stripy states [Fig. 9(b)]. If
we see a unit consisting of one site on the $A$ or $B$ sublattice and its
three neighbors, in the zigzag states three spins are ferromagnetic (parallel)
and the other is antiferromagnetic (antiparallel); in the stripy states two
spins are ferromagnetic and the others are antiferromagnetic [Fig. 9(b)]. In
our analyses we consider two kinds of $\bm{Q}=\bm{K}/2$, i.e., one is
$\bm{Q}_{\textrm{ZZ-}X}=\bm{K}_{01}/2$ and the other is
$\bm{Q}_{\textrm{ZZ-}Z}=\bm{K}_{10}/2$ (for the details of $\bm{K}_{01}$ and
$\bm{K}_{10}$ see Appendix F). In the zigzag state with
$\bm{Q}=\bm{Q}_{\textrm{ZZ-}X}$ or $\bm{Q}=\bm{Q}_{\textrm{ZZ-}Z}$ the spins
on the $X$ or the $Z$ bonds, respectively, are antiferromagnetic; in the
stripy state with $\bm{Q}=\bm{Q}_{\textrm{ZZ-}X}$ or
$\bm{Q}=\bm{Q}_{\textrm{ZZ-}Z}$ those are ferromagnetic. This difference
between the zigzag (or stripy) states with $\bm{Q}=\bm{Q}_{\textrm{ZZ-}X}$ and
$\bm{Q}_{\textrm{ZZ-}Z}$ is partly due to the momentum dependences of
$K_{\mu}(\bm{q})$’s. Namely, since Eqs. (94)–(96) show

$\displaystyle
K_{x}(\bm{Q}_{\textrm{ZZ-}X})=\frac{K_{X}}{2}e^{i\pi}e^{-i\pi/3}=\frac{-K_{X}}{2}e^{-i\pi/3},$
(103) $\displaystyle
K_{y}(\bm{Q}_{\textrm{ZZ-}X})=\frac{K_{Y}}{2}e^{-i\pi/3},$ (104)
$\displaystyle K_{z}(\bm{Q}_{\textrm{ZZ-}X})=\frac{K_{Z}}{2}e^{-i\pi/3},$
(105)

and

$\displaystyle K_{x}(\bm{Q}_{\textrm{ZZ-}Z})=\frac{K_{X}}{2}e^{-i\pi/3},$
(106) $\displaystyle
K_{y}(\bm{Q}_{\textrm{ZZ-}Z})=\frac{K_{Y}}{2}e^{-i\pi/3},$ (107)
$\displaystyle
K_{z}(\bm{Q}_{\textrm{ZZ-}Z})=\frac{K_{Z}}{2}e^{i\pi}e^{-i\pi/3}=\frac{-K_{Z}}{2}e^{-i\pi/3},$
(108)

the effective Kitaev interaction of the $X$ bonds for
$\bm{Q}=\bm{Q}_{\textrm{ZZ-}X}$ or of the $Z$ bonds for
$\bm{Q}=\bm{Q}_{\textrm{ZZ-}Z}$ has the opposite sign and, as a result, the
spins on the $X$ bonds for $\bm{Q}=\bm{Q}_{\textrm{ZZ-}X}$ or the $Z$ bonds
for $\bm{Q}=\bm{Q}_{\textrm{ZZ-}Z}$ are aligned in the opposite direction
compared with the spins on the other bonds. The third group consists of the
magnetic state with $\bm{Q}=\bm{Q}_{120}$, the $120^{\circ}$ order state [Fig.
9(c)]. The magnetic states explained above are realized in Mott insulators
with strong SOC on the honeycomb lattice Rau-PRL ; PD1 .

###  IV.2 Results

Figure 10:  The $|u_{ij}|$ dependences of the energies of the magnetic states
for $\bm{Q}=\bm{0}$, $\bm{Q}_{\textrm{ZZ-}X}(=\bm{K}_{01}/2)$,
$\bm{Q}_{\textrm{ZZ-}Z}(=\bm{K}_{10}/2)$, and $\bm{Q}_{120}$ within the MFA in
the first case of our model with $\bm{E}_{\textrm{circ}}(t)$ [(a) and (d)],
$\bm{E}_{\textrm{linear-}b}(t)$ [(b) and (e)], and
$\bm{E}_{\textrm{linear-}a}(t)$ [(c) and (f)]. The value of $\omega$ is $3$ eV
in (a)–(c) and $2.2$ eV in (d)–(f). In this case the bond-averaged nearest-
neighbor hopping integrals are considered.  Figure 11:  The $|u_{ij}|$
dependences of the energies of the magnetic states for $\bm{Q}=\bm{0}$,
$\bm{Q}_{\textrm{ZZ-}X}(=\bm{K}_{01}/2)$,
$\bm{Q}_{\textrm{ZZ-}Z}(=\bm{K}_{10}/2)$, and $\bm{Q}_{120}$ within the MFA in
the second case of our model with $\bm{E}_{\textrm{circ}}(t)$ [(a) and (d)],
$\bm{E}_{\textrm{linear-}b}(t)$ [(b) and (e)], and
$\bm{E}_{\textrm{linear-}a}(t)$ [(c) and (f)]. The value of $\omega$ is $3$ eV
in (a)–(c) and $2.2$ eV in (d)–(f). In this case the bond-anisotropic nearest-
neighbor hopping integrals are considered. In contrast to the first case (Fig.
10), the degeneracy of the magnetic states for $\bm{Q}=\bm{K}_{01}/2$ and
$\bm{K}_{10}/2$ is lifted even for $\bm{E}_{\textrm{circ}}(t)$ due to the bond
anisotropy of the nearest-neighbor hopping integrals.  Figure 12:  The
$|u_{ij}|$ dependences of the energies of the magnetic states for
$\bm{Q}=\bm{0}$, $\bm{Q}_{\textrm{ZZ-}X}(=\bm{K}_{01}/2)$,
$\bm{Q}_{\textrm{ZZ-}Z}(=\bm{K}_{10}/2)$, and $\bm{Q}_{120}$ within the MFA in
the third case of our model with $\bm{E}_{\textrm{circ}}(t)$ [(a) and (d)],
$\bm{E}_{\textrm{linear-}b}(t)$ [(b) and (e)], and
$\bm{E}_{\textrm{linear-}a}(t)$ [(c) and (f)]. The value of $\omega$ is $3$ eV
in (a)–(c) and $2.2$ eV in (d)–(f). In this case the bond-averaged nearest-
neighbor hopping integrals and the third-neighbor one are considered. In
contrast to the first case (Fig. 10), the magnetic state for
$\bm{Q}=\bm{K}_{01}/2$ or $\bm{K}_{10}/2$ has the lower energy than that for
$\bm{Q}=\bm{0}$ owing to finite $J^{\textrm{3rd}}_{\delta}$ induced by the
third-neighbor hopping integral.  Figure 13:  The $|u_{ij}|$ dependences of
$J_{\delta}^{\textrm{3rd}}$ in the third or fourth case of our model, in which
$t_{\textrm{3rd}}=-40$ meV, with $\bm{E}_{\textrm{circ}}(t)$ [(a) and (b)],
$\bm{E}_{\textrm{linear-}b}(t)$ [(c) and (d)], and
$\bm{E}_{\textrm{linear-}a}(t)$ [(e) and (f)]. The value of $\omega$ is $3$ eV
in (a)–(c) and $2.2$ eV in (d)–(f). In (a) and (b)
$J_{X}^{\textrm{3rd}}=J_{Y}^{\textrm{3rd}}=J_{Z}^{\textrm{3rd}}=J_{\textrm{3rd}}$.
The similar $|u_{ij}|$ dependences are obtained in the fifth case.  Figure 14:
The $|u_{ij}|$ dependences of the energies of the magnetic states for
$\bm{Q}=\bm{0}$, $\bm{Q}_{\textrm{ZZ-}X}(=\bm{K}_{01}/2)$,
$\bm{Q}_{\textrm{ZZ-}Z}(=\bm{K}_{10}/2)$, and $\bm{Q}_{120}$ within the MFA in
the fourth case of our model with $\bm{E}_{\textrm{circ}}(t)$ [(a) and (d)],
$\bm{E}_{\textrm{linear-}b}(t)$ [(b) and (e)], and
$\bm{E}_{\textrm{linear-}a}(t)$ [(c) and (f)]. The value of $\omega$ is $3$ eV
in (a)–(c) and $2.2$ eV in (d)–(f). In this case the bond-anisotropic nearest-
neighbor hopping integrals and the third-neighbor one are considered.  Figure
15:  The $|u_{ij}|$ dependences of the energies of the magnetic states for
$\bm{Q}=\bm{0}$, $\bm{Q}_{\textrm{ZZ-}X}(=\bm{K}_{01}/2)$,
$\bm{Q}_{\textrm{ZZ-}Z}(=\bm{K}_{10}/2)$, and $\bm{Q}_{120}$ within the MFA in
the fifth case of our model with $\bm{E}_{\textrm{circ}}(t)$ [(a) and (d)],
$\bm{E}_{\textrm{linear-}b}(t)$ [(b) and (e)], and
$\bm{E}_{\textrm{linear-}a}(t)$ [(c) and (f)]. The value of $\omega$ is $3$ eV
in (a)–(c) and $2.2$ eV in (d)–(f). In this case the bond-anisotropic nearest-
neighbor hopping integrals and the third-neighbor one are considered; the
difference between the parameters in the fourth and fifth cases is in the
values of $t_{\textrm{3rd}}$.

To study how the light fields affect the magnetic states, we numerically
calculate their energies within the MFA in the five cases of our model at
$\omega=3$ and $2.2$ eV. (The results at $\omega=1.8$ eV, which are not shown,
are qualitatively the same as those at $\omega=2.2$ eV.) As explained in Sec.
IV A, the energy of a magnetic state characterized by $\bm{Q}$ corresponds to
the lowest eigenvalue obtained by diagonalizing Eq. (102) at $\bm{q}=\bm{Q}$.
Furthermore, the $\bm{Q}$’s considered in this study are $\bm{0}$,
$\bm{Q}_{\textrm{ZZ-}X}(=\bm{K}_{01}/2)$,
$\bm{Q}_{\textrm{ZZ-}Z}(=\bm{K}_{10}/2)$, and $\bm{Q}_{120}$. In the numerical
calculations we impose the periodic boundary condition and set
$N_{1}=N_{2}=120$ in Eq. (161). The parameters of our model except
$t_{\textrm{3rd}}$ are chosen in the way described in Sec. III B. We set
$t_{\textrm{3rd}}=-40$ meV in the third or fourth case and
$t_{\textrm{3rd}}=-60$ meV in the fifth case. Note that the former value
corresponds to the average of the intraorbital hopping integral of the
$d_{xy}$ orbital for the $Z_{3}$ bonds, that of the $d_{yz}$ orbital for the
$X_{3}$ bonds, and that of the $d_{xz}$ orbital for the $Y_{3}$ bonds which
are obtained in the first-principles calculations Valenti-PRB ; and that the
latter value is also considered to clarify the effects of
$J_{\delta}^{\textrm{3rd}}$ in detail.

Before discussing the properties for $|u_{ij}|\neq 0$, we comment on the
properties at $u_{ij}=0$ in the five cases of our model. First, the
ferromagnetic state, a magnetic state with $\bm{Q}=\bm{0}$, has the lowest
energy in the first, second, and fourth cases, whereas the zigzag state for
$\bm{Q}=\bm{Q}_{\textrm{ZZ-}Z}$, a magnetic state with $\bm{Q}=\bm{K}_{10}/2$,
is the lowest-energy state in the third and fifth cases; the $120^{\circ}$
order state has the highest energy in all the cases. This result indicates
that the stability of the zigzag or ferromagnetic state is sensitive to the
value of $t_{\textrm{3rd}}$ and the degree of the bond anisotropy of the
hopping integrals. Since the main effect of $t_{\textrm{3rd}}$ is to induce
$J_{\delta}^{\textrm{3rd}}$, our result is consistent with the result obtained
in a minimal model of $\alpha$-RuCl3 Valenti-PRB . Furthermore, the
competition between these magnetic states agrees with the experimental result
al-RuCl3-FM-compete . The sensitivity to $J_{\delta}^{\textrm{3rd}}$ can be
understood by estimating its energy in the MFA: since the expectation value of
Eq. (79) can be written in the MFA as $\sum_{\langle\langle\langle
i,j\rangle\rangle\rangle}J_{\delta}^{\textrm{3rd}}\langle\bm{S}_{i}\rangle\cdot\langle\bm{S}_{j}\rangle$,
the corresponding energies per spin surrounded by three third neighbors for
the ferromagnetic, zigzag, and $120^{\circ}$ order states are
$3J_{\textrm{3rd}}S^{2}$, $-3J_{\textrm{3rd}}S^{2}$, and $0$, respectively,
where
$J_{X}^{\textrm{3rd}}=J_{Y}^{\textrm{3rd}}=J_{Z}^{\textrm{3rd}}=J_{\textrm{3rd}}$.
Then, the zigzag states for $\bm{Q}=\bm{Q}_{\textrm{ZZ-}Z}$ and
$\bm{Q}_{\textrm{ZZ-}X}$, ones of the magnetic states with
$\bm{Q}=\bm{K}_{10}/2$ and $\bm{K}_{01}/2$, are degenerate in the first and
third cases; in the second, fourth, and fifth cases this degeneracy is lifted
and the former state is of lower energy. This lifting is due to the bond
anisotropy of the hopping integrals; in addition, the lower energy of the
zigzag state for $\bm{Q}=\bm{Q}_{\textrm{ZZ-}Z}$ is mainly due to
$|J_{X}+K_{X}|>|J_{Z}+K_{Z}|$, which makes the ferromagnetic spin alignment of
the $X$ bonds more stable than that of the $Z$ bonds. (As explained in Sec. IV
A, two spins on a $X$ or $Z$ bond are antiferromagnetic in the zigzag state
for $\bm{Q}=\bm{Q}_{\textrm{ZZ-}X}$ or $\bm{Q}_{\textrm{ZZ-}Z}$,
respectively.) Note that the energy difference between these zigzag states is
about $0.1$ meV per spin.

We now present the $|u_{ij}|$ dependences of the energies of the magnetic
states in the first case of our model with $\bm{E}_{\textrm{circ}}(t)$,
$\bm{E}_{\textrm{linear-}b}(t)$, or $\bm{E}_{\textrm{linear-}a}(t)$. These
dependences at $\omega=3$ and $2.2$ eV are shown in Fig. 10. The results with
$\bm{E}_{\textrm{circ}}(t)$ [Figs. 10(a) and 10(d)] show that the energies of
the magnetic states with $\bm{Q}=\bm{0}$, $\bm{Q}_{\textrm{ZZ-}X}$, and
$\bm{Q}_{\textrm{ZZ-}Z}$ are close even in the range of $0<|u_{ij}|\leq 1$ and
that the energy of the $120^{\circ}$ order state is much higher than them
except near the $|u_{ij}|$’s at which the exchange interactions are very small
in magnitude [see Figs. 2(b) and 10(d)]. The similar properties hold even in
the results with $\bm{E}_{\textrm{linear-}b}(t)$ or
$\bm{E}_{\textrm{linear-}a}(t)$ at $\omega=3$ eV [Fig. 10(b) or 10(c)].

There are several properties characteristic of the linearly polarized light
fields. One is the lifting of the degeneracy of the magnetic states with
$\bm{Q}=\bm{Q}_{\textrm{ZZ-}X}$ and $\bm{Q}_{\textrm{ZZ-}Z}$. For example, in
the case with $\bm{E}_{\textrm{linear-}b}(t)$ at $\omega=3$ eV the magnetic
state with $\bm{Q}=\bm{Q}_{\textrm{ZZ-}X}$ is of lower energy than that with
$\bm{Q}=\bm{Q}_{\textrm{ZZ-}Z}$ [Fig. 10(b)], whereas in the case with
$\bm{E}_{\textrm{linear-}a}(t)$ at $\omega=3$ eV the latter is of lower energy
[Fig. 10(c)]. This lifting is due to the light-induced bond anisotropy of the
exchange interactions, one of the characteristics of linearly polarized light.
The property that the magnetic state with $\bm{Q}=\bm{Q}_{\textrm{ZZ-}X}$ is
of lower energy than that with $\bm{Q}=\bm{Q}_{\textrm{ZZ-}Z}$ in the case
with $\bm{E}_{\textrm{linear-}b}(t)$ at $\omega=3$ eV comes from the facts
that the exchange interactions for the $Z$ bonds are larger in magnitude than
those for the $X$ or $Y$ bonds and that $J_{Z}$ and $K_{Z}$ are ferromagnetic
[see Figs. 4(a)–(c) in the range of $0\leq|u_{ij}|\leq 1$]; the property in
the case with $\bm{E}_{\textrm{linear-}a}(t)$ at $\omega=3$ eV can be
similarly understood. The other characteristic properties are the changes in
the competing magnetic states. From Fig. 10(e) we see the competing magnetic
states in the range of $0<|u_{ij}|\leq 0.4$ are the magnetic states with
$\bm{Q}=\bm{0}$ and $\bm{Q}_{\textrm{ZZ-}Z}$, whereas those in the range of
$0.6\leq|u_{ij}|\leq 0.8$ become the magnetic states with
$\bm{Q}=\bm{Q}_{\textrm{ZZ-}Z}$, $\bm{Q}_{\textrm{ZZ-}X}$, and $\bm{Q}_{120}$.
In addition, from Fig. 10(f) we see a similar change in the competing magnetic
states. These results are related to the light-induced bond anisotropy of the
exchange interactions because in the range of $0.6\leq|u_{ij}|\leq 0.8$ with
$\bm{E}_{\textrm{linear-}b}(t)$ or $\bm{E}_{\textrm{linear-}a}(t)$ the signs
of the exchange interactions only for the $Z$ bonds or for the $X$ and $Y$
bonds, respectively, are changed and their magnitudes become larger than those
for the other bonds [see Figs. 4(d)–(f) for $\bm{E}_{\textrm{linear-}b}(t)$
and Figs. 5(d)–(f) for $\bm{E}_{\textrm{linear-}a}(t)$].

We turn to the results in the second case of our model. Figure 11 shows the
$|u_{ij}|$ dependences of the energies of the magnetic states in this case
with $\bm{E}_{\textrm{circ}}(t)$, $\bm{E}_{\textrm{linear-}b}(t)$, or
$\bm{E}_{\textrm{linear-}a}(t)$. These results are qualitatively the same as
those in the first case except that the degeneracy of the magnetic states with
$\bm{Q}=\bm{Q}_{\textrm{ZZ-}X}$ and $\bm{Q}_{\textrm{ZZ-}Z}$ is lifted even
for $\bm{E}_{\textrm{circ}}(t)$. (As described above, this lifting results
from the bond anisotropy of the hopping integrals, which is absent in the
first case and present in the second case.) Combining this result with the
results at $u_{ij}=0$, we find that the main effects of the bond anisotropy of
the hopping integrals are to lift the degeneracy of the magnetic states with
$\bm{Q}=\bm{Q}_{\textrm{ZZ-}X}$ and $\bm{Q}_{\textrm{ZZ-}Z}$ and to decrease
the energy of the ferromagnetic state compared with those of the zigzag states
at $u_{ij}=0$.

As well as the results in the second cases, the results in the third case
(Fig. 12) are similar to those in the first case except for two differences.
The two differences are that at $u_{ij}=0$ the energies of the zigzag states
are lower than that of the ferromagnetic state and that the magnetic state
with $\bm{q}=\bm{0}$ becomes of lower energy than those of the magnetic states
with $\bm{Q}=\bm{Q}_{\textrm{ZZ-}X}$ and $\bm{Q}_{\textrm{ZZ-}Z}$ above a
certain value of $|u_{ij}|$ [e.g., see the values at $|u_{ij}|=0.46$ in Fig.
12(a)]. These differences result from the effects of
$J_{\delta}^{\textrm{3rd}}$ because the latter difference can be understood
from the $|u_{ij}|$ dependences of $J_{\delta}^{\textrm{3rd}}$ (Fig. 13); for
example, in the case with $\bm{E}_{\textrm{circ}}(t)$ at $\omega=3$ eV the
blue and red lines of Fig. 12(a) cross at the value of $|u_{ij}|$ at which
$J_{\textrm{3rd}}$ is small in magnitude [Fig. 13(a)]. Then, as described
above, the former difference can be understood from the difference in the
energies due to $J_{\delta}^{\textrm{3rd}}$.

The similar effects of $J_{\delta}^{\textrm{3rd}}$ appear in the fourth and
fifth cases (Figs. 14 and 15). However, in the presence of the bond anisotropy
of the nearest-neighbor hopping integrals $t_{\textrm{3rd}}=-40$ meV is not
sufficient for making the energies of the zigzag states lower at $u_{ij}=0$
than that of the ferromagnetic state (Fig. 14), although the zigzag states are
of lower energy at $u_{ij}=0$ for $t_{\textrm{3rd}}=-60$ meV (Fig. 15).
Nevertheless, we believe this does not contradict the experimental result that
the zigzag state is stabilized in $\alpha$-RuCl3 because
$J_{\delta}^{\textrm{3rd}}$ is underestimated in our simplified treatment (see
Sec. III A and Appendix D); a more accurate calculation is beyond the scope of
this paper.

The most important thing seen from Figs. 11, 12, 14, and 15 is that the
characteristic properties found in the first case remain qualitatively
unchanged even in the other four cases. Thus, we believe the lifting of the
magnetic states with $\bm{Q}=\bm{Q}_{\textrm{ZZ-}X}$ and
$\bm{Q}_{\textrm{ZZ-}Z}$ and the change in the competing magnetic states are
characteristic of linearly polarized light. Although the former can be
realized by using the bond-anisotropic nearest-neighbor hopping integrals, the
latter is a unique effect of linearly polarized light.

##  V Discussion

First, we discuss the validity of our model. Our model has the three nearest-
neighbor hopping integrals, including their bond anisotropy, and the third-
neighbor one. For the Mott insulating state the effective Hamiltonian consists
of the nearest-neighbor Heisenberg, Kitaev, and off-diagonal symmetric
exchange interactions and the third-neighbor Heisenberg interaction (i.e.,
$J_{\delta}$, $K_{\delta}$, $\Gamma_{\delta}$, and
$J_{\delta}^{\textrm{3rd}}$). We believe our model is sufficient for
describing the magnetic properties of $\alpha$-RuCl3 because its first-
principles calculations Valenti-PRB  showed that the leading hopping integrals
are $t_{2}$ and $t_{3}$, which are at least an order of magnitude larger than
the others and suggested that its minimal spin model consists of the bond-
averaged $J_{\delta}$, $K_{\delta}$, $\Gamma_{\delta}(=-K_{\delta})$, and
$J_{\delta}^{\textrm{3rd}}$. In addition, the signs of the exchange
interactions of our model are consistent with those estimated by fitting a
magnetization measurement Exp-FM-Kitaev . Note that according to this
measurement, another off-diagonal symmetric exchange interaction
$\Gamma^{\prime}$, which does not appear in our model, is smaller than these
four exchange interactions. The similar results are obtained in other studies
al-RuCl3-FM-compete ; Valenti-PRB . This is reasonable because
$\Gamma^{\prime}$ is proportional to the other nearest-neighbor interorbital
hopping integral Rau-PRL , which is small in the case of $\alpha$-RuCl3
Valenti-PRB .

We should note that although the effect of SOC on the coefficients of the
exchange interactions is necessary for discussing their values quantitatively
Valenti-PRB , the following arguments indicate that its effect may be not
large in the case of $\alpha$-RuCl3. In general, $H_{\textrm{SOC}}$ affects
the coefficients of the exchange interactions of Mott insulators, as we can
see from Eq. (55). Since the $LS$-type SOC induces the onsite interorbital
excitations, it can connect $|i;\Gamma,g_{\Gamma}\rangle$ with
$|i;\Gamma,g_{\Gamma}^{\prime}\rangle$ for $g_{\Gamma}^{\prime}\neq
g_{\Gamma}$; such off-diagonal terms result in the degeneracy lifting of the
states $|i;\Gamma,g_{\Gamma}\rangle$’s for given $\Gamma$. Thus, the main
effect of $H_{\textrm{SOC}}$ on the coefficients of the exchange interactions
is to change the energies of the intermediate states in the second-order
perturbation processes considered to derive the exchange interactions; the
energy for a certain $\Gamma$, including this effect of $H_{\textrm{SOC}}$,
could be expressed as $E_{\Gamma}\pm c_{\Gamma}\lambda$, where $\lambda$ is
the coupling constant of $H_{\textrm{SOC}}$ and $c_{\Gamma}=O(1)$ or $O(0.1)$.
Accordingly, the modulations of the exchange interactions which come from the
intermediate states for $\Gamma$ are roughly given by
$E_{\Gamma}/(E_{\Gamma}\pm c_{\Gamma}\lambda)\sim 1\mp
c_{\Gamma}(\lambda/E_{\Gamma})$. Since $\lambda=O(0.1\textrm{eV})$ in
$\alpha$-RuCl3, we have $(\lambda/E_{\Gamma})=O(0.1)$, and thus the
corrections due to SOC are small for $\alpha$-RuCl3. From the above arguments,
we conclude that our treatment, in which the effect of SOC on the coefficients
of the exchange interactions is neglected, is sufficient for discussing the
exchange interactions of $\alpha$-RuCl3 qualitatively.

Next, we make some remarks about heating effects. Since the periodic driving
field causes the system to heat up, it eventually approaches an infinite-
temperature state Heat-InfT1 . However, there are intermediate times
$t\lesssim\tau$ at which the periodically driven system can be approximately
described by the Floquet Hamiltonian Heat-Floq1 ; Heat-Floq2 . Since $\tau$ is
roughly given by $\tau\approx T\exp(\omega/J_{\textrm{ex}})$ Heat-Floq2  and
our parameters satisfy $\omega=O(1\textrm{eV})$ and
$J_{\textrm{ex}}=\textrm{max}(J_{\delta},K_{\delta},\Gamma_{\delta},J_{\delta}^{\textrm{3rd}})=O(10\textrm{meV})$,
the intermediate times of our system may be sufficiently large. Note that
because of these values of $\omega$ and $J_{\textrm{ex}}$, the correction to
the Floquet Hamiltonian, the second term of Eq. (58), is negligible. Then,
since our $\omega$ is non-resonant, the heating effect due to the doublon
creation induced by the driving field is also negligible Floquet-MultiMott .
(Since this heating effect is non-negligible for resonant or nearly resonant
$\omega$, we have shown the results for some non-resonant $\omega$’s in Secs.
III B and IV B.) It should be noted that the effects of the doublon creation
due to the driving field can be described by $\bar{H}_{\textrm{KE}}$ in Eq.
(55), and that they become negligible if $\omega$ is non-resonant in the sense
that the denominator of Eq. (55) does not diverge Floquet-MultiMott . (If
$\omega$ is resonant or nearly resonant, $\bar{H}_{\textrm{KE}}$ induces a
non-negligible imaginary part, resulting in the heating effects Floquet-
MultiMott .) Thus, we believe our results based on the Floquet Hamiltonian are
meaningful as the properties of periodically driven $\alpha$-RuCl3.

We also remark on a property induced by the field of circularly polarized
light. It has been shown for a single-orbital Hubbard model driven by
circularly polarized light Chirality1 ; Chirality2  that when the light
frequency is comparable with the Hubbard interaction, the effective
Hamiltonian of the Mott insulator could acquire a spin scalar chirality, which
is non-negligible compared with the antiferromagnetic Heisenberg interaction.
Note that the spin scalar chirality terms come from the fourth-order
perturbation processes in which the kinetic terms are treated as perturbation,
whereas the exchange interactions come from the second-order ones. Although a
similar mechanism might work in more complicated models, our rough estimate
shown below indicates that such contributions may be negligible in our cases.
Since the dominant Bessel functions appearing in the exchange interactions and
the spin scalar chirality terms are $\mathcal{J}_{0}(u_{ij})$ and
$\mathcal{J}_{1}(u_{ij})$, a ratio of $J_{\textrm{chi}}$, one of the spin
scalar chirality terms, to $J_{\textrm{exch}}$, one of the exchange
interactions, may be roughly given by
$|J_{\textrm{chi}}/J_{\textrm{exch}}|\sim
t^{2}\mathcal{J}_{1}(u_{ij})^{2}/(U_{\textrm{int}}|U_{\textrm{int}}-\omega|)$,
where $t$ is of the order of the hopping integrals and $U_{\textrm{int}}$ is
of the order of the onsite Coulomb interactions; in some cases of our analyses
our parameters correspond to $t=O(0.1\textrm{eV})$,
$U_{\textrm{int}}=O(1\textrm{eV})$, and
$|U_{\textrm{int}}-\omega|=O(0.1\textrm{eV})$. (Precisely speaking, our
$|U_{\Gamma}-n\omega|$’s for any $\Gamma$ and allowed $n$ are larger than the
hopping integrals.) In addition, since $\mathcal{J}_{1}(u_{ij})^{2}=O(0.1)$,
we have $|J_{\textrm{chi}}/J_{\textrm{ex}}|\sim 10^{-2}$. Although there is
another contribution to the spin scalar chirality terms,
$J_{\textrm{chi}}^{\prime}$, it may be smaller than the above contribution
because of an additional factor $|U_{\textrm{int}}-\omega|/\omega=O(10^{-1})$
[i.e.,
$|J_{\textrm{chi}}^{\prime}/J_{\textrm{chi}}|\sim|U_{\textrm{int}}-\omega|/\omega=O(10^{-1})$].
Thus, we believe the spin scalar chirality terms are negligible and the
results shown in Sec. IV B remain qualitatively unchanged.

We now address an experimental observation of our results. Controlling the
exchange interactions via a light field, in principle, is achievable by
performing pump-probe measurements. However, a strong light field is necessary
because the $\omega$’s considered in our study are high. For example, to
realize the periodically driven $\alpha$-RuCl3 at $\omega=2.2$ eV for
$|u_{ij}|\sim 0.4$ or $0.6$, the amplitude of a light field, $E_{0}$, should
be $E_{0}\sim 26$ or $39$ MVcm-1, respectively [for the relation between
$u_{ij}$ and $E_{0}$ see Eq. (45)]; in the cases at $\omega=1.8$ eV for
$|u_{ij}|\sim 0.4$ or $0.6$, we have $E_{0}\sim 21$ or $32$ MVcm-1,
respectively. In these estimates we have used $a_{\textrm{NN}}\sim 3.4\times
10^{-8}$cm, which corresponds to the average of the experimentally observed
lengths of two nearest-neighbor bonds al-RuCl3-lattice . Then, the properties
of the magnetic states could be observed by using, for example, neutron
scattering measurements. Specifically, the changes in the competing magnetic
states due to linearly polarized light could be detected as the evolution of
short-range correlations characterized by the corresponding ordering vectors.
Since a light field of the order of $10$ MVcm-1 can be realized experimentally
Iwai-review , we hope our main results, i.e., the changes in the exchange
interactions, their bond anisotropy, and the competing magnetic states via
linearly polarized light, will be observed by experiments.

Finally, we comment on several directions for further relevant research.
First, our theory can be extended to the cases of $\alpha$-RuCl3 with both a
light field and an external magnetic field and of other periodically driven
Mott insulators on the honeycomb lattice (e.g., some Ir oxides). Our theory
may be also useful for formulating a theory in the case on another lattice
with strong SOC. The studies in these cases are contained in the possible
research directions. Another research direction is to study a possibility of
Kitaev spin liquids in a more elaborate method than the MFA. Such recent
studies include Ref. Floquet-StrongSOC, . In addition, our results may be
useful for realizing a gapped spin liqud, a toric code phase Kitaev ; PRR-
Fujimoto , because it could be stabilized in the presence of strong bond
anisotropy of the exchange interactions. One of the possible situations might
be $\alpha$-RuCl3 with both an inplane magnetic field $H_{ab}$ and a linearly
polarized light field because a gapless spin liquid could be stabilized for
$\alpha$-RuCl3 in the range of $7.5\textrm{T}<H_{ab}<16\textrm{T}$ Hind-SL
and a linearly polarized light field could induce the strong bond anisotropy.
Then, our results about the exchange interactions could be used to control
magnetization dynamics and spintronics phenomena of Mott insulators with
strong SOC because the key quantities to describe them are the exchange
interactions review-opt ; review-spintronics . Thus, the extensions to
dynamical or transport properties are ones of the important future research
directions. Another important future research direction is an extension to the
case for small $\omega$, at which the heating effects are no longer
negligible.

##  VI Conclusion

We have studied the magnetic properties of $\alpha$-RuCl3 driven by circularly
or linearly polarized light. We showed that, as well as the magnitudes and
signs of the exchange interactions, their bond anisotropy can be changed by
tuning the amplitude and frequency of one of the fields of linearly polarized
light. This is one of the characteristics of linearly polarized light because
the bond anisotropy is not induced by circularly polarized light. Since the
light-induced bond anisotropy can be used to change the ratios of $J$, $K$,
and $\Gamma$ for the $Z$ bonds to those for the $X$ or $Y$ bonds, the
honeycomb-network spin system could be transformed into weakly coupled zigzag
or step chains for $\bm{E}_{\textrm{linear-}b}(t)$ or
$\bm{E}_{\textrm{linear-}a}(t)$, respectively. We also showed that the
$|u_{ij}|$ dependences of the exchange interactions obtained in the first case
remain qualitatively unchanged except for the degeneracy lifting of the
exchange interactions for the $Z$ bonds and for the $X$ and $Y$ bonds. Then,
we showed that the competing magnetic states can be changed only for linearly
polarized light. Such a situation could be realized by using the strong field
of the order of $10$MVcm-1 in pump-probe measurements. We also showed that the
bond anisotropy of the nearest-neighbor hopping integrals and the third-
neighbor hopping integral do not change qualitatively the results obtained in
the first case of our model except for the stability of the zigzag states at
$|u_{ij}|=0$ and the degeneracy lifting of the magnetic states with
$\bm{Q}=\bm{Q}_{\textrm{ZZ-}X}$ and $\bm{Q}_{\textrm{ZZ-}Z}$. We believe this
paper is useful for further research of $\alpha$-RuCl3 and the relevant
materials such as the honeycomb iridates and provides an important step
towards a comprehensive understanding of magnetic properties of periodically
driven Mott insulators with strong SOC.

###### Acknowledgements.

This work was supported by JST CREST Grant No. JPMJCR1901, JSPS KAKENHI Grants
No. JP19K14664 and No. JP16K05459, and MEXT Q-LEAP Grant No. JP-MXS0118067426.

##  Appendix A Derivation of Eq. (44)

We calculate $\bar{H}_{\textrm{KE}}$’s for
$\bm{A}(t)=\bm{A}_{\textrm{circ}}(t)$, $\bm{A}_{\textrm{linear-}b}(t)$, and
$\bm{A}_{\textrm{linear-}a}(t)$. $\bar{H}_{\textrm{KE}}$ is given by

$\displaystyle\bar{H}_{\textrm{KE}}=$
$\displaystyle\mathcal{P}_{1}\frac{\omega}{2\pi}\int_{0}^{2\pi/\omega}dtH_{\textrm{KE}}\mathcal{P}_{1}$
$\displaystyle=$
$\displaystyle\mathcal{P}_{1}\frac{\omega}{2\pi}\int_{0}^{2\pi/\omega}dt\sum_{i,j}\sum_{a,b}\sum_{\sigma}t_{iajb}e^{-ie(\bm{R}_{i}-\bm{R}_{j})\cdot\bm{A}(t)}$
$\displaystyle\times c_{ia\sigma}^{\dagger}c_{jb\sigma}\mathcal{P}_{1}.$ (109)

To perform the time integral, we rewrite the Peierls phase factor using Eqs.
(6)–(8). For $\bm{A}(t)=\bm{A}_{\textrm{circ}}(t)$, we write it as follows
Floquet-NA :

$\displaystyle e^{-ie(\bm{R}_{i}-\bm{R}_{j})\cdot\bm{A}(t)}=\begin{cases}\
e^{iu_{ij}\sin(\omega t+\frac{5\pi}{3})}\ \ (X\ \textrm{bonds}),\\\ \
e^{iu_{ij}\sin(\omega t+\frac{\pi}{3})}\ \ \ (Y\ \textrm{bonds}),\\\ \
e^{iu_{ij}\sin(\omega t+\pi)}\ \ \ \ (Z\ \textrm{bonds}),\end{cases}$ (110)

where $u_{ij}$ is defined in Eq. (45). Similarly, we obtain

$\displaystyle e^{-ie(\bm{R}_{i}-\bm{R}_{j})\cdot\bm{A}(t)}=\begin{cases}\
e^{i\frac{1}{2}u_{ij}\sin(\omega t-\frac{\pi}{2})}\ \ (X\ \textrm{or}\ Y\
\textrm{bonds}),\\\ \ e^{iu_{ij}\sin(\omega t+\frac{\pi}{2})}\ \ \ \ (Z\
\textrm{bonds})\end{cases}$ (111)

for $\bm{A}(t)=\bm{A}_{\textrm{linear-}b}(t)$, and

$\displaystyle e^{-ie(\bm{R}_{i}-\bm{R}_{j})\cdot\bm{A}(t)}=\begin{cases}\
e^{i\frac{\sqrt{3}}{2}u_{ij}\sin(\omega t+\frac{\pi}{2})}\ \ (X\
\textrm{bonds}),\\\ \ e^{i\frac{\sqrt{3}}{2}u_{ij}\sin(\omega
t-\frac{\pi}{2})}\ \ (Y\ \textrm{bonds}),\\\ \ 1\ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ (Z\ \textrm{bonds})\end{cases}$ (112)

for $\bm{A}(t)=\bm{A}_{\textrm{linear-}a}(t)$. By combining Eqs. (110)–(112)
with Eq. (109) and using the relations

$\displaystyle e^{ix\sin\theta}=$
$\displaystyle\sum_{n=-\infty}^{\infty}\mathcal{J}_{n}(x)e^{in\theta}$ (113)

and

$\displaystyle\frac{\omega}{2\pi}\int_{0}^{2\pi/\omega}dte^{in\omega
t}=\delta_{n,0},$ (114)

we obtain Eq. (44).

##  Appendix B Derivation of Eq. (50)

We solve Eq. (49). By integrating both sides and choosing for the initial
state the state in which the lower-limit contributions of the integrals cancel
each other out, we have

$\displaystyle
ie^{i(\bar{H}_{\textrm{KE}}+\tilde{H}_{\textrm{int}})t}|\Psi_{1}\rangle_{t}=$
$\displaystyle\int^{t}dt^{\prime}e^{i(\bar{H}_{\textrm{KE}}+\tilde{H}_{\textrm{int}})t^{\prime}}H_{\textrm{KE}}|\Psi_{0}\rangle_{t^{\prime}}.$
(115)

Furthermore, since the time variation of $|\Psi_{0}\rangle_{t^{\prime}}$ is
slow, we could write Eq. (115) as

$\displaystyle
ie^{i(\bar{H}_{\textrm{KE}}+\tilde{H}_{\textrm{int}})t}|\Psi_{1}\rangle_{t}\approx$
$\displaystyle\int^{t}dt^{\prime}e^{i(\bar{H}_{\textrm{KE}}+\tilde{H}_{\textrm{int}})t^{\prime}}H_{\textrm{KE}}|\Psi_{0}\rangle_{t}.$
(116)

By combining Eq. (116) with Eq. (5) and Eqs. (110)–(113) and performing the
time integral, we obtain Eq. (50).

##  Appendix C Derivation of Eq. (63)

We calculate the possible terms of Eq. (60) for the $Z$ bonds. The
calculations consist of two steps.

First, we calculate $\langle i;\Gamma,g_{\Gamma}|T_{ij}|\textrm{i}\rangle$ for
$T_{ij}=T_{ij}^{Z}$. Since

$\displaystyle|\textrm{i}\rangle=\\{|+\rangle_{1}|+\rangle_{2},|+\rangle_{1}|-\rangle_{2},|-\rangle_{1}|+\rangle_{2},|-\rangle_{1}|-\rangle_{2}\\},$
(117)

we calculate the finite terms of $\langle
i;\Gamma,g_{\Gamma}|T_{12}^{Z}|\textrm{i}\rangle$; the contributions from
$\langle i;\Gamma,g_{\Gamma}|T_{21}^{Z}|\textrm{i}\rangle$ to
$\bar{H}_{\textrm{eff}}$ can be taken into account by multiplying those from
$\langle i;\Gamma,g_{\Gamma}|T_{12}^{Z}|\textrm{i}\rangle$ by two. Using Eqs.
(18), (19), and (61), we have

$\displaystyle T_{12}^{Z}|+,+\rangle=$
$\displaystyle\frac{1}{3}\Bigl{[}(t_{3}-t_{1})c_{1d_{yz}\downarrow}^{\dagger}c_{1d_{xy}\uparrow}^{\dagger}-i(t_{1}-t_{3})c_{1d_{zx}\downarrow}^{\dagger}c_{1d_{xy}\uparrow}^{\dagger}$
$\displaystyle+2t_{2}c_{1d_{yz}\downarrow}^{\dagger}c_{1d_{zx}\downarrow}^{\dagger}-it_{2}c_{1d_{yz}\downarrow}^{\dagger}c_{1d_{xy}\uparrow}^{\dagger}$
$\displaystyle-
t_{2}c_{1d_{zx}\downarrow}^{\dagger}c_{1d_{xy}\uparrow}^{\dagger}\Bigr{]}|0\rangle,$
(118) $\displaystyle T_{12}^{Z}|-,-\rangle=$
$\displaystyle\frac{1}{3}\Bigl{[}(t_{1}-t_{3})c_{1d_{yz}\uparrow}^{\dagger}c_{1d_{xy}\downarrow}^{\dagger}-i(t_{1}-t_{3})c_{1d_{zx}\uparrow}^{\dagger}c_{1d_{xy}\downarrow}^{\dagger}$
$\displaystyle+2t_{2}c_{1d_{yz}\uparrow}^{\dagger}c_{1d_{zx}\uparrow}^{\dagger}-it_{2}c_{1d_{yz}\uparrow}^{\dagger}c_{1d_{xy}\downarrow}^{\dagger}$
$\displaystyle+t_{2}c_{1d_{zx}\uparrow}^{\dagger}c_{1d_{xy}\downarrow}^{\dagger}\Bigr{]}|0\rangle,$
(119) $\displaystyle T_{12}^{Z}|+,-\rangle=$
$\displaystyle\frac{1}{3}\Bigl{[}-(t_{1}-it_{2})c_{1d_{yz}\uparrow}^{\dagger}c_{1d_{yz}\downarrow}^{\dagger}-(t_{1}+it_{2})c_{1d_{zx}\uparrow}^{\dagger}c_{1d_{zx}\downarrow}^{\dagger}$
$\displaystyle-(it_{1}+t_{2})c_{1d_{yz}\uparrow}^{\dagger}c_{1d_{zx}\downarrow}^{\dagger}+(it_{1}-t_{2})c_{1d_{zx}\uparrow}^{\dagger}c_{1d_{yz}\downarrow}^{\dagger}$
$\displaystyle-(t_{1}-it_{2})c_{1d_{yz}\uparrow}^{\dagger}c_{1d_{xy}\uparrow}^{\dagger}+(it_{1}-t_{2})c_{1d_{zx}\uparrow}^{\dagger}c_{1d_{xy}\uparrow}^{\dagger}$
$\displaystyle-
t_{3}c_{1d_{yz}\downarrow}^{\dagger}c_{1d_{xy}\downarrow}^{\dagger}-it_{3}c_{1d_{zx\downarrow}}^{\dagger}c_{1d_{xy}\downarrow}^{\dagger}$
$\displaystyle-
t_{3}c_{1d_{xy}\uparrow}^{\dagger}c_{1d_{xy}\downarrow}^{\dagger}\Bigr{]}|0\rangle,$
(120) $\displaystyle T_{12}^{Z}|-,+\rangle=$
$\displaystyle\frac{1}{3}\Bigl{[}(t_{1}+it_{2})c_{1d_{yz}\uparrow}^{\dagger}c_{1d_{yz}\downarrow}^{\dagger}+(t_{1}-it_{2})c_{1d_{zx}\uparrow}^{\dagger}c_{1d_{zx}\downarrow}^{\dagger}$
$\displaystyle+(it_{1}+t_{2})c_{1d_{yz}\uparrow}^{\dagger}c_{1d_{zx}\downarrow}^{\dagger}-(it_{1}-t_{2})c_{1d_{zx}\uparrow}^{\dagger}c_{1d_{yz}\downarrow}^{\dagger}$
$\displaystyle+(t_{1}+it_{2})c_{1d_{yz}\downarrow}^{\dagger}c_{1d_{xy}\downarrow}^{\dagger}+(it_{1}+t_{2})c_{1d_{zx}\downarrow}^{\dagger}c_{1d_{xy}\downarrow}^{\dagger}$
$\displaystyle+t_{3}c_{1d_{yz}\uparrow}^{\dagger}c_{1d_{xy}\uparrow}^{\dagger}-it_{3}c_{1d_{zx\uparrow}}^{\dagger}c_{1d_{xy}\uparrow}^{\dagger}$
$\displaystyle+t_{3}c_{1d_{xy}\uparrow}^{\dagger}c_{1d_{xy}\downarrow}^{\dagger}\Bigr{]}|0\rangle.$
(121)

By using Eqs. (118)–(121) and Eqs. (26)–(40), we can calculate $\langle
i;\Gamma,g_{\Gamma}|T_{12}^{Z}|\textrm{i}\rangle$’s; as a result, the finite
terms are given by

$\displaystyle\langle i;A_{1}|T_{12}^{Z}|+,-\rangle$ $\displaystyle=-\langle
i;A_{1}|T_{12}^{Z}|-,+\rangle$
$\displaystyle=-\frac{1}{3\sqrt{3}}(2t_{1}+t_{3}),$ (122)
$\displaystyle\langle i;E,u|T_{12}^{Z}|+,-\rangle$ $\displaystyle=-\langle
i;E,u|T_{12}^{Z}|-,+\rangle$
$\displaystyle=-\frac{2}{3\sqrt{6}}(t_{1}-t_{3}),$ (123) $\displaystyle\langle
i;E,v|T_{12}^{Z}|+,-\rangle$ $\displaystyle=\langle
i;E,v|T_{12}^{Z}|-,+\rangle=\frac{2}{3\sqrt{2}}it_{2},$ (124)
$\displaystyle\langle i;T_{1},\alpha_{+}|T_{12}^{Z}|-,-\rangle$
$\displaystyle=\langle
i;T_{1},\alpha_{-}|T_{12}^{Z}|+,+\rangle=\frac{2}{3}t_{2},$ (125)
$\displaystyle\langle i;T_{1},\alpha|T_{12}^{Z}|+,-\rangle$
$\displaystyle=-\langle
i;T_{1},\alpha|T_{12}^{Z}|-,+\rangle=-\frac{2}{3\sqrt{2}}it_{1},$ (126)
$\displaystyle\langle i;T_{2},\alpha|T_{12}^{Z}|+,-\rangle$
$\displaystyle=-\langle
i;T_{2},\alpha|T_{12}^{Z}|-,+\rangle=-\frac{2}{3\sqrt{2}}t_{2},$ (127)
$\displaystyle\langle i;T_{1},\beta_{+}|T_{12}^{Z}|+,-\rangle$
$\displaystyle=\frac{1}{3}(it_{1}-t_{2}),$ (128) $\displaystyle\langle
i;T_{1},\beta_{+}|T_{12}^{Z}|-,+\rangle$ $\displaystyle=\langle
i;T_{1},\beta_{-}|T_{12}^{Z}|+,-\rangle=-\frac{i}{3}t_{3},$ (129)
$\displaystyle\langle i;T_{1},\beta_{-}|T_{12}^{Z}|-,+\rangle$
$\displaystyle=\frac{1}{3}(it_{1}+t_{2}),$ (130) $\displaystyle\langle
i;T_{1},\beta|T_{12}^{Z}|+,+\rangle$
$\displaystyle=-\frac{1}{3\sqrt{2}}[i(t_{1}-t_{3})+t_{2}],$ (131)
$\displaystyle\langle i;T_{1},\beta|T_{12}^{Z}|-,-\rangle$
$\displaystyle=-\frac{1}{3\sqrt{2}}[i(t_{1}-t_{3})-t_{2}],$ (132)
$\displaystyle\langle i;T_{2},\beta|T_{12}^{Z}|+,+\rangle$
$\displaystyle=\frac{1}{3\sqrt{2}}[i(t_{1}-t_{3})+t_{2}],$ (133)
$\displaystyle\langle i;T_{2},\beta|T_{12}^{Z}|-,-\rangle$
$\displaystyle=-\frac{1}{3\sqrt{2}}[i(t_{1}-t_{3})-t_{2}],$ (134)
$\displaystyle\langle i;T_{1},\gamma_{+}|T_{12}^{Z}|+,-\rangle$
$\displaystyle=\frac{1}{3}(t_{1}-it_{2}),$ (135) $\displaystyle\langle
i;T_{1},\gamma_{+}|T_{12}^{Z}|-,+\rangle$ $\displaystyle=-\langle
i;T_{1},\gamma_{-}|T_{12}^{Z}|+,-\rangle$ $\displaystyle=-\frac{1}{3}t_{3},$
(136) $\displaystyle\langle i;T_{1},\gamma_{-}|T_{12}^{Z}|-,+\rangle$
$\displaystyle=-\frac{1}{3}(t_{1}+it_{2}),$ (137) $\displaystyle\langle
i;T_{1},\gamma|T_{12}^{Z}|+,+\rangle$
$\displaystyle=\frac{1}{3\sqrt{2}}[(t_{1}-t_{3})+it_{2}],$ (138)
$\displaystyle\langle i;T_{1},\gamma|T_{12}^{Z}|-,-\rangle$
$\displaystyle=-\frac{1}{3\sqrt{2}}[(t_{1}-t_{3})-it_{2}],$ (139)
$\displaystyle\langle i;T_{2},\gamma|T_{12}^{Z}|+,+\rangle$
$\displaystyle=\frac{1}{3\sqrt{2}}[(t_{1}-t_{3})+it_{2}],$ (140)
$\displaystyle\langle i;T_{2},\gamma|T_{12}^{Z}|-,-\rangle$
$\displaystyle=\frac{1}{3\sqrt{2}}[(t_{1}-t_{3})-it_{2}],$ (141)

and the others are zero.

Then, we express Eq. (60) for the $Z$ bonds in terms of the exchange
interactions. By combining Eqs. (122)–(141) and Eqs. (22)–(25) with Eq. (60),
we can write $\bar{H}_{\textrm{eff}}$ for the $Z$ bonds,
$\bar{H}_{\textrm{eff}}^{Z}$, as follows:

$\displaystyle\bar{H}_{\textrm{eff}}^{Z}=H_{A_{1}}+H_{E}+H_{T_{1}}+H_{T_{2}},$
(142)

where

$\displaystyle H_{A_{1}}$ $\displaystyle=\sum_{\langle
i,j\rangle_{Z}}\sum_{n=-\infty}^{\infty}\frac{4(2t_{1}+t_{3})^{2}\mathcal{J}_{n}(u_{ij}^{Z})^{2}}{27(U+2J^{\prime}-n\omega)}\bm{S}_{i}\cdot\bm{S}_{j},$
(143) $\displaystyle H_{E}$ $\displaystyle=\sum_{\langle
i,j\rangle_{Z}}\sum_{n=-\infty}^{\infty}\frac{8\mathcal{J}_{n}(u_{ij}^{Z})^{2}}{27(U-J^{\prime}-n\omega)}$
$\displaystyle\times\\{[(t_{1}-t_{3})^{2}-3t_{2}^{2}]\bm{S}_{i}\cdot\bm{S}_{j}+6t_{2}^{2}S_{i}^{z}S_{j}^{z}\\},$
(144) $\displaystyle H_{T_{1}}$ $\displaystyle=\sum_{\langle
i,j\rangle_{Z}}\sum_{n=-\infty}^{\infty}\frac{4\mathcal{J}_{n}(u_{ij}^{Z})^{2}}{9(U^{\prime}-J_{\textrm{H}}-n\omega)}$
$\displaystyle\times\\{-[(t_{1}-t_{3})^{2}+t_{2}^{2}-2t_{1}^{2}-4t_{1}t_{3}]\bm{S}_{i}\cdot\bm{S}_{j}$
$\displaystyle+[(t_{1}-t_{3})^{2}+t_{2}^{2}](S_{i}^{y}S_{j}^{y}+S_{i}^{x}S_{j}^{x})$
$\displaystyle+2t_{2}(t_{1}-t_{3})(S_{i}^{x}S_{j}^{y}+S_{i}^{y}S_{j}^{x})$
$\displaystyle+2[(t_{1}^{2}+t_{2}^{2}+t_{3}^{2})-2t_{2}^{2}-2t_{1}t_{3}]S_{i}^{z}S_{j}^{z}\\},$
(145) $\displaystyle H_{T_{2}}$ $\displaystyle=\sum_{\langle
i,j\rangle_{Z}}\sum_{n=-\infty}^{\infty}\frac{4\mathcal{J}_{n}(u_{ij}^{Z})^{2}}{9(U^{\prime}+J_{\textrm{H}}-n\omega)}$
$\displaystyle\times\\{[(t_{1}-t_{3})^{2}+t_{2}^{2}](S_{i}^{x}S_{j}^{x}+S_{i}^{y}S_{j}^{y})$
$\displaystyle-2t_{2}(t_{1}-t_{3})(S_{i}^{x}S_{j}^{y}+S_{i}^{y}S_{j}^{x})$
$\displaystyle-[(t_{1}-t_{3})^{2}-t_{2}^{2}]\bm{S}_{i}\cdot\bm{S}_{j}\\}.$
(146)

In deriving these equations we have used the relations of operators in the
$j_{\textrm{eff}}=1/2$ subspace [e.g.,
$|-,+\rangle\langle-,+|=(\frac{1}{2}-S_{1}^{z})(\frac{1}{2}+S_{2}^{z})$ and
$|-,+\rangle\langle+,-|=S_{1}^{-}S_{2}^{+}$] and omitted the constant terms.
By using the identity
$S_{i}^{x}S_{j}^{x}+S_{i}^{y}S_{j}^{y}=\bm{S}_{i}\cdot\bm{S}_{j}-S_{i}^{z}S_{j}^{z}$,
we can rewrite Eqs. (145) and (146). A combination of the resultant equations
and Eqs. (142)–(144) gives Eq. (63).

##  Appendix D Estimates of the values of $J_{\delta}^{\textrm{3rd}}$ and
$K_{\delta}^{\textrm{3rd}}$

We estimate $J_{\delta}^{\textrm{3rd}}$ and $K_{\delta}^{\textrm{3rd}}$ at
$E_{0}=0$ in two cases. In the following estimation we calculate their values
for the $Z_{3}$ bonds (Fig. 1) because the bond anisotropy of the third-
neighbor hopping integrals is weak Valenti-PRB ; the values of the third-
neighbor hopping integrals used below are consistent with those of Ref.
Valenti-PRB, . Then, we set $J^{\prime}=J_{\textrm{H}}$,
$U^{\prime}=U-2J_{\textrm{H}}$, $U=3$ eV, and $J_{\textrm{H}}=0.5$ eV.

If we consider only $t_{\textrm{3rd}}$, the intraorbital hopping integral of
the $d_{xy}$ orbital, among the third-neighbor hopping integrals,
$J_{Z}^{\textrm{3rd}}$ and $K_{Z}^{\textrm{3rd}}$ at $E_{0}=0$ are given by

$\displaystyle
J_{Z}^{\textrm{3rd}}=\frac{4t_{\textrm{3rd}}^{2}}{27(U+2J_{\textrm{H}})}+\frac{8t_{\textrm{3rd}}^{2}}{27(U-J_{\textrm{H}})},$
(147)

and

$\displaystyle
K_{Z}^{\textrm{3rd}}=\frac{4}{9}t_{\textrm{3rd}}^{2}\Bigl{(}\frac{1}{U-3J_{\textrm{H}}}-\frac{1}{U-J_{\textrm{H}}}\Bigr{)},$
(148)

respectively. If we set $t_{\textrm{3rd}}\sim-40$ meV, we have
$J_{Z}^{\textrm{3rd}}\sim 0.25$ meV and $K_{Z}^{\textrm{3rd}}\sim 0.19$ meV.

Next, we consider not only $t_{\textrm{3rd}}$, but also two additional terms
of third-neighbor hopping integrals, $t_{\textrm{3rd}}^{\prime}$ and
$t_{\textrm{3rd}}^{\prime\prime}$. Here $t_{\textrm{3rd}}^{\prime}$ represents
the intraorbital hopping integral of the $d_{yz}$ or $d_{zx}$ orbital and
$t_{\textrm{3rd}}^{\prime\prime}$ represents the interorbital hopping integral
between these orbitals. Then, we can write $J_{Z}^{\textrm{3rd}}$ and
$K_{Z}^{\textrm{3rd}}$ at $E_{0}=0$ as follows:

$\displaystyle J_{Z}^{\textrm{3rd}}=$
$\displaystyle\frac{4(2t_{\textrm{3rd}}^{\prime}+t_{\textrm{3rd}})^{2}}{27(U+2J_{\textrm{H}})}+\frac{8(t_{\textrm{3rd}}^{\prime}-t_{\textrm{3rd}})^{2}}{27(U-J_{\textrm{H}})}$
$\displaystyle+\frac{8t_{\textrm{3rd}}^{\prime}(t_{\textrm{3rd}}^{\prime}+2t_{\textrm{3rd}})}{9(U-3J_{\textrm{H}})},$
(149) $\displaystyle K_{Z}^{\textrm{3rd}}=$
$\displaystyle\frac{4}{9}\Bigl{[}(t_{\textrm{3rd}}^{\prime}-t_{\textrm{3rd}})^{2}-3(t_{\textrm{3rd}}^{\prime\prime})^{2}\Bigr{]}$
$\displaystyle\times\Bigl{(}\frac{1}{U-3J_{\textrm{H}}}-\frac{1}{U-J_{\textrm{H}}}\Bigr{)}.$
(150)

Setting $t_{\textrm{3rd}}\sim-40$ meV, $t_{\textrm{3rd}}^{\prime}\sim-8$ meV,
and $t_{\textrm{3rd}}^{\prime\prime}\sim-7$ meV, we have
$J_{Z}^{\textrm{3rd}}\sim 0.65$ meV and $K_{Z}^{\textrm{3rd}}\sim 0.1$ meV.

Comparing the estimated values in the above two cases, we see
$J_{Z}^{\textrm{3rd}}$ and $K_{Z}^{\textrm{3rd}}$ are underestimated and
overestimated, respectively, if we consider only $t_{\textrm{3rd}}$ among the
third-neighbor hopping integrals. Since $K_{\delta}^{\textrm{3rd}}$ is much
smaller than $J_{\delta}^{\textrm{3rd}}$ even in a more realistic situation
Valenti-PRB , we consider only $J_{\delta}^{\textrm{3rd}}$ and neglect
$K_{\delta}^{\textrm{3rd}}$ in our analyses.

##  Appendix E Derivation of Eq. (88)

We rewrite Eq. (83) using Eqs. (86) and (87). First, we can rewrite the terms
of $J_{\delta}$’s as

$\displaystyle\sum_{\langle
i,j\rangle}J_{\delta}\langle\bm{S}_{i}\rangle\cdot\langle\bm{S}_{j}\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{i=1}^{N/2}\sum_{j=1}^{z_{\textrm{NN}}}J_{\delta}\langle\bm{S}_{i}\rangle\cdot\langle\bm{S}_{j}\rangle+\frac{1}{2}\sum_{j=1}^{N/2}\sum_{i=1}^{z_{\textrm{NN}}}J_{\delta}\langle\bm{S}_{i}\rangle\cdot\langle\bm{S}_{j}\rangle$
$\displaystyle=$
$\displaystyle\sum_{\bm{q}}\sum_{\mu=x,y,z}\Bigl{[}J_{1}(\bm{q})\langle
S_{-\bm{q}A}^{\mu}\rangle\langle
S_{\bm{q}B}^{\mu}\rangle+J_{1}(\bm{q})^{\ast}\langle
S_{-\bm{q}B}^{\mu}\rangle\langle S_{\bm{q}A}^{\mu}\rangle\Bigr{]},$ (151)

where

$\displaystyle J_{1}(\bm{q})=$
$\displaystyle\sum_{j=1}^{z_{\textrm{NN}}}\frac{J_{\delta}}{2}e^{-i\bm{q}\cdot(\bm{R}_{i}-\bm{R}_{j})}$
$\displaystyle=$
$\displaystyle\frac{J_{X}}{2}e^{-i\frac{q_{x}}{2}+i\frac{\sqrt{3}}{2}q_{y}}+\frac{J_{Y}}{2}e^{-i\frac{q_{x}}{2}-i\frac{\sqrt{3}}{2}q_{y}}+\frac{J_{Z}}{2}e^{iq_{x}},$
(152)

and $z_{\textrm{NN}}$ denotes the number of nearest-neighbor sites at a
certain cite on the honeycomb lattice. In Eq. (152) the first, second, and
third terms correspond to the contributions from the $X$, $Y$, and $Z$ bonds
(Fig. 1), respectively. Similarly, we can express the other terms as follows:

$\displaystyle\sum_{\langle i,j\rangle}K_{\delta}\langle
S_{i}^{\gamma}\rangle\langle S_{j}^{\gamma}\rangle$ $\displaystyle=$
$\displaystyle\sum_{\bm{q}}\Bigl{[}K_{x}(\bm{q})\langle
S_{-\bm{q}A}^{x}\rangle\langle S_{\bm{q}B}^{x}\rangle+K_{y}(\bm{q})\langle
S_{-\bm{q}A}^{y}\rangle\langle S_{\bm{q}B}^{y}\rangle$ $\displaystyle\ \ \ \
+K_{z}(\bm{q})\langle S_{-\bm{q}A}^{z}\rangle\langle
S_{\bm{q}B}^{z}\rangle\Bigr{]}$
$\displaystyle+\sum_{\bm{q}}\Bigl{[}K_{x}(\bm{q})^{\ast}\langle
S_{-\bm{q}B}^{x}\rangle\langle
S_{\bm{q}A}^{x}\rangle+K_{y}(\bm{q})^{\ast}\langle
S_{-\bm{q}B}^{y}\rangle\langle S_{\bm{q}A}^{y}\rangle$ $\displaystyle\ \ \ \
+K_{z}(\bm{q})^{\ast}\langle S_{-\bm{q}B}^{z}\rangle\langle
S_{\bm{q}A}^{z}\rangle\Bigr{]},$ (153) $\displaystyle\sum_{\langle
i,j\rangle}\Gamma_{\delta}(\langle S_{i}^{\alpha}\rangle\langle
S_{j}^{\beta}\rangle+\langle S_{i}^{\beta}\rangle\langle
S_{j}^{\alpha}\rangle)$ $\displaystyle=$
$\displaystyle\sum_{\bm{q}}\Bigl{[}\Gamma_{x}(\bm{q})(\langle
S_{-\bm{q}A}^{y}\rangle\langle S_{\bm{q}B}^{z}\rangle+\langle
S_{-\bm{q}A}^{z}\rangle\langle S_{\bm{q}B}^{y}\rangle)$ $\displaystyle\ \ \ \
+\Gamma_{y}(\bm{q})(\langle S_{-\bm{q}A}^{z}\rangle\langle
S_{\bm{q}B}^{x}\rangle+\langle S_{-\bm{q}A}^{x}\rangle\langle
S_{\bm{q}B}^{z}\rangle)$ $\displaystyle\ \ \ \ +\Gamma_{z}(\bm{q})(\langle
S_{-\bm{q}A}^{x}\rangle\langle S_{\bm{q}B}^{y}\rangle+\langle
S_{-\bm{q}A}^{y}\rangle\langle S_{\bm{q}B}^{x}\rangle)\Bigr{]}$
$\displaystyle+\sum_{\bm{q}}\Bigl{[}\Gamma_{x}(\bm{q})^{\ast}(\langle
S_{-\bm{q}B}^{y}\rangle\langle S_{\bm{q}A}^{z}\rangle+\langle
S_{-\bm{q}B}^{z}\rangle\langle S_{\bm{q}A}^{y}\rangle)$ $\displaystyle\ \ \ \
+\Gamma_{y}(\bm{q})^{\ast}(\langle S_{-\bm{q}B}^{z}\rangle\langle
S_{\bm{q}A}^{x}\rangle+\langle S_{-\bm{q}B}^{x}\rangle\langle
S_{\bm{q}A}^{z}\rangle)$ $\displaystyle\ \ \ \
+\Gamma_{z}(\bm{q})^{\ast}(\langle S_{-\bm{q}B}^{x}\rangle\langle
S_{\bm{q}A}^{y}\rangle+\langle S_{-\bm{q}B}^{y}\rangle\langle
S_{\bm{q}A}^{x}\rangle)\Bigr{]},$ (154)
$\displaystyle\sum_{\langle\langle\langle
i,j\rangle\rangle\rangle}J_{\delta}^{\textrm{3rd}}\langle\bm{S}_{i}\rangle\cdot\langle\bm{S}_{j}\rangle$
$\displaystyle=$
$\displaystyle\sum_{\bm{q}}\sum_{\mu=x,y,z}\Bigl{[}J_{3}(\bm{q})\langle
S_{-\bm{q}A}^{\mu}\rangle\langle
S_{\bm{q}B}^{\mu}\rangle+J_{3}(\bm{q})^{\ast}\langle
S_{-\bm{q}B}^{\mu}\rangle\langle S_{\bm{q}A}^{\mu}\rangle\Bigr{]},$ (155)

where $K_{\mu}(\bm{q})$’s and $\Gamma_{\mu}(\bm{q})$’s have been defined in
Eq. (94)–(99), and $J_{3}(\bm{q})$ is given by

$\displaystyle J_{3}(\bm{q})=$
$\displaystyle\frac{J_{X}^{\textrm{3rd}}}{2}e^{iq_{x}-i\sqrt{3}q_{y}}+\frac{J_{Y}^{\textrm{3rd}}}{2}e^{iq_{x}+i\sqrt{3}q_{y}}+\frac{J_{Z}^{\textrm{3rd}}}{2}e^{-2iq_{x}}.$
(156)

By combining Eqs. (151)–(156) and setting
$J(\bm{q})=J_{1}(\bm{q})+J_{3}(\bm{q})$, we obtain Eq. (88).

##  Appendix F Some remarks on momentum

We remark on momentum in the case of the honeycomb lattice. The remarks are
about its expression and the reciprocal lattice vector with the periodic
boundary condition. In the case of the honeycomb lattice (Fig. 1) a set of
primitive vectors, $\bm{a}_{1}$ and $\bm{a}_{2}$, can be written as

$\displaystyle\bm{a}_{1}={}^{t}\Bigl{(}\frac{\sqrt{3}}{2}a_{\textrm{2nd}}\ \
\frac{1}{2}a_{\textrm{2nd}}\Bigr{)}={}^{t}\Bigl{(}\frac{3}{2}\ \
\frac{\sqrt{3}}{2}\Bigr{)},$ (157) $\displaystyle\bm{a}_{2}={}^{t}\Bigl{(}0\ \
a_{\textrm{2nd}}\Bigr{)}={}^{t}\Bigl{(}0\ \ \sqrt{3}\Bigr{)},$ (158)

and thus the primitive vectors for the reciprocal lattice, $\bm{b}_{1}$ and
$\bm{b}_{2}$, are given by

$\displaystyle\bm{b}_{1}={}^{t}\Bigl{(}\frac{4\pi}{3}\ \ 0\Bigr{)},$ (159)
$\displaystyle\bm{b}_{2}={}^{t}\Bigl{(}-\frac{2\pi}{3}\ \
\frac{2\pi}{\sqrt{3}}\Bigr{)}.$ (160)

(As described in the caption of Fig. 1, we can represent the honeycomb lattice
as a triangular Bravais lattice with a two-sublattice structure Ashcro-Merm .)
By imposing the periodic boundary condition Ashcro-Merm , we can express
momentum $\bm{q}$ as

$\displaystyle\bm{q}=\frac{l_{1}}{N_{1}}\bm{b}_{1}+\frac{l_{2}}{N_{2}}\bm{b}_{2},$
(161)

where the integers $l_{1}$ and $l_{2}$ satisfy $0\leq l_{1}<N_{1}$ and $0\leq
l_{2}<N_{2}$ with $N_{1}N_{2}=N$. As described in Sec. IV B, we set
$N_{1}=N_{2}=120$ in our analyses. The values of $l_{1}$ and $l_{2}$ for the
magnetic states considered in Sec. IV B are given as follows: for
$\bm{q}=\bm{Q}_{\textrm{ZZ-}X}$ $l_{1}=N_{1}/4$ and $l_{2}=N_{1}/2$; for
$\bm{q}=\bm{Q}_{\textrm{ZZ-}Z}$ $l_{1}=N_{1}/2$ and $l_{2}=0$; and for
$\bm{q}=\bm{Q}_{120}$ $l_{1}=l_{2}=(2N_{1}/3)$. Then, the reciprocal lattice
vector $\bm{K}$ can be written as

$\displaystyle\bm{K}={}^{t}(K_{x}\ \ K_{y})=m_{1}\bm{b}_{1}+m_{2}\bm{b}_{2},$
(162)

where

$\displaystyle K_{x}=\frac{2\pi}{3}(2m_{1}-m_{2}),\
K_{y}=\frac{2\pi}{\sqrt{3}}m_{2},$ (163)

and $m_{1}$ and $m_{2}$ are integers. From these equations we see that
$\bm{K}$ for $(m_{1},m_{2})=(0,1)$, $\bm{K}_{01}$, equals
$2\bm{Q}_{\textrm{ZZ-}X}$, and that $\bm{K}$ for $(m_{1},m_{2})=(1,0)$,
$\bm{K}_{10}$, equals $2\bm{Q}_{\textrm{ZZ-}Z}$. Namely, we have

$\displaystyle\bm{Q}_{\textrm{ZZ-}X}=\frac{\bm{K}_{01}}{2},$ (164)
$\displaystyle\bm{Q}_{\textrm{ZZ-}Z}=\frac{\bm{K}_{10}}{2}.$ (165)

In contrast to $\bm{Q}_{\textrm{ZZ-}X}$ and $\bm{Q}_{\textrm{ZZ-}Z}$,
$\bm{Q}_{120}$ is not equal to $\bm{K}/2$ for any allowed $m_{1}$ and $m_{2}$.

## References

  * (1) J. H. Shirley, Solution of the Schrödinger Equation with a Hamiltonian Periodic in Time, Phys. Rev. 138, B979 (1965).
  * (2) H. Sambe, Steady States and Quasienergies of a Quantum-Mechanical System in an Oscillating Field, Phys. Rev. A 7, 2203 (1973).
  * (3) M. Bukov, L. D’Alessio, and A. Polkovnikov, Universal High-Frequency Behavior of Periodically Driven Systems: from Dynamical Stabilization to Floquet Engineering, Advances in Physics 64, 139 (2015).
  * (4) D. A. Abanin, W. De Roeck, and F. Huveneers, Exponentially Slow Heating in Periodically Driven Many-Body Systems, Phys. Rev. Lett. 115, 256803 (2015).
  * (5) T. Kuwahara, T. Mori, and K. Saito, Floquet-Magnus Theory and Generic Transient Dynamics in Periodically Driven Many-Body Quantum Systems, Ann. Phys. (Amsterdam, Neth.) 367, 96 (2016).
  * (6) J. H. Mentink, K. Balzer, and M. Eckstein, Ultrafast and reversible control of the exchange interaction in Mott insulators, Nat. Commun. 6, 6708 (2015).
  * (7) B. J. Kim, H. Jin, S. J. Moon, J.-Y. Kim, B.-G. Park, C. S. Leem, J. Yu, T. W. Noh, C. Kim, S.-J. Oh, J.-H. Park, V. Durairaj, G. Cao, and E. Rotenberg, Novel $J_{\textrm{eff}}=1/2$ Mott State Induced by Relativistic Spin-Orbit Coupling in Sr2IrO4, Phys. Rev. Lett. 101, 076402 (2008).
  * (8) W. Witczak-Krempa, G. Chen, Y. B. Kim, and L. Balents, Correlated Quantum Phenomena in the Strong Spin-Orbit Regime, Annu. Rev. Condens. Matter Phys. 5, 57 (2014).
  * (9) T. Takayama, J. Chaloupka, A. Smerald, G. Khaliullin, and H. Takagi, Spin-orbit-entangled electronic phases in $4d$ and $5d$ transition-metal compounds, J. Phys. Soc. Jpn. 90, 062001 (2021).
  * (10) G. Jackeli and G. Khaliullin, Mott Insulators in the Strong Spin-Orbit Coupling Limit: From Heisenberg to a Quantum Compass and Kitaev Models, Phys. Rev. Lett. 102, 017205 (2009).
  * (11) J. Chaloupka, G. Jackeli, and G. Khaliullin, Kitaev-Heisenberg Model on a Honeycomb Lattice: Possible Exotic Phases in Iridium Oxides A2IrO3, Phys. Rev. Lett. 105, 027204 (2010).
  * (12) J. G. Rau, E. K.-H. Lee, and H.-Y. Kee, Generic Spin Model for the Honeycomb Iridates beyond the Kitaev Limit, Phys. Rev. Lett. 112, 077204 (2014).
  * (13) S. M. Winter, Y. Li, H. O. Jeschke, and R. Valentí, Challenges in design of Kitaev materials: Magnetic interactions from competing energy scales, Phys. Rev. B 93, 214431 (2016).
  * (14) N. Arakawa, Vector chirality for effective total momentum $J_{\textrm{eff}}$ in a nonfrustrated Mott insulator: Effects of strong spin-orbit coupling and broken inversion symmetry, Phys. Rev. B 94, 174416 (2016).
  * (15) J. G. Rau and H.-Y. Kee, Trigonal distortion in the honeycomb iridates: Proximity of zigzag and spiral phases in Na2IrO3, arXiv:1408.4811 (unpublished).
  * (16) J. Chaloupka and G. Khaliullin, Hidden symmetries of the extended Kitaev-Heisenberg model: Implications for the honeycomb-lattice iridates A2IrO3, Phys. Rev. B 92, 024413 (2015).
  * (17) L. Janssen, E. C. Andrade, and M. Vojta, Magnetization processes of zigzag states on the honeycomb lattice: Identifying spin models for $\alpha$-RuCl3 and Na2IrO3, Phys. Rev. B 96, 064430 (2017).
  * (18) J. Wang, B. Normand, and Z.-X. Liu, One Proximate Kitaev Spin Liquid in the $K$-$J$-$\Gamma$ Model on the Honeycomb Lattice, Phys. Rev. Lett. 123, 197201 (2019).
  * (19) J. Rusnačko, D. Gotfryd, and J. Chaloupka, Kitaev-like honeycomb magnets: Global phase behavior and emergent effective models, Phys. Rev. B 99, 064425 (2019).
  * (20) N. Arakawa and K. Yonemitsu, Floquet engineering of Mott insulators with strong spin-orbit coupling, Phys. Rev. B 103, L100408 (2021).
  * (21) A. Kirilyuk, A. V. Kimel, and T. Rasing, Ultrafast optical manipulation of magnetic order, Rev. Mod. Phys. 82, 2731 (2010).
  * (22) K. Hejazi, J. Liu, and L. Balents, Floquet spin and spin-orbital Hamiltonians and doublon-holon generations in periodically driven Mott insulators, Phys. Rev. B 99, 205111 (2019).
  * (23) N. W. Ashcroft and N. D. Mermin, Solid State Physics (Thomson Learning Inc., USA, 1976).
  * (24) J. Kanamori, Electron Correlation and Ferromagnetism of Transition Metals, Prog. Theor. Phys. 30, 275 (1963).
  * (25) K. W. Plumb, J. P. Clancy, L. J. Sandilands, V. V. Shankar, Y. F. Hu, K. S. Burch, H.-Y. Kee, and Y.-J. Kim, $\alpha$-RuCl3: A spin-orbit assisted Mott insulator on a honeycomb lattice, Phys. Rev. B 90, 041112(R) (2014).
  * (26) H.-S. Kim, V. Shankar V., A. Catuneanu, and H.-Y. Kee, Kitaev magnetism in honeycomb RuCl3 with intermediate spin-orbit coupling, Phys. Rev. B 91, 241110(R) (2015).
  * (27) S. Ishihara, T. Hatakeyama, and S. Maekawa, Magnetic Ordering, Orbital Ordering and Resonant X-ray Scattering in Perovskite Titanates, Phys. Rev. B 65, 064442 (2002).
  * (28) J. M. Luttinger and L. Tisza, Theory of Dipole Interaction in Crystals, Phys. Rev. 70, 954 (1946).
  * (29) J. M. Luttinger, A Note on the Ground State in Antiferromagnetics, Phys. Rev. 81, 1015 (1951).
  * (30) For the analysis in the region surronded by the dashed white line in Fig. 2(a) of Ref. Rau-PRL, , the authors imposed an additional restriction.
  * (31) H. Suzuki, H. Liu, J. Bertinshaw, K. Ueda, H. Kim, S. Laha, D. Weber, Z. Yang, L. Wang, H. Takahashi, K. Fur̈sich, M. Minola, B. V. Lotsch, B. J. Kim, H. Yavasş, M. Daghofer, J. Chaloupka, G. Khaliullin, H. Gretarsson, and B. Keimer, Proximate ferromagnetic state in the Kitaev model material $\alpha$-RuCl3, Nat. Comm. 12, 4512 (2021).
  * (32) J. A. Sears, L. E. Chern, S. Kim, P. J. Bereciartua, S. Francoual, Y. B. Kim, and Y.-J. Kim, Ferromagnetic Kitaev interaction and the origin of large magnetic anisotropy in $\alpha$-RuCl3, Nat. Phys. 16, 837 (2020).
  * (33) L. D’Alessio and M. Rigol, Long-time Behavior of Isolated Periodically Driven Interacting Lattice Systems, Phys. Rev. X 4, 041048 (2014).
  * (34) M. Claassen, H.-C. Jiang, B. Moritz, and T. P. Devereaux, Dynamical time-reversal symmetry breaking and photo-induced chiral spin liquids in frustrated Mott insulators, Nat. Commun. 8, 1192 (2017).
  * (35) S. Kitamura, T. Oka, and H. Aoki, Probing and controlling spin chirality in Mott insulators by circularly polarized laser, Phys. Rev. B 96, 014406 (2017).
  * (36) R. D. Johnson, S. C. Williams, A. A. Haghighirad, J. Singleton, V. Zapf, P. Manuel, I. I. Mazin, Y. Li, H. O. Jeschke, R. Valentí, and R. Coldea, Monoclinic crystal structure of $\alpha$-RuCl3 and the zigzag antiferromagnetic ground state, Phys. Rev. B 92, 235119 (2015).
  * (37) Y. Kawakami, H. Itoh, K. Yonemitsu, and S. Iwai, Strong light-field effects driven by nearly single-cycle 7fs light-field in correlated organic conductors, J. Phys. B: At. Mol. Opt. Phys. 51, 174005 (2018).
  * (38) A. Sriram and M. Claassen, Light-Induced Control of Magnetic Phases in Kitaev Quantum Magnets, arXiv:2105.01062 (unpublished).
  * (39) A. Kitaev, Anyons in an exactly solved model and beyond, Annals of Physics 321, 2 (2006).
  * (40) M. O. Takahashi, M. G. Yamada, D. Takikawa, T. Mizushima, and S. Fujimoto, Topological nematic phase transition in Kitaev magnets under applied magnetic fields, Phys. Rev. Research 3, 023189 (2021).
  * (41) J. Zheng, K. Ran, T. Li, J. Wang, P. Wang, B. Liu, Z.-X. Liu, B. Normand, J. Wen, and W. Yu, Gapless Spin Excitations in the Field-Induced Quantum Spin Liquid Phase of $\alpha$-RuCl3, Phys. Rev. Lett. 119, 227208 (2017).
  * (42) I. Žutić, J. Fabian, and S. D. Sarma, Spintronics: Fundamentals and applications, Rev. Mod. Phys. 76, 323 (2004).