#  Bulk photovoltaic effect in antiferromagnet: Role of collective spin
dynamics

Junta Iguchi<EMAIL_ADDRESS>Department of Applied
Physics, The University of Tokyo, Bunkyo, Tokyo 113-8656, Japan  Hikaru
Watanabe<EMAIL_ADDRESS>Research Center for Advanced
Science and Technology, The University of Tokyo, Meguro, Tokyo 153-8904, Japan
Yuta Murakami  Center for Emergent Matter Science, RIKEN, Wako, Saitama
351-0198, Japan  Takuya Nomoto  Research Center for Advanced Science and
Technology, The University of Tokyo, Meguro, Tokyo 153-8904, Japan  Ryotaro
Arita  Research Center for Advanced Science and Technology, The University of
Tokyo, Meguro, Tokyo 153-8904, Japan  Center for Emergent Matter Science,
RIKEN, Wako, Saitama 351-0198, Japan

###### Abstract

Inspired by recent advancements in the bulk photovoltaic effect which can
extend beyond the independent particle approximation (IPA), this study delves
into the influence of collective spin dynamics in an antiferromagnetic on
photocurrent generation using a time domain calculation. In the linear and
photocurrent conductivity spectra, we observe peaks below the bandgap regime,
attributed to the resonant contributions of collective modes, alongside
broadband modifications resulting from off-resonant spin dynamics. Notably,
the emergence of spin dynamics allows various types of photocurrent, which are
absent in the IPA framework. Furthermore, we emphasize the importance of
energy scale proximity between electronic and spin degrees of freedom in
enabling efficient feedback between them. These findings offer new avenues for
efficient energy harvesting and optoelectronic applications.

##  I introduction

The photovoltaic effect is a phenomenon in which light is converted into an
electric current known as a photocurrent. It can be understood as an even-
order optical response, hence requires inversion breaking in the system. This
effect is commonly observed in heterostructures such as p-n junctions [1],
where inversion symmetry is artificially broken. Recently, the growing demand
for energy harvesting has provoked the exploration for photocurrent generation
in bulk single-phase crystals without inversion symmetry, termed the bulk
photovoltaic effect (BPVE) [2, 3]. BPVE has attracted much attention due to
its distinct features. For instance, unlike conventional PVE, photovoltage is
not limited by the band gap energy, and the Shockley-Queisser limit [4] could
be overcome [5].

Extensive studies have revealed the complete classification of BPVE in non-
magnetic and magnetic systems based on the independent particle approximation
(IPA) [6, 7, 8]. For a representative example, BPVE in insulators can be
classified into two contributions: shift current and injection current. The
shift current [9, 10, 11, 12] arises from the difference in intracell
coordinates of excited carriers and is related to the geometric quantity
called Berry connection. On the other hand, injection current results from the
group velocity difference of the excited electron-hole pair. Based on this
classification, extensive searches for efficient BPVE have been carried out
across a wide range of materials, including Weyl semimetals [13], topological
insulators [14], and transition metal dichalcogenides [15].

This classification is based on the IPA with the rigid band picture. However,
recent progress in theoretical studies of BPVE has revealed that electron
correlation effects in solids also play an important role in photocurrent
generation [16, 17]. In particular, photocurrent can originate from elementary
excitations below the bandgap, such as excitons [18, 19], soft phonons in
ferroelectrics [20], electromagnons in multiferroics [21], and collective
modes in excitonic insulators [22]. Understanding these phenomena requires
explicit consideration of collective dynamics beyond the IPA.

In this paper, we focus on the BPVE in an antiferromagnetic system and discuss
the effects of the collective spin dynamics, exploring BPVE beyond IPA. BPVE
in magnetic materials is particularly intriguing for several reasons. Firstly,
there are various kinds of magnetic textures, in some of which the magnetic
order breaks the inversion symmetry and thereby ties the charge degree of
freedom with spins. This strong spin-charge coupling is also attracting much
attention from the field of multiferroics [23, 24, 25] and antiferromagnetic
spintronics [26, 27, 28, 29]. Secondly, the spin textures accompany the
breaking of time-reversal symmetry, which may allow the nonzero contribution
from injection current under linearly polarized and unpolarized light. This is
advantageous because the amplitude of the linear injection can dominate that
of the shift currents, particularly in clean systems, potentially enabling
highly efficient photovoltaic devices. Thirdly, not limited to the
antiferromagnetic system, the external magnetic field can tune the magnetic
structures, and we can expect a significant modulation of photocurrent
response. In fact, a recent study on the BPVE effect in magnetic systems
within the IPA shows that photocurrent direction can be controlled by external
field [30, 31], demonstrating the high tunability of the photocurrent in
magnets [32]. As illustrated, BPVE phenomena in magnetic systems exhibit rich
characteristics, even within the IPA framework. In such systems, the role of
spin dynamics on BPVE remains largely unexplored. Especially in
antiferromagnetic systems, the energy scale of the collective spin dynamics
typically lies in the THz regime, which is proximate to that of the electronic
system. Hence, we expect that the spin dynamics significantly affects the
response of the electronic systems. Moreover, the emergence of fictitious
fields by the spin dynamics may lead to richer optoelectronic properties,
which are absent within IPA [22]. By investigating the effects of spin
dynamics on PVE, we aim to provide a more comprehensive understanding of BPVE,
moving beyond the IPA.

In this paper, we investigate the role of collective spin dynamics in BPVE in
a simple model representing a one-dimensional antiferromagnet. We employ a
unified framework for real-time simulations of conduction electrons and
localized spin moments based on the methods introduced in previous studies
[22, 33]. This approach allows us to investigate the influence of collective
modes on BPVE and gives a comprehensive understanding of their impact on the
optical response. By employing the symmetry analysis, we identify an optically
active collective mode and investigate its role in both linear and nonlinear
optical responses. In the conductivity spectrum, we observe a peak below the
bandgap regime stemming from resonant contributions of collective modes, as
well as broadband modifications resulting from off-resonant contributions of
collective spin dynamics. Additionally, we found that photocurrent arising
from the spin dynamics can be classified into several types, which are absent
in the IPA. Here, the symmetry and the phase degrees of freedom of the driven
forces play crucial roles. Furthermore, we demonstrate the significance of
proximity in energy scales between collective modes and electronic
excitations. This insight provides a guiding principle for harnessing
substantial photocurrents from spin dynamics.

##  II method

In this section, we elaborate on the model and explain how we incorporate the
effect of collective spin dynamics on optical responses. Firstly, we introduce
the model used in this study, then move on to our time-dependent simulation
scheme.

###  II.1 Model

Figure 1: (a) One-dimensional chain model with canted antiferromagnetic order.
(b) band dispersion of the system with the parameters
$t_{h}=1,\lambda=0.8,J=0.6,K_{z}=0.2,h_{y}=0.2$.

This study focuses on a locally noncentrosymmetric system in one dimension,
where conduction electrons are coupled to localized spin moments with a canted
antiferromagnetic structure (see FIG. 1). This model is a simple
representative example of magnetoelectric materials [34], where we expect
large spin and charge coupling.

The Hamiltonian of the model is expressed as

$\displaystyle\hat{\mathcal{H}}=\hat{\mathcal{H}}_{\mathrm{ele}}+\hat{\mathcal{H}}_{\mathrm{exc}}+\mathcal{H}_{\mathrm{spin}}+\hat{\mathcal{H}}_{E}.$
(1)

The first term,

$\displaystyle\hat{\mathcal{H}}_{\mathrm{ele}}=-2t_{h}\sum_{k}\sum_{\sigma}\cos
k\quantity[\hat{c}_{A\sigma}^{\dagger}(k)\hat{c}_{B\sigma}(k)+\hat{c}_{B\sigma}^{\dagger}(k)\hat{c}_{A\sigma}(k)]$
$\displaystyle-\lambda\sum_{k}\sum_{\sigma\sigma^{\prime}}\quantity(\mathbf{g}(k)\cdot{\bf\it\sigma})^{\sigma\sigma^{\prime}}\quantity[\hat{c}_{A\sigma}^{\dagger}(k)\hat{c}_{A\sigma^{\prime}}(k)-\hat{c}_{B\sigma}^{\dagger}(k)\hat{c}_{B\sigma^{\prime}}(k)]$
(2)

is the Hamiltonian of the electronic system, consisting of the nearest
neighbor hopping $t_{h}$, and the sublattice dependent antisymmetric spin-
orbit coupling $\lambda$ (sASOC). $\hat{c}_{\alpha\sigma}^{\dagger}(k)$
($\hat{c}_{\alpha\sigma}(k)$) is the creation (annihilation) operator of the
electron on sublattice $\alpha$ ($\alpha=A,B$) having spin $\sigma$
($\sigma=\uparrow,\downarrow$), and $k$ is the wave vector. $\mathbf{g}(k)$ is
called g-vector, and satisfies $\mathbf{g}(-k)=-\mathbf{g}(k)$. We set
$\mathbf{g}(k)=(0,0,\sin 2k)$, implicitly assuming the presence of ligand
structure around each sublattice which breaks the local inversion symmetry.
This sASOC is the essential ingredient to enhance magnetoelectric coupling
[34, 35, 36].

The second term

$\displaystyle\hat{\mathcal{H}}_{\mathrm{exc}}=-J\sum_{k}\sum_{\alpha}\sum_{\sigma\sigma^{\prime}}\hat{c}_{\alpha\sigma}^{\dagger}(k)\quantity({\bf\it\sigma}\cdot\mathbf{S}_{\alpha}(t))^{\sigma\sigma^{\prime}}\hat{c}_{\alpha\sigma^{\prime}}(k)$
(3)

corresponds to the interaction between electronic degrees of freedom and
localized spin moments with coupling constant $J$. Here we assume that the
localized spins $\mathrm{S}_{\alpha}$ are classical spins and set
$\absolutevalue{\mathbf{S}_{\alpha}}=1$. We note that, in Eq. (3), we assume
that the spin configuration is commensurate with the lattice structure as in
FIG. 1 and the local spins can be classified only by the sublattices.

The third term

$\displaystyle\mathcal{H}_{\mathrm{spin}}=-\sum_{\alpha}\quantity[K_{z}\quantity(\mathrm{S}_{\alpha}^{z})^{2}-h_{y}\mathrm{S}_{\alpha}^{y}]$
(4)

is the Hamiltonian for the localized spins. To stabilize the antiferromagnetic
order along the chain direction ($z$ direction), we consider the uniaxial
anisotropy $K_{z}$. We also apply the external magnetic field perpendicular to
the chain direction ($y$ direction) to induce the canted configurations. This
is because, without this canted moment, the light-induced spin accumulation on
each sublattice, which acts as spin transfer torque for the localized spin
moments, would be parallel to $\mathbf{S}_{\alpha}$. Therefore, there would be
no spin dynamics induced by light irradiation.

The last term

$\displaystyle\hat{\mathcal{H}}_{E}=-E^{z}(t)\sum_{kk^{\prime}}\sum_{\alpha}\sum_{\sigma}\quantity[i\partialderivative{k}\delta(k-k^{\prime})]\hat{c}_{k\alpha\sigma}^{\dagger}\hat{c}_{k^{\prime}\alpha\sigma}$
(5)

describes the light-matter coupling in the length gauge, where $E^{z}(t)$ is a
time-dependent electric field along the chain direction, where we set lattice
constant $a=1$, and elementary charge $e=1$. In the expression Eq. (5), we
assume that the Wannier state of a conducting electron is well-localized at a
given site and ignore the light-matter coupling originating from its spatial
extension. This light-matter coupling is the same as the celebrated Peierls
substitution in the dipolar gauge [37]. The difficulty of dealing with the
$k$-derivative of the delta function in Eq. (5) can be resolved in the next
subsection.

###  II.2 Calculation Scheme

Here, we introduce two coupled equations that describe the dynamics of
electronic degrees of freedom and localized spin moment. Firstly, the time
evolution of the conduction electrons can be described by the single-particle
density matrix (SPDM)
$\rho_{\alpha\beta}^{\sigma\sigma^{\prime}}(k)=\left\langle{\hat{c}_{\beta\sigma^{\prime}}^{\dagger}(k)\hat{c}_{\alpha\sigma}(k)}\right\rangle$,
and SPDM satisfies the following equation called the von Neumann equation
[38],

$\displaystyle\begin{split}\partialderivative{{\bf\it\rho}(k,t)}{t}=-i\quantity[{\bf\it
H}(k,t),\,{\bf\it\rho}(k,t)]-E^{z}(t)\partialderivative{{\bf\it\rho}(k,t)}{k}\\\
-\gamma({\bf\it\rho}(k,t)-{\bf\it\rho}_{\text{eq}}(k)).\end{split}$ (6)

Secondly, the time evolution of the localized spin system is described by the
Landau-Lifshitz-Gilbert (LLG) equation.

$\displaystyle\dfrac{d\bm{\mathrm{S}}_{\alpha}}{dt}$
$\displaystyle=\dfrac{1}{1+\alpha_{G}^{2}}\left(\bm{\mathrm{H}}_{\alpha}^{\mathrm{eff}}\times\bm{\mathrm{S}}_{\alpha}+\alpha_{G}\bm{\mathrm{S}}_{\alpha}\times\left(\bm{\mathrm{S}}_{\alpha}\times\bm{\mathrm{H}}_{\alpha}^{\mathrm{eff}}\right)\right),$
(7) $\displaystyle\bm{\mathrm{H}}_{\alpha}^{\mathrm{eff}}$
$\displaystyle=-J\langle\bm{\sigma}_{\alpha}\rangle+\dfrac{\delta\mathcal{H}_{\text{spin}}}{\delta\bm{\mathrm{S}}_{\alpha}}.$
(8)

In Eq. (6)$,{\bf\it H}(k,t)$ is the time-dependent electronic Hamiltonian at
each $k$ point defined as follows

$\displaystyle\hat{\mathcal{H}}_{\mathrm{ele}}+\hat{\mathcal{H}}_{\mathrm{exc}}=\sum_{k}\sum_{\alpha\beta}\sum_{\sigma\sigma^{\prime}}\quantity({\bf\it
H}(k,t))_{\alpha\beta}^{\sigma\sigma^{\prime}}\hat{c}_{\alpha\sigma}^{\dagger}(k)\hat{c}_{\beta\sigma^{\prime}}(k).$
(9)

The $k$-derivative of the delta function in Eq. (5) becomes the $k$-derivative
of the SPDM, which is computationally manageable. In the LLG equation Eq. (7),
$\left\langle{{\bf\it\sigma}_{\alpha}}\right\rangle$ is the sublattice-
dependent spin density of itinerant electrons, and this can be calculated from
SPDM.

Additionally, real materials exhibit relaxation of excited carriers due to
electron-electron correlations, electron-phonon interactions, and impurity
scattering. To ensure a physically reasonable response to light, we account
for these effects phenomenologically by using the relaxation time
approximation in the von Neumann equation as
$\gamma({\bf\it\rho}(k,t)-{\bf\it\rho}_{\text{eq}}(k))$ in Eq. (6), and the
Gilbert damping $\alpha_{G}$ in Eq. (7). Here ${\bf\it\rho}_{\text{eq}}(k)$ is
the SPDM in the equilibrium at the temperature $T=0$ as shown below.

The SPDM at equilibrium ${\bf\it\rho}_{\text{eq}}(k)$ represents the SPDM in
the initial state at temperature $T=0$. SPDM in the band basis
$\tilde{{\bf\it{\rho}}}_{\text{eq}}(k)$ is diagonal, and the diagonal
component represents the occupation number $\Theta(\mu-\epsilon_{nk})$, with
$\epsilon_{nk}$ being the eigenvalue of the Hamiltonian ${\bf\it H}(k)$.
Therefore $\tilde{{\bf\it\rho}}_{\text{eq}}(k)$ is written as,

