# Fractal and multifractal properties of a family of fractal networks

Bao-Gen Li1, Zu-Guo Yu1,2 and Yu Zhou1
1 Hunan Key Laboratory for Computation and Simulation in Science and
Engineering and
Key Laboratory of Intelligent Computing and Information Processing of Ministry
of Education,
Xiangtan University, Xiangtan, Hunan 411105, China.
2School of Mathematical Sciences, Queensland University of Technology,
GPO Box 2434, Brisbane, Q4001, Australia. Corresponding author, email:
<EMAIL_ADDRESS>
###### Abstract

In this work, we study the fractal and multifractal properties of a family of
fractal networks introduced by Gallos et al. ( Proc. Natl. Acad. Sci. U.S.A.,
2007, 104: 7746). In this fractal network model, there is a parameter $e$
which is between $0$ and $1$, and allows for tuning the level of fractality in
the network. Here we examine the multifractal behavior of these networks,
dependence relationship of fractal dimension and the multifractal parameters
on the parameter $e$. First, we find that the empirical fractal dimensions of
these networks obtained by our program coincide with the theoretical formula
given by Song et al. ( Nat. Phys, 2006, 2: 275). Then from the shape of the
$\tau(q)$ and $D(q)$ curves, we find the existence of multifractality in these
networks. Last, we find that there exists a linear relationship between the
average information dimension $<D(1)>$ and the parameter $e$.

Key words: Complex network, scale-free, multifractality, box covering
algorithm.

PACS: 89.75.Hc, 05.45.Df, 47.53.+n

##  1 Introduction

Complex networks have caused extensive attention due to their close connection
with so many real-world systems, such as the world-wide web, the internet,
energy landscapes, and biological and social systems [1]. The fractality and
percolation transition [2], fractal transition [3] in complex networks, and
properties of a scale-free Koch networks [6, 5, 4] have turned to be hot
topics in recent years.

Fractal analysis (using the fractal dimension) is a useful method to describe
global properties of complex fractal sets [7, 8, 9]. Song et al. [1, 10]
proposed an algorithm to calculate the fractal dimension of complex networks
which can unfold their self-similar property. They mentioned that the box
counting fractal analysis is an effective tool for the further study of
complex networks. But the fractal analysis is not enough when the object
studied can not be described by a single fractal dimension. It has been found
that the multifractal analysis (MFA) is a powerful tool in both the theory and
practice to describe the spatial heterogeneity of fractal object
systematically [11, 12]. The MFA was originally raised to handle turbulence
data, and now it has been successfully applied in many fields, such as
financial modelling [13, 14], biological systems [15, 16, 17, 18, 19, 20, 21,
22, 23, 24] and geophysical systems [25, 26, 27, 28, 29, 30, 31]. Lee and Jung
[32] found that MFA is the best tool to describe the probability distribution
of the clustering coefficient of a complex network. Furuya and Yakubo [33]
analytically and numerically demonstrated the possibility that the fractal
property of a scale-free network cannot be characterized by a single fractal
dimension when the network takes a multifractal structure. Almost at the same
time, Wang et al. [34] proposed a modified fixed-size box-counting algorithm
to study the multifractal property of complex networks.

In this paper, we study the fractal and multifractal properties of a family of
complex networks introduced by Gallos et al. [35]. In order to imitate the
fractal property of many scale-free networks found in nature, Song et al. [36]
developed a network model to describe the fractality of networks. The main
characteristic of this model is the introduction of a parameter $e$ which
could be used to control the original hubs whether continue to form
connections between the nodes in the process of the growth of complex
networks. The authors [36] pointed out that the parameter $e$ can be regarded
as a level of fractality of the network. The network corresponds to a pure
fractal network which is a pure fractal set (defined by Mandelbrot [7]) when
$e=0$, and a pure small world network when $e=1$ [2]. Later on, Gallos et al.
[35] proposed a generalized version of this network model.

In Section 2, we introduce the generalized version of the network model in
Ref. [35] and some of its topological properties. In Section 3, we examine the
fractal dimension of these networks . A new fixed-size box-counting algorithm
for MFA of networks modified from the one proposed in Ref. [34] is given in
Section 4. The multifractal properties of the model networks and their results
are also given in this Section. Some conclusions are presented in Section 5.

##  2 Network model

A graph (or network) is a collection of nodes which denote the elements of a
system, and links or edges which identify the relations or interactions among
these elements. In this section, the algorithm of the generalized version of
the network model in Ref. [35] is presented. The network could be obtained by
a method described as follows.

First we give a real number $0\leq e\leq 1$, and two positive integers $m$ and
$x$ ($x\leq m$). In the generation $n=0$, we start with only two nodes and one
edge between them. In order to get the network of the generation $n+1$, every
endpoint of each edge $L$ in the network of the generation $n$ is attached to
$m$ new nodes. Then we generate a random number $p$ which obeys the uniform
distribution between 0 and 1. If $0\leq p<e$, each edge $L$ of the generation
$n$ is kept and $x-1$ new edges are appended to connect pairs of the new nodes
attached to the endpoints of $L$; otherwise, for each edge $L$ of the
generation $n$, we add $x$ new edges matching new nodes at the ends of $L$ and
remove $L$ (see Fig. 1). As shown in Fig. 2, if we take $e=0,m=2,x=2$, we add
$m=2$ new nodes to the two endpoints of the sole edge in the generation $n=0$
to get the network of the generation $n=1$. Due to $e=0$, we take away the
edge of $n=0$ and add $x=2$ edges between the new nodes. Notice that when
$x=1$, we get tree structure without loop for any value of $e$.

Figure 1: Construction of network. The link between hub remains with
probability $e$, otherwise, it is replaced by another link between new nodes
with probability $1-e$.

Figure 2: Construction of a pure fractal network. Example of the network
modelof generations $n=0,1,2$ with parameters $m=2,x=2,e=0$.

According to above description, if we denote $M_{n}$ the number of edges in
the network of the generation $n$, we can have $M_{n+1}=(2m+x)*M_{n}$ in the
generation $n+1$. Hence ${{M_{n}}}={(2m+x)^{n}}$. Meanwhile in the growth of
the network from the generation $n-1$ to the generation $n$, each edge in the
network of the generation $n-1$ produces $2m$ new nodes. Hence we have
$N_{n}=2mM_{n-1}+N_{n-1}$, where $N_{n}$ is the number of nodes in the network
of the generation $n$. Therefore,

${N_{n}}=\frac{2m}{2m+x-1}(2m+x)^{n}+2-\frac{2m}{2m+x-1}$ (1)

It was proved in Song et al. [36] that the degree distribution $P(k)$ of the
network model satisfies a power law relationship $P(k)\approx{k^{-\gamma}}$
with $\gamma=1+\ln b/\ln s$, where $b$ is the scaling of the node number and
$s$ is the scaling of the node degrees between two adjacent generations in the
process of the growth of the network. From the algorithm described above, we
know that if the degree of a node in the network of the generation $n$ is
$K_{n}$, then it should be $mK_{n}+K_{n}$ with probability $e$, or $mK_{n}$
with probability $1-e$ in generation $n+1$, therefore
$K_{n+1}=e(mK_{n}+K_{n})+(1-e)mK_{n}=(m+e)K_{n}$. It is easy to know that
$b=2m+x$ from Eq. (1). So for the above model, we have

$\gamma=1+\frac{{\ln(2m+x)}}{{\ln(m+e)}}$ (2)

In addition, a simple analysis shows that the clustering coefficient of the
network model is $0$ for any value of $e$.

##  3 The fractal dimension

We find that if the distance between two nodes of the generation $n$ is
$L_{n}$, then in the network of the generation $n+1$, the distance $L_{n+1}$
would be $L_{n}$ with probability $e$, or $3L_{n}$ with probability $1-e$.
Hence $L_{n+1}=eL_{n}+3(1-e)L_{n}=(3-2e)L_{n}$. From Ref. [36], we know that
the theoretical fractal dimension of the model networks is ${d_{f}^{T}}=\ln
b/\ln a$, where $a=L_{n+1}/L_{n}$. Therefore,

${d_{f}^{T}}=\frac{{\ln(2m+x)}}{{\ln(3-2e)}}$ (3)

If $x=2,m=2$, we have $d_{f}^{T}=\ln(6)/\ln(3-2e)$.

We can also numerically calculate the fractal dimension of the model networks
using some algorithms (e.g. [10, 37]). Here we adopt the random sequential
box-counting algorithm proposed by Kim et al. [37] to estimate the fractal
dimension of networks (two examples for estimating fractal dimension are shown
in Fig. 3). We denote $d_{f}^{N}$ the fractal dimension of the network
obtained numerically. First we want to check whether values of $d_{f}^{N}$
coincide with the theoretical values of fractal dimension $d_{f}^{T}$. If the
numerical and theoretical fractal dimensions coincide with each other, we will
have confidence on our process and program to estimate the multifractal curves
$\tau(q)$ and $D(q)$ of these networks. Due to the limit of computational
capacity of our computers, we only generate the networks up to the generation
$n=5$.

For each value of $e$, we generate 100 networks (100 realizations) and
calculate the average value of $d_{f}^{N}$ over the 100 realizations. The $e$
vs $<d_{f}^{N}>$ plot is presented in Fig. 4. From Fig. 4, we can see that the
numerical $<d_{f}^{N}>$ coincides with the theoretical $<d_{f}^{T}>$
perfectly.

Figure 3:  Two examples to estimate the fractal dimension of networks for
$e=0.1$ and $0.5$, here parameters $n=5,m=2,x=2$. We can see the estimated
fractal dimension is very close to the theoretical result 2.5850 (for $e=0.5$)
and 1.7402 (for $e=0.1$) respectively.

Figure 4: Numerical result of the relationship between the fractal dimension
$<d_{f}^{N}>$ of the networks and $e$ with parameters $n=5,m=2,x=2$. Here
$<\cdot>$ means the average over 100 realizations.

##  4 Multifractal analysis

In this section, we first introduce a new algorithm for MFA of networks
modified from the one proposed in Ref. [34] and then apply it to the model
networks presented in Section 2.

Two networks which have the same fractal dimension may look completely
different. In addition, when the networks have rich scale and self-similar
structures, they exhibit different dimensions in different scales. MFA is a
powerful method to study the networks with such characteristics.

At present, the fixed-size box-counting algorithm is the most common algorithm
for MFA [12]. For a given probability measure $0\leq\mu\leq 1$ with support
set $E$ in metric space, we consider the partition function

${Z_{r}}(q)=\sum\limits_{\mu(B)\neq 0}{{{[\mu(B)]}^{q}}}$ (4)

where $q\in R$ , the result is the sum of all the different non-overlapping
boxes $B$ with a given size $r$ in the covering of the support set $E$. It is
easy to know that ${Z_{r}}(q)\geq 0$ and ${Z_{r}}(0)=1$. We define the mass
exponent function $\tau(q)$ of the measure $\mu$ as

$\tau(q)=\mathop{\lim}\limits_{r\to 0}\frac{{\ln{Z_{r}}(q)}}{{\ln r}}$ (5)

Then we get the generalized fractal dimensions of the measure $\mu$ by

${D(q)}=\frac{{\tau(q)}}{{q-1}},q\neq 1$ (6)

and

${D(1)}=\mathop{\lim}\limits_{r\to 0}\frac{{{Z_{(1,r)}}}}{{\ln r}},q=1$ (7)

where${Z_{(1,r)}}=\sum\limits_{\mu(B)\neq 0}{\mu(B)\ln\mu(B)}$.

In practice, the generalized fractal dimensions are usually obtained by linear
regression. Specifically, $D(0)$ is the fractal dimension of the support set
of the measure $\mu$, ${D(1)}$ and ${D(2)}$ are called the information
dimension and the correlation dimension respectively.

For a network, the measure $\mu$ of each box can be defined as the ratio of
the number of nodes covered by the box to the total number of nodes in the
network [34, 37]. We need to complete the following two steps before we
proceed MFA.

i) Map a network to an adjacent matrix ${A_{N\times N}}$, where $N$ is the
total number of nodes in the network. It is easy to know that ${A_{N\times
N}}$ is a symmetric matrix where the elements ${a_{ij}}=1$ when there is an
edge between the nodes $i$ and $j$, otherwise ${a_{ij}}=0$. Here, the edge
from node $i$ to node $i$ is not considered, so ${a_{ii}}=0$.

ii) Use ${A_{N\times N}}$ to calculate the shortest distance between any two
nodes in the network and store them into another matrix ${B_{N\times N}}$.
Here, in our study, we use Dijkstra s algorithm of MatLab toolbox to calculate
the shortest distance between two nodes in the network.

After finishing the two steps presented above, we can use matrix ${B_{N\times
N}}$ as the input of MFA of the network model described in Section 2 based on
a new algorithm for MFA of networks modified from the one proposed in Ref.
[34] as follows.

(I) Ensure that all nodes in the network are not covered, and no node has been
selected as the center of a box.

(II) According to the size of our networks $N=6222$ (with the parameters
$n=5;m=2,x=2$), we set $t=1,2,\ldots,T$. Here we take $T=1000$, then we
rearrange the nodes number into $T=1000$ different random orders. That is to
ensure that the nodes of a network are randomly chosen as center nodes.

(III) Set the radius $r$ of boxes which will be used to cover the nodes in the
range $[1,d]$, where $d$ is the diameter of the network (i.e. the longest
distance between nodes in the network).

(IV) Treat the nodes of the $t$th kind of random orders that we have got in
(II) as the center of a box successively, then search all the other nodes. If
a node has a distance to the center node within $r$ and has not been covered
yet, then cover it.

(V) If no more new nodes can be covered by this box, then we abandon this box.

(VI) Repeat (IV) - (VI) until all the nodes are covered by the corresponding
boxes. We denote the number of boxes in this box covering as $N(t,r)$.

(VII) Repeat steps (III) and (VI) for all the random orders to find a box
covering with minimal number of boxes $N(t,r)$.

(VIII) For each nonempty box $B$ in the first box covering with minimal number
of boxes, we define its measure as ${\mu(B)={N_{B}}/6222}$, where ${N_{B}}$ is
the number of nodes covered by the box $B$. For each $r$, we calculate the
partition sum $Z_{r}(q)=\sum\limits_{\mu(B)\neq 0}{{{[\mu(B)]}^{q}}}$.

(IX) For different $r$, we repeat (III)-(VIII). Then we use $Z_{r}(q)$ for
linear regression.

Remark 1: In the algorithm of MFA of networks proposed in Ref. [34], we use
$\overline{Z}_{r}(q)$ (the average of $Z_{r}(q)$ for all $T=1000$ different
random orders of the nodes) for linear regression to get $\tau(q)$ (hence
$D(q)$). But when $q=0$, $D(0)$ got in this way is not the box-counting
dimension of the network because there requires minimum number of boxes which
cover the fractal set (network here) [9]. Here we modify to use $Z_{r}(q)$ of
a covering with minimum number of boxes for linear regression to get $\tau(q)$
(hence $D(q)$). So when $q=0$, $D(0)$ is exact the box-count dimension of the
network. It is a more reasonable extension from the traditional MFA.

In order to get the range $r\in[{r_{\min,}}{r_{\max}}]$ in which the networks
obey the power law and then to get the mass exponents $\tau(q)$ and the
generalized fractal dimensions $D_{q}$, linear regression is an important
step. In our calculation, we run the linear regression of $\ln{Z_{r}}(q)$
against $\ln r$ to get $\tau(q)$, and then get $D(q)$ through formula
$D(q)=\tau(q)/(q-1)$ for $q\neq 1$ and $D(1)$ through the linear regression of
$Z_{(1,r)}=\sum\limits_{\mu(B)\neq 0}{\mu(B)\ln\mu(B)}$ against $\ln r$ for
$q=1$.

By applying the new fixed-size box-counting algorithm described above on the
model networks, we get the following results:

First, for each value of $e$ (here we take $e=0.1,0.2,...,0.8$), we generate
100 networks (we take 100 realizations because the MFA for networks is very
time consuming when the network is large), and calculate the $\tau(q)$ and
$D(q)$ curves for each network using the new fixed-size box-counting
algorithm. Then we take average for these $\tau(q)$ and $D(q)$ curves over the
100 realizations. The shape of the $<\tau(q)>$ curves shown in Fig. 5 and the
$<D_{q}>$ curves shown in Fig. 6 are all nonlinear, which indicate that all
the networks we studied have multifractal property. We also find that the
value of $\Delta(<D(q)>)$ defined by $\max(<D(q)>)-\min(<D(q)>)$ increases
with the increase of the parameter $e$, which indicates that the multifractal
property of the model networks becomes more obvious when the value of the
parameter $e$ becomes larger.

The multifractal property of the model networks revealed by our work indicates
that the model networks are very complicated which cannot be characterized by
a single fractal dimension. The MFA algorithm proposed here can be used to
provide a more accurate characterization for the model networks, even for some
other complicated networks.

Second, we find that the average information dimension $<D(1)>$ has a linear
relation with the parameter $e$, i.e. $<D(1)>=1.5053*e+1.4735$ as shown in
Fig. 7, which is different from that of $D(0)$ shown in Eq. (3).

Figure 5: The $\tau(q)$ curves of the network model, here $<\cdot>$ means the
average over 100 realizations.

Figure 6: The $D(q)$ curves of the network model, here $<\cdot>$ means the
average over 100 realizations.

Figure 7: The relationship of $<D(1)>$ against parameter $e$ respectively,
here $<\cdot>$ means the average over 100 realizations.

Remark 2: Furuya and Yakubo [33] also proposed an algorithm for MFA of complex
networks. The difference between the algorithm in Ref. [33] and our algorithm
is the definition of the measure $\mu$. In the algorithm in Ref. [33], it
allows that any two boxes in the box covering have overlap and defines the
measure $\mu_{i}$ by counting the times of overlaps of each node, hence it is
not a natural extension of the traditional MFA (see Eq. (4)). In our
algorithm, overlap of any two boxes in the box covering is not allowed, so it
is a natural extension of the traditional MFA. Our network model with $e>0$ is
different from the ($u,v$)-flower network model. Only the network model with
$e=0$ corresponds to the deterministic ($u,v$)-flower network model with
$u=v=3$. Furuya and Yakubo [33] also gave a theoretical formula for the
$\tau(q)$ function of ($u,v$)-flower network model (Eq. (11) of Ref. [33]).
When $u=v=3$, $\tau(q)$ has the formula [33]: $\tau(q)=q$ if
$q\geq\ln(6)/\ln(2)=2.5850$, and $\tau(q)=(q-1)\frac{\ln(6)}{\ln(3)}$ if
$q<\ln(6)/\ln(2)=2.5850$. We compared this formula with our numerical result
for $e=0$ in Fig.5, and found that we also have $\tau(q)\simeq q$ if $q\geq
2.5850$, but the $\tau(q)$ values are different from
$(q-1)\frac{\ln(6)}{\ln(3)}$ if $q<2.5850$.

##  5 Conclusion

We have studied the fractal and multifractal properties of a family of model
networks that were originally proposed to explain the origin of fractality in
complex networks. This model introduces a parameter $e$, which can be used to
tune the fractality level of the network. One can get a pure fractal network
when $e=0$ and obtain a small-world network when $e=1$. We investigated the
fractal and multifractal properties through numerical calculation. To make the
calculation feasible and accurate, we calculated the model with parameters
${n=5;m=2,x=2}$; ${e=0}$, ${0.1,0.2,\ldots,0.8}$; ${q=-10,\ldots,10}$. The
result of the $\tau(q),D(q)$ (including $D(0)$ and $D(1)$) are averaged over
$100$ realizations (networks). The shape of $\tau(q)$ and $D(q)$ curves are
all nonlinear, which indicates that all the networks we studied have
multifractal property. We also found that the value of
$\Delta(<D(q)>)=\max(<D(q)>)-\min(<D(q)>)$ increases with the increase of the
parameter $e$, which indicates that the multifractal property of the model
becomes more obvious when the value of the parameter $e$ becomes larger.

We also found that the average information dimension $<D(1)>$ has a linear
relation with the parameter $e$, i.e. $<D(1)>=1.5053*e+1.4735$.

The MFA algorithm proposed here can be used to provide a more accurate
characterization for the model networks, even for some other complicated
networks.

## Acknowledgments

This project was supported by the Natural Science Foundation of China (Grant
Nos. 11071282 and 11371016), the Chinese Program for Changjiang Scholars and
Innovative Research Team in University (PCSIRT) (Grant No. IRT1179); the
Research Foundation of Education Commission of Hunan Province of China (Grant
No. 11A122); the Lotus Scholars Program of Hunan province of China. The
authors would like to thank the editor and the reviewers for their insights,
comments and suggestions to improve this paper.