$\displaystyle\quantity(\tilde{{\bf\it{\rho}}}_{\text{eq}}(k))_{nn^{\prime}}=\delta_{nn^{\prime}}\Theta(\mu-\epsilon_{nk}).$
(10)

By using this expression, the SPDM in the original basis
${\bf\it\rho}_{\text{eq}}(k)$ can be calculated by

$\displaystyle{\bf\it\rho}_{\text{eq}}(k)={\bf\it
U}(k)\tilde{{\bf\it\rho}}_{\text{eq}}(k){\bf\it U}^{\dagger}(k)$ (11)

with ${\bf\it U}$ is a unitary matrix which diagonalizes the Hamiltonian as,

$\displaystyle{\bf\it U}^{\dagger}(k){\bf\it H}(k){\bf\it
U}(k)={\bf\it\mathcal{E}}(k)$ (12)
$\displaystyle\quantity({\bf\it\mathcal{E}}(k))_{nn^{\prime}}=\delta_{nn^{\prime}}\epsilon_{nk}.$
(13)

To obtain initial spin configuration, where the localized spin moment
$\mathbf{S}_{\alpha}(t=0)$ is parallel to the effective field
$\mathbf{H}_{\alpha}^{\mathrm{eff}}$ acting on $\mathbf{S}_{\alpha}(t=0)$, we
perform the self-consistent calculation.

In each time step, the current density is evaluated from SPDM as follows. The
current operator in the length gauge is expressed as

$\displaystyle\hat{J}(t)$
$\displaystyle=\sum_{k}\sum_{\alpha\beta}\sum_{\sigma\sigma^{\prime}}\partialderivative{\quantity({\bf\it
H}(k,t))_{\alpha\beta}^{\sigma\sigma^{\prime}}}{k}\hat{c}_{\alpha\sigma}^{\dagger}(k)\hat{c}_{\beta\sigma^{\prime}}(k)$
(14)
$\displaystyle\equiv\sum_{k}\sum_{\alpha\beta}\sum_{\sigma\sigma^{\prime}}\quantity({\bf\it
J}(k))_{\alpha\beta}^{\sigma\sigma^{\prime}}\hat{c}_{\alpha\sigma}^{\dagger}(k)\hat{c}_{\beta\sigma^{\prime}}(k).$
(15)

It is noteworthy that since the $\hat{\mathcal{H}}_{\text{exc}}$ term in the
Hamiltonian does not depend on momentum $k$, the current operator is
independent of local spin dynamics.

We solve the coupled equations Eq. (6) and Eq. (7) by the fourth-order Runge-
Kutta method. In each time step, we calculate the sublattice-dependent spin
density using the SPDM and update the exchange Hamiltonian with the newly
obtained spin configurations. Additionally, the expectation value of the
current operator is obtained by using the density matrix as

$\displaystyle J^{z}(t)=\sum_{k}\operatorname{Tr}\quantity[{\bf\it
J}(k){\bf\it\rho}(k,t)].$ (16)

To evaluate the $\partial{\bf\it\rho}(k,t)/\partial k$ in von Neumann equation
Eq. (6), we use symmetric derivative

$\displaystyle\partialderivative{{\bf\it\rho}(k,t)}{k}=\dfrac{{\bf\it\rho}(k+dk,t)-{\bf\it\rho}(k-dk,t)}{2dk},$
(17)

where $dk=\pi/N$, with $N=5000$.

Based on this scheme, we calculate the linear and nonlinear susceptibility to
the electric field from the obtained current response, which will be shown in
the next section. In the following calculations, we use the parameters
$t_{h}=1,\lambda=0.8,J=0.6,K_{z}=0.2,h_{y}=0.2,\alpha_{G}=0.01,\gamma=0.01$,
otherwise explicitly mentioned. In FIG. 1 (b), we show the band structure of
the system in equilibrium. Due to the inversion breaking and time reversal
breaking in the system, the band structure is asymmetric about $k$, namely
$\epsilon(k)\neq\epsilon(-k)$. We set the chemical potential $\mu$ to be $0$
within the bandgap.

##  III Results

This section presents the result of the light-induced dynamics of the system.
Starting from the symmetry analysis, we clarify the optically active
collective excitations in III.1. Secondly, in III.2, we investigate the effect
of spin dynamics on optical response in the linear response regime. In III.3
and III.4, we study the effect of collective mode on photocurrent
conductivity. Here, we decompose the photocurrent into different processes,
where the linear susceptibility of spins to electric field obtained in III.2
plays a crucial role in photocurrent generation. The last subsection III.5 is
dedicated to demonstrating the tunability of photocurrent generation by
changing the canted angle of the local spin moments.

###  III.1 Symmetry analysis

Figure 2: Collecive mode of canted antiferromagnetic moment. The red (blue) arrow indicates the spin moment in the $A$ ($B$) sublattice, respectively. The black arrow denotes the summation of the sublattice spin moment. Table 1: Sign under symmetry operation of magnetic point group $\mathcal{G}$. $\mathrm{M}^{a}$ means ferroic configurations of spin moment along the $a$ direction, while $\mathrm{L}^{a}$ means staggered spin moments along $a$ direction. $+$ means the observable does not change its sign, while $-$ means the observable flip its sign under the operation. | $1$ | $\theta 2_{x}$ | $\bar{2}_{y}$ | $\theta\bar{2}_{z}$
---|---|---|---|---
$E^{z}$ | $+$ | $-$ | $+$ | $-$
$J^{z}$ | $+$ | $+$ | $+$ | $+$
$\mathrm{M}^{x}$ | $+$ | $+$ | $-$ | $-$
$\mathrm{L}^{x}$ | $+$ | $-$ | $+$ | $-$
$\mathrm{M}^{y}$ | $+$ | $+$ | $+$ | $+$
$\mathrm{L}^{y}$ | $+$ | $+$ | $-$ | $-$
$\mathrm{M}^{z}$ | $+$ | $+$ | $-$ | $-$
$\mathrm{L}^{z}$ | $+$ | $+$ | $+$ | $+$

In this subsection, we analyze the collective spin dynamics that is linearly
coupled to the external electric field based on the symmetry analysis. Since
photocurrent response is largely restricted by the symmetry of the system and
the external field, symmetry analysis of collective dynamics lays the
foundation of our study.

The system has the following symmetries and belongs to the magnetic point
group $\mathcal{G}=2^{\prime}mm^{\prime}$, explicitly given by

$\displaystyle\mathcal{G}=\quantity{1,\theta
2_{x},\bar{2}_{y},\theta\bar{2}_{z}}.$ (18)

Here, $1,\theta,2_{a}$, and $\bar{2}_{a}$ are the identity operator, time
reversal operator, twofold rotation around the $a$ axis, rotatory inversion
about the $a$ axis, respectively.

Let us identify the collective modes in the canted antiferromagnetic moments.
Since our method is based on the momentum-space formulation, we only consider
$k=0$ magnon, or antiferromagnetic resonance. As we show in FIG. 2, two
collective modes exist in our system, named $\alpha$-mode and $\beta$-mode
[39]. In the $\alpha$-mode, the total spin moment oscillates along the
$y$-direction, whereas the $\beta$-mode denotes the precession of net spin
moment around the $y$-axis.

To determine which mode is optically active, we analyze the sign of
observables under the symmetry operation of $\mathcal{G}$. For example, the
$\bar{2}_{y}$ operation can be explicitly expressed as

$\displaystyle\bar{2}_{y}=\quantity(-i\sigma_{y})\otimes\tau_{x},$ (19)

where $\sigma$ and $\tau$ are Pauli matrices representing the spin and
sublattice degrees of freedom. Under this operation, staggered spin moment
along the $x$ direction $\sigma_{x}\otimes\tau_{z}$ can be transformed as

$\displaystyle\bar{2}_{y}\quantity(\sigma_{x}\otimes\tau_{z})\bar{2}_{y}^{-1}=\sigma_{x}\otimes\tau_{z}.$
(20)

Therefore, $\sigma_{x}\otimes\tau_{z}$ does not change its sign under the
$\bar{2}_{y}$ operation. In TABLE 1, we summarize how the physical quantities
change their sign under the operation in $\mathcal{G}$. Here $E^{z},J^{z}$ are
electric field and electric current along the $z$ direction, respectively.
$\mathrm{M}^{a}$ means ferroic configurations of spin moment along the $a$
direction, while $\mathrm{L}^{a}$ means staggered spin moments along the $a$
direction. Owing to the incompatibility of the external electric field $E^{z}$
and $\mathrm{M}^{x},\mathrm{L}^{y},\mathrm{M}^{z}$, under $\bar{2}_{y}$
operation, the $\beta$-mode is not linearly coupled to the electric field.
Although the symmetry of $\mathrm{M}^{y},\mathrm{L}^{z}$ and $E^{z}$ are
different under $\theta 2_{x},\theta\bar{2}_{z}$ operation, these can be
linearly excited by electric current $J^{z}$. This is because, in
nonequilibrium phenomena, time-reversal symmetry is effectively broken by the
dissipation process. This can be viewed as an extension of the frequency
response of the magnetoelectric effect and Edelstein effect [40, 41].
Consequently, only components related to $\mathrm{L}^{x}$, $\mathrm{M}^{y}$,
and $\mathrm{L}^{z}$ can be linearly excited by light, indicating that the
$\alpha$-mode can be linearly coupled to the electric field. Since the
symmetry argument does not depend on frequency profile, it is valid for both
the resonant and off-resonant conditions.

There is no symmetry operation inverting the electric current $J^{z}$ since
the model shows the polar symmetry along the $z$ axis similar to that of the
typical photodiodes. Notably, the electric current changes its sign under the
time reversal operation, coinciding with the symmetry of the toroidal moment.
At the same time, the electric polarization is forbidden to be in the chain
direction due to the anti-unitary symmetry. This unique combination of
magnetic and polar symmetries leads to the absence of shift current and the
presence of injection current within the IPA, as proved in the Appendix B.

###  III.2 Linear Response Functions

Figure 3: (a) Linear electromagnetic susceptibility to the external electric
field. (b) Linear optical conductivity of the system. The red solid line and
the blue dashed line indicate the calculation with and without updating the
spin configurations, respectively. The black dashed line indicates the band
gap frequency.

Before we move on to the nonlinear optical response, we examine the optical
response in the linear response regime, which provides valuable insights into
the mutual interaction between charge and spin degrees of freedom. Here, we
investigate the two types of linear response functions. The first one is the
electromagnetic susceptibility, which has direct information on the spin
collective mode induced by light. The second one is the linear optical
conductivity of the system, by which we can understand how the collective
excitation of the local spin affects the optical response of the electronic
systems. Electromagnetic susceptibility is defined as

$\displaystyle\Delta\mathrm{M}^{a}(\omega)$
$\displaystyle=\chi_{\text{M}^{a}E}(\omega)E^{z}(\omega),$ (21)
$\displaystyle\Delta\mathrm{L}^{a}(\omega)$
$\displaystyle=\chi_{\text{L}^{a}E}(\omega)E^{z}(\omega).$ (22)

Here, $\Delta\mathrm{M}^{a}(\omega),\Delta\mathrm{L}^{a}(\omega)$ are the
Fourier components of the deviation of the spin configurations from the
equilibrium states, which are defined as

$\displaystyle\Delta\mathrm{M}^{a}(t)$
$\displaystyle=\frac{1}{2}\quantity(\Delta\mathrm{S}_{A}^{a}(t)+\Delta\mathrm{S}_{B}^{a}(t)),$
(23) $\displaystyle\Delta\mathrm{L}^{a}(t)$
$\displaystyle=\frac{1}{2}\quantity(\Delta\mathrm{S}_{A}^{a}(t)-\Delta\mathrm{S}_{B}^{a}(t)).$
(24)

We define the spin dynamics and electromagnetic susceptibilities related to
the $\alpha$-mode component as

$\displaystyle\Delta\mathbf{S}(\omega)$
$\displaystyle=\quantity(\Delta\mathrm{L}^{x}(\omega),\Delta\mathrm{M}^{y}(\omega),\Delta\mathrm{L}^{z}(\omega))^{T},$
(25) $\displaystyle{\bf\it\chi}_{\mathbf{S}E}(\omega)$
$\displaystyle=\quantity(\chi_{\mathrm{L}^{x}E}(\omega),\chi_{\mathrm{M}^{y}E}(\omega),\chi_{\mathrm{L}^{z}E}(\omega))^{T}.$
(26)

By using these quantities, spin dynamics linearly coupled to the electric
field are characterized by

$\displaystyle\Delta\mathbf{S}(\omega)={\bf\it\chi}_{\mathbf{S}E}(\omega)E^{z}(\omega).$
(27)

Additionally, linear conductivity is defined as

$\displaystyle J^{z}(\omega)=\sigma(\omega)E^{z}(\omega).$ (28)

Here, $J^{z}(\omega)$ is the Fourier transform of the $J^{z}(t)$.

To calculate these response functions, we applied a Gaussian pulse electric
field

$\displaystyle
E^{z}(t)=\dfrac{E_{0}}{\sqrt{2\pi\sigma^{2}}}\exp(-\dfrac{t^{2}}{2\sigma^{2}})$
(29)

at $t=0$. In the limit of $\sigma\to 0$, the gaussian pulse Eq. (29) becomes
$E^{z}(t)=E_{0}\delta(t)$, containing all the frequencies. The linear response
functions can be calculated by the Fourier transform of the resulting response
in the time domain as follows,

$\displaystyle\sigma(\omega)=\dfrac{1}{E_{0}}e^{\sigma^{2}\omega^{2}}\int_{0}^{\infty}e^{i\omega
t}J^{z}(t)dt.$ (30)

We plot the frequncy dependence of the linear response functions
${\bf\it\chi}_{\mathbf{S}E}(\omega)$, and $\sigma(\omega)$ in FIG. 3. In FIG.
3 (a-1)-(a-3), we show the imaginary part of the electromagnetic
susceptibility $\imaginary\chi_{\mathbf{S}E}(\omega)$. Electromagnetic
susceptibility exhibits the resonance structure below the optical gap of the
electronic system. This resonant structure corresponds to the optically active
$\alpha$-mode. Additionally, there are off-resonant contributions in the above
bandgap regime, coming from electronic excitation. Moreover, in the higher
frequency around $\omega\sim 1.5$, ${\bf\it\chi}_{\mathbf{S}E}$ shows small
amplitude, indicating that the local spins are unaffected by the electric
field. This is because the local spin moments with small resonance frequency
cannot follow the fast oscillations by the external field. In the linear
response calculations, we confirmed that only the $\alpha$-mode is linearly
coupled to the electric field, in agreement with the symmetry analysis.

In FIG. 3(b), we show the linear optical conductivity
$\operatorname{Re}\sigma^{zz}(\omega)$. The blue dashed line, denoted as
without LLG, shows the linear conductivity spectra without updating the spin
configuration, which corresponds to the IPA. The red solid line, denoted as
with LLG, indicates the conductivity spectrum with the effect of spin
dynamics. There are three characteristics of the effect of spin dynamics.
Firstly, $\operatorname{Re}\sigma(\omega)$ with LLG simulation shows the
resonant peak in the in-gap regime, which is absent in the IPA. The frequency
of this peak corresponds to that of the $\alpha$-mode we saw in the
electromagnetic susceptibility. Secondly, the presence of spin dynamics
introduces a suppression of optical conductivity within the regime above the
band gap. This suppression arises due to off-resonant contribution from
collective spin dynamics. It is noteworthy that due to the sum rule
$\int\operatorname{Re}\sigma(\omega)d\omega=\mathrm{const.}$, the emergence of
the in-gap peak results in the suppression of conductivity in the above-gap
regime. Note that this does not apply to the photocurrent conductivity
discussed below. Lastly, deviations from the IPA calculation become negligible
in the higher frequency region around $\omega\sim 1.5$, as the
$\operatorname{Im}{\bf\it\chi}_{\mathrm{S}E}$ lose its amplitude.