## References

  * [1]
  * [1] Song C, Havlin S and Makse H A, Self-similarity of complex networks, Nature 433 (2005) 392-395.
  * [2] Rozenfeld H D and Makse H A, Fractality and the percolation transition in complex networks, Chem. Eng. Sci. 64 (2009) 4572-4575.
  * [3] Rozenfeld H D, Song C and Makse H A, Small-World to Fractal Transition in Complex Networks: A Renormalization Group Approach, Phys. Rev. Lett. 104 (2010) 025701.
  * [4] Zhang Z., Zhou S., Zou T. and Chen G.,Fractal scale-free networks resistant to disease spread,  J. Stat. Mech.: Theor. Exp. 9 (2008) P09008.
  * [5] Zhang Z, Zhou S, Xie W, Chen L, Lin Y, and Guan J, Standard random walks and trapping on the Koch network with scale-free behavior and small-world effect, Phys. Rev. E 79 (2009) 061113.
  * [6] Liu J X and Kong X M, Establishment and structure properties of the scale-free Koch network (in Chinese), Acta Phy. Sin. 59 (2010) 2244-2249.
  * [7] Mandelbrot B B, 1983 The Fractal Geometry of Nature (New York: Academic Press).
  * [8] Feder J, 1988 Fractals (New York: Plenum).
  * [9] Falconer K, 1997 Techniques in Fractal Geometry (New York: Wiley).
  * [10] Song C, Gallos L K, Havlin S and Makse H A, How to calculate the fractal dimension of a complex network: the box covering algorithm,  J. Stat. Mech.: Theor. Exp. 3 (2007) P03006.
  * [11] Grassberger P and Procaccia I, Characterization of Strange Attractors, Phys. Rev. Lett 50 (1983) 346-349.
  * [12] Halsey T C, Jensen M H, Kadanoff L P, Procaccia I and Shraiman B I, Fractal measures and their singularities: The characterization of strange sets, Phys. Rev. A 33 (1986) 1141-1151.
  * [13] Canessa E, Multifractality in time series, J. Phys. A 33 (2000) 3637-3651.
  * [14] Anh V V, Tieng Q M and Tse Y K, Cointegration of stochastic multifractals with application to foreign exchange rates, Int. Trans. Oper. Res. 7 (2000) 349-363.
  * [15] Yu Z G, Anh V and Lau K S, Multifractal characterisation of length sequences of coding and noncoding segments in a complete genome, Physica A 31 (2001) 351-361.
  * [16] Yu Z G, Anh V and Lau K S, Measure representation and multifractal analysis of complete genomes, Phys. Rev. E 64 (2001) 031903\.
  * [17] Yu Z G, Anh V and Lau K S, Multifractal and correlation analyses of protein sequences from complete genomes, Phys. Rev. E 68 (2003) 021913\.
  * [18] Yu Z G, Anh V and Lau K S, Chaos game representation of protein sequences based on the detailed HP model and their multifractal and correlation analyses, J. Theor. Biol. 226 (2004) 341-348.
  * [19] Yu Z G, Anh V V, Lau K S and Zhou L Q, Clustering of protein structures using hydrophobic free energy and solvent accessibility of proteins, Phys. Rev. E 73 (2006) 031920.
  * [20] Yu Z G, Xiao Q J, Shi L, Yu J W, and Anh V, Chaos game representation of functional protein sequences, and simulation and multifractal analysis of induced measures, Chin. Phys. B 19 (2010) 068701.
  * [21] Anh V V, Lau K S and Yu Z G, Recognition of an organism from fragments of its complete genome, Phys. Rev. E 66 (2002) 031910\.
  * [22] Zhou L Q,Yu Z G, Deng J Q, Anh V and Long S C, A fractal method to distinguish coding and non-coding sequences in a complete genome based on a number sequence representation, J. Theor. Biol. 232 (2005) 559-567.
  * [23] Han J J and Fu W J, Wavelet-based multifractal analysis of DNA sequences by using chaos-game representation, Chin. Phys. B 19 (2010) 010205.
  * [24] Zhu S M, Yu Z G, Ahn V, Protein structural classification and family identification by multifractal analysis and wavelet spectrum, Chin. Phys. B 20 (2011) 010505\.
  * [25] Kantelhardt J W, Koscielny-Bunde E, Rybski D, Braun P, Bunde A and Havlin S, Long-term persistence and multifractality of precipitation and river runoff records, J. Geophys. Res. 111 (2006) D01106.
  * [26] Veneziano D, Langousis A and Furcolo P, Multifractality and rainfall extremes: A review, Water Resour. Res. 42 (2006) W06D15.
  * [27] Venugopal V, Roux S G, Foufoula-Georgiou E and Arneodo A, Revisiting multifractality of high-resolution temporal rainfall using a wavelet-based formalism, Water Resour. Res. 42 (2006) W06D14.
  * [28] Yu Z G, Anh V V, Wanliss J A, and Watson S M, Chaos game representation of the Dst index and prediction of geomagnetic storm events, Chaos, Solitons and Fractals 31 (2007) 736-746.
  * [29] Yu Z G, Anh V and Eastes R, 2009 J. Geophys. Res. 114 A05214
  * [30] Yu Z G, Anh V, Wang Y, Mao D and Wanliss J, Multifractal analysis of geomagnetic storm and solar flare indices and their class dependence, J. Geophys. Res. 115 (2010) A10219
  * [31] Zang B J and Shang P J, Multifractal analysis of the Yellow River flow, Chin. Phys. 16 (2007) 565-569.
  * [32] Lee C Y and Jung S H, Statistical self-similar properties of complex networks, Phys. Rev. E 73 (2006) 066102\.
  * [33] Furuya S and Yakubo K, Multifractality of complex networks, Phys. Rev. E 84 (2011) 036118
  * [34] Wang D L, Yu Z G and Anh V, Multifractal analysis of complex networks, Chin. Phys. B 21 (2012) 080504\.
  * [35] Gallos L K, Song C, Havlin S and Makse H A, Proc. Natl. Acad. Sci. U.S.A. 104 (2007) 7746-7751.
  * [36] Song C, Havlin S and Makse H A, Origins of fractality in the growth of complex networks, Nat. Phys. 2 (2006) 275-281.
  * [37] Kim J S, Goh K I, Salvi G, Oh E, Kahng B and Kim D, Fractality in complex networks: Critical and supercritical skeletons, Phys. Rev. E 75 (2007) 016110\.
  * [2]


# Direct and alternating magnon spin currents across a junction interface
irradiated by linearly polarized laser

Kouki Nakata  Advanced Science Research Center, Japan Atomic Energy Agency,
Tokai, Ibaraki 319-1195, Japan  Yuichi Ohnuma  Research Center for Advanced
Science and Technology (RCAST), The University of Tokyo, Meguro, Tokyo
153-8904, Japan

###### Abstract

The developments in the field of quantum optics raise expectations that laser-
matter coupling is a promising building block for magnonics. Here, we propose
a method for the generation of direct and alternating spin currents of magnons
across the junction interface irradiated by linearly polarized laser. In a
junction of ferromagnetic insulators with a large electronic gap, the spin
angular momentum is exchanged during the tunneling process of magnons across
the junction interface. The advanced technology in the field of plasmonics and
metamaterials realizes that spins irradiated by the laser field interact only
with the magnetic component of the laser through the Zeeman coupling. Using an
analytic perturbation theory, we provide a general formula for magnon
transport induced by the inversion symmetry breaking across the junction
interface. Then, we show that those spin currents are enhanced by the
ferromagnetic resonance, and the period of the ac spin current is one-half of
that of the laser magnetic field. Finally, we estimate the magnitude of the
spin current, and find that it will be within experimental reach.

##  I Introduction

Toward efficient transmission of information that goes beyond what is offered
by conventional electronics, the last decade has seen a rapid development of
magnon-based spintronics, dubbed magnonics Chumak  _et al._ (2015); Nakata
_et al._ (2017a), aiming at utilizing the quantized spin waves, magnons, as a
carrier of information in units of the Bohr magneton $\mu_{\text{B}}$. A most
promising strategy for this goal is to use magnetic insulators, which are free
from drawbacks of conventional electronics such as significant energy loss due
to Joule heating. Thus, inventing methods to handle a flexible manipulation of
magnon transport in magnetic insulators, in the complete absence of any
conducting metallic elements, is a central task in the field of magnonics.

The recent developments in the field of plasmonics Ciappina  _et al._ (2017)
and metamaterials Mukai  _et al._ (2014) raise expectations that laser-matter
coupling Mukai  _et al._ (2014); Bossini  _et al._ (2016); Ciappina  _et al._
(2017); Arikawa  _et al._ (2017) is an important building block for this holy
grail of magnonics. Using the optical method Kirilyuk  _et al._ (2010); Kimel
_et al._ (2004, 2005); Hansteen  _et al._ (2005); Stanciu  _et al._ (2007a);
N$\check{\rm{e}}$mec  _et al._ (2018), the reversal of magnetization was
achieved experimentally Stanciu  _et al._ (2007b); Radu  _et al._ (2011);
Ostler  _et al._ (2012); Vahaplar  _et al._ (2009). An optical analog of the
conventional Barnett effect Barnett (1915, 1935); Matsuo  _et al._ (2017),
dubbed the optical Barnett effect Rebei and Hohlfeld (2008a, b) (i.e., laser-
induced magnetization Takayoshi  _et al._ (2014a, b)) 111 See Refs. Takayoshi
_et al._ (2014a, b) for the importance of modulating laser frequency
adiabatically by the chirping technique Sato  _et al._ (2013); Kamada  _et
al._ (2013). See also Ref. Nakata and Takayoshi (2020) for the difference
between the optical Barnett effect and the inverse Faraday effect Kirilyuk
_et al._ (2010); Kimel  _et al._ (2005)., and the resulting magnon Josephson
effect Nakata (2021) were proposed theoretically by means of circularly
polarized laser.

In this paper, using linearly polarized laser, we propose a method for the
generation of direct (dc) and alternating (ac) spin currents of magnons across
a junction interface due to the inversion symmetry breaking. We consider a
junction formed by ferromagnetic insulators with a large electronic gap, where
the spin angular momentum is exchanged during the tunneling process of magnons
across the junction interface. In addition to the large electronic gap, owing
to the state-of-the-art technology in the field of plasmonics and
metamaterials that realizes strong magnetic field pulses with small electric
field Ciappina  _et al._ (2017); Mukai  _et al._ (2014), spins in the laser
field interact only with the magnetic component of the laser through the
Zeeman coupling. By means of an analytic perturbation theory, we provide a
general formula for magnon transport induced by the inversion symmetry
breaking across the junction interface. Then, we find that the ferromagnetic
resonance enhances those spin currents, and the period of the ac spin current
is one-half of that of the laser magnetic field. Finally, we estimate the
magnitude of the spin current, and show that it will be within experimental
reach.

This paper is organized as follows. We explain the setup of our study in Sec.
II, and investigate magnon transport across the junction interface in Sec.
III. Then, estimating the values of the magnon spin current, we propose a
method for observation, and provide a few remarks in Sec. IV. Finally, we
summarize in Sec. V. Technical details are described in the Appendix.