###  III.3 Photocurrent Spectra

Figure 4: (a) Photocurrent spectra. The blue dashed and red solid lines
indicate the IPA calculation and calculation incorporated the spin dynamics,
respectively. (b,c) Dependence of photocurrent conductivity at the above
bandgap frequency ($\omega=0.8$) on (b) electric-field strength $E_{0}$ and
(c) relaxation time $\tau$. The blue and red lines indicated the IPA
calculation and calculation incorporated the spin dynamics, respectively. The
black solid line in (b) indicates the line which is proportional to
$E_{0}^{2}$.

Now, we show the result of photocurrent spectra, revealing the effect of
collective spin dynamics on the BPVE. The second-order BPVE can be described
as

$\displaystyle\begin{split}J^{z}(\omega=0,\omega_{p})=\sigma^{z;zz}(0;\omega_{p},-\omega_{p})E^{z}(\omega_{p})E^{z}(-\omega_{p})\\\
+\sigma^{z;zz}(0;-\omega_{p},\omega_{p})E^{z}(-\omega_{p})E^{z}(\omega_{p}).\end{split}$
(31)

Here $J^{z}(\omega=0,\omega_{p})$ is the DC component of output current along
chain direction $z$ induced by the electric field with frequency $\omega_{p}$.
To calculate the photocurrent spectra, we apply the continuous electric field
$E^{z}(t)=E_{0}\sin\omega_{p}t$. Here we define the photocurrent conductivity
$\sigma(0;\omega_{p})$ as

$\displaystyle\sigma(\omega=0;\omega_{p})$
$\displaystyle=\dfrac{1}{2}\quantity(\sigma^{z;zz}(0;\omega_{p},-\omega_{p})+\sigma^{z;zz}(0;-\omega_{p},\omega_{p})).$
(32)

In this definition, the photocurrent conductivity is calculated by the
following formula [22]

$\displaystyle\sigma(\omega=0,\omega_{p})=\dfrac{2}{E_{0}^{2}NT_{p}}\int_{t_{\text{sat}}}^{t_{\text{sat}}+NT_{p}}J^{z}(t)dt.$
(33)

Here $T_{p}=2\pi/\omega_{p}$ is the period of the external light. We set
$t_{\text{sat}}$ large enough so that the system reaches the time-periodic
steady state and use $N>10$ to get the average.

In FIG. 4(a), we show the photocurrent spectra, with and without the
contributions from spin dynamics. The blue dashed line indicates the IPA
calculation, and the red solid line indicates the photocurrent spectra with
the effect of spin dynamics.

Firstly, we observed the sharp resonant contribution from collective spin
dynamics in the in-gap regime. This corresponds to the $\alpha$-mode resonance
we observed in the linear response functions in FIG. 3. Secondly, we confirmed
the substantial enhancement above the band gap regime, which comes from off-
resonant spin dynamics. Moreover, we observed the dip structure at the band
edge frequency, which comes from the interference effect of the external field
and internal spin dynamics, as shown below. In FIG. 4 (b), (c), we show the
field amplitude and the relaxation time $\tau=\gamma^{-1}$ dependence of the
photocurrent at the above-band gap frequency $\omega_{p}=0.8$. Although the
photocurrent response obtained in our method includes a higher-order nonlinear
optical response to the electric field in principle, we confirmed that the
second-order response is dominant. Additionally, the relaxation time
dependence of the photocurrent spectra scaling linearly to the relaxation time
indicates that the injection current contribution is dominant, which agrees
with symmetry analysis. It is noteworthy that our time-dependent calculation
captures photocurrent not only from resonant contributions but also from off-
resonant contributions of collective spin dynamics, which were not discussed
in the earlier studies [21, 42] of BPVE from magnetic excitation.

###  III.4 Decomposition of Photocurrent

Figure 5: (a) Schematic picture of three different processes of the
photocurrent in the presence of collective spin dynamics. The wavy and dashed
lines indicate the electric field and interaction $J$, respectively. $\otimes$
represents the output photocurrent. The solid triangles describe the
photocurrent susceptibility evaluated with the IPA. $\Delta\mathbf{S}$ is the
light induced spin dynamics. (b) Photocurrent spectra originating from the
three different processes. The inset shows the magnified view of
$\sigma_{\text{col-E}}$ component around the collective mode frequency
$\omega_{p}\sim 0.35$. Figure 6: Relaxation time dependence of each
photocurrent conductivity component. The insets in (b) and (c) show the
magnified view of photocurrent spectra around the collective mode frequency
$\omega_{p}\sim 0.35$.

To obtain further insight into the role of spin dynamics on the PVE, we
decompose the photocurrent contributions into three distinct processes. Our
calculation is essentially equivalent to the random phase approximation [43,
44], wherein the light-induced internal spin dynamics can be considered as an
additional external field applied to the conduction electrons (see Eq. (3) and
Eq. (9)). Consequently, as illustrated in FIG. 5(a), three distinct processes
contribute to the photocurrent response, described as

$\displaystyle J^{z}(\omega=0,\omega_{p})$
$\displaystyle=J^{z}_{0}(\omega_{p})+J^{z}_{\text{col-E}}(\omega_{p})+J^{z}_{\text{col-
col}}(\omega_{p}),$ (34)

where

$\displaystyle J_{0}^{z}(\omega_{p})$
$\displaystyle=\sigma_{EE}^{z;zz}(0;-\omega_{p},\omega_{p})E^{z}(-\omega_{p})E^{z}(\omega_{p})$
$\displaystyle+\sigma_{EE}^{z;zz}(0;\omega_{p},-\omega_{p})E^{z}(\omega_{p})E^{z}(-\omega_{p}),$
(35) $\displaystyle J_{\text{col-E}}^{z}(\omega_{p})$
$\displaystyle=\sum_{\lambda}\sigma_{E\mathrm{S}}^{z;z\lambda}(0;-\omega_{p},\omega_{p})E^{z}(-\omega_{p})\Delta\mathrm{S}^{\lambda}(\omega_{p})$
$\displaystyle+\sum_{\lambda}\sigma_{E\mathrm{S}}^{z;z\lambda}(0;\omega_{p},-\omega_{p})E^{z}(\omega_{p})\Delta\mathrm{S}^{\lambda}(-\omega_{p}),$
(36) $\displaystyle J_{\text{col-col}}^{z}(\omega_{p})$
$\displaystyle=\sum_{\nu\lambda}\sigma_{\mathrm{S}\mathrm{S}}^{z;\nu\lambda}(0;-\omega_{p},\omega_{p})\Delta\mathrm{S}^{\nu}(-\omega_{p})\Delta\mathrm{S}^{\lambda}(\omega_{p})$
$\displaystyle+\sum_{\nu\lambda}\sigma_{\mathrm{S}\mathrm{S}}^{z;\nu\lambda}(0;\omega_{p},-\omega_{p})\Delta\mathrm{S}^{\nu}(\omega_{p})\Delta\mathrm{S}^{\lambda}(-\omega_{p}).$
(37)

Here we use the spin dynamics $\Delta\mathbf{S}(\omega)$ defined in Eq. (25),
and indices $\nu,\lambda$ represent the components of
$\Delta\mathbf{S}(\omega)$. Firstly, $J_{0}$ depicts the photocurrent in the
absense of collective spin dynamics, described by the photocurrent
conductivity $\sigma_{EE}$. This contribution corresponds to the IPA. The
second contribution is $J_{\text{col-E}}$, where the external electric field
and light-induced spin dynamics synergistically generate the photocurrent, and
is characterized by $\sigma_{E\mathrm{S}}$. This can be understood as
interference between the external electric field and collective spin dynamics.
The third contribution $J_{\text{col-col}}$ involves solely the light-induced
internal spin fields, to generate the photocurrent, and is characterized by
$\sigma_{\mathrm{SS}}$. Given that we have confirmed the linear coupling of
spin dynamics to the electric field, all these components yield a second-order
response to the electric field. We emphasize that
$\sigma_{EE},\sigma_{E\mathrm{S}},\sigma_{\mathrm{SS}}$ are calculated within
the IPA framework. The detailed expressions of them are given in Appendix A.
In linear response, the total output current with multiple fields is simply
the superposition of the individual outputs for each field. However, in
nonlinear response, the total output is influenced by interference effects
that result from the interaction of different fields. Based on the
photocurrent classification in Eq. (35), Eq. (36), Eq. (37), we can define
photocurrent conductivity in the presence of collective spin dynamics as

$\displaystyle J_{0}^{z}(\omega_{p})$
$\displaystyle=2\sigma_{0}(0;\omega_{p})E^{z}(\omega_{p})E^{z}(-\omega_{p}),$
(38) $\displaystyle J_{\text{col-E}}^{z}(\omega_{p})$
$\displaystyle=2\sigma_{\text{col-E}}(0;\omega_{p})E^{z}(\omega_{p})E^{z}(-\omega_{p}),$
(39) $\displaystyle J_{\text{col-col}}^{z}(\omega_{p})$
$\displaystyle=2\sigma_{\text{col-
col}}(0;\omega_{p})E^{z}(\omega_{p})E^{z}(-\omega_{p}),$ (40)

where

$\displaystyle\sigma_{0}(0;\omega_{p})$
$\displaystyle=\dfrac{1}{2}\quantity(\sigma_{EE}^{z;zz}(0;-\omega_{p},\omega_{p})+\sigma_{EE}^{z;zz}(0;\omega_{p},-\omega_{p})),$
(41)
$\displaystyle\begin{split}\sigma_{\text{col-E}}(0;\omega_{p})&=\dfrac{1}{2}\sum_{\lambda}\sigma_{E\mathrm{S}}^{z;z\lambda}(0;-\omega_{p},\omega_{p})\chi_{\mathrm{S}^{\lambda}E}(\omega_{p})\\\
&+\dfrac{1}{2}\sum_{\lambda}\sigma_{E\mathrm{S}}^{z;z\lambda}(0;\omega_{p},\omega_{p})\chi_{\mathrm{S}^{\lambda}E}(-\omega_{p}),\end{split}$
(42) $\displaystyle\begin{split}\sigma_{\text{col-
col}}(0;\omega_{p})&=\dfrac{1}{2}\sum_{\nu\lambda}\sigma_{\mathrm{SS}}^{z;\nu\lambda}(0;-\omega_{p},\omega_{p})\chi_{\mathrm{S}^{\nu}E}(-\omega_{p})\chi_{\mathrm{S}^{\lambda}E}(\omega_{p})\\\
&+\dfrac{1}{2}\sum_{\nu\lambda}\sigma_{\mathrm{SS}}^{z;\nu\lambda}(0;\omega_{p},-\omega_{p})\chi_{\mathrm{S}^{\nu}E}(\omega_{p})\chi_{\mathrm{S}^{\lambda}E}(-\omega_{p}).\end{split}$
(43)

Here we use electromagnetic susceptibilities defined in Eq. (26). Our time-
dependent calculations naturally give the decomposition as follows. Firstly,
$\sigma_{0}(\omega_{p})$ corresponds to the calculation without updating spin
configurations. Secondly, $\sigma_{\text{col-col}}(\omega_{p})$ can be
calculated by switching off the external electric field $E^{z}(t)$ and
updating the Hamiltonian with $\mathbf{S}(t)$ obtained in the calculations
with the LLG equation. The $\sigma_{\text{col-E}}(\omega_{p})$ spectrum can be
calculated by subtracting $\sigma_{0}(\omega_{p})$ and $\sigma_{\text{col-
col}}(\omega_{p})$ from $\sigma(\omega_{p})$.

In FIG. 5(b), we show the photocurrent spectra corresponding to each
component. $\sigma_{0}(\omega_{p})$ is the same as the blue dashed line in
FIG. 4(a). $\sigma_{\text{col-E}}(\omega_{p})$ component undergoes sign change
near the bandgap frequency. This sign change may come from the interference
effect between the external electric field and field-induced spin dynamics,
which might give a new way to control the photocurrent direction. In addition,
we confirmed that the strong resonance peak mainly comes from
$\sigma_{\text{col-col}}(\omega_{p})$ component. Since the resonant structure
of $\sigma(\omega_{p})$ in the in-gap regime originates from the resonant
structure in $\chi_{\mathbf{S}E}(\omega_{p})$, $\sigma_{\text{col-
col}}\propto\quantity(\chi_{\mathbf{S}E})^{2}$ gives a shaper resonance peak
than $\sigma_{\text{col-E}}\propto\chi_{\mathbf{S}E}$ does.

Table 2: Photocurrent generation induced by the electric field and spin dynamics. For instance, Electric field $E$ and uniform spin dynamics along $y$ direction $\mathrm{M}^{y}$ induces $\sigma_{\text{shift}}$ and $\sigma_{\text{Enj}}$. $\sigma_{\text{shift}}$ is induced by in-phase component of $E$ and $\mathrm{M}^{y}$, while the $\sigma_{\text{Enj}}$ is induced by out-of-phase component of $E$ and $\mathrm{M}^{y}$. The detailed expression for $\sigma_{\text{shift}},\sigma_{\text{gyro}},\sigma_{\text{Enj}},\sigma_{\text{Mnj}}$ are given in the Appendix A | $E$ | $\mathrm{L}^{x}$ | $\mathrm{M}^{y}$ | $\mathrm{L}^{z}$
---|---|---|---|---
$E$ | $\sigma_{\text{Mnj}}$ |  $\sigma_{\text{Mnj}}$, $\sigma_{\text{gyro}}$ |  $\sigma_{\text{shift}}$, $\sigma_{\text{Enj}}$ |  $\sigma_{\text{shift}}$, $\sigma_{\text{Enj}}$
$\mathrm{L}^{x}$ | $\sigma_{\text{Mnj}},\sigma_{\text{gyro}}$ | $\sigma_{\text{Mnj}}$ | $\sigma_{\text{Enj}},\sigma_{\text{shift}}$ | $\sigma_{\text{Enj}},\sigma_{\text{shift}}$
$\mathrm{M}^{y}$ | $\sigma_{\text{shift}},\sigma_{\text{Enj}}$ | $\sigma_{\text{shift}},\sigma_{\text{Enj}}$ | $\sigma_{\text{Mnj}}$ | $\sigma_{\text{Mnj}}$
$\mathrm{L}^{z}$ | $\sigma_{\text{shift}},\sigma_{\text{Enj}}$ | $\sigma_{\text{shift}},\sigma_{\text{Enj}}$ | $\sigma_{\text{Mnj}}$ | $\sigma_{\text{Mnj}}$

It is worth noting that symmetry relations and the phase degrees of freedom of
driving fields constrain the types of photocurrent generated. As a general
case, let us focus on noninteracting electrons and consider the photocurrent
induced by two external fields $X(\omega_{p})$ and $Y(\omega_{p})$ within the
IPA framework. This process can be written by using the photocurrent
conductivity $\sigma_{XY}(0;\omega_{p},-\omega_{p})$ as

$\displaystyle\begin{split}J_{XY}(\omega_{p})&=\sigma_{XY}(0;-\omega_{p},\omega_{p})X(-\omega_{p})Y(\omega_{p})\\\
&+\sigma_{XY}(0;\omega_{p},-\omega_{p})X(\omega_{p})Y(-\omega_{p}).\end{split}$
(44)

As we show in Appendix A, $\sigma_{XY}$ can be classified into the following
four components,

$\displaystyle\sigma_{XY}=\sigma_{XY,\text{shift}}+\sigma_{XY,\text{gyro}}+\sigma_{XY,\text{Mnj}}+\sigma_{XY,\text{Enj}}.$
(45)

Here, $\sigma_{XY,\text{shift}}$ and $\sigma_{XY,\text{gyro}}$ correspond to
shift current and gyration current, while $\sigma_{XY,\text{Mnj}}$ and
$\sigma_{XY,\text{Enj}}$ correspond to magnetic injection current and electric
injection current, as discussed in the context of normal photocurrent [8, 7].
These contributions are different in terms of relaxation time dependence.
$\sigma_{XY,\text{shift}}$ and $\sigma_{XY,\text{gyro}}$ are independent of
relaxation time, whereas $\sigma_{XY,\text{Mnj}}$ and $\sigma_{XY,\text{Enj}}$
are proportional to relaxation time. Additionally,
$\sigma_{XY,\text{shift}},\sigma_{XY,\text{Mnj}}$
($\sigma_{XY,\text{gyro}},\sigma_{XY,\text{Enj}}$) is the real (pure-
imaginary) quantities, which describes the counterpart of the photocurrent
induced by linearly (circularly) polarized light. Therefore,
$\sigma_{XY,\text{shift}},\sigma_{XY,\text{Mnj}}$ are finite only when two
fields $X(\omega_{p}),Y(\omega_{p})$ are in-phase to each other, while
$\sigma_{XY,\text{gyro}},\sigma_{XY,\text{Enj}}$ are finite only when
$X(\omega_{p}),Y(\omega_{p})$ are out-of-phase to each other. Besides this
restriction from phase degrees of freedom, the symmetry of the two fields also
restricts the type of photocurrent.

In our case, the field $X(\omega_{p}),Y(\omega_{p})$ corresponds to the
electric field $E(\omega_{p})$ and/or the spin dynamics
$\mathrm{S}^{\nu}(\omega_{p})$. For instance, we confirmed that the spin
dynamics $\mathrm{L}^{x}(\omega_{p})$ and $\mathrm{M}^{y}(\omega)$ has both
in-phase and out-of-phase components, and their sign under $\theta 2_{x}$ are
different. This fact restricts the photocurrent generation from these fields
to shift current and electric injection current similar to the photocurrent
response to time-reversal-symmetric systems [8]. A similar argument can be
applied to the case of the photocurrent originating from the interference of
electric field and collective spin dynamics described by
$\sigma_{E\mathrm{S}}$. We summarize the photocurrent classification in TABLE
2.

From this analysis, we reveal that unlike $\sigma_{0}$, $\sigma_{\mathrm{SS}}$
and $\sigma_{E\mathrm{S}}$ can include the component of shift current-like
contribution. Moreover, although we consider the photocurrent induced by the
linearly polarized light, $\sigma_{\mathrm{SS}}$ and $\sigma_{E\mathrm{S}}$
has $\sigma_{\text{Enj}},\sigma_{\text{gyro}}$ component which is the
counterpart of circularly polarized induced photocurrent. Practically, we can
observe the in-gap resonant structure in $\sigma_{\text{col-E}}(\omega_{p})$
is mainly coming from the injection current-like contribution, since
$\operatorname{Re}\chi_{\mathrm{S}^{\nu}E}$
($\operatorname{Im}\chi_{\mathrm{S}^{\nu}E}$) represents the in-phse (out-of-
phase) component of spin dynamics to the electric field (Appendix B for a more
detailed discussion). Here, we reveal that the fictitious fields arising from
the spin collective dynamics result in various types of photocurrent, which is
not allowed in the IPA.

In FIG. 6, we show the relaxation time $\tau$ dependence of photocurrent
spectra. Photocurrent spectra increase linearly with relaxation time,
indicating the injection current-like contribution is dominant. Although we
have revealed the presence of shift current-like contribution in
$\sigma_{\text{col-E}}$ and $\sigma_{\text{col-col}}$, the primary
contributions to $\sigma_{\text{col-E}}$ and $\sigma_{\text{col-col}}$ are
injection current-like one, namely $\sigma_{\text{Enj}},\sigma_{\text{Mnj}}$.

###  III.5 Tuning of photocurrent by external magnetic field

Figure 7: (a) Magnetic field dependence of linear electromagnetic
susceptibility $\operatorname{Im}\chi_{\mathrm{L}^{x}E}(\omega)$ to the
external electric field. (b) Photocurrent spectra
$\sigma(\omega=0;\omega_{p})-\sigma_{0}(\omega=0;\omega_{p})$ with modulating
the canted angle of antiferromagnetic moments by changing the magnetic field
along the $y$ direction

Here, we demonstrate the tunability of photocurrent by applying the external
magnetic field perpendicular to the chain direction, changing the canted angle
of antiferromagnetic moments. In FIG. 7, we show
${\bf\it\chi}_{\mathbf{S}E}(\omega_{p})$ and the photocurrent spectra
originating from spin dynamics by changing the amplitude of the external
magnetic field.

In FIG. 7 (a), we show the linear electromagnetic susceptibility of
$\operatorname{Im}\chi_{\mathrm{L}^{x}E}$. As we increase the magnetic field
along the $y$ axis, the amplitude of $\operatorname{Im}\chi_{\mathrm{L}^{x}E}$
gradually increases, showing the strong modulation of feedback between charge
and spin degrees of freedom by the external field. This modulation arises from
the fact that the alteration of the canted angle induces changes in the
optical gap of the system and the energy scale of the collective spin mode.
Consequently, this modulation allows us to tune the strength of the feedback
from spin dynamics.

Consistent with the modulation of the electromagnetic susceptibility, the
photocurrent conductivity spectra also experience large modulation due to the
external field, as shown in FIG. 7 (b). As explained in the previous
subsection, ${\bf\it\chi}_{\mathbf{S}E}(\omega_{p})$ plays a critical role in
photocurrent generation in the presence of collective spin dynamics.
Therefore, larger ${\bf\it\chi}_{\mathbf{S}E}(\omega_{p})$ results in a more
significant modulation in photocurrent generation.

As these results demonstrate the importance of energy scale matching between
collective modes and electronic excitations, we can expect substantial
modulation of photocurrent from magnetic excitation in the narrow gap system
with the higher collective mode frequencies.

##  IV Summary and discussion

We investigated the effect of antiferromagnetic excitations on linear and
nonlinear optical responses by a real-time calculation method. Firstly, we
observed the optical response arising from an electrically active
antiferromagnetic resonance mode, which aligned with our symmetry analysis. We
further found that collective spin dynamics, originating from both resonant
and off-resonant contributions, can significantly enhance the photocurrent
response. In addition, we delved into a comprehensive analysis by classifying
it into different processes. The component includes the unique photocurrent
arising from the interference between the electric field and the internal spin
field, which is inherent to the nonlinear response. Moreover, based on the
symmetry analysis, we revealed that the emergence of the spin dynamics allows
the various types of photocurrent, which are totally absent in the IPA.
Additionally, we showcased the tunability of the photocurrent in our system by
varying the canted angles and emphasized the significance of matching in the
energy scales of collective modes and electronic excitations.

Now, we estimate the strength of modulation by spin dynamics. We asuume the
energy scale of hopping $t_{h}=$0.1\text{\,}\mathrm{e}\mathrm{V}$$ and lattice
constant $a=$0.5\text{\,}\mathrm{n}\mathrm{m}$$. Additionally, to compare the
obtained photocurrent conductivity in a one-dimensional system to that of a
three-dimensional system, we assume the three-dimensional stack of a one-
dimensional chain with a lattice spacing of $1.0\text{\,}\mathrm{n}\mathrm{m}$
[45]. In this condition, the modulation of the photocurrent conductivity by
spin dynamics $\sigma_{\text{col-E}}+\sigma_{\text{col-col}}$ with relaxation
time $\tau\sim$60\text{\,}\mathrm{f}\mathrm{s}$$ reaches the order of
$\sim$10\text{\,}\mathrm{m}\mathrm{A}\mathrm{/}\mathrm{V}^{2}$$, which can be
larger than that of Weyl semimetals [13], demonstrating the significance of
the effect of spin dynamics on photocurrent response.

While our calculations are limited to the simple two-sublattice model, we can
extend these to more complex systems, like skyrmions [46, 47, 48], and
hedgehogs [49]. The essential ingredient is the optically active collective
spin dynamics, whose energy scale is close to electronic degrees of freedom.
Therefore, searching for a real material should be in the realm of inversion
broken narrow gap semiconductors with optically active magnons, such as
$\ce{MnBi2Te4}$ [50, 51].

## acknowledgement

This work is supported by Grant-in-Aid for Scientific Research from JSPS,
KAKENHI Grant No. JP23K13058 (H.W.), No. JP21H05017 (Y.M.), No. 21H04990
(R.A.), No. JP21H04437 (R.A.), No. JP20K14412 (Y.M.), No. 19H05825 (R.A.),
JST-PRESTO No. JPMJPR20L7 (T.N.), and JST-CREST No. JPMJCR1901 (Y.M.), No.
JPMJCR18T3 (R.A.).

##  Appendix A Nonlinear optical response functions based on density matrix
formalism

This section is dedicated to showing the formalism of perturbative expansion
of the von Neumann equation Eq. (6) in the presence of both the electric field
and light-induced spin dynamics. The Hamiltonian for electrons under an
external electric field and effective magnetic fields from the spin dynamics
is expressed as

$\displaystyle\mathcal{H}(t)=\mathcal{H}_{0}+\Delta\mathcal{H}(t).$ (46)

$\Delta\mathcal{H}(t)$ is the perturbative Hamiltonian and consists of an
external electric field and light-induced spin dynamics defined as

$\displaystyle\Delta\mathcal{H}(t)$
$\displaystyle=\Delta\mathcal{H}_{E}(t)+\Delta\mathcal{H}_{\text{spin}}(t),$
(47) $\displaystyle\Delta\mathcal{H}_{E}$ $\displaystyle=-r^{\mu}E^{\mu},$
(48) $\displaystyle\Delta\mathcal{H}_{\text{spin}}(t)$
$\displaystyle=-JA^{\mu}\Delta\mathrm{S}^{\mu}(t),$ (49)
$\displaystyle\Delta\mathbf{S}(t)$
$\displaystyle=\quantity(\Delta\mathrm{L}^{x}(t),\Delta\mathrm{M}^{y}(t),\Delta\mathrm{L}^{z}(t))^{T},$
(50) $\displaystyle{\bf\it A}$
$\displaystyle=\quantity(\sigma_{x}\tau_{z},\sigma_{y}\tau_{0},\sigma_{z}\tau_{z})^{T}.$
(51)

Here, $\sigma^{\mu},\tau^{\mu}$ are the Pauli matrices that describe the spin
and sublattice degrees of freedom, respectively. $\mathcal{H}_{0}$ is the
nonperturbative Hamiltonian given by

$\displaystyle\mathcal{H}_{0}=\int\dfrac{d{\bf\it
k}}{(2\pi)^{d}}\sum_{a}\epsilon_{{\bf\it k}a}c_{{\bf\it
k}a}^{\dagger}c_{{\bf\it k}a}.$ (52)

Here $c_{{\bf\it k}a}^{\dagger}(c_{{\bf\it k}a})$ is the creation
(annihilation) operators of the Bloch state $\ket{\psi_{{\bf\it
k}a}}=\exp(i{\bf\it k}\cdot\mathbf{r})\ket{u_{{\bf\it k}a}}$. $\ket{u_{{\bf\it
k}a}}$ is the cell-periodic part of the Bloch state and satisfies the
following equation.

$\displaystyle H_{0}({\bf\it k})\ket{u_{{\bf\it k}a}}=\epsilon_{{\bf\it
k}a}\ket{u_{{\bf\it k}a}}.$ (53)

According to [52], in the infinite volume limit the position operator in Eq.
(48) is written in the Bloch representation as

$\displaystyle\quantity[{\bf\it
r_{k}}]_{ab}=i\nabla_{\mathbf{k}}\delta_{ab}+{\bf\it\xi}_{ab}.$ (54)

Here ${\bf\it\xi}_{ab}$ is the Berry connection defined as
${\bf\it\xi}_{ab}=i\matrixelement{u_{{\bf\it k}a}}{\nabla_{{\bf\it
k}}}{u_{{\bf\it k}b}}$. In the following calculation, $k$ dependence is
omitted, otherwise explicitly mentioned.

As we explained in the main text, the von Neumann equation of SPDM
$\rho_{ab}=\langle c_{{\bf\it k}b}^{\dagger}c_{{\bf\it k}a}\rangle$ is given
by,

$\displaystyle i\partialderivative{\rho_{ab}(t)}{t}$
$\displaystyle=\quantity[\mathcal{H},\rho]_{ab}=\quantity[\mathcal{H}_{0},\rho]_{ab}+\quantity[\Delta\mathcal{H}(t),\rho]_{ab}.$
(55)

We define the Fourier transformation as

$\displaystyle f(t)=\int\dfrac{d\omega}{2\pi}f(\omega)e^{-i(\omega+i\eta)t},$
(56)

and we get the frequency-space representation of the von Neumann equation,

$\displaystyle(\omega+i\eta-\epsilon_{ab})\rho_{ab}(\omega)=\int\dfrac{d\omega_{1}}{2\pi}\quantity[\Delta\mathcal{H}(\omega_{1}),\rho(\omega-\omega_{1})]_{ab}.$
(57)

Here, we introduced the infinitesimal parameter $\eta$ to describe the
adiabatic application of the external field. We can perturbatively expand this
equation by introducing the $\rho^{(n)}$, which is the $n$-th order with
respect to the external field,

$\displaystyle(\omega+i\eta-\epsilon_{ab})\rho_{ab}^{(n+1)}(\omega)=\int\dfrac{d\omega_{1}}{2\pi}\quantity[\Delta\mathcal{H}(\omega_{1}),\rho^{(n)}(\omega-\omega_{1})]_{ab}.$
(58)

We introduce matrix $d(\omega)$ as

$\displaystyle
d_{ab}(\omega)=\dfrac{1}{\omega+i\eta-\epsilon_{ab}}=\mathcal{P}\dfrac{1}{\omega-\epsilon_{ab}}-i\pi\delta(\omega-\epsilon_{ab}),$
(59)

and by using the Hadamard product $(A\odot B)_{ab}=A_{ab}B_{ab}$ [53], the von
Neumann equation Eq. (58) can be recast as

$\displaystyle\rho_{ab}^{(n+1)}(\omega)=\int\dfrac{d\omega_{1}}{2\pi}\quantity(d(\omega)\odot[\Delta\mathcal{H}(\omega_{1}),\rho^{(n)}(\omega-\omega_{1})])_{ab}.$
(60)

In the following calculation, we solve this equation with the boundary
condition of $\rho_{ab}^{(0)}(\omega)=2\pi\delta(\omega)f_{ab}\delta_{ab}$.
Since we focus on the second-order response to the external field, we solve
the equation