Figure 1: Schematic picture of the junction. The junction consists of three-
dimensional ferromagnetic insulators where localized spins are polarized along
the $z$ direction. The two ferromagnetic insulators are separated by a thin
film of a nonmagnetic insulator and the boundary spins in the left and the
right ferromagnetic insulators are weakly exchange coupled, which results in
the tunneling process of magnons with the amplitude $J_{\text{ex}}$ across the
junction interface. The linearly polarized laser, i.e., time-periodic
transverse magnetic field, with the field strength of being spatially uniform
is applied to the two ferromagnetic insulators.

##  II System

We consider a junction of three-dimensional ferromagnetic insulators with a
large electronic gap. The two ferromagnetic insulators are aligned along the
$x$ direction and the localized spins are polarized along the $z$ direction,
see Fig. 1 Nakata  _et al._ (2015a). A finite overlap of the wave function of
the boundary spins in the left and the right ferromagnetic insulators results
in an exchange interaction, which induces a tunneling process of magnons
across the junction interface. During the tunneling process, the spin angular
momentum is exchanged and it may be described in general by the Hamiltonian
Nakata  _et al._ (2015a, 2018)
${\mathcal{H}}_{\text{ex}}=-J_{\text{ex}}\sum_{{\mathbf{k}}_{\perp}}\sum_{k_{x},k_{x}^{\prime}}a_{{\text{L}},{\mathbf{k}}}a^{\dagger}_{{\text{R}},{\mathbf{k}}^{\prime}}+{\text{H.
c.}}$, where the bosonic operator $a^{(\dagger)}_{\text{L/R}}$ annihilates
(creates) magnons at the boundary of the left/right insulator with the wave
numbers ${\mathbf{k}}\coloneqq(k_{x},k_{y},k_{z})$,
${\mathbf{k}}^{\prime}\coloneqq(k_{x}^{\prime},k_{y},k_{z})$, and
${\mathbf{k}}_{\perp}\coloneqq(0,k_{y},k_{z})$, and where $J_{\text{ex}}>0$
represents the tunneling amplitude of magnons across the junction interface
and the value depends on the spin quantum number of the boundary spins in each
insulator. The ferromagnetic insulators are separated by a thin film of a
nonmagnetic insulator, and the boundary spins in the left and the right
ferromagnetic insulators are weakly exchange coupled in that the magnitude of
$J_{\text{ex}}$ is small compared with the spin exchange interaction between
the nearest-neighbor spins in each insulator. See Refs. Nakata  _et al._
(2015a, 2018) for the details of the relation between the magnon Hamiltonian
${\mathcal{H}}_{\text{ex}}$ and the microscopic spin model through a Holstein-
Primakoff expansion to leading order Holstein and Primakoff (1940).

Then, we apply the linearly polarized laser of the frequency $\Omega$ to the
junction of the two ferromagnetic insulators with a large electronic gap. Due
to the large electronic gap, spins in the laser interact only with the
magnetic component of the laser field through the Zeeman coupling. Since we
assume that the magnetic field amplitude is spatially uniform, the laser field
(i.e., time-periodic transverse magnetic field) excites the zero mode of
magnons in each insulator 222Note that in Ref. Nakata  _et al._ (2015b), the
laser field is applied only to one side. and the process in the left (right)
insulator is described Nakata  _et al._ (2015b) by the time-periodic
Hamiltonian $V_{\text{L(R)}}(t)$:

$\displaystyle V_{\text{L}}(t)$
$\displaystyle=\Gamma_{\text{L}}(t)(a_{{\text{L}},{\mathbf{k}}=0}+a^{\dagger}_{{\text{L}},{\mathbf{k}}=0}),$
(1a) $\displaystyle V_{\text{R}}(t)$
$\displaystyle=\Gamma_{\text{R}}(t)(a_{{\text{R}},{\mathbf{k}}=0}+a^{\dagger}_{{\text{R}},{\mathbf{k}}=0}),$
(1b)

with

$\displaystyle\Gamma_{\text{L}}(t)$
$\displaystyle=2\Gamma_{0({\text{L}})}{{\text{cos}}(\Omega t)},$ (2a)
$\displaystyle\Gamma_{\text{R}}(t)$
$\displaystyle=2\Gamma_{0({\text{R}})}{{\text{cos}}(\Omega t)},$ (2b)

where $2\Gamma_{0({\text{L/R}})}$ represents the magnitude of the interaction
between the laser magnetic field and magnons in the left/right insulator.
Since the magnitude depends on the spin quantum number, it is
$\Gamma_{0(\text{L})}\neq\Gamma_{0(\text{R})}$ in general. For simplicity,
assuming that the $g$-factor of spins in the left insulator is identical to
that of the right one, we denote it as $g$. Since we apply a weak laser field
in that the magnitude of $2\Gamma_{0({\text{L/R}})}$ is small compared with
the spin exchange interaction between the nearest-neighbor spins in each
insulator, the description by means of magnon degrees of freedom through the
Holstein-Primakoff expansion to leading order [Eqs. (1a) and (1b)] remains
valid and the laser magnetic field $V_{\text{L(R)}}(t)$ can be treated as a
perturbation term.

We remark that since the laser excites the zero mode of magnons, the result
does not depend on the details of the energy dispersion relation, e.g., the
spin stiffness constant, except the magnon energy gap $\Delta_{\text{L(R)}}$
for the zero mode in the left (right) insulator. The magnon energy gap is
formed by, e.g., a spin anisotropy Nakata  _et al._ (2017b). Later, for
convenience, we phenomenologically introduce the lifetime Nakata and Ohnuma
(2021) of magnons $\tau_{\text{L(R)}}$ in the left (right) insulator through
the Green’s function Tatara  _et al._ (2008).

##  III Magnon transport in junction

###  III.1 Spin currents of magnons

Through the Heisenberg equation of motion, the expectation value of the magnon
spin current across the junction interface $\langle I(t)\rangle$, from the
right to the left insulator, is given as See

$\langle
I(t)\rangle=2g\mu_{\text{B}}\frac{J_{\text{ex}}}{\hbar}{\text{Im}}\langle{\mathcal{A}}(t)\rangle+O({J_{\text{ex}}}^{2}),$
(3)

where $\hbar$ denotes the reduced Planck constant and

${\mathcal{A}}(t)\coloneqq\sum_{{\mathbf{k}},k_{x}^{\prime}}a_{{\text{L}},{\mathbf{k}}}(t)a^{\dagger}_{{\text{R}},{\mathbf{k}}^{\prime}}(t).$
(4)

Thus, first, we calculate
$\langle{\mathcal{A}}(t)\rangle=O({J_{\text{ex}}}^{0}{\Gamma_{0({\text{L/R}})}}^{2})$
using a perturbation theory. Then, taking the imaginary component of
$\langle{\mathcal{A}}(t)\rangle$, we evaluate the magnon spin current $\langle
I(t)\rangle=O(J_{\text{ex}}{\Gamma_{0({\text{L/R}})}}^{2})$ in the junction
subjected to the laser magnetic field.

A straightforward perturbative calculation in $\Gamma_{0({\text{L/R}})}$ based
on the Schwinger-Keldysh formalism Schwinger (1961); Martin and Schwinger
(1959); Keldysh (1965), and up to ${\Gamma_{0({\text{L/R}})}}^{2}$, provides
See

$\displaystyle\langle{\mathcal{A}}(t)\rangle=\frac{1}{\hbar^{2}}$
$\displaystyle\int_{{\text{c}}}d\tau_{1}\int_{{\text{c}}}d\tau_{2}\Gamma_{\text{L}}(\tau_{1})\Gamma_{\text{R}}(\tau_{2})$
$\displaystyle\cdot
G_{{\text{L}},{\mathbf{k}}=0}(\tau,\tau_{1})G_{{\text{R}},{\mathbf{k}}=0}(\tau_{2},\tau)|_{\tau\rightarrow
t},$ (5)

where $G_{{\text{L(R)}}}$ is the contour-ordered Green’s function Haug and
Jauho (2007); Kita (2010); Tatara  _et al._ (2008) of magnons in the left
(right) insulator, and $\tau$ and $\tau_{1(2)}$ are the contour variables
defined on the Schwinger-Keldysh closed time path Rammer (2007) c,
${\text{c}}\coloneqq{\text{c}}_{\rightarrow}+{\text{c}}_{\leftarrow}$,
consisting of the forward path ${\text{c}}_{\rightarrow}$ and the backward
path ${\text{c}}_{\leftarrow}$. Here, for convenience, we take $\tau$ on the
forward path. Even when $\tau$ is located on the backward path, the result
remains unchanged. Using the Langreth theorem Langreth and Wilkins (1972), it
becomes

$\displaystyle\langle{\mathcal{A}}(t)\rangle=\frac{1}{\hbar^{2}}$
$\displaystyle\int_{-\infty}^{\infty}dt_{1}\int_{-\infty}^{\infty}dt_{2}\Gamma_{\text{L}}(t_{1})\Gamma_{\text{R}}(t_{2})$
$\displaystyle\cdot[G^{\text{t}}_{{\text{L}},{\mathbf{k}}=0}(t,t_{1})G^{\text{t}}_{{\text{R}},{\mathbf{k}}=0}(t_{2},t)$
$\displaystyle-G^{\text{t}}_{{\text{L}},{\mathbf{k}}=0}(t,t_{1})G^{\text{>}}_{{\text{R}},{\mathbf{k}}=0}(t_{2},t)$
$\displaystyle-G^{\text{<}}_{{\text{L}},{\mathbf{k}}=0}(t,t_{1})G^{\text{t}}_{{\text{R}},{\mathbf{k}}=0}(t_{2},t)$
$\displaystyle+G^{\text{<}}_{{\text{L}},{\mathbf{k}}=0}(t,t_{1})G^{\text{>}}_{{\text{R}},{\mathbf{k}}=0}(t_{2},t)],$
(6)

where $G^{\text{t}}$, $G^{<}$, and $G^{>}$ are the time-ordered, lesser, and
greater Green’s functions, respectively. Through the identities,
$G^{\text{t}}_{{\mathbf{k}}}(t,t_{1})=G^{\text{r}}_{{\mathbf{k}}}(t,t_{1})+G^{<}_{{\mathbf{k}}}(t,t_{1})$
and
$G^{<}_{{\mathbf{k}}}(t,t_{1})-G^{>}_{{\mathbf{k}}}(t,t_{1})=G^{\text{a}}_{{\mathbf{k}}}(t,t_{1})-G^{\text{r}}_{{\mathbf{k}}}(t,t_{1})$,
it reduces to
$\langle{\mathcal{A}}(t)\rangle=(1/\hbar^{2})\int_{-\infty}^{\infty}dt_{1}\int_{-\infty}^{\infty}dt_{2}\Gamma_{\text{L}}(t_{1})\Gamma_{\text{R}}(t_{2})G^{\text{r}}_{{\text{L}}}(t,t_{1})G^{\text{a}}_{{\text{R}}}(t_{2},t)$,
where $G^{\text{r(a)}}$ is the retarded (advanced) Green’s function. After the
Fourier transform of the Green’s functions, we perform the time integral and
the resulting Dirac delta function for the frequency. Finally, we obtain