$\displaystyle\begin{split}\rho_{ab}^{(2)}(\omega)&=\int\dfrac{d\omega_{1}}{2\pi}\quantity(d(\omega)\odot[\Delta\mathcal{H}(\omega_{1}),\rho^{(1)}(\omega-\omega_{1})])_{ab}\\\
&=\int\dfrac{d\omega_{1}d\omega_{2}}{(2\pi)^{2}}\quantity(d(\omega)\odot[\Delta\mathcal{H}(\omega_{1}),d(\omega-\omega_{1})\odot[\Delta\mathcal{H}(\omega_{2}),\rho^{(0)}(\omega-\omega_{1}-\omega_{2})]])_{ab}.\end{split}$
(61)

Considering both the electric field and the spin dynamics,
$\rho_{ab}^{(2)}(\omega)$ can be divided into three components as,

$\displaystyle\rho_{ab}^{(2)}(\omega)=\rho_{EE,ab}(\omega)+\rho_{E\mathrm{S},ab}(\omega)+\rho_{\mathrm{S}\mathrm{S},ab}(\omega),$
(62)

where

$\displaystyle\rho_{EE,ab}(\omega)$
$\displaystyle=\dfrac{1}{2}\int\dfrac{d\omega_{1}d\omega_{2}}{(2\pi)^{2}}E^{\nu}(\omega_{1})E^{\lambda}(\omega_{2})\quantity(d(\omega)\odot[r^{\nu},d(\omega-\omega_{1})\odot[r^{\lambda},\rho^{(0)}(\omega-\omega_{1}-\omega_{2})]])_{ab}+[(\nu,\omega_{1})\leftrightarrow(\lambda,\omega_{2})],$
(63)
$\displaystyle\begin{split}\rho_{E\mathrm{S},ab}(\omega)&=\dfrac{J}{2}\int\dfrac{d\omega_{1}d\omega_{2}}{(2\pi)^{2}}E^{\nu}(\omega_{1})\mathrm{S}^{\lambda}(\omega_{2})\quantity(d(\omega)\odot[r^{\nu},d(\omega-\omega_{1})\odot[A^{\lambda},\rho^{(0)}(\omega-\omega_{1}-\omega_{2})]])_{ab}\\\
&+\dfrac{J}{2}\int\dfrac{d\omega_{1}d\omega_{2}}{(2\pi)^{2}}\mathrm{S}^{\lambda}(\omega_{1})E^{\nu}(\omega_{2})\quantity(d(\omega)\odot[A^{\lambda},d(\omega-\omega_{1})\odot[r^{\nu},\rho^{(0)}(\omega-\omega_{1}-\omega_{2})]])_{ab},\end{split}$
(64) $\displaystyle\rho_{\mathrm{S}\mathrm{S},ab}(\omega)$
$\displaystyle=\dfrac{J^{2}}{2}\int\dfrac{d\omega_{1}d\omega_{2}}{(2\pi)^{2}}\mathrm{S}^{\nu}(\omega_{1})\mathrm{S}^{\lambda}(\omega_{2})\quantity(d(\omega)\odot[A^{\nu},d(\omega-\omega_{1})\odot[A^{\lambda},\rho^{(0)}(\omega-\omega_{1}-\omega_{2})]])_{ab}+[(\nu,\omega_{1})\leftrightarrow(\lambda,\omega_{2})].$
(65)

Here, we symmetrize the expression of $\rho$ since the perturbative
calculation should be invariant under the arbitrary exchange of external
field. Based on this decomposition, we formulate the photocurrent formula in
the presence of spin dynamics.

###  A.1 Electric field induced photocurrent

Here we derive the electric field-induced photocurrent based on the
$\rho_{EE}$. It is convenient to decompose the position operator in Eq. (54)
into inter- (${\bf\it r}_{e}$) and intra-band (${\bf\it r}_{i}$) component as

$\displaystyle\begin{split}\quantity({\bf\it
r}_{i})_{ab}&=\delta_{ab}\quantity(i\nabla_{\mathbf{k}}+\xi_{aa}),\\\
\quantity({\bf\it r}_{e})_{ab}&=\quantity(1-\delta_{ab})\xi_{ab}.\end{split}$
(66)

Based on this decomposition, $\rho_{EE}(\omega)$ can be classified into the
following four component,

$\displaystyle\rho_{EE,ab}(\omega)=\rho_{EE,ab}^{(ii)}+\rho_{EE,ab}^{(ie)}+\rho_{EE,ab}^{(ei)}+\rho_{EE,ab}^{(ee)}.$
(67)

Each term is explicitly written as,

$\displaystyle\begin{split}\rho_{EE,ab}^{(ii)}&=\dfrac{1}{2}\int\dfrac{d\omega_{1}d\omega_{2}}{(2\pi)^{2}}E^{\nu}(\omega_{1})E^{\lambda}(\omega_{2})\quantity(d(\omega)\odot[r_{i}^{\nu},d(\omega-\omega_{1})\odot[r_{i}^{\lambda},\rho^{(0)}(\omega-\omega_{1}-\omega_{2})]])_{ab}+[(\nu,\omega_{1})\leftrightarrow(\lambda,\omega_{2})]\\\
&=-\int\dfrac{d\omega_{1}d\omega_{2}}{(2\pi)^{2}}E^{\nu}(\omega_{1})E^{\lambda}(\omega_{2})d_{ab}(\omega)d_{ab}(\omega-\omega_{1})\partial_{\nu}\partial_{\lambda}f(\epsilon_{{\bf\it
k}a})\delta_{ab}2\pi\delta(\omega-\omega_{1}-\omega_{2})+[(\nu,\omega_{1})\leftrightarrow(\lambda,\omega_{2})],\end{split}$
(68)
$\displaystyle\begin{split}\rho_{EE,ab}^{(ie)}&=\dfrac{1}{2}\int\dfrac{d\omega_{1}d\omega_{2}}{(2\pi)^{2}}E^{\nu}(\omega_{1})E^{\lambda}(\omega_{2})\quantity(d(\omega)\odot[r_{i}^{\nu},d(\omega-\omega_{1})\odot[r_{e}^{\lambda},\rho^{(0)}(\omega-\omega_{1}-\omega_{2})]])_{ab}+[(\nu,\omega_{1})\leftrightarrow(\lambda,\omega_{2})]\\\
&=-\dfrac{i}{2}\int\dfrac{d\omega_{1}\omega_{1}}{(2\pi)^{2}}E^{\nu}(\omega_{1})E^{\lambda}(\omega_{2})d_{ab}(\omega)\quantity[\partial_{\nu}\quantity(d_{ab}(\omega-\omega_{1})\xi_{ab}^{\lambda}f_{ab})-i\quantity(\xi_{aa}^{\nu}-\xi_{bb}^{\nu})d_{ab}(\omega-\omega_{1})\xi_{ab}^{\lambda}f_{ab}]\end{split},$
(69)
$\displaystyle\begin{split}\rho_{EE,ab}^{(ei)}&=\dfrac{1}{2}\int\dfrac{d\omega_{1}d\omega_{2}}{(2\pi)^{2}}E^{\nu}(\omega_{1})E^{\lambda}(\omega_{2})\quantity(d(\omega)\odot[r_{e}^{\nu},d(\omega-\omega_{1})\odot[r_{i}^{\lambda},\rho^{(0)}(\omega-\omega_{1}-\omega_{2})]])_{ab}+[(\nu,\omega_{1})\leftrightarrow(\lambda,\omega_{2})]\\\
&=-\dfrac{i}{2}\int\dfrac{d\omega_{1}d\omega_{2}}{(2\pi)^{2}}E^{\nu}(\omega_{1})E^{\lambda}(\omega_{2})d_{ab}(\omega)d_{aa}(\omega-\omega_{1})\xi_{ab}^{\nu}\partial_{\lambda}f_{ab}2\pi\delta(\omega-\omega_{1}-\omega_{2}),\end{split},$
(70)
$\displaystyle\begin{split}\rho_{EE,ab}^{(ee)}&=\dfrac{1}{2}\int\dfrac{d\omega_{1}d\omega_{2}}{(2\pi)^{2}}E^{\nu}(\omega_{1})E^{\lambda}(\omega_{2})\quantity(d(\omega)\odot[r_{e}^{\nu},d(\omega-\omega_{1})\odot[r_{e}^{\lambda},\rho^{(0)}(\omega-\omega_{1}-\omega_{2})]])_{ab}+[(\nu,\omega_{1})\leftrightarrow(\lambda,\omega_{2})]\\\
&=\dfrac{1}{2}\sum_{c}\int\dfrac{d\omega_{1}d\omega_{2}}{(2\pi)^{2}}E^{\nu}(\omega_{1})E^{\lambda}(\omega_{2})d_{ab}(\omega)\quantity[d_{cb}(\omega-\omega_{1})\xi_{ac}^{\mu}\xi_{cb}^{\mu}f_{bc}-d_{ac}(\omega-\omega_{1})\xi_{cb}^{\mu}\xi_{ac}^{\nu}f_{ca}]2\pi\delta(\omega-\omega_{1}-\omega_{2}).\end{split}$
(71)

Due to the Fermi surface factor $\partial_{\nu}f$,
$\rho^{(ii)}_{EE},\rho^{(ei)}_{EE}$ vanishes in the limit of cold insulators.
In the following calculation, we neglect these components and derive the shift
and injection current formula.

Using $\rho_{EE,ab}$, we can get the following second-order current response,

$\displaystyle\begin{split}J^{\mu}_{EE}(\omega)&=\int\dfrac{d{\bf\it
k}}{(2\pi)^{d}}\sum_{abc}J_{ab}^{\mu}\rho_{EE,ba}(\omega)\\\
&\eqqcolon\int\dfrac{d\omega_{1}d\omega_{2}}{(2\pi)^{2}}{\sigma}^{\mu;\nu\lambda}_{EE}(\omega,\omega_{1},\omega_{2})E^{\nu}(\omega_{1})E^{\lambda}(\omega_{2})2\pi\delta(\omega-\omega_{1}-\omega_{2}).\end{split}$
(72)

As we divide the SPDM into the four different contributions,
$\sigma_{EE}^{\mu;\nu\lambda}(\omega,\omega_{1},\omega_{2})$ is also
decomposed into the corresponding the four components defined as,

$\displaystyle\sigma_{EE}^{\mu;\nu\lambda}=\sigma_{EE,(ii)}^{\mu;\nu\lambda}+\sigma_{EE,(ie)}^{\mu;\nu\lambda}+\sigma_{EE,(ei)}^{\mu;\nu\lambda}+\sigma_{EE,(ee)}^{\mu;\nu\lambda}.$
(73)

Photocurrent conductivity can be obtained by setting

$\displaystyle\omega=0,\omega_{1}=-\Omega,\omega_{2}=\Omega.$ (74)

In the following calculation, we focus on $\sigma_{EE,(ie)}^{\mu;\nu\lambda},$
$\sigma_{EE,(ee)}^{\mu;\nu\lambda}$, from which we will derive shift current
and injection current conductivity.

###  A.2 Injection current

Firstly, we focus on $\sigma_{EE,(ee)}^{\mu;\nu\lambda}$ with the diagonal
component of the current matrix denoted as
$\sigma_{EE,(ee;d)}^{\mu;\nu\lambda}$;

$\displaystyle\begin{split}&\sigma^{\mu;\nu\lambda}_{EE,(ee;d)}(\omega,\omega_{1},\omega_{2})\\\
&=\dfrac{1}{2}\int\dfrac{d{\bf\it k}}{(2\pi)^{d}}\sum_{a\neq
c}J_{aa}^{\mu}d_{aa}(\omega)\quantity[d_{ca}(\omega-\omega_{1})\xi^{\nu}_{ac}\xi^{\lambda}_{ca}f_{ac}-d_{ac}(\omega-\omega_{1})\xi_{ca}^{\nu}\xi_{ac}^{\lambda}f_{cb}]+[(\nu,\omega_{1})\leftrightarrow(\lambda,\omega_{2})]\\\
&=\dfrac{1}{2}\dfrac{1}{\omega+i\eta}\int\dfrac{d{\bf\it
k}}{(2\pi)^{d}}\sum_{a\neq
c}(J_{aa}^{\mu}-J_{cc}^{\mu})\quantity[d_{ca}(\omega-\omega_{1})\xi^{\nu}_{ac}\xi^{\lambda}_{ca}f_{ac}]+[(\nu,\omega_{1})\leftrightarrow(\lambda,\omega_{2})]\\\
&=\dfrac{1}{2}\dfrac{1}{\omega+i\eta}\int\dfrac{d{\bf\it
k}}{(2\pi)^{d}}\sum_{a\neq
c}(J_{aa}^{\mu}-J_{cc}^{\mu})\xi^{\nu}_{ac}\xi^{\lambda}_{ca}f_{ac}\quantity[d_{ca}(\omega-\omega_{1})+d_{ac}(\omega-\omega_{2})]\\\
&=\dfrac{1}{2}\dfrac{1}{\omega+i\eta}\int\dfrac{d{\bf\it
k}}{(2\pi)^{d}}\sum_{a\neq
c}\Delta_{ac}^{\mu}\xi^{\nu}_{ac}\xi^{\lambda}_{ca}f_{ac}\quantity[d_{ca}(\omega-\omega_{1})+d_{ac}(\omega-\omega_{2})].\end{split}$
(75)

Here we introduced velocity difference matrix
$\Delta_{ac}^{\mu}=J_{aa}^{\mu}-J_{cc}^{\mu}=\partial_{\mu}\epsilon_{a}-\partial_{\mu}\epsilon_{b}$.
In the $\omega\to 0$ limit, the expression diverges due to the factor of
$1/\omega+i\eta$. We can eliminate this unphysical divergence by considering
the infinitesimal parameter $\eta$ as finite system relaxation $\gamma$. By
using this method and relation Eq. (59), we can obtain the injection current
formula as

$\displaystyle\sigma_{EE,\text{inj}}^{\mu;\nu\lambda}$
$\displaystyle=\dfrac{\pi}{\gamma}\int\dfrac{d{\bf\it
k}}{(2\pi)^{d}}\sum_{a\neq
c}\Delta_{ac}^{\mu}\xi^{\nu}_{ac}\xi^{\lambda}_{ca}f_{ac}\delta(\Omega-\epsilon_{ca}).$
(76)

Here we decompose the $\sigma_{EE,\text{inj}}^{\mu;\nu\lambda}$ into real and
imaginary part as

$\displaystyle\sigma_{EE,\text{Mnj}}^{\mu;\nu\lambda}$
$\displaystyle=\dfrac{\pi}{\gamma}\int\dfrac{d{\bf\it
k}}{(2\pi)^{d}}\sum_{a\neq
c}\Delta_{ac}^{\mu}\operatorname{Re}\quantity[\xi^{\nu}_{ac}\xi^{\lambda}_{ca}]f_{ac}\delta(\Omega-\epsilon_{ca}),$
(77) $\displaystyle\sigma_{EE,\text{Enj}}^{\mu;\nu\lambda}$
$\displaystyle=i\dfrac{\pi}{\gamma}\int\dfrac{d{\bf\it
k}}{(2\pi)^{d}}\sum_{a\neq
c}\Delta_{ac}^{\mu}\operatorname{Im}\quantity[\xi^{\nu}_{ac}\xi^{\lambda}_{ca}]f_{ac}\delta(\Omega-\epsilon_{ca}).$
(78)

$\sigma_{\text{Mnj}}$ ($\sigma_{\text{Enj}}$) are called magnetic (electric)
injection current, which is induced by linearly (circularly) polarized light.

###  A.3 shift current