$\displaystyle\langle{\mathcal{A}}(t)\rangle=\Gamma_{0(\text{L})}\Gamma_{0(\text{R})}$
$\displaystyle[(G^{\text{r}}_{{\text{L}},{\mathbf{k}}=0,-\Omega}G^{\text{a}}_{{\text{R}},{\mathbf{k}}=0,-\Omega}$
$\displaystyle+G^{\text{r}}_{{\text{L}},{\mathbf{k}}=0,\Omega}G^{\text{a}}_{{\text{R}},{\mathbf{k}}=0,\Omega})$
$\displaystyle+G^{\text{r}}_{{\text{L}},{\mathbf{k}}=0,-\Omega}G^{\text{a}}_{{\text{R}},{\mathbf{k}}=0,\Omega}{\text{e}}^{2i\Omega
t}$
$\displaystyle+G^{\text{r}}_{{\text{L}},{\mathbf{k}}=0,\Omega}G^{\text{a}}_{{\text{R}},{\mathbf{k}}=0,-\Omega}{\text{e}}^{-2i\Omega
t}].$ (7)

Thus, we find from Eqs. (7) and (3) that the tunneling spin current of magnons
is described essentially as the product of the magnetic susceptibility for
each insulator. Those are summarized as follows:

$\langle{\mathcal{A}}(t)\rangle=\Gamma_{0(\text{L})}\Gamma_{0(\text{R})}\sum_{n=0,\pm
2}A_{n}{\text{e}}^{ni\Omega t},$ (8)

with the coefficients

$\displaystyle A_{2}$
$\displaystyle=G^{\text{r}}_{{\text{L}},{\mathbf{k}}=0,-\Omega}G^{\text{a}}_{{\text{R}},{\mathbf{k}}=0,\Omega},$
(9a) $\displaystyle A_{-2}$
$\displaystyle=G^{\text{r}}_{{\text{L}},{\mathbf{k}}=0,\Omega}G^{\text{a}}_{{\text{R}},{\mathbf{k}}=0,-\Omega},$
(9b) $\displaystyle A_{0}$
$\displaystyle=G^{\text{r}}_{{\text{L}},{\mathbf{k}}=0,-\Omega}G^{\text{a}}_{{\text{R}},{\mathbf{k}}=0,-\Omega}+G^{\text{r}}_{{\text{L}},{\mathbf{k}}=0,\Omega}G^{\text{a}}_{{\text{R}},{\mathbf{k}}=0,\Omega}.$
(9c)

###  III.2 dc and ac components

We evaluate ${\text{Im}}{\langle{\mathcal{A}}(t)\rangle}$ and provide an
analytical formula for the magnon spin current across the junction interface
[Eq. (3)]. To this end, defining

$A_{\Omega}\coloneqq A_{2}\in{\mathbb{C}},$ (10)

we introduce the real variables $A_{\Omega}^{\text{R}}\in{\mathbb{R}}$ and
$A_{\Omega}^{\text{I}}\in{\mathbb{R}}$ as
$A_{\Omega}=A_{\Omega}^{\text{R}}+iA_{\Omega}^{\text{I}}\in{\mathbb{C}}$ for
convenience, and describe the magnon spin current in terms of
$A_{\Omega}^{\text{R}}={\text{Re}}(A_{2})$ and
$A_{\Omega}^{\text{I}}={\text{Im}}(A_{2})$. Using the identity

$A_{-\Omega}=A_{-2},$ (11)

a straightforward calculation starting from Eq. (8) provides

$\displaystyle{\text{Im}}{\langle{\mathcal{A}}(t)\rangle}=$
$\displaystyle\Gamma_{0(\text{L})}\Gamma_{0(\text{R})}[(A_{\Omega}^{\text{I}}+A_{-\Omega}^{\text{I}}){\text{cos}}({2\Omega
t})$
$\displaystyle+(A_{\Omega}^{\text{R}}-A_{-\Omega}^{\text{R}}){\text{sin}}({2\Omega
t})+{\text{Im}}(A_{0})].$ (12)

We find that the magnon spin current [Eq. (3)] consists of two parts: The dc
component, ${\text{Im}}(A_{0})$, and the ac one. Finally, we obtain the magnon
spin current across the junction interface as

$\langle I(t)\rangle=\langle I_{\text{{ac}}}(t)\rangle+\langle
I_{\text{dc}}\rangle+O({J_{\text{ex}}}^{2}),$ (13)

where the ac component $\langle
I_{\text{{ac}}}(t)\rangle=O(J_{\text{ex}}\Gamma_{0(\text{L})}\Gamma_{0(\text{R})})$
and the dc one $\langle
I_{\text{dc}}\rangle=O(J_{\text{ex}}\Gamma_{0(\text{L})}\Gamma_{0(\text{R})})$
are described as

$\displaystyle\langle I_{\text{{ac}}}(t)\rangle\coloneqq$ $\displaystyle
2g\mu_{\text{B}}\frac{J_{\text{ex}}}{\hbar}\sqrt{(A_{\Omega}^{\text{I}}+A_{-\Omega}^{\text{I}})^{2}+(A_{\Omega}^{\text{R}}-A_{-\Omega}^{\text{R}})^{2}}$
$\displaystyle\cdot\Gamma_{0(\text{L})}\Gamma_{0(\text{R})}{\text{cos}}({2\Omega
t}+\theta_{0}),$ (14a) $\displaystyle\langle I_{\text{dc}}\rangle\coloneqq$
$\displaystyle
2g\mu_{\text{B}}\frac{J_{\text{ex}}}{\hbar}\Gamma_{0(\text{L})}\Gamma_{0(\text{R})}{\text{Im}}(A_{0}),$
(14b)

with the phase $\theta_{0}$ characterized as

$\displaystyle{\text{sin}}\theta_{0}$
$\displaystyle=-\frac{A_{\Omega}^{\text{R}}-A_{-\Omega}^{\text{R}}}{\sqrt{(A_{\Omega}^{\text{I}}+A_{-\Omega}^{\text{I}})^{2}+(A_{\Omega}^{\text{R}}-A_{-\Omega}^{\text{R}})^{2}}},$
(15a) $\displaystyle{\text{cos}}\theta_{0}$
$\displaystyle=\frac{A_{\Omega}^{\text{I}}+A_{-\Omega}^{\text{I}}}{\sqrt{(A_{\Omega}^{\text{I}}+A_{-\Omega}^{\text{I}})^{2}+(A_{\Omega}^{\text{R}}-A_{-\Omega}^{\text{R}})^{2}}}.$
(15b)

This is the general formula for magnon transport in the junction irradiated by
the laser field: The main result of this paper. For convenience, introducing
the variables $C_{\text{I}}:=A_{\Omega}^{\text{I}}+A_{-\Omega}^{\text{I}}$ and
$C_{\text{R}}:=A_{\Omega}^{\text{R}}-A_{-\Omega}^{\text{R}}$, the phase of the
ac component $\theta_{0}$ is characterized, as an example, as follows:

$\displaystyle\theta_{0}$ $\displaystyle=0\ \ \ \ \ \ \ \text{for}\
C_{\text{I}}>0\ \text{and}\ C_{\text{R}}=0.$ (16a) $\displaystyle\theta_{0}$
$\displaystyle=-\pi/2\ \text{for}\ C_{\text{I}}=0\ \text{and}\
C_{\text{R}}>0.$ (16b) $\displaystyle\theta_{0}$ $\displaystyle=-\pi/4\
\text{for}\ C_{\text{I}}=C_{\text{R}}>0.$ (16c) $\displaystyle\theta_{0}$
$\displaystyle=\pi/4\ \ \ \text{for}\ C_{\text{I}}=-C_{\text{R}}>0.$ (16d)

We find from Eq. (14a) that the ac spin current of the period $\pi/\Omega$,

$\langle I_{\text{{ac}}}(t)\rangle=\langle
I_{\text{{ac}}}(t+\pi/\Omega)\rangle,$ (17)

arises from the laser magnetic field of the period $2\pi/\Omega$ [Eqs. (2a)
and (2b)]

$\Gamma_{\text{L(R)}}(t)=\Gamma_{\text{L(R)}}(t+2\pi/\Omega).$ (18)

The period of the ac spin current is one-half of that of the laser magnetic
field. We remark that this spin current of magnons is analogous to the
Josephson spin current of magnons Nakata  _et al._ (2014); Nakata (2021) in
that those spin currents across the junction interface arise as the process of
$O(J_{\text{ex}})$ in $J_{\text{ex}}$. This is in contrast to the spin current
of thermal magnons, which arises as the process of $O({J_{\text{ex}}}^{2})$ in
$J_{\text{ex}}$, induced by the applied temperature difference Nakata  _et
al._ (2015a).

###  III.3 Inversion symmetry breaking

From a practical viewpoint, we describe $A_{\Omega}^{\text{R(I)}}$ and $A_{0}$
in terms of the magnon lifetime by phenomenologically introducing it 333 Here,
we phenomenologically assume that the finite lifetime of magnons is caused by
impurity scattering Nakata and Ohnuma (2021). through the Green’s function
Tatara  _et al._ (2008), e.g.,
$G^{\text{r}}_{{\text{L}},{\mathbf{k}}=0,\Omega}=1/[\hbar\Omega-\Delta_{\text{L}}+i\hbar/(2\tau_{\text{L}})]=(G^{\text{a}}_{{\text{L}},{\mathbf{k}}=0,\Omega})^{*}$,
where $\tau_{\text{L(R)}}$ and $\Delta_{\text{L(R)}}$ are the lifetime and the
energy gap of magnons for the zero mode in the left (right) insulator,
respectively, and $(G^{\text{a}}_{{\text{L}},{\mathbf{k}}=0,\Omega})^{*}$
represents the complex conjugate of
$G^{\text{a}}_{{\text{L}},{\mathbf{k}}=0,\Omega}$. Then, those are
characterized as See

$\displaystyle A_{\Omega}^{\text{R}}=$
$\displaystyle-\frac{(\hbar\Omega+\Delta_{\text{L}})(\hbar\Omega-\Delta_{\text{R}})-(\hbar/2\tau_{\text{L}})(\hbar/2\tau_{\text{R}})}{[(\hbar\Omega+\Delta_{\text{L}})^{2}+(\hbar/2\tau_{\text{L}})^{2}][(\hbar\Omega-\Delta_{\text{R}})^{2}+(\hbar/2\tau_{\text{R}})^{2}]},$
(19a) $\displaystyle A_{\Omega}^{\text{I}}=$
$\displaystyle-\frac{(\hbar/2\tau_{\text{L}})(\hbar\Omega-\Delta_{\text{R}})+(\hbar/2\tau_{\text{R}})(\hbar\Omega+\Delta_{\text{L}})}{[(\hbar\Omega+\Delta_{\text{L}})^{2}+(\hbar/2\tau_{\text{L}})^{2}][(\hbar\Omega-\Delta_{\text{R}})^{2}+(\hbar/2\tau_{\text{R}})^{2}]},$
(19b) $\displaystyle{\text{Im}}(A_{0})$
$\displaystyle=\frac{(\hbar/2\tau_{\text{L}})(\hbar\Omega+\Delta_{\text{R}})-(\hbar/2\tau_{\text{R}})(\hbar\Omega+\Delta_{\text{L}})}{[(\hbar\Omega+\Delta_{\text{L}})^{2}+(\hbar/2\tau_{\text{L}})^{2}][(\hbar\Omega+\Delta_{\text{R}})^{2}+(\hbar/2\tau_{\text{R}})^{2}]}$
$\displaystyle+\frac{(\hbar/2\tau_{\text{L}})(-\hbar\Omega+\Delta_{\text{R}})-(\hbar/2\tau_{\text{R}})(-\hbar\Omega+\Delta_{\text{L}})}{[(\hbar\Omega-\Delta_{\text{L}})^{2}+(\hbar/2\tau_{\text{L}})^{2}][(\hbar\Omega-\Delta_{\text{R}})^{2}+(\hbar/2\tau_{\text{R}})^{2}]}.$
(19c)

These show that tuning the laser frequency,

$\hbar\Omega=\Delta_{\text{L(R)}},$ (20)

the ferromagnetic resonance enhances those magnon spin currents across the
junction interface.

We remark that when the left insulator is identical to the right one, i.e.,
$\Delta_{\text{L}}=\Delta_{\text{R}}$ and $\tau_{\text{L}}=\tau_{\text{R}}$,
the coefficients become

$\displaystyle A_{\Omega}^{\text{I}}+A_{-\Omega}^{\text{I}}$
$\displaystyle=0,$ (21a) $\displaystyle
A_{\Omega}^{\text{R}}-A_{-\Omega}^{\text{R}}$ $\displaystyle=0,$ (21b)
$\displaystyle{\text{Im}}(A_{0})$ $\displaystyle=0.$ (21c)

This results in

$\displaystyle\langle I_{\text{{ac}}}(t)\rangle$ $\displaystyle=0,$ (22a)
$\displaystyle\langle I_{\text{dc}}\rangle$ $\displaystyle=0,$ (22b)
$\displaystyle\langle I(t)\rangle$ $\displaystyle=0,$ (22c)

and there are no magnon spin currents across the junction interface. This
shows that the dc component of the magnon spin current arises from the
difference of the magnon energy gap in the junction subjected to the laser
field: In the junction out of equilibrium, the magnon energy gap works as a
nonequilibrium spin chemical potential Basso  _et al._ (2016); Cornelissen
_et al._ (2016a, b); Du  _et al._ (2017); Demokritov  _et al._ (2006) (i.e., a
potential in an effective magnetic field Johnson and Silsbee (1987)) and the
difference induces magnon transport across the junction interface Haldane and
Arovas (1995); Meier and Loss (2003); Fujimoto (2009); Nakata  _et al._
(2017b, 2015a, 2018). When $\tau_{\text{L}}=\tau_{\text{R}}=:\tau$ while
$\Delta_{\text{R}}\neq\Delta_{\text{L}}$, the inversion symmetry across the
junction interface violates. Then, the dc component is proportional to the
difference of the magnon energy gap, $\Delta_{\text{R}}-\Delta_{\text{L}}$, as
follows:

$\displaystyle\langle I_{\text{dc}}\rangle$
$\displaystyle=2g\mu_{\text{B}}\frac{J_{\text{ex}}}{\hbar}\Gamma_{0(\text{L})}\Gamma_{0(\text{R})}\frac{\hbar}{2\tau}(\Delta_{\text{R}}-\Delta_{\text{L}})$
(23)
$\displaystyle\times\Bigg{[}\frac{1}{[(\hbar\Omega+\Delta_{\text{L}})^{2}+(\hbar/2\tau)^{2}][(\hbar\Omega+\Delta_{\text{R}})^{2}+(\hbar/2\tau)^{2}]}$
$\displaystyle+\frac{1}{[(\hbar\Omega-\Delta_{\text{L}})^{2}+(\hbar/2\tau)^{2}][(\hbar\Omega-\Delta_{\text{R}})^{2}+(\hbar/2\tau)^{2}]}\Bigg{]}.$

Hence, when $\Delta_{\text{R}}>\Delta_{\text{L}}$, the dc spin current flows
from the right to the left insulator due to the inversion symmetry breaking.

##  IV Discussion

For an estimate, we assume the following parameter values Heinrich  _et al._
(2011); Tserkovnyak  _et al._ (2005); Tupitsyn  _et al._ (2008):
$\Delta_{\text{L}}=\hbar\Omega=5\ \mu$eV, $\Delta_{\text{R}}=6\ \mu$eV,
$\Gamma_{0(\text{L})}=\Gamma_{0(\text{R})}=0.5\ \mu$eV, $\tau=10\ \mu$s
Heinrich  _et al._ (2011); Tserkovnyak  _et al._ (2005); Nakata and Ohnuma
(2021), $J_{\text{ex}}=35\ \mu$eV Tupitsyn  _et al._ (2008). Then, following
Ref. Proskurin  _et al._ (2018), we estimate the dc spin current $\langle
I_{\text{dc}}\rangle$ in electric units $e/g\mu_{\text{B}}$, and find that it
amounts to $\langle I_{\text{dc}}\rangle e/g\mu_{\text{B}}\sim 0.05$ mA. We
believe, while being small compared with the one ($0.1$ mA) Proskurin  _et
al._ (2018) for antiferromagnets, still it will be within experimental reach.
We expect, to the best of our knowledge, that to use the inverse spin Hall
effect Saitoh  _et al._ (2006) by attaching a metal to the insulating magnet
will be one of the most promising strategies for observation.

In this paper focusing on the junction and assuming the weak exchange coupling
between the left and right ferromagnetic insulator, we have studied the
tunneling spin current of magnons across the junction interface irradiated by
linearly polarized laser. It will be of significance to develop this work into
bulk insulators by taking into account higher-order terms in the Holstein-
Primakoff expansion and thus studying the nonlinear terms in the magnetic
susceptibility. We leave the advanced study for future work. We remark that
Ref. Proskurin  _et al._ (2018) proposed the optically induced spin current of
magnons in the bulk of antiferromagnetic insulators irradiated by a circularly
polarized electromagnetic field, where, analogous to our present work, the
field-induced spin current is described by the product of the susceptibility,
and the current is enhanced by (antiferromagnetic) resonance.

##  V Conclusion

Using the junction of the ferromagnetic insulators irradiated by linearly
polarized laser, we have proposed the method for the generation of dc and ac
spin currents of magnons across the junction interface, and provided the
analytical formula for magnon transport induced by the inversion symmetry
breaking. We have shown that the applied laser field excites the zero mode of
magnons, and tuning the laser frequency to the magnon energy gap for the zero
mode, the ferromagnetic resonance enhances those spin currents across the
junction interface. The period of the ac spin current is one-half of that of
the laser magnetic field. We hope that our proposal serves as a key ingredient
for magnonics.

###### Acknowledgements.

We would like to thank S. Takayoshi for fruitful discussions and feedback on
this manuscript. The author (K. N.) is grateful also to S. Uchino and Y. Araki
for useful discussions, and H. Chudo for helpful feedback on the experimental
feasibility. We also sincerely thank anonymous referees for the valuable
comments that helped us drastically improve the manuscript. We acknowledge
support by JSPS KAKENHI Grant Number JP20K14420 (K. N.) and JP22K03519 (K.
N.), by Leading Initiative for Excellent Young Researchers, MEXT, Japan (K.
N.), and by JST ERATO Grant No. JPMJER1601 (Y. O.). *

##  Appendix A Magnon transport across the junction interface

In this Appendix, we provide some details of the straightforward calculation
for the magnon spin current across the junction interface $\langle
I(t)\rangle=O(J_{\text{ex}}{\Gamma_{0({\text{L/R}})}}^{2})$,

$\langle
I(t)\rangle=2g\mu_{\text{B}}\frac{J_{\text{ex}}}{\hbar}{\text{Im}}\langle{\mathcal{A}}(t)\rangle+O({J_{\text{ex}}}^{2}),$
(24)

by evaluating
$\langle{\mathcal{A}}(t)\rangle=O({J_{\text{ex}}}^{0}{\Gamma_{0({\text{L/R}})}}^{2})$
of

${\mathcal{A}}(t)=\sum_{{\mathbf{k}},k_{x}^{\prime}}a_{{\text{L}},{\mathbf{k}}}(t)a^{\dagger}_{{\text{R}},{\mathbf{k}}^{\prime}}(t).$
(25)

The number operator of magnons in the left insulator is
$N_{\text{L}}=\sum_{{\mathbf{q}}}a^{\dagger}_{{\text{L}},{\mathbf{q}}}a_{{\text{L}},{\mathbf{q}}}$.
Through the Heisenberg equation of motion, we define the operator for the
magnon current across the junction interface $I_{\text{m}}(t)$, from the right
to the left insulator, as

$I_{\text{m}}(t)\coloneqq\frac{1}{i\hbar}[N_{\text{L}},{\mathcal{H}}_{\text{ex}}].$
(26)

Then, the magnon spin current, i.e., the spin current carried by magnons,
across the junction interface $I(t)\coloneqq g\mu_{\text{B}}I_{\text{m}}(t)$
Meier and Loss (2003); Nakata  _et al._ (2015a, 2018) is given as

$I(t)=-ig\mu_{\text{B}}\frac{J_{\text{ex}}}{\hbar}\sum_{{\mathbf{k}},k_{x}^{\prime}}a_{{\text{L}},{\mathbf{k}}}(t)a^{\dagger}_{{\text{R}},{\mathbf{k}}^{\prime}}(t)+{\text{H.
c.}},$ (27)

which results in Eqs. (24) and (25).

Using the contour-ordered Green’s function Haug and Jauho (2007); Kita (2010);
Tatara  _et al._ (2008) we perform the perturbation calculation in
$\Gamma_{0({\text{L/R}})}$ based on the Schwinger-Keldysh formalism Schwinger
(1961); Martin and Schwinger (1959); Keldysh (1965), and up to
${\Gamma_{0({\text{L/R}})}}^{2}$, which provides