Here we derive the shift current formula based on
$\sigma_{EE,(ee)}^{\mu;\nu\lambda}$ and $\sigma_{EE,(ie)}^{\mu;\nu\lambda}$.
First, we consider $\sigma_{EE,(ee)}^{\mu;\nu\lambda}$ with the off-diagonal
component of the velocity operator denoted as
$\sigma^{\mu;\nu\lambda}_{EE,(ee;o)}$;

$\displaystyle\begin{split}&\sigma^{\mu;\nu\lambda}_{EE,(ee;o)}(0,-\Omega,\Omega)\\\
&=\dfrac{1}{2}\int\dfrac{d{\bf\it k}}{(2\pi)^{d}}\sum_{a\neq b\neq
c}J_{ab}^{\mu}d_{ba}(0)\quantity[d_{ca}(\Omega)\xi^{\nu}_{bc}\xi^{\lambda}_{ca}f_{ac}-d_{bc}(\Omega)\xi_{ca}^{\nu}\xi_{bc}^{\lambda}f_{cb}]+[(\nu,-\Omega)\leftrightarrow(\lambda,\Omega)]\\\
&=\dfrac{1}{2}\int\dfrac{d{\bf\it k}}{(2\pi)^{d}}\sum_{a\neq b\neq
c}i\epsilon_{ab}\xi_{ab}^{\mu}d_{ba}(0)\quantity[d_{ca}(\Omega)\xi^{\nu}_{bc}\xi^{\lambda}_{ca}f_{ac}-d_{bc}(\Omega)\xi_{ca}^{\nu}\xi_{bc}^{\lambda}f_{cb}]+[(\nu,-\Omega)\leftrightarrow(\lambda,\Omega)]\\\
&=\dfrac{1}{2}\int\dfrac{d{\bf\it k}}{(2\pi)^{d}}\sum_{a\neq b\neq
c}i\xi_{ab}^{\mu}\quantity[d_{ca}(\Omega)\xi^{\nu}_{bc}\xi^{\lambda}_{ca}f_{ac}-d_{bc}(\Omega)\xi_{ca}^{\nu}\xi_{bc}^{\lambda}f_{cb}]+[(\nu,-\Omega)\leftrightarrow(\lambda,\Omega)]\\\
&=\dfrac{1}{2}\int\dfrac{d{\bf\it k}}{(2\pi)^{d}}\sum_{a\neq b\neq
c}i(\xi_{ab}^{\mu}\xi_{bc}^{\nu}-\xi_{bc}\xi_{ab}^{\nu})\xi_{ca}^{\lambda}f_{ac}d_{ac}(\Omega)+[(\nu,-\Omega)\leftrightarrow(\lambda,\Omega)]\\\
&=\dfrac{1}{2}\int\dfrac{d{\bf\it k}}{(2\pi)^{d}}\sum_{a\neq
c}\quantity(\quantity[D^{\mu}\xi^{\nu}]_{ac}-\quantity[D^{\nu}\xi^{\mu}]_{ac})\xi_{ca}^{\lambda}f_{ac}d_{ac}(\Omega)+[(\nu,-\Omega)\leftrightarrow(\lambda,\Omega)].\end{split}$
(79)

In the last line, we defined $U(1)$-covariant derivative

$\displaystyle\quantity[D^{\mu}O]_{ab}=\partialderivative{O_{ab}}{k^{\mu}}-i\quantity(\xi_{aa}^{\mu}-\xi_{bb}^{\mu})O_{ab},$
(80)

and applied the following formula under the condition of $a\neq c$;

$\displaystyle\quantity[D^{\mu}\xi^{\nu}]_{ac}-\quantity[D^{\nu}\xi^{\mu}]_{ac}=\sum_{b}i(\xi_{ab}^{\mu}\xi_{bc}^{\nu}-\xi_{bc}\xi_{ab}^{\nu}).$
(81)

Next, we consider the $\sigma_{EE,(ie)}^{\mu;\nu\lambda}$ component;

$\displaystyle\begin{split}&\sigma^{\mu;\nu\lambda}_{EE,(ie)}(0;-\Omega,\Omega)\\\
&=-\dfrac{i}{2}\int\dfrac{d{\bf\it k}}{(2\pi)^{d}}\sum_{a\neq
b}J_{ab}^{\mu}d_{ba}(\omega)\quantity[\partialderivative{k^{\nu}}\quantity(d_{ba}(\Omega)f_{ba}\xi_{ba}^{\lambda}-i\quantity(\xi_{aa}^{\nu}-\xi_{bb}^{\mu})d_{ba}(\Omega)f_{ba}\xi_{ba}^{\lambda})]+[(\nu,-\Omega)\leftrightarrow(\lambda,\Omega)]\\\
&=-\dfrac{i}{2}\int\dfrac{d{\bf\it k}}{(2\pi)^{d}}\sum_{a\neq
b}i\epsilon_{ab}\xi_{ab}^{\mu}d_{ba}(\omega)\quantity[\partialderivative{k^{\nu}}\quantity(d_{ba}(\Omega)f_{ba}\xi_{ba}^{\lambda}-i\quantity(\xi_{aa}^{\nu}-\xi_{bb}^{\mu})d_{ba}(\Omega)f_{ba}\xi_{ba}^{\lambda})]+[(\nu,-\Omega)\leftrightarrow(\lambda,\Omega)]\\\
&=\dfrac{1}{2}\int\dfrac{d{\bf\it k}}{(2\pi)^{d}}\sum_{a\neq
b}\xi_{ab}^{\mu}\quantity[\partialderivative{k^{\nu}}\quantity(d_{ba}(\Omega)f_{ba}\xi_{ba}^{\lambda})-i\quantity(\xi_{aa}^{\nu}-\xi_{bb}^{\mu})d_{ba}(\Omega)f_{ba}\xi^{\lambda}_{ba}]+[(\nu,-\Omega)\leftrightarrow(\lambda,\Omega)]\\\
&=-\dfrac{1}{2}\int\dfrac{d{\bf\it k}}{(2\pi)^{d}}\sum_{a\neq
b}\quantity[D^{\nu}\xi^{\mu}]_{ab}d_{ba}(\Omega)f_{ba}\xi_{ba}^{\lambda}+[(\nu,-\Omega)\leftrightarrow(\lambda,\Omega)].\end{split}$
(82)

By summing up these contributions, we obtain

$\displaystyle\begin{split}\sigma^{\mu;\nu\lambda}_{EE,(ee;o)+(ie)}&=\sigma^{\mu;\nu\lambda}_{EE,(ee;o)}(0,-\Omega,\Omega)+\sigma^{\mu;\nu\lambda}_{EE,(ie)}(0;-\Omega,\Omega)\\\
&=\dfrac{1}{2}\int\dfrac{d{\bf\it k}}{(2\pi)^{d}}\sum_{a\neq
b}\quantity[D^{\mu}\xi^{\nu}]_{ab}\xi_{ba}^{\lambda}f_{ab}d_{ba}(\Omega)+[(\nu,-\Omega)\leftrightarrow(\lambda,\Omega)]\\\
&=\dfrac{1}{2}\int\dfrac{d{\bf\it k}}{(2\pi)^{d}}\sum_{a\neq
b}\quantity[D^{\mu}\xi^{\nu}]_{ab}\xi_{ba}^{\lambda}f_{ab}d_{ba}(\Omega)+\quantity[D^{\mu}\xi^{\lambda}]_{ab}\xi_{ba}^{\nu}f_{ab}d_{ba}(-\Omega)\\\
&=\dfrac{1}{2}\int\dfrac{d{\bf\it k}}{(2\pi)^{d}}\sum_{a\neq
b}\quantity[\quantity[D^{\mu}\xi^{\nu}]_{ab}\xi_{ba}^{\lambda}+\quantity[D^{\mu}\xi^{\lambda}]_{ba}\xi_{ab}^{\nu}]f_{ab}\mathcal{P}\dfrac{1}{\Omega-\epsilon_{ba}}\\\
&-\dfrac{i\pi}{2}\int\dfrac{d{\bf\it k}}{(2\pi)^{d}}\sum_{a\neq
b}\quantity[\quantity[D^{\mu}\xi^{\nu}]_{ab}\xi_{ba}^{\lambda}-\quantity[D^{\mu}\xi^{\lambda}]_{ab}\xi_{ba}^{\nu}]f_{ab}\delta(\Omega-\epsilon_{ba}).\end{split}$
(83)

Absorptive part $\delta(\Omega-\epsilon_{ba})$ corresponds to the shift
current conductivity, and we divide the real and imaginary parts as

$\displaystyle\sigma_{EE,\text{shift}}^{\mu;\nu\lambda}=\dfrac{\pi}{2}\int\dfrac{d{\bf\it
k}}{(2\pi)^{d}}\sum_{a\neq
b}\operatorname{Im}\quantity[\quantity[D^{\mu}\xi^{\nu}]_{ab}\xi_{ba}^{\lambda}-\quantity[D^{\mu}\xi^{\lambda}]_{ba}\xi_{ab}^{\nu}]f_{ab}\delta(\Omega-\epsilon_{ba}),$
(84)
$\displaystyle\sigma_{EE,\text{gyro}}^{\mu;\nu\lambda}=-\dfrac{i\pi}{2}\int\dfrac{d{\bf\it
k}}{(2\pi)^{d}}\sum_{a\neq
b}\operatorname{Re}\quantity[\quantity[D^{\mu}\xi^{\nu}]_{ab}\xi_{ba}^{\lambda}-\quantity[D^{\mu}\xi^{\lambda}]_{ba}\xi_{ab}^{\nu}]f_{ab}\delta(\Omega-\epsilon_{ba}).$
(85)

$\sigma_{\text{shift}}$ ($\sigma_{\text{gyro}}$) are called shift (gyration)
currents, which are induced by linearly (circularly) polarized light.

###  A.4 Spin dynamics induced photocurrent

Here, we derive the photocurrent formula induced by spin dynamics. Using the
SPDM $\rho_{\mathrm{SS}}$, current response to the spin field can be expressed
as

$\displaystyle J^{\mu}_{\mathrm{SS}}$ $\displaystyle=\int\dfrac{d{\bf\it
k}}{(2\pi)^{d}}\sum_{abc}J_{ab}^{\mu}\rho_{\mathrm{S}\mathrm{S},ba}^{(2)}(\omega)$
(86)
$\displaystyle\eqqcolon\int\dfrac{d\omega_{1}d\omega_{2}}{(2\pi)^{2}}\sigma^{\mu;\nu\lambda}_{\mathrm{SS}}(\omega,\omega_{1},\omega_{2})\Delta\mathrm{S}^{\nu}(\omega_{1})\Delta\mathrm{S}^{\lambda}(\omega_{2})2\pi\delta(\omega-\omega_{1}-\omega_{2})$
(87)

In the previous subsection, we derive the photocurrent conductivity based on
the perturbative expansion of the von Neumann equation and obtain the shift
current and injection current formula. In the same manner, we can get the
shift current and injection current-like formula, using $\rho_{\mathrm{SS}}$,
just by replacing the Berry connection term $\xi$ related to the external
electric field appearing in the $\sigma_{EE}$ by spin operator $A$. As a
result, we get the shift current and injection-like current formula induced by
spin dynamics;

$\displaystyle\sigma_{\mathrm{S}\mathrm{S},\text{shift}}^{\mu;\nu\lambda}$
$\displaystyle=J^{2}\dfrac{\pi}{2}\int\dfrac{d{\bf\it
k}}{(2\pi)^{d}}\sum_{a\neq
b}\operatorname{Im}\quantity[\quantity[D^{\mu}A^{\nu}]_{ab}A_{ba}^{\lambda}-\quantity[D^{\mu}A^{\lambda}]_{ba}A_{ab}^{\nu}]f_{ab}\delta(\Omega-\epsilon_{ba}),$
(88) $\displaystyle\sigma_{\mathrm{S}\mathrm{S},\text{gyro}}^{\mu;\nu\lambda}$
$\displaystyle=-J^{2}\dfrac{i\pi}{2}\int\dfrac{d{\bf\it
k}}{(2\pi)^{d}}\sum_{a\neq
b}\operatorname{Re}\quantity[\quantity[D^{\mu}A^{\nu}]_{ab}A_{ba}^{\lambda}-\quantity[D^{\mu}A^{\lambda}]_{ba}A_{ab}^{\nu}]f_{ab}\delta(\Omega-\epsilon_{ba}),$
(89) $\displaystyle\sigma_{\mathrm{S}\mathrm{S},\text{Mnj}}^{\mu;\nu\lambda}$
$\displaystyle=J^{2}\dfrac{\pi}{\gamma}\int\dfrac{d{\bf\it
k}}{(2\pi)^{d}}\sum_{a\neq
c}\Delta_{ac}^{\mu}\operatorname{Re}\quantity[A^{\nu}_{ac}A^{\lambda}_{ca}]f_{ac}\delta(\Omega-\epsilon_{ca}),$
(90) $\displaystyle\sigma_{\mathrm{S}\mathrm{S},\text{Enj}}^{\mu;\nu\lambda}$
$\displaystyle=J^{2}\dfrac{i\pi}{\gamma}\int\dfrac{d{\bf\it
k}}{(2\pi)^{d}}\sum_{a\neq
c}\Delta_{ac}^{\mu}\operatorname{Im}\quantity[A^{\nu}_{ac}A^{\lambda}_{ca}]f_{ac}\delta(\Omega-\epsilon_{ca}).$
(91)

###  A.5 Interference of Electric field and Spin dynamics

Following the previous subsection, we derive the photocurrent formula
generated by the interference of the electric field and spin dynamics. Using
$\rho_{E\mathrm{S}}$, photocurrent formula can be written as

$\displaystyle J_{E\mathrm{S}}^{\mu}(\omega)$
$\displaystyle=\int\dfrac{d{\bf\it
k}}{(2\pi)^{d}}\sum_{abc}J_{ab}^{\mu}\rho_{E\mathrm{S},ba}^{(2)}(\omega)$ (92)
$\displaystyle\eqqcolon\int\dfrac{d\omega_{1}d\omega_{2}}{(2\pi)^{2}}\tilde{\sigma}_{E\mathrm{S}}^{\mu;\nu\lambda}(\omega,\omega_{1},\omega_{2})E^{\nu}(\omega_{1})\Delta\mathrm{S}^{\lambda}(\omega_{2}),$
(93)

where $\tilde{\sigma}_{E\mathrm{S}}^{\mu;\nu\lambda}$ can be classified into
the following four component.

$\displaystyle\sigma_{E\mathrm{S},\text{shift}}^{\mu;\nu\lambda}$
$\displaystyle=J\dfrac{\pi}{2}\int\dfrac{d{\bf\it k}}{(2\pi)^{d}}\sum_{a\neq
b}\operatorname{Im}\quantity[\quantity[D^{\mu}\xi^{\nu}]_{ab}A_{ba}^{\lambda}-\quantity[D^{\mu}A^{\lambda}]_{ba}\xi_{ab}^{\nu}]f_{ab}\delta(\Omega-\epsilon_{ba}),$
(94) $\displaystyle\sigma_{E\mathrm{S},\text{gyro}}^{\mu;\nu\lambda}$
$\displaystyle=-J\dfrac{i\pi}{2}\int\dfrac{d{\bf\it k}}{(2\pi)^{d}}\sum_{a\neq
b}\operatorname{Re}\quantity[\quantity[D^{\mu}\xi^{\nu}]_{ab}A_{ba}^{\lambda}-\quantity[D^{\mu}A^{\lambda}]_{ba}\xi_{ab}^{\nu}]f_{ab}\delta(\Omega-\epsilon_{ba}),$
(95) $\displaystyle\sigma_{E\mathrm{S},\text{Mnj}}^{\mu;\nu\lambda}$
$\displaystyle=J\dfrac{\pi}{\gamma}\int\dfrac{d{\bf\it
k}}{(2\pi)^{d}}\sum_{a\neq
c}\Delta_{ac}^{\mu}\operatorname{Re}\quantity[\xi^{\nu}_{ac}A^{\lambda}_{ca}]f_{ac}\delta(\Omega-\epsilon_{ca}),$
(96) $\displaystyle\sigma_{E\mathrm{S},\text{Enj}}^{\mu;\nu\lambda}$
$\displaystyle=J\dfrac{i\pi}{\gamma}\int\dfrac{d{\bf\it
k}}{(2\pi)^{d}}\sum_{a\neq
c}\Delta_{ac}^{\mu}\operatorname{Im}\quantity[\xi^{\nu}_{ac}A^{\lambda}_{ca}]f_{ac}\delta(\Omega-\epsilon_{ca}).$
(97)