$\langle{\mathcal{A}}(t)\rangle=\frac{1}{\hbar^{2}}\int_{{\text{c}}}d\tau_{1}\int_{{\text{c}}}d\tau_{2}\Gamma_{\text{L}}(\tau_{1})\Gamma_{\text{R}}(\tau_{2})G_{{\text{L}},{\mathbf{k}}=0}(\tau,\tau_{1})G_{{\text{R}},{\mathbf{k}}=0}(\tau_{2},\tau)|_{\tau\rightarrow
t},$ (28)

where the applied laser magnetic field excites the zero mode of magnons.
Taking the Schwinger-Keldysh closed time path Rammer (2007)
${\text{c}}={\text{c}}_{\rightarrow}+{\text{c}}_{\leftarrow}$,

$\int_{{\text{c}}={\text{c}}_{\rightarrow}+{\text{c}}_{\leftarrow}}d\tau_{1}\int_{{\text{c}}={\text{c}}_{\rightarrow}+{\text{c}}_{\leftarrow}}d\tau_{2}=\int_{{\text{c}}_{\rightarrow}}d\tau_{1}\int_{{\text{c}}_{\rightarrow}}d\tau_{2}+\int_{{\text{c}}_{\rightarrow}}d\tau_{1}\int_{{\text{c}}_{\leftarrow}}d\tau_{2}+\int_{{\text{c}}_{\leftarrow}}d\tau_{1}\int_{{\text{c}}_{\rightarrow}}d\tau_{2}+\int_{{\text{c}}_{\leftarrow}}d\tau_{1}\int_{{\text{c}}_{\leftarrow}}d\tau_{2},$
(29)

and using the Langreth theorem Langreth and Wilkins (1972); Tatara  _et al._
(2008); Haug and Jauho (2007), it becomes

$\displaystyle\langle{\mathcal{A}}(t)\rangle=$
$\displaystyle\frac{1}{\hbar^{2}}\int_{-\infty}^{\infty}dt_{1}\int_{-\infty}^{\infty}dt_{2}\Gamma_{\text{L}}(t_{1})\Gamma_{\text{R}}(t_{2})$
(30)
$\displaystyle\cdot[G^{\text{t}}_{{\text{L}},{\mathbf{k}}=0}(t,t_{1})G^{\text{t}}_{{\text{R}},{\mathbf{k}}=0}(t_{2},t)-G^{\text{t}}_{{\text{L}},{\mathbf{k}}=0}(t,t_{1})G^{\text{>}}_{{\text{R}},{\mathbf{k}}=0}(t_{2},t)-G^{\text{<}}_{{\text{L}},{\mathbf{k}}=0}(t,t_{1})G^{\text{t}}_{{\text{R}},{\mathbf{k}}=0}(t_{2},t)+G^{\text{<}}_{{\text{L}},{\mathbf{k}}=0}(t,t_{1})G^{\text{>}}_{{\text{R}},{\mathbf{k}}=0}(t_{2},t)].$

Since Haug and Jauho (2007); Kita (2010); Tatara  _et al._ (2008)
$G^{\text{t}}_{{\mathbf{k}}}(t,t_{1})=G^{\text{r}}_{{\mathbf{k}}}(t,t_{1})+G^{<}_{{\mathbf{k}}}(t,t_{1})$
and
$G^{<}_{{\mathbf{k}}}(t,t_{1})-G^{>}_{{\mathbf{k}}}(t,t_{1})=G^{\text{a}}_{{\mathbf{k}}}(t,t_{1})-G^{\text{r}}_{{\mathbf{k}}}(t,t_{1})$,
it reduces to

$\langle{\mathcal{A}}(t)\rangle=\frac{1}{\hbar^{2}}\int_{-\infty}^{\infty}dt_{1}\int_{-\infty}^{\infty}dt_{2}\Gamma_{\text{L}}(t_{1})\Gamma_{\text{R}}(t_{2})G^{\text{r}}_{{\text{L}},{\mathbf{k}}=0}(t,t_{1})G^{\text{a}}_{{\text{R}},{\mathbf{k}}=0}(t_{2},t),$
(31)

where

$\Gamma_{\text{L}}(t_{1})\Gamma_{\text{R}}(t_{2})=\Gamma_{0(\text{L})}\Gamma_{0(\text{R})}[{\text{e}}^{i\Omega(t_{1}+t_{2})}+{\text{e}}^{i\Omega(t_{1}-t_{2})}+{\text{e}}^{i\Omega(-t_{1}+t_{2})}+{\text{e}}^{i\Omega(-t_{1}-t_{2})}].$
(32)

The Fourier transform of the Green’s functions results in

$\displaystyle\langle{\mathcal{A}}(t)\rangle=$
$\displaystyle\Gamma_{0(\text{L})}\Gamma_{0(\text{R})}\int_{-\infty}^{\infty}dt_{1}\int_{-\infty}^{\infty}dt_{2}\int\frac{d\omega_{1}}{2\pi}\int\frac{d\omega_{2}}{2\pi}G^{\text{r}}_{{\text{L}},{\mathbf{k}}=0,\omega_{1}}G^{\text{a}}_{{\text{R}},{\mathbf{k}}=0,\omega_{2}}{\text{e}}^{i(-\omega_{1}+\omega_{2})t}$
$\displaystyle{\cdot}[{\text{e}}^{i(\omega_{1}+\Omega)t_{1}}{\text{e}}^{i(-\omega_{2}+\Omega)t_{2}}+{\text{e}}^{i(\omega_{1}+\Omega)t_{1}}{\text{e}}^{i(-\omega_{2}-\Omega)t_{2}}+{\text{e}}^{i(\omega_{1}-\Omega)t_{1}}{\text{e}}^{i(-\omega_{2}+\Omega)t_{2}}+{\text{e}}^{i(\omega_{1}-\Omega)t_{1}}{\text{e}}^{i(-\omega_{2}-\Omega)t_{2}}].$
(33)

Performing the time integral and the resulting Dirac delta function for the
frequency, we obtain

$\displaystyle\langle{\mathcal{A}}(t)\rangle=\Gamma_{0(\text{L})}\Gamma_{0(\text{R})}[$
$\displaystyle
G^{\text{r}}_{{\text{L}},{\mathbf{k}}=0,-\Omega}G^{\text{a}}_{{\text{R}},{\mathbf{k}}=0,\Omega}{\text{e}}^{2i\Omega
t}+G^{\text{r}}_{{\text{L}},{\mathbf{k}}=0,\Omega}G^{\text{a}}_{{\text{R}},{\mathbf{k}}=0,-\Omega}{\text{e}}^{-2i\Omega
t}$
$\displaystyle+(G^{\text{r}}_{{\text{L}},{\mathbf{k}}=0,-\Omega}G^{\text{a}}_{{\text{R}},{\mathbf{k}}=0,-\Omega}+G^{\text{r}}_{{\text{L}},{\mathbf{k}}=0,\Omega}G^{\text{a}}_{{\text{R}},{\mathbf{k}}=0,\Omega})].$
(34)

This is summarized as follows:

$\langle{\mathcal{A}}(t)\rangle=\Gamma_{0(\text{L})}\Gamma_{0(\text{R})}\sum_{n=0,\pm
2}A_{n}{\text{e}}^{ni\Omega t},$ (35)

with the coefficients

$\displaystyle A_{2}$
$\displaystyle=G^{\text{r}}_{{\text{L}},{\mathbf{k}}=0,-\Omega}G^{\text{a}}_{{\text{R}},{\mathbf{k}}=0,\Omega},$
(36a) $\displaystyle A_{-2}$
$\displaystyle=G^{\text{r}}_{{\text{L}},{\mathbf{k}}=0,\Omega}G^{\text{a}}_{{\text{R}},{\mathbf{k}}=0,-\Omega},$
(36b) $\displaystyle A_{0}$
$\displaystyle=G^{\text{r}}_{{\text{L}},{\mathbf{k}}=0,-\Omega}G^{\text{a}}_{{\text{R}},{\mathbf{k}}=0,-\Omega}+G^{\text{r}}_{{\text{L}},{\mathbf{k}}=0,\Omega}G^{\text{a}}_{{\text{R}},{\mathbf{k}}=0,\Omega}.$
(36c)

The coefficients are described in terms of the magnon lifetime and the energy
gap for the zero mode as

$\displaystyle A_{2}$ $\displaystyle\eqqcolon
A_{\Omega}=-\frac{[\hbar\Omega+\Delta_{\text{L}}+i\hbar/(2\tau_{\text{L}})][\hbar\Omega-\Delta_{\text{R}}+i\hbar/(2\tau_{\text{R}})]}{[(\hbar\Omega+\Delta_{\text{L}})^{2}+(\hbar/2\tau_{\text{L}})^{2}][(\hbar\Omega-\Delta_{\text{R}})^{2}+(\hbar/2\tau_{\text{R}})^{2}]},$
(37a) $\displaystyle A_{-2}$ $\displaystyle=A_{-\Omega},$ (37b) $\displaystyle
A_{0}$
$\displaystyle=\frac{[\hbar\Omega+\Delta_{\text{L}}+i\hbar/(2\tau_{\text{L}})][\hbar\Omega+\Delta_{\text{R}}-i\hbar/(2\tau_{\text{R}})]}{[(\hbar\Omega+\Delta_{\text{L}})^{2}+(\hbar/2\tau_{\text{L}})^{2}][(\hbar\Omega+\Delta_{\text{R}})^{2}+(\hbar/2\tau_{\text{R}})^{2}]}+\frac{[-\hbar\Omega+\Delta_{\text{L}}+i\hbar/(2\tau_{\text{L}})][-\hbar\Omega+\Delta_{\text{R}}-i\hbar/(2\tau_{\text{R}})]}{[(\hbar\Omega-\Delta_{\text{L}})^{2}+(\hbar/2\tau_{\text{L}})^{2}][(\hbar\Omega-\Delta_{\text{R}})^{2}+(\hbar/2\tau_{\text{R}})^{2}]}.$
(37c)