##  Appendix B Symmetry analysis of the photocurrent

In this section, we show the symmetry analysis of photocurrent in the presence
of spin dynamics. Due to the $\theta 2_{x}$ symmetry of the system and phase
matching between the different fields, the photocurrent generation is
restricted. Here, we consider a one-dimensional system, applying the electric
field along the $z$ direction and observing the photocurrent along the $z$
direction, as in the case of the main text. In the following calculation, we
assume the wave vector ${\bf\it k}$ to be scaler.

Under $\theta 2_{x}$ operation, the wave vector along $z$ direction transforms
as $k\rightarrow k$. Therefore, Bloch states that momentum $k$ is related to
those with the same momentum $k$.

$\displaystyle\theta 2_{x}\ket{u_{a}(k)}=\ket{{u_{\tilde{a}}}(k)}.$ (98)

Since $(\theta 2_{x})^{2}=1$, there is no Kramers degeneracy. Therefore, by
choosing the phase factor, we can assume

$\displaystyle\theta 2_{x}\ket{u_{a}(k)}=\ket{{u_{a}}(k)}.$ (99)

Due to the $\theta 2_{x}$ symmetry, Berry connection $\xi_{ab}^{z}(k)$
satisfies the following relation

$\displaystyle\begin{split}\xi_{ab}^{z}(k)&=\matrixelement{u_{a}(k)}{i\partial_{k}}{u_{b}(k)}\\\
&=\matrixelement{u_{b}(k)}{(\theta 2_{x})i\partial_{k}(\theta
2_{x})^{-1}}{u_{a}(k)}\\\
&=\matrixelement{u_{b}(k)}{-i\partial_{k}}{u_{a}(k)}\\\
&=-\xi_{ba}^{z}(k).\end{split}$ (100)

Moreover, the $A^{\nu}$ under the $\theta 2_{x}$ symmetry satisfies the
following relation.

$\displaystyle\begin{split}A_{ab}^{\nu}&=\matrixelement{u_{a}(k)}{A^{\nu}}{u_{b}(k)}\\\
&=\matrixelement{u_{b}(k)}{(\theta 2_{x})(A^{\nu})^{\dagger}(\theta
2_{x})}{u_{a}(k)}\\\ &=\matrixelement{u_{b}(k)}{(\theta 2_{x})(A^{\nu})(\theta
2_{x})}{u_{a}(k)}\\\
&=\sigma_{A^{\nu}}\matrixelement{u_{b}(k)}{A^{\nu}}{u_{a}(k)}\\\
&=\sigma_{A^{\nu}}A_{ba}^{\nu}.\end{split}$ (101)

Here $\sigma_{A^{\nu}}$ is the sign under $\theta 2_{x}$ operation. As for the
spin operator coupled to the $\alpha$-mode, $\sigma_{A^{\nu}}$ satisfies the
following relation.

$\displaystyle\sigma_{\sigma_{x}\tau_{z}}=-1,\sigma_{\sigma_{y}\tau_{0}}=1,\sigma_{\sigma_{z}\tau_{0}}=1.$
(102)

###  B.1 Elecctric field induced photocurrent

Drawing on the previous subsection, the photocurrent along the $z$ direction
induced by the electric field along the $z$ direction can be written as

$\displaystyle
J_{0}^{z}=\int\dfrac{d\Omega}{2\pi}\sigma_{EE}^{z;zz}E^{z}(-\Omega)E^{z}(\Omega),$
(103)

where $\sigma^{zzz}(\Omega)$ can be classified into following four
contributions,

$\displaystyle\sigma_{EE,\text{shift}}^{z;zz}$
$\displaystyle=\dfrac{\pi}{2}\int\dfrac{dk}{2\pi}\sum_{a\neq
b}\operatorname{Im}\quantity[\quantity[D^{z}\xi^{z}]_{ab}\xi_{ba}^{z}+\quantity[D^{z}\xi^{z}]_{ab}\xi_{ba}^{z}]f_{ab}\delta(\Omega-\epsilon_{ab}),$
(104) $\displaystyle\sigma_{EE,\text{gyro}}^{z;zz}$
$\displaystyle=-\dfrac{i\pi}{2}\int\dfrac{dk}{2\pi}\sum_{a\neq
b}\operatorname{Re}\quantity[\quantity[D^{\mu}\xi^{\nu}]_{ab}\xi_{ba}^{\lambda}-\quantity[D^{\mu}\xi^{\lambda}]_{ac}\xi_{ba}^{\nu}]f_{ab}\delta(\Omega-\epsilon_{ab}),$
(105) $\displaystyle\sigma_{EE,\text{Mnj}}^{z;zz}$
$\displaystyle=\dfrac{\pi}{\gamma}\int\dfrac{dk}{2\pi}\sum_{a\neq
c}\Delta_{ab}^{z}\operatorname{Re}\quantity[\xi^{z}_{ab}\xi^{z}_{ba}]f_{ab}\delta(\Omega-\epsilon_{ba}),$
(106) $\displaystyle\sigma_{EE,\text{Enj}}^{z;zz}$
$\displaystyle=i\dfrac{\pi}{\gamma}\int\dfrac{dk}{2\pi}\sum_{a\neq
b}\Delta_{ab}^{z}\operatorname{Im}\quantity[\xi^{z}_{ab}\xi^{z}_{ba}]f_{ab}\delta(\Omega-\epsilon_{ba}).$
(107)

In Eq. (103), since the left hand side is real, and
$E^{z}(\Omega)\quantity[E^{z}(\Omega)]^{\ast}$ is also real,
$\sigma^{zzz}_{EE}(\Omega)$ should be real. Therefore only
$\sigma_{\text{shift}}^{z;zz},\sigma_{\text{Mnj}}^{z;zz}$ can contribute to
the photocurrent generation. Moreover due to the $\theta 2_{x}$ in the system,
covariant derivative $\quantity[D^{z}\xi^{z}]_{ab}\xi_{ba}$ satisfies
following relation

$\displaystyle\begin{split}\quantity[D^{z}\xi^{z}]_{ab}\xi^{z}_{ba}&=\quantity[\partialderivative{\xi^{z}_{ab}}{k^{z}}-i\quantity(\xi^{z}_{aa}-\xi_{bb}^{z})\xi_{ab}^{z}]\xi_{ba}^{z}\\\
&=-\quantity[-\partialderivative{\xi^{z}_{ba}(k)}{k}-i\quantity(\xi^{z}_{aa}(k)-\xi^{z}_{bb}(k))\xi^{z}_{ba}(k)]\xi^{z}_{ab}(k)\\\
&=\quantity[\partialderivative{\xi^{z}_{ba}(k)}{k}-i\quantity(\xi^{z}_{bb}(k)-\xi^{z}_{aa}(k))\xi^{z}_{ba}(k)]\xi^{z}_{ab}(k)\\\
&=\quantity[D^{z}(k)\xi^{z}(k)]_{ba}\xi^{z}_{ab}(k)\\\
&=\quantity(\quantity[D^{z}(k)\xi^{z}(k)]_{ab}\xi^{z}_{ba}(k))^{\ast}.\end{split}$
(108)

Therefore $\sigma_{\text{shift}}$ contribution vanishes under the $\theta
2_{x}$ symmetry, while $\sigma_{\text{Mnj}}$ remain intact.

###  B.2 Spin dynamics induced photocurrent

Photocurrent induced solely by collective spin dynamics can be expressed as

$\displaystyle\begin{split}J^{\mu}_{SS}&=\int\dfrac{d\Omega}{2\pi}\sigma_{\mathrm{SS}}^{\mu;\nu\lambda}(0;-\Omega,\Omega)\Delta\mathrm{S}^{\nu}(-\Omega)\Delta\mathrm{S}^{\lambda}(\Omega)\\\
&=\int\dfrac{d\Omega}{2\pi}\sigma_{\mathrm{SS}}^{\mu;\nu\lambda}(0;-\Omega,\Omega)\quantity[\Delta\mathrm{S}^{\nu}(\Omega)]^{\ast}\Delta\mathrm{S}^{\lambda}(\Omega).\end{split}$
(109)

Here $\sigma_{\mathrm{SS}}$ can be classified into the following four
components

$\displaystyle\sigma_{\mathrm{S}\mathrm{S},\text{shift}}^{\mu;\nu\lambda}$
$\displaystyle=J^{2}\dfrac{\pi}{2}\int\dfrac{d{\bf\it
k}}{(2\pi)^{d}}\sum_{a\neq
b}\operatorname{Im}\quantity[\quantity[D^{\mu}A^{\nu}]_{ab}A_{ba}^{\lambda}-\quantity[D^{\mu}A^{\lambda}]_{ba}A_{ab}^{\nu}]f_{ab}\delta(\Omega-\epsilon_{ba}),$
(110)
$\displaystyle\sigma_{\mathrm{S}\mathrm{S},\text{gyro}}^{\mu;\nu\lambda}$
$\displaystyle=-J^{2}\dfrac{i\pi}{2}\int\dfrac{d{\bf\it
k}}{(2\pi)^{d}}\sum_{a\neq
b}\operatorname{Re}\quantity[\quantity[D^{z}A^{z}]_{ab}A_{ba}^{\lambda}-\quantity[D^{z}A^{z}]_{ba}A_{ab}^{z}]f_{ab}\delta(\Omega-\epsilon_{ba}),$
(111) $\displaystyle\sigma_{\mathrm{S}\mathrm{S},\text{Mnj}}^{\mu;\nu\lambda}$
$\displaystyle=J^{2}\dfrac{\pi}{\gamma}\int\dfrac{d{\bf\it
k}}{(2\pi)^{d}}\sum_{a\neq
c}\Delta_{ac}^{\mu}\operatorname{Re}\quantity[A^{\nu}_{ac}A^{\lambda}_{ca}]f_{ac}\delta(\Omega-\epsilon_{ca}),$
(112) $\displaystyle\sigma_{\mathrm{S}\mathrm{S},\text{Enj}}^{\mu;\nu\lambda}$
$\displaystyle=J^{2}\dfrac{i\pi}{\gamma}\int\dfrac{d{\bf\it
k}}{(2\pi)^{d}}\sum_{a\neq
c}\Delta_{ac}^{\mu}\operatorname{Im}\quantity[A^{\nu}_{ac}A^{\lambda}_{ca}]f_{ac}\delta(\Omega-\epsilon_{ca}).$
(113)

By considering the $\theta 2_{x}$ operation, the matrix elements in each
photocurrent conductivity satisfy the following relations.

$\displaystyle\operatorname{Re}\quantity[\quantity[D^{z}A^{\nu}]_{ab}A_{ba}^{\lambda}-\quantity[D^{z}A^{\lambda}]_{ba}A_{ab}^{\nu}]$
$\displaystyle=\dfrac{1+\sigma_{A^{\nu}}\sigma_{A^{\lambda}}}{2}\quantity(\quantity[D^{z}A^{\nu}]_{ab}A_{ba}^{\lambda}-\quantity[D^{z}A^{\lambda}]_{ba}A_{ab}^{\nu}),$
(114)
$\displaystyle\operatorname{Im}\quantity[\quantity[D^{z}A^{\nu}]_{ab}A_{ba}^{\lambda}-\quantity[D^{z}A^{\lambda}]_{ba}A_{ab}^{\nu}]$
$\displaystyle=\dfrac{1-\sigma_{A^{\nu}}\sigma_{A^{\lambda}}}{2i}\quantity(\quantity[D^{z}A^{\nu}]_{ab}A_{ba}^{\lambda}-\quantity[D^{z}A^{\lambda}]_{ba}A_{ab}^{\nu}),$
(115) $\displaystyle\operatorname{Re}\quantity[A_{ab}^{\nu}A_{ba}^{\lambda}]$
$\displaystyle=\dfrac{1+\sigma_{A^{\nu}}\sigma_{A^{\lambda}}}{2}A_{ab}^{\nu}A_{ba}^{\lambda},$
(116) $\displaystyle\operatorname{Im}\quantity[A_{ab}^{\nu}A_{ba}^{\lambda}]$
$\displaystyle=\dfrac{1-\sigma_{A^{\nu}}\sigma_{A^{\lambda}}}{2}A_{ab}^{\nu}A_{ba}^{\lambda}.$
(117)

In the equation Eq. (109), when two fields $\mathrm{S}^{\nu}$ and
$\mathrm{S}^{\lambda}$ are in-phase (out-of-phase), the quantity
$\quantity[\Delta\mathrm{S}^{\nu}(\Omega)]^{\ast}\Delta\mathrm{S}^{\lambda}(\Omega)$
becomes real (pure-imaginary). Consequently, the photocurrent induced by two
fields that are in phase with each other can be characterized by the real part
of $\sigma_{\mathrm{SS}}$, namely $\sigma_{\mathrm{SS,shift}}$ and
$\sigma_{\mathrm{SS,Mnj}}$. Conversely, the photocurrent induced by two fields
that are out-of-phase with each other can be characterized by the imaginary
part of $\sigma_{\mathrm{SS}}$, namely $\sigma_{\mathrm{SS,gyro}}$ and
$\sigma_{\mathrm{SS,Enj}}$. In our calculation, we confirmed that
$\mathrm{M}^{y}$ and $\mathrm{L}^{z}$ are in-phase each other, while
$\mathrm{L}^{x}$ is not perfectly in or out of phase with $\mathrm{M}^{y}$ and
$\mathrm{L}^{z}$. Combining this phase analysis with the $\theta 2_{x}$
symmetry constraints in Eq. (102), we can identify the photocurrent
contribution as summarized in TABLE 2. Although we consider the photocurrent
response to the linearly polarized light, $\sigma_{\mathrm{SS}}$ contains the
$\sigma_{\text{Enj}}$ which is the counterpart of circularly polarized induced
photocurrent. Moreover, the photocurrent induced by spin dynamics can contain
the shift current contribution, namely $\sigma_{\text{shift}}$ and
$\sigma_{\text{gyro}}$, which may survive in the disordered system [54].

###  B.3 Interference of Electric Field and Spin Dynamics

Photocurrent generated by the interference of electric field and spin dynamics
can be described as follows.

$\displaystyle
J^{z}_{E\mathrm{S}}=\int\dfrac{d\Omega}{2\pi}\sigma_{E\mathrm{S}}^{z;\nu\lambda}(0;-\Omega,\Omega)E^{z}(-\Omega)\Delta\mathrm{S}^{\lambda}(\Omega).$
(118)

Here $\sigma_{E\mathrm{S}}^{z;z\lambda}$ have four contributions and are
explicitly written as follows.

$\displaystyle\sigma_{E\mathrm{S},\text{shift}}^{z;z\lambda}$
$\displaystyle=J\dfrac{\pi}{2}\int\dfrac{d{\bf\it k}}{(2\pi)^{d}}\sum_{a\neq
b}\operatorname{Im}\quantity[\quantity[D^{z}\xi^{z}]_{ab}A_{ba}^{\lambda}-\quantity[D^{z}A^{\lambda}]_{ba}\xi_{ab}^{z}]f_{ab}\delta(\Omega-\epsilon_{ba}),$
(119) $\displaystyle\sigma_{E\mathrm{S},\text{gyro}}^{z;z\lambda}$
$\displaystyle=-J\dfrac{i\pi}{2}\int\dfrac{d{\bf\it k}}{(2\pi)^{d}}\sum_{a\neq
b}\operatorname{Re}\quantity[\quantity[D^{z}\xi^{z}]_{ab}A_{ba}^{\lambda}-\quantity[D^{z}A^{\lambda}]_{ba}\xi_{ab}^{z}]f_{ab}\delta(\Omega-\epsilon_{ba}),$
(120) $\displaystyle\sigma_{E\mathrm{S},\text{Mnj}}^{z;z\lambda}$
$\displaystyle=J\dfrac{\pi}{\gamma}\int\dfrac{d{\bf\it
k}}{(2\pi)^{d}}\sum_{a\neq
c}\Delta_{ac}^{z}\operatorname{Re}\quantity[\xi^{z}_{ac}A^{\lambda}_{ca}]f_{ac}\delta(\Omega-\epsilon_{ca}),$
(121) $\displaystyle\sigma_{E\mathrm{S},\text{Enj}}^{z;z\lambda}$
$\displaystyle=J\dfrac{i\pi}{\gamma}\int\dfrac{d{\bf\it
k}}{(2\pi)^{d}}\sum_{a\neq
c}\Delta_{ac}^{z}\operatorname{Im}\quantity[\xi^{z}_{ac}A^{\lambda}_{ca}]f_{ac}\delta(\Omega-\epsilon_{ca}).$
(122)

Under the $\theta 2_{x}$ operation, the matrix element related to photocurrent
generation satisfies the following relations.

$\displaystyle\operatorname{Im}\quantity[\quantity[D^{z}\xi^{\nu}]_{ab}A_{ba}^{\lambda}-\quantity[D^{z}A^{\lambda}]_{ab}\xi_{ba}^{\lambda}]$
$\displaystyle=\dfrac{1+\sigma_{A^{\lambda}}}{2i}\quantity[\quantity[D^{z}\xi^{\nu}]_{ab}A_{ba}^{\lambda}-\quantity[D^{z}A^{\lambda}]_{ab}\xi_{ba}^{\lambda}],$
(123)
$\displaystyle\operatorname{Re}\quantity[\quantity[D^{z}\xi^{\nu}]_{ab}A_{ba}^{\lambda}-\quantity[D^{z}A^{\lambda}]_{ab}\xi_{ba}^{\lambda}]$
$\displaystyle=\dfrac{1-\sigma_{A^{\lambda}}}{2}\quantity[\quantity[D^{z}\xi^{\nu}]_{ab}A_{ba}^{\lambda}-\quantity[D^{z}A^{\lambda}]_{ab}\xi_{ba}^{\lambda}],$
(124)
$\displaystyle\operatorname{Re}\quantity[\xi_{ac}^{\nu}A_{ca}^{\lambda}]$
$\displaystyle=\dfrac{1-\sigma_{A}^{\lambda}}{2}\xi_{ac}^{\nu}A_{ca}^{\lambda},$
(125)
$\displaystyle\operatorname{Im}\quantity[\xi_{ac}^{\nu}A_{ca}^{\lambda}]$
$\displaystyle=\dfrac{1+\sigma_{A}^{\lambda}}{2}\xi_{ac}^{\nu}A_{ca}^{\lambda}.$
(126)

Besides the $\theta 2_{x}$ symmetry restriction, phase degrees of freedom
between the electric field and fictitious spin field play an important role in
determining the photocurrent contribution. Rewriting the equation Eq. (118) by
using the electromagnetic susceptibility $\chi_{\mathrm{S}^{\lambda}E}$ to the
electric field

$\displaystyle\begin{split}J_{ES}^{z}&=\int\dfrac{d\Omega}{2\pi}\sigma_{E\mathrm{S}}^{\mu;\nu\lambda}(0;-\Omega,\Omega)\chi_{\mathrm{S}^{\lambda}E}(\omega)E^{z}(-\Omega)E^{z}(\Omega)\\\
&=\int\dfrac{d\Omega}{2\pi}\sigma_{E\mathrm{S}}^{\mu;\nu\lambda}(0;-\Omega,\Omega)\quantity[\operatorname{Re}\chi_{\mathrm{S}^{\lambda}E^{z}}(\Omega)+i\operatorname{Im}\chi_{\mathrm{S}^{\lambda}(\Omega)E^{z}}]\quantity[E^{z}(\Omega)]^{\ast}E^{z}(\Omega)\\\
&=\int\dfrac{d\Omega}{2\pi}\quantity[\operatorname{Re}\sigma_{E\mathrm{S}}^{\mu;\nu\lambda}(0;-\Omega,\Omega)\operatorname{Re}\chi_{\mathrm{S}^{\lambda}E^{z}}(\Omega)-\operatorname{Im}\sigma_{E\mathrm{S}}^{\mu;\nu\lambda}(0;-\Omega,\Omega)\operatorname{Im}\chi_{\mathrm{S}^{\lambda}E^{z}}(\Omega)]\quantity[E^{z}(\Omega)]^{\ast}E^{z}(\Omega).\end{split}$
(127)

Here, we used the fact that the left-hand side is real, and
$\quantity[E^{z}(\Omega)]^{\ast}E^{z}(\Omega)$ is also real. By considering
the symmetry constraints, the photocurrent generation by the interference of
electric field and spin dynamics are summarized in the table 2. Unlike the
case of IPA, the generation of various kinds of photocurrent can be allowed in
the $\sigma_{\text{ES}}$.

## References

  * Williams [1960] R. Williams, The Journal of Chemical Physics 32, 1505 (1960).
  * Fridkin [2001] V. M. Fridkin, Crystallography Reports 46, 654 (2001).
  * von Baltz and Kraut [1981] R. von Baltz and W. Kraut, Physical Review B 23, 5590 (1981).
  * Shockley and Queisser [1961] W. Shockley and H. J. Queisser, Journal of Applied Physics 32, 510 (1961).
  * Spanier  _et al._ [2016] J. E. Spanier, V. M. Fridkin, A. M. Rappe, A. R. Akbashev, A. Polemi, Y. Qi, Z. Gu, S. M. Young, C. J. Hawley, D. Imbrenda, G. Xiao, A. L. Bennett-Jackson, and C. L. Johnson, Nature Photonics 10, 611 (2016).
  * Sipe and Shkrebtii [2000] J. E. Sipe and A. I. Shkrebtii, Physical Review B 61, 5337 (2000).
  * Ahn  _et al._ [2020] J. Ahn, G.-Y. Guo, and N. Nagaosa, Physical Review X 10, 041041 (2020).
  * Watanabe and Yanase [2021] H. Watanabe and Y. Yanase, Physical Review X 11, 011001 (2021).
  * Young and Rappe [2012] S. M. Young and A. M. Rappe, Physical Review Letters 109, 116601 (2012).
  * Sotome  _et al._ [2019] M. Sotome, M. Nakamura, J. Fujioka, M. Ogino, Y. Kaneko, T. Morimoto, Y. Zhang, M. Kawasaki, N. Nagaosa, Y. Tokura, and N. Ogawa, Proceedings of the National Academy of Sciences 116, 1929 (2019).
  * Sotome  _et al._ [2021] M. Sotome, M. Nakamura, T. Morimoto, Y. Zhang, G.-Y. Guo, M. Kawasaki, N. Nagaosa, Y. Tokura, and N. Ogawa, Physical Review B 103, L241111 (2021).
  * Tan  _et al._ [2016] L. Z. Tan, F. Zheng, S. M. Young, F. Wang, S. Liu, and A. M. Rappe, npj Computational Materials 2, 16026 (2016).
  * Osterhoudt  _et al._ [2019] G. B. Osterhoudt, L. K. Diebel, M. J. Gray, X. Yang, J. Stanco, X. Huang, B. Shen, N. Ni, P. J. W. Moll, Y. Ran, and K. S. Burch, Nature Materials 18, 471 (2019).
  * Braun  _et al._ [2016] L. Braun, G. Mussler, A. Hruban, M. Konczykowski, T. Schumann, M. Wolf, M. MÃ¼nzenberg, L. Perfetti, and T. Kampfrath, Nature Communications 7 (2016), 10.1038/ncomms13259.
  * Yu  _et al._ [2016] W. J. Yu, Q. A. Vu, H. Oh, H. G. Nam, H. Zhou, S. Cha, J.-Y. Kim, A. Carvalho, M. Jeong, H. Choi, A. H. C. Neto, Y. H. Lee, and X. Duan, Nature Communications 7, 13278 (2016).
  * Morimoto and Nagaosa [2018] T. Morimoto and N. Nagaosa, Scientific Reports 8, 2973 (2018).
  * Michishita and Peters [2021] Y. Michishita and R. Peters, Physical Review B 103, 195133 (2021).
  * Morimoto and Nagaosa [2016] T. Morimoto and N. Nagaosa, Physical Review B 94, 035117 (2016).
  * Chan  _et al._ [2021] Y.-H. Chan, D. Y. Qiu, F. H. D. Jornada, S. G. Louie, and F. H. D. J. Performed,  118 (2021), 10.1073/pnas.1906938118/-/DCSupplemental.
  * Okamura  _et al._ [2022] Y. Okamura, T. Morimoto, N. Ogawa, Y. Kaneko, G.-Y. Guo, M. Nakamura, M. Kawasaki, N. Nagaosa, Y. Tokura, and Y. Takahashi, Proceedings of the National Academy of Sciences 119 (2022), 10.1073/pnas.2122313119.
  * Morimoto and Nagaosa [2019] T. Morimoto and N. Nagaosa, Physical Review B 100, 235138 (2019).
  * Kaneko  _et al._ [2021] T. Kaneko, Z. Sun, Y. Murakami, D. GoleÅ¾, and A. J. Millis, Physical Review Letters 127, 127402 (2021).
  * Tokura  _et al._ [2014] Y. Tokura, S. Seki, and N. Nagaosa, Reports on Progress in Physics 77, 076501 (2014).
  * Spaldin and Ramesh [2019] N. A. Spaldin and R. Ramesh, Nature Materials 18, 203 (2019).
  * Fiebig [2005] M. Fiebig, Journal of Physics D: Applied Physics 38, R123 (2005).
  * Jungwirth  _et al._ [2016] T. Jungwirth, X. Marti, P. Wadley, and J. Wunderlich, Nature Nanotechnology 11, 231 (2016).
  * Baltz  _et al._ [2018] V. Baltz, A. Manchon, M. Tsoi, T. Moriyama, T. Ono, and Y. Tserkovnyak, Rev. Mod. Phys. 90, 015005 (2018).
  * Manchon  _et al._ [2019] A. Manchon, J. Železný, I. M. Miron, T. Jungwirth, J. Sinova, A. Thiaville, K. Garello, and P. Gambardella, Rev. Mod. Phys. 91, 035004 (2019).
  * Amin  _et al._ [2020] O. J. Amin, K. W. Edmonds, and P. Wadley, Applied Physics Letters 117 (2020), 10.1063/5.0013917.
  * Okumura  _et al._ [2021] S. Okumura, T. Morimoto, Y. Kato, and Y. Motome, Physical Review B 104, L180407 (2021).
  * Zhang  _et al._ [2019] Y. Zhang, T. Holder, H. Ishizuka, F. de Juan, N. Nagaosa, C. Felser, and B. Yan, Nature Communications 10, 3783 (2019).
  * Watanabe and Yanase [2020] H. Watanabe and Y. Yanase, Physical Review Research 2, 043081 (2020).
  * Ono and Ishihara [2021] A. Ono and S. Ishihara, npj Computational Materials 7, 171 (2021).
  * Yanase [2014] Y. Yanase, Journal of the Physical Society of Japan 83, 014703 (2014).
  * Å½eleznÃ½  _et al._ [2014] J. Å½eleznÃ½, H. Gao, K. VÃ½bornÃ½, J. Zemen, J. MaÅ¡ek, A. Manchon, J. Wunderlich, J. Sinova, and T. Jungwirth, Physical Review Letters 113, 157201 (2014).
  * Hayami  _et al._ [2014] S. Hayami, H. Kusunose, and Y. Motome, Physical Review B 90, 024432 (2014).
  * Murakami and SchÃ¼ler [2022] Y. Murakami and M. SchÃ¼ler, Physical Review B 106, 035204 (2022).
  * Yue and Gaarde [2022] L. Yue and M. B. Gaarde, Journal of the Optical Society of America B 39, 535 (2022).
  * Rezende  _et al._ [2019] S. M. Rezende, A. Azevedo, and R. L. RodrÃguez-SuÃ¡rez, Journal of Applied Physics 126, 151101 (2019).
  * Watanabe and Yanase [2017] H. Watanabe and Y. Yanase, Physical Review B 96, 064432 (2017).
  * Hayami  _et al._ [2018] S. Hayami, M. Yatsushiro, Y. Yanagi, and H. Kusunose, Physical Review B 98, 165110 (2018).
  * Morimoto  _et al._ [2021] T. Morimoto, S. Kitamura, and S. Okumura, Physical Review B 104, 075139 (2021).
  * Murakami  _et al._ [2020] Y. Murakami, D. GoleÅ¾, T. Kaneko, A. Koga, A. J. Millis, and P. Werner, Physical Review B 101, 195118 (2020).
  * Kurebayashi and Nagaosa [2021] D. Kurebayashi and N. Nagaosa, Communications Physics 4, 260 (2021).
  * Cook  _et al._ [2017] A. M. Cook, B. M. Fregoso, F. de Juan, S. Coh, and J. E. Moore, Nature Communications 8, 14176 (2017).
  * Nagaosa and Tokura [2013] N. Nagaosa and Y. Tokura, Nature Nanotechnology 8, 899 (2013).
  * Tokura and Kanazawa [2021] Y. Tokura and N. Kanazawa, Chemical Reviews 121, 2857 (2021).
  * Zhang  _et al._ [2016] X. Zhang, Y. Zhou, and M. Ezawa, Scientific Reports 6, 24795 (2016).
  * Okumura  _et al._ [2020] S. Okumura, S. Hayami, Y. Kato, and Y. Motome, Physical Review B 101, 144416 (2020).
  * Li  _et al._ [2020] H. Li, S. Liu, C. Liu, J. Zhang, Y. Xu, R. Yu, Y. Wu, Y. Zhang, and S. Fan, Physical Chemistry Chemical Physics 22, 556 (2020).
  * Fei  _et al._ [2020] R. Fei, W. Song, and L. Yang, Physical Review B 102, 035440 (2020).
  * Adams and Blount [1959] E. Adams and E. Blount, Journal of Physics and Chemistry of Solids 10, 286 (1959).
  * Ventura  _et al._ [2017] G. B. Ventura, D. J. Passos, J. M. B. L. dos Santos, J. M. V. P. Lopes, and N. M. R. Peres, Physical Review B 96, 035431 (2017).
  * Hatada  _et al._ [2020] H. Hatada, M. Nakamura, M. Sotome, Y. Kaneko, N. Ogawa, T. Morimoto, Y. Tokura, and M. Kawasaki, Proceedings of the National Academy of Sciences 117, 20411 (2020).