## References

  * Chumak  _et al._ (2015) A. V. Chumak, V. I. Vasyuchka, A. A. Serga, and B. Hillebrands, Nat. Phys. 11, 453 (2015).
  * Nakata  _et al._ (2017a) K. Nakata, P. Simon, and D. Loss, J. Phys. D: Appl. Phys. 50, 114004 (2017a).
  * Ciappina  _et al._ (2017) M. F. Ciappina, J. A. P.-Hernandez, A. S. Landsman, W. A. Okell, S. Zherebtsov, B. Forg, J. Schotz, L. Seiffert, T. Fennel, T. Shaaran, T. Zimmermann, A. Chacon, R. Guichard, A. Zair, J. W. G. Tisch, J. P. Marangos, T. Witting, A. Braun, S. A. Maier, L. Roso, M. Kruger, P. Hommelhoff, M. F. Kling, F. Krausz, and M. Lewenstein, Rep. Prog. Phys. 80, 054401 (2017).
  * Mukai  _et al._ (2014) Y. Mukai, H. Hirori, T. Yamamoto, H. Kageyama, and K. Tanaka, Appl. Phys. Lett. 105, 022410 (2014).
  * Bossini  _et al._ (2016) D. Bossini, V. I. Belotelov, A. K. Zvezdin, A. N. Kalish, and A. V. Kimel, ACS Photon. 3, 1385 (2016).
  * Arikawa  _et al._ (2017) T. Arikawa, S. Morimoto, and K. Tanaka, Opt. Express 25, 13728 (2017).
  * Kirilyuk  _et al._ (2010) A. Kirilyuk, A. V. Kimel, and T. Rasing, Rev. Mod. Phys. 82, 2731 (2010).
  * Kimel  _et al._ (2004) A. V. Kimel, A. Kirilyuk, A. Tsvetkov, R. V. Pisarev, and T. Rasing, Nature 429, 850 (2004).
  * Kimel  _et al._ (2005) A. V. Kimel, A. Kirilyuk, P. A. Usachev, R. V. Pisarev, A. M. Balbashov, and T. Rasing, Nature 435, 655 (2005).
  * Hansteen  _et al._ (2005) F. Hansteen, A. Kimel, A. Kirilyuk, and T. Rasing, Phys. Rev. Lett. 95, 047402 (2005).
  * Stanciu  _et al._ (2007a) C. D. Stanciu, F. Hansteen, A. V. Kimel, A. Tsukamoto, A. Itoh, A. Kirilyuk, and T. Rasing, Phys. Rev. Lett. 98, 207401 (2007a).
  * N$\check{\rm{e}}$mec  _et al._ (2018) P. N$\check{\rm{e}}$mec, M. Fiebig, T. Kampfrath, and A. V. Kimel, Nat. Phys. 14, 229 (2018).
  * Stanciu  _et al._ (2007b) C. D. Stanciu, F. Hansteen, A. V. Kimel, A. Kirilyuk, A. Tsukamoto, A. Itoh, and T. Rasing, Phys. Rev. Lett. 99, 047601 (2007b).
  * Radu  _et al._ (2011) I. Radu, K. Vahaplar, C. Stamm, T. Kachel, N. Pontius, H. A. D$\ddot{{\mathrm{u}}}$rr, T. A. Ostler, J. Barker, R. F. L. Evans, R. W. Chantrell, A. Tsukamoto, A. Itoh, A. Kirilyuk, T. Rasing, and A. V. Kimel, Nature 472, 205 (2011).
  * Ostler  _et al._ (2012) T. A. Ostler, J. Barker, R. F. L. Evans, R. W. Chantrell, U. Atxitia, O. C.-Fesenko, S. E. Moussaoui, L. L. Guyader, E. Mengotti, L. J. Heyderman, F. Nolting, A. Tsukamoto, A. Itoh, D. Afanasiev, B. A. Ivanov, A. M. Kalashnikova, K. Vahaplar, J. Mentink, A. Kirilyuk, T. Rasing, and A. Kimel, Nat. Commun. 3, 666 (2012).
  * Vahaplar  _et al._ (2009) K. Vahaplar, A. M. Kalashnikova, A. V. Kimel, D. Hinzke, U. Nowak, R. Chantrell, A. Tsukamoto, A. Itoh, A. Kirilyuk, and T. Rasing, Phys. Rev. Lett. 103, 117201 (2009).
  * Barnett (1915) S. J. Barnett, Phys. Rev. 6, 239 (1915).
  * Barnett (1935) S. J. Barnett, Rev. Mod. Phys. 7, 129 (1935).
  * Matsuo  _et al._ (2017) M. Matsuo, E. Saitoh, and S. Maekawa, J. Phys. Soc. Jpn. 86, 011011 (2017).
  * Rebei and Hohlfeld (2008a) A. Rebei and J. Hohlfeld, Phys. Lett. A 372, 1915 (2008a).
  * Rebei and Hohlfeld (2008b) A. Rebei and J. Hohlfeld, J. Appl. Phys. 103, 07B118 (2008b).
  * Takayoshi  _et al._ (2014a) S. Takayoshi, M. Sato, and T. Oka, Phys. Rev. B 90, 214413 (2014a).
  * Takayoshi  _et al._ (2014b) S. Takayoshi, H. Aoki, and T. Oka, Phys. Rev. B 90, 085150 (2014b).
  * Note (1) See Refs. Takayoshi  _et al._ (2014a, b) for the importance of modulating laser frequency adiabatically by the chirping technique Sato  _et al._ (2013); Kamada  _et al._ (2013). See also Ref. Nakata and Takayoshi (2020) for the difference between the optical Barnett effect and the inverse Faraday effect Kirilyuk  _et al._ (2010); Kimel  _et al._ (2005).
  * Nakata (2021) K. Nakata, Phys. Rev. B 104, 104402 (2021).
  * Nakata  _et al._ (2015a) K. Nakata, P. Simon, and D. Loss, Phys. Rev. B 92, 134425 (2015a).
  * Nakata  _et al._ (2018) K. Nakata, Y. Ohnuma, and M. Matsuo, Phys. Rev. B 98, 094430 (2018).
  * Holstein and Primakoff (1940) T. Holstein and H. Primakoff, Phys. Rev. 58, 1098 (1940).
  * Note (2) Note that in Ref. Nakata  _et al._ (2015b), the laser field is applied only to one side.
  * Nakata  _et al._ (2015b) K. Nakata, P. Simon, and D. Loss, Phys. Rev. B 92, 014422 (2015b).
  * Nakata  _et al._ (2017b) K. Nakata, S. K. Kim, J. Klinovaja, and D. Loss, Phys. Rev. B 96, 224414 (2017b).
  * Nakata and Ohnuma (2021) K. Nakata and Y. Ohnuma, Phys. Rev. B 104, 064408 (2021).
  * Tatara  _et al._ (2008) G. Tatara, H. Kohno, and J. Shibata, Phys. Rep. 468, 213 (2008).
  * (34) See the Appendix for details.
  * Schwinger (1961) J. Schwinger, J. Math. Phys. 2, 407 (1961).
  * Martin and Schwinger (1959) P. C. Martin and J. Schwinger, Phys. Rev. 115, 1342 (1959).
  * Keldysh (1965) L. V. Keldysh, Sov. Phys. JETP 20, 1018 (1965).
  * Haug and Jauho (2007) H. Haug and A. Jauho,  _Quantum Kinetics in Transport and Optics of Semiconductors_ (Springer New York, 2007).
  * Kita (2010) T. Kita, Prog. Theor. Phys. 123, 581 (2010).
  * Rammer (2007) J. Rammer,  _Quantum Field Theory of Non-equilibrium States_ (Cambridge University Press, Cambridge, 2007).
  * Langreth and Wilkins (1972) D. C. Langreth and J. W. Wilkins, Phys. Rev. B 6, 3189 (1972).
  * Nakata  _et al._ (2014) K. Nakata, K. A. van Hoogdalem, P. Simon, and D. Loss, Phys. Rev. B 90, 144419 (2014).
  * Note (3) Here, we phenomenologically assume that the finite lifetime of magnons is caused by impurity scattering Nakata and Ohnuma (2021).
  * Basso  _et al._ (2016) V. Basso, E. Ferraro, and M. Piazzi, Phys. Rev. B 94, 144422 (2016).
  * Cornelissen  _et al._ (2016a) L. J. Cornelissen, K. J. H. Peters, G. E. W. Bauer, R. A. Duine, and B. J. van Wees, Phys. Rev. B 94, 014412 (2016a).
  * Cornelissen  _et al._ (2016b) L. J. Cornelissen, J. Shan, and B. J. van Wees, Phys. Rev. B 94, 180402 (2016b).
  * Du  _et al._ (2017) C. Du, T. V. der Sar, T. X. Zhou, P. Upadhyaya, F. Casola, H. Zhang, M. C. Onbasli, C. A. Ross, R. L. Walsworth, Y. Tserkovnyak, and A. Yacoby, Science 357, 195 (2017).
  * Demokritov  _et al._ (2006) S. O. Demokritov, V. E. Demidov, O. Dzyapko, G. A. Melkov, A. A. Serga, B. Hillebrands, and A. N. Slavin, Nature (London) 443, 430 (2006).
  * Johnson and Silsbee (1987) M. Johnson and R. H. Silsbee, Phys. Rev. B 35, 4959 (1987).
  * Haldane and Arovas (1995) F. D. M. Haldane and D. P. Arovas, Phys. Rev. B 52, 4223 (1995).
  * Meier and Loss (2003) F. Meier and D. Loss, Phys. Rev. Lett. 90, 167204 (2003).
  * Fujimoto (2009) S. Fujimoto, Phys. Rev. Lett. 103, 047203 (2009).
  * Heinrich  _et al._ (2011) B. Heinrich, C. Burrowes, E. Montoya, B. Kardasz, E. Girt, Y.-Y. Song, Y. Sun, and M. Wu, Phys. Rev. Lett. 107, 066604 (2011).
  * Tserkovnyak  _et al._ (2005) Y. Tserkovnyak, A. Brataas, G. E. W. Bauer, and B. I. Halperin, Rev. Mod. Phys. 77, 1375 (2005).
  * Tupitsyn  _et al._ (2008) I. S. Tupitsyn, P. C. E. Stamp, and A. L. Burin, Phys. Rev. Lett. 100, 257202 (2008).
  * Proskurin  _et al._ (2018) I. Proskurin, A. S. Ovchinnikov, J.-i. Kishine, and R. L. Stamps, Phys. Rev. B 98, 134422 (2018).
  * Saitoh  _et al._ (2006) E. Saitoh, M. Ueda, H. Miyajima, and G. Tatara, Appl. Phys. Lett. 88, 182509 (2006).
  * Sato  _et al._ (2013) M. Sato, T. Higuchi, N. Kanda, K. Konishi, K. Yoshioka, T. Suzuki, K. Misawa, and M. K.-Gonokami, Nat. Photonics 7, 724 (2013).
  * Kamada  _et al._ (2013) S. Kamada, S. Murata, and T. Aoki, Appl. Phys. Express 6, 032701 (2013).
  * Nakata and Takayoshi (2020) K. Nakata and S. Takayoshi, Phys. Rev. B 102, 094417 (2020).