#  Triplon current generation in solids

Yao Chen1 Masahiro Sato2,† Yifei Tang1 Yuki Shiomi3 Koichi Oyanagi1,4
Takatsugu Masuda5 Yusuke Nambu1,6,7 Masaki Fujita1 & Eiji Saitoh1,8,9,10,11,∗

###### Abstract

## Abstract

A triplon refers to a fictitious particle that carries angular momentum $S$ =
1 corresponding to the elementary excitation in a broad class of quantum
dimerized spin systems. Such systems without magnetic order have long been
studied as a testing ground for quantum properties of spins. Although triplons
have been found to play a central role in thermal and magnetic properties in
dimerized magnets with singlet correlation, a spin angular momentum flow
carried by triplons, a triplon current, has not been detected yet. Here we
report spin Seebeck effects induced by a triplon current: triplon spin Seebeck
effect, using a spin-Peierls system CuGeO3. The result shows that the heating-
driven triplon transport induces spin current whose sign is positive, opposite
to the spin-wave cases in magnets. The triplon spin Seebeck effect persists
far below the spin-Peierls transition temperature, being consistent with a
theoretical calculation for triplon spin Seebeck effects.

Institute for Materials Research, Tohoku University, Sendai 980-8577, Japan

Department of Physics, Ibaraki University, Mito, Ibaraki 310-8512, Japan

Department of Basic Science, The University of Tokyo, Tokyo 153-8902, Japan

Faculty of Science and Engineering, Iwate University, Morioka 020-8551, Japan

Institute of Solid State Physics, The University of Tokyo, Kashiwa, Chiba
277-8581 Japan

FOREST, Japan Science and Technology Agency, Kawaguchi, Saitama 332-0012,
Japan

Organization for Advanced Studies, Tohoku University, Sendai 980-8577, Japan

Department of Applied Physics, The University of Tokyo, Tokyo 113-8656, Japan

Institute for AI and Beyond, The University of Tokyo, Tokyo 113-8656, Japan

Advanced Institute for Materials Research, Tohoku University, Sendai 980-8577,
Japan

Advanced Science Research Center, Japan Atomic Energy Agency, Tokai 319-1195,
Japan

* e-mail address<EMAIL_ADDRESS>
$\dagger$ e-mail address<EMAIL_ADDRESS>
## Introduction

Spin Seebeck effects[1] (SSEs) refer to the generation of a spin current, a
flow of spin angular momentum of electrons, from a temperature gradient
applied to a bilayer system comprising a magnet and heavy metal such as Pt.
Spin current in the magnet propagates along the temperature gradient and
reaches the heavy metal[2]. When the randomly polarized electrons in the heavy
metal scatter at the interface with magnet, the electron spin experience an
exchange interaction form the spin current carrier of the magnet. Such
interaction aligns spins in the heavy metal along the spin direction of spin
current carrier in the magnet and induces spin accumulation at the heavy metal
side of interface. The diffusion of such spin accumulation causes a pure spin
current in the heavy metal, in which spin current can be detected as a voltage
signal via the inverse spin-Hall effect (ISHE)[3, 4]. The SSE has been
observed in various insulating systems, in which parasitic thermal effects
originating from itinerant electrons can be eliminated. Some types of spin
carriers including magnons in magnets[5, 6, 7, 8], paramagnons[9],
antiferromagnetic magnons[7, 8], and spinons in a spin liquid[10, 11] have
been demonstrated for different mechanisms of the SSE.

Among quantum spin systems without magnetic order, dimerized magnets spin
systems occupy an important position, in which two neighboring spins are
frozen as $S=0$ singlets in the ground state. The elementary spin excitations
are $S=1$ triplet states, called a triplon. A triplon has been expected to
carry spin angular momentum.

In this work, we report the observation of spin current in CuGeO3, carried by
triplons in terms of longitudinal SSE measurements. A typical spin-dimer
system is a spin-Peierls material CuGeO3, which contains one-dimensonal
spin-1/2 chains with antiferromagnetic exchange interaction for nearest-
neighbor spins. In the spin chain, with lowering the temperature, neighboring
spins dimerizes to form a spin-gapped phase; the transition is called a spin
Peierls (SP) transition[12]. Below the SP transition temperature
$T_{\mathrm{SP}}$, the chain distorts so that the distances between
neighboring spins change alternately. The bond-alternating exchange
interaction causes neighboring spins to dimerize to reduce the total energy,
creating a gap in the spin excitation energy spectrum. CuGeO3 is the firstly
discovered inorganic system that exhibits a SP transition[13]. High-quality
single crystals of CuGeO3 are easier to obtain than organic SP materials and
CuGeO3 has typically been used to study spin excitations in SP systems[14].

## Results

###  0.1 Sample structure and the concept of the study.

Figure 1a shows the crystal structure of CuGeO3. One dimensional Cu2+ spin
($S=1/2$) chains are aligned along the $c$-axis, as illustrated as dotted
lines. The nearest-neighbor exchange interaction along the $c$-axis is
$J_{c}\sim 120$ K, much greater than the interchain exchange coupling[15],
$J_{b}\sim 0.1J_{c}$ and $J_{a}\sim-0.01J_{c}$ along the $b$ and $a$-axis,
respectively. CuGeO3 is thus considered as a quasi-one-dimensional spin
system. As the temperature $T$ decreases below the SP transition temperature
$T_{\mathrm{SP}}\sim 14.5$ K, the lattice of CuGeO3 is spontaneously distorted
to form a bond alternating configuration[16]. The spin configuration in the
ground state is shown schematically in Fig. 1b, where all spins are frozen as
dimerized $S=0$ singlets. The elementary spin excitation from the ground state
is an $S=1$ triplon [17], also shown in Fig. 1b. The excitation gap of a
triplon is estimated to be $\sim 23$ K from ESR[18] and neutron scattering[19,
20, 15] experiments.

Based on the previous neutron scattering results[19] and theoretical
calculations[21], the dispersion relation of triplons along the $c$-axis is
sketched in Fig. 1c. In the absence of a magnetic field, $\mu_{0}H=0$, the
triplon states are threefold degenerated:
$\left|\uparrow\uparrow\right\rangle$,
$\Bigl{(}\left|\uparrow\downarrow\right\rangle+\left|\downarrow\uparrow\right\rangle\Bigr{)}/\sqrt{2}$,
and $\left|\downarrow\downarrow\right\rangle$ corresponding to states with
different spin quantum numbers, $S_{z}=1,0$, and $-1$, respectively. When a
magnetic field is applied, the degenerated triplet bands split into three due
to the Zeeman effect, where the triplet state with the parallel spin direction
to $H$, $\left|\uparrow\uparrow\right\rangle$, has the lowest energy[18, 19,
20]. The $\left|\uparrow\uparrow\right\rangle$ state thus exhibits the highest
probability for thermal excitation, while
$\left|\downarrow\downarrow\right\rangle$ the lowest. Driven by a thermal
gradient, different occupancy among the three states can result in a flow of a
net spin angular momentum in an external magnetic field: triplon SSE. One of
the key features of the triplon SSE is that it should exhibit the opposite
sign of ISHE signal to the magnon SSE in ferri- and ferro- magnets. This is
because triplon spin currents are carried by excitation with the parallel spin
direction to the external field; while magnons carry antiparallel spins to the
external field, as shown in Fig. 1d. Due to the opposite polarization
direction of spin current carrier in these two systems, the electron spin in
Pt feels reverse direction of exchange interaction, resulting in spin
accumulation with reverse polarization direction and different sign in ISHE.
Our theoretical calculation based on the tunnel spin current (see
Supplementary Note G3 for details) suggests that the opposite sign of triplon
SSE to magnon SSE is a universal property for gapped triplet-spin systems.

###  0.2 Material characterization and measurement setup.

First, we measured temperature dependence of magnetization, $M(T)$, for the
CuGeO3 sample. Figure 1e shows the magnetic susceptibility $\chi(T)=M(T)/H$
measured at $\mu_{0}H=1$ T $\parallel$ $b$-axis. When $T$ is decreased from 20
K, $\chi$ exhibits a sudden drop at the SP transition temperature
$T_{\mathrm{SP}}\sim 14.5$ K and rapidly decreases towards zero. This
indicates that Cu spins along the spin chain ($\parallel c$-axis) form dimers
and a spin gap develops below $T_{\mathrm{SP}}$. By measuring $\chi(T)$ at
different values of $H$ (see Supplementary Fig. 1 for details), we obtain the
$H$ \- $T$ spin-phase diagram of the CuGeO3 sample in Fig. 1f. In the phase
diagram, for large $H$ in the SP phase, CuGeO3 undergoes another phase
transition into a magnetic phase at a transition field, $H_{\mathrm{m}}$,
where magnetization is recovered by forming a spin soliton lattice in the
chain[22]. The transition can be confirmed as a sudden jump of magnetization
as a function of fields[23] (also see Supplementary Fig. 1). The obtained
phase diagram is in good agreement with a previous study[24].

The experimental setup for the SSE measurement and the sample used in the
present study are shown in Figs. 2a and 2b. A trilayer structure comprising a
Pt wire for ISHE, an insulating SiO2 layer, and an Au wire as a heater was
fabricated on the top of the $(001)$ plane of a CuGeO3 single crystal. The
spin Seebeck voltage signal $V_{\mathrm{SSE}}$ in the Pt wire is measured by a
lock-in method, where an a.c. current is applied to the Au layer under an in-
plane magnetic field, $H$, applied perpendicular to the Pt wire. We normalize
$V_{\mathrm{SSE}}$ by the heater power and the detector resistance as
$\tilde{V}_{\mathrm{SSE}}=V_{\mathrm{SSE}}/(P_{\mathrm{Au}}\cdot
R_{\mathrm{Pt}})$.

###  0.3 Observation of the triplon SSE.

The magnetic field dependence of $\tilde{V}_{\mathrm{SSE}}$ at some selected
temperatures is shown in Figs. 2c and 2d. At $15$ K, just above
$T_{\mathrm{SP}}\sim 14.5$ K, no $H$-dependent signal is recognized. Most
notably, a clear voltage signal appears when $T$ is decreased down to $2$ K.
The signal is an odd function of $H$, reflecting the symmetry of the ISHE[25].
The sign of the $\tilde{V}_{\mathrm{SSE}}$ signal is opposite to that of the
magnon-mediated SSE[26, 27], which is confirmed in a similar setup for a
ferrimagnetic insulator Y3Fe5O12/Pt sample in Figs. 2e and 2f. The result
indicates that spin current carriers in these two systems have different spin
polarization directions with respect to $H$, consistent with the triplon-
current senario. Note that the sign of the SSE signals in spin nematic
LiCuVO4/Pt (ref. [28]) is also opposite to the present CuGeO3/Pt case. We also
note that, unlike the magnon-mediated SSE in ferro/ferrimagnets, the triplon-
mediated SSE dose not scale with the $M(H)$ curve. For the magnon-mediated
SSE[27], $M(H)$ represents the increase in the spin-current polarization. In
contrast, the magnetization increase in the SP phase means larger density of
broken dimers and the scattering probability between triplons and broken
dimers becomes lager. The unusual $H$-dependence of the triplon SSE originates
from the different properties between ferro/ferrimagnets and SP phases, and it
is one of the characteristics of triplon SSE. (see Supplementary Fig. 2 for
details).

As shown in Figs. 3a and 3b, the temperature dependence of
$\tilde{V}_{\mathrm{SSE}}$ shows that the signal magnitude increases
monotonically as $T$ decreases from $T_{\mathrm{SP}}$ to $3$ K. When $T<3$ K,
the SSE voltage remains almost unchanged down to $2$ K. At higher
temperatures, the dimerized ground state is weakened because the thermally
excited phonon disturbs the spontaneous displacement of the Cu2+ ion. In other
words, the singlet ground state and triplon excitation is no longer a good
picture for the low-energy excitation, and the triplon lifetime shrinks toward
higher temperatures and completely disappears at $T_{\mathrm{SP}}$. This
effect has been observed as an increase in the peak width of the triplon
excitation in neutron scattering experiments[29]. We note that, by
phenomenologically incorporating the triplon-phonon scattering term in a
calculation, the $T$ dependence of $\tilde{V}_{\mathrm{SSE}}$ at moderate
temperatures can be explained qualitatively (see Supplementary Note G4 for
details).

###  0.4 Influence of crystal orientation and Zn-doping on the triplon SSE.

Let us turn to the $H$ dependence of $\tilde{V}_{\mathrm{SSE}}$ in CuGeO3/Pt.
In the entire temperature range, as $H$ increases from zero,
$\tilde{V}_{\mathrm{SSE}}$ increases first and takes its maximum at
$\mu_{0}H\sim 2$ T. As $H$ further increases, $\tilde{V}_{\mathrm{SSE}}$ is
then suppressed. As shown theoretically in the following, triplon scattering
by inevitable impurities and defects in the nominally pure CuGeO3 sample can
explain the observed $H$ dependence of $\tilde{V}_{\mathrm{SSE}}$.

To show that the observed voltage is related to the triplon in the SP phase,
we performed another control experiment with $\nabla T$ applied perpendicular
to the spin chain configuration. Figure 3c shows $\tilde{V}_{\mathrm{SSE}}(H)$
for $\nabla T\perp c$ in a reference sample prepared with the same procedure
as $\nabla T\parallel c$. The observed SSE signal in $\nabla T\perp c$ is only
$1/7$ in magnitude as compared to the $\nabla T\parallel c$ configuration,
consistent with the highly anisotropic transport of triplons. This result
corroborates that the origin of the observed signal is the 1-D triplon
transport. Moreover, additional reference experiments are performed for two
nonmagnetic Zn-doped CuGeO3 samples, in which the triplon spin current should
be suppressed by Zn-induced scattering. Figures 3d and 3e show the
$\tilde{V}_{\mathrm{SSE}}(H)$ for Cu0.99Zn0.01GeO3/Pt and Cu0.97Zn0.03GeO3/Pt
in the same temperature range as the CuGeO3/Pt measurements. Clearly,
$\tilde{V}_{\mathrm{SSE}}$ is suppressed in the $1\%$ Zn-doped sample and
disappeared in the $3\%$ Zn-doped sample.

Figures 3f and 3g show $\chi$ for Cu0.99Zn0.01GeO3 and Cu0.97Zn0.03GeO3,
respectively. For Cu0.99Zn0.01GeO3, the SP transition is observed at 13.2 K.
The decrease of $T_{\mathrm{SP}}$ in doped CuGeO3 agrees with the result from
the previous reports[30]. On the other hand, the SP transition is not visible
in $\chi(T)$ in Cu0.97Zn0.03GeO3. Instead, a clear cusp is observed in
$\chi_{c}(T)$ (with $H\parallel c$-axis) at $T_{\mathrm{N}}\sim 4.3$ K,
indicating phase transition to an antiferromagnetic phase. The
$T_{\mathrm{N}}$ is consistent with the previously reported value for Zn-$3\%$
doped CuGeO3 (ref. [30]). The doping effect implies that non-magnetic Zn
breaks the Cu-spin chains and introduces unpaired Cu free spins, called
solitons, into the spin chain. Solitons disturb the dimerization, and local
antiferromagnetic correlation may lead to long-range antiferromagnetic
order[31] at low temperatures in Cu0.97Zn0.03GeO3. As shown in Fig. 3h,
$\tilde{V}_{\mathrm{SSE}}$ is significantly suppressed in the $\nabla T\perp$
spin-chain configuration and the Zn-doped CuGeO3 samples in all $T$ range. In
the doped samples, the spin-excitation gap is still present in the SP phase,
but their transport may be blocked by Zn. The result suggests that the triplon
spin current is strongly influenced by Zn-induced scattering. Curie-Weiss
fitting of $\chi(T)$ at low $T$ ($<3$ K)[13, 32](see Supplementary Note F for
details) reveals that the density of free $S=1/2$ spins is estimated to be
$0.02\%$ in the nominally pure CuGeO3 sample, so the effect of triplon
scattering by impurities cannot be ignored. Therefore, we attempted to
formulate the triplon spin current in the manner of Boltzmann equation, which
is capable of incorporating scattering effects.

###  0.5 Theoretical model for the triplon SSE.

Figure 4 compares the $H$ dependence of the observed SSE voltage signal with
theoretical results. As we already discussed, $|\tilde{V}_{\mathrm{SSE}}|$ is
expected to monotonically grow with increasing $H$ because the triplon density
of $S_{z}=1$, i.e., the dominant carrier of spin current, increases with
lowering the spin gap up to the critical field $H_{\mathrm{m}}$. However, the
observed broad peak (see Fig. 4a) and the Zn-doped results shown in Fig. 3
strongly suggest the importance of triplon scattering processes. To
theoretically analyze triplon currents in $\rm CuGeO_{3}$, we here apply the
Boltzmann equation[33]. The spin current $J_{s}$ computed by the Boltzmann
equation (See Supplementary Note G2) is given by

$J_{s}=\sum_{S_{z}}\int\frac{k}{2\pi}\hbar
S_{z}(v_{S_{z}}(k))^{2}\tau_{k}k_{\mathrm{B}}\Delta
T\frac{\epsilon_{S_{z}}(k)}{(k_{\mathrm{B}}T)^{2}}\frac{\exp(\epsilon_{S_{z}}(k)/k_{\mathrm{B}}T)}{(\exp\left[\epsilon_{S_{z}}(k)/k_{\mathrm{B}}T\right]-1)^{2}},$
(1)

where $\Delta T$ is the temperature difference along the spin chain direction.
$v_{S_{z}}(k)$, $\tau_{k}$ and $\epsilon_{S_{z}}(k)$ are the group velocity,
relaxation time, and energy dispersion for the triplons[17], respectively. We
assume that the triplon can be viewed as a bosonic particle and the final part
of Equation (1) stems from the Bose distribution function. The $S_{z}=0$
branch of the triplon does not contribute to the spin current. At $H=0$, the
distribution functions of $S_{z}=\pm 1$ branches are identical, which cancel
out the total spin current. By varying $T$ and $H$, the thermal excitation and
the Zeeman effect alter the triplon distributions. Assuming a constant
$\tau_{k}$, the population difference between $S_{z}=\pm 1$ increases
monotonically at a fixed $T(<T_{\rm SP})$ due to the Zeeman effect, resulting
in an almost $H$-linear increase in $J_{s}$. The $H$-linear behavior can also
be reproduced by using a tunnel spin-current formula,[34, 10, 28, 35] which is
another microscopic theory for thermal spin current (see Supplementary Note
G3).

As shown in Fig. 4a, $|\tilde{V}_{\mathrm{SSE}}(H)|$ takes a maximum at around
3 T. To explain the non-monotonic behavior of $\tilde{V}_{\mathrm{SSE}}$, we
considered two mechanisms of triplon scattering processes. One is the
scattering by nonmagnetic impurities, where the scattering probability
$1/\tau_{k,\mathrm{non-mag}}$ is constant with respect to $H$. Another source
of scattering is free spins of broken dimers and unpaired Cu2+ (S = $1/2$).
The density of scatterers and scattering potential are assumed to depend on
$H$ as a Brillouin function. The total scattering probability is then the sum
of the two scattering mechanisms,
$1/\tau_{k,\mathrm{total}}=1/\tau_{k,\mathrm{non-
mag}}+1/\tau_{k,\mathrm{mag}}$. Assuming that the scattering effect by
magnetic impurities is dominant, we can reproduce the experimentally obtained
$H$ dependence of $\tilde{V}_{\mathrm{SSE}}$ shown in Fig. 4b. In the low $H$
regime, where the spin of magnetic impurities is not aligned and the
scattering effect is weak, the Zeeman effect dominates the generation and
transport of the triplon current. With increasing $H$, the triplon mode of
$S_{z}=+1$ ($S_{z}=-1$) shifts downward (upward) due to the Zeeman effect, and
the imbalance between $S_{z}=\pm 1$ triplons becomes larger, resulting in an
increase in the total triplon current. As $H$ is further increased, the
impurities are magnetized and magnetic scattering dominates the triplon
transport, resulting in a broad peak as shown in Fig. 4b. The result agrees
with the experimentally obtained data shown in Fig. 4a in a semi-quantitative
level.

## Discussion

Finally, we here examine other possible origins for the observed
thermoelectric voltage. One candidate could be the paramagnetic SSE due to
free Cu spins. Several papers have addressed the spin current transport in
paramagnets, including Gd3Ga5O12 (refs. [9, 36]), DyScO3 (ref. [9]) and
La2NiMnO6 (ref. [37]). Paramagnon is believed to be responsible for spin
transport in paramagnets as a result of short-range magnetic correlation or
long-range dipole interactions. The spin diffusion length of Gd3Ga5O12 is
estimated to be about 1.8 $\mathrm{\mu m}$ at $5$ K[36]. However, in our
CuGeO3 sample with the free spin density of $0.02\%$, the average distance
between free spins is estimated to be around $0.29$ nm$/0.02\%\sim\ 1.5\
\mathrm{\mu m}$ (the distance between neighboring Cu2+ ions along the $c$-axis
is $0.29$ nm[38]). This means that the paramagnetic spins are too dilute for
spin correlation and spin current transport. Furthermore, the suppression of
$V_{\mathrm{SSE}}$ in the Zn-doped samples, which exhibit larger paramagnetic
moments than CuGeO3, also rules out the possibility of paramagnetic SSE.

The normal and anomalous Nernst effects may also contribute to an electric
voltage under a temperature gradient, but these effects contradict with the
suppression of the signal in the high temperature spin-liquid phase; the
vanishing signals in the $\nabla T\perp c$ configuration and the Zn-doped
samples (see Figs. 3a-3e). Furthermore, the normal Nernst effect should be
linear with respect to the applied magnetic field and can not explain the
observed magnetic field dependence of $V_{\mathrm{SSE}}$ shown in Figs. 3a and
3b.

Note that the voltage observed at 15 K in the CuGeO3/Pt shows almost no
$H$-dependence, as shown in Fig. 2c. Above $T_{\mathrm{SP}}$, the elementary
excitation of spin chains are gapless spinons[21]. A spinon spin current, a
spin current carried by spinons, was found in Sr2CuO3 (ref. [10]). The spinon
spin current may contribute to a small $H$-linear voltage signal above
$T_{\mathrm{SP}}$ (see Supplementary Note C for details). In the spin-Peierls
phase, the gapless spinon is replaced by the gapped triplon[21].

In summary, we observed the triplon spin-Seebeck effect in CuGeO3/Pt. Due to
the triplon excitation in the system, the sign of the observed SSE is opposite
to that of the conventional magnon SSE. The triplon SSE signal persists down
to $2$ K in the spin-Peierls phase and the magnetic field dependence of the
triplon SSE is consistent with microscopic calculation. Our result shows that
the spin-Seebeck effect also acts as a probe for spin excitations in gapped
spin systems, and can also be applied to other materials with exotic spin
excitation, such as a spin ladder system SrCu2O3 (ref. [39]) and a spin dimer
system SrCu2(BO3)2 (ref. [40]).

###  0.6 Sample fabrication.

CuGeO3 single crystals were grown by a traveling solvent floating zone (TSFZ)
method. Samples were cut into size of $\sim$ 7 mm $\times$ 3 mm for further
measurements. 40nm-thick ferrimagnetic insulator Y3Fe5O12 single crystals were
grown on Gd3Ga5O${}_{12}(111)$ substrates by magnetron sputtering. Prior to
the deposition, Gd3Ga5O12 substrates were annealed in the air at 825 ∘C for 30
min in a face-to-face configuration. To crystalize the as-grown amorphous
Y3Fe5O12, the samples were post-annealed in air at 825 ∘C for 200 s. For SSE
measurements, on-chip devices with Pt detector and SiO2/Au heater were
fabricated on top of samples (both pure/doped CuGeO3 and Y3Fe5O12) by an
electron-beam lithography, magnetron sputtering and lift-off technique (see
Supplementary Note A).

###  0.7  SSE measurement.

We adapted the on-chip heating method for SSE, as illustrated in Fig.1b. The
resistance between the Au and Pt layers is much greater than
$\mathrm{M\Omega}$ at room temperature. A longitudinal temperature gradient
was generated by applying an a.c. current with the frequency of $f=13$ Hz to
the Au heater. Since the SSE is driven by Joule heating, the frequency of the
SSE signal on the Pt wire is $2f$. The SSE signal was detected via a lock-in
method.

###  0.8 Magnetization measurement.

The magnetization of CuGeO3 was measured in a Quantum Design Physical
Properties Measurement System (PPMS) with a vibrating sample magnetometer
(VSM) option.

This is a post-peer-review, pre-copyedit version of an article published in
Nature Communications. The final authenticated version is available online at:
https://doi.org/10.1038/s41467-021-25494-7. We thank D. Hirobe, T. Kikkawa, G.
E. W. Bauer and Y. Hirayama for valuable discussions. This research was
supported by JST ERATO Spin Quantum Rectification Project (No. JPMJER1402),
JPSP KAKENHI ( the Grant-in-Aid for Scientific Research (S) No. JP19H05600 and
JP17H06137; the Grant-in-Aid for Scientific Research (A) No. JP16H02125; the
Grant-in-Aid for Scientific Research (B) No. JP20H01830 and No. JP19H02424;
the Grant-in-Aid for Scientific Research (C) No. JP17K05513; the Grant-in-Aid
for Challenging Research (Exploratory) No. JP19K22124 and No. JP20K20896; the
Grant-in-Aid for Scientific Research on Innovative Areas (Research in a
proposed research area) No. JP20H04631, No. JP20H05153, No. JP19H05825 and No.
JP19H04683; the Grant-in-Aid for Research Activity start-up No. JP20K22476),
JST CREST (No. JPMJCR20C1 and No. JPMJCR20T2) and MEXT (Innovative Area “Nano
Spin Conversion Science” (No. 26103005)). Y.C., Y.T. and K.O. are supported by
GP-Spin at Tohoku University. Y.C. is also supported by the Japan Society for
Promotion of Science through a research fellowship for young scientists (No.
JP18J21304).

Y.C. designed the experiments in discuss with Y.S. and E.S. Y.T., Y.N., T.M.
and M.F. grew single crystals for this study. Y.C. and K.O. nano-fabricated
devices. Y.C. collected and analyzed the experimental data. M.S. preformed the
theoretical calculations. E.S. supervised this work. Y.C., M.S., Y.S. and E.S.
wrote the manuscript. All authors discussed the results and commented on the
manuscript.

The authors declare that they have no competing interests.

The data that support the findings of this study are available from the
corresponding authors upon reasonable request.

The codes used in theoretical simulations and calculations are available from
the corresponding authors upon reasonable request.

Correspondence and requests for materials should be addressed to Y.C (e-mail:
<EMAIL_ADDRESS>and M.S (e-mail:masahiro.sato.phys@vc.ibaraki.ac.jp).

Figure 1: Spin chain and spin-Peierls (SP) phase of CuGeO3. a, Crystalline
structure of CuGeO3. The Cu spin chains are shown as dotted lines. Spin-$1/2$
spin chains along the $c$-axis in CuGeO3 are formed by staking CuO2 chains.
$J_{c}$ and $J_{b}$ denote the exchange interactions along the $c$ and
$b$-axis, respectively. b, When temperature ($T$) falls down to the SP
transition temperature $T_{\mathrm{SP}}$, CuGeO3 undergoes a SP transition.
After the transition, the Cu spins dimerize and the ground state becomes a
spin singlet state
($\Bigl{(}\left|\uparrow\downarrow\right\rangle-\left|\downarrow\uparrow\right\rangle\Bigr{)}/\sqrt{2}$
with $S=0$). Elementary excitation in the SP phase is a triplon with $S=1$. c,
Schematic band structure of triplon in CuGeO3[19, 21]. Three triplet states
($\left|\downarrow\downarrow\right\rangle$,
$\Bigl{(}\left|\uparrow\downarrow\right\rangle+\left|\downarrow\uparrow\right\rangle\Bigr{)}/\sqrt{2}$
and $\left|\uparrow\uparrow\right\rangle$) are threefold degenerated at zero
field. By applying a magnetic field ($H$), the Zeeman effect lift the
degeneracy. d, The ground state and excitation state of a ferromagnet.
Elementary excitation is a magnon with $S=1$. The magnons reduce the
magnetization ($\mathbf{M}$) along the magnetic field ($\mathbf{H}$) as
$\delta S^{z}$. e, $T$ dependence of magnetic susceptibility ($\chi$) of
CuGeO3 under an applied magnetic field of $1$ T along the $b$-axis. The SP
transition temperature $T_{\mathrm{SP}}$ is illustrated as the arrow in the
figure. f, Magnetic field - temperature ($H$ \- $T$) spin phase diagram for
the CuGeO3 sample obtained from the $M(H)$ and $\chi(T)$ measurement results.
The diamond symbols represent data taken from a previous study[24]. The dotted
curves are guide for the eyes.  Figure 2:  SSE measurement results for
CuGeO3/Pt. a, Experimental setup for the SSE measurement. Temperature gradient
($\nabla T$) is applied by using an Au heater electrically insulated from the
Pt layer by a SiO2 film. An a.c. current is applied to the Au heater and SSE
voltage in the Pt is measured with a lock-in amplifier as the second harmonic
voltage $V_{2f}$. b, An optical micrograph for a on-chip SSE device made on
the top of CuGeO3. c and d, The $H$ dependence of $\tilde{V}_{\mathrm{SSE}}$
(the SSE signal normalized to the heating power and the detector resistance)
for the CuGeO3/Pt sample at $15$ K and $2$ K, respectively. e and f,
$\tilde{V}_{\mathrm{SSE}}(H)$ measured for Y3Fe5O12/Pt at 2 K and 15 K.
Figure 3:  $\mathbfit{T}$ dependence of $\mathbfit{V}_{\mathbf{SSE}}$. a and
b, $H$ dependence of $\tilde{V}_{\mathrm{SSE}}$ at selected temperatures for
CuGeO3/Pt, with $\nabla T$ applied along the spin chain. c, $H$ dependence of
$\tilde{V}_{\mathrm{SSE}}$ for CuGeO3/Pt, with $\nabla T$ applied
perpendicular to the spin chain. d and e, $H$ dependence of
$\tilde{V}_{\mathrm{SSE}}$ at selected temperatures for Cu0.99Zn0.01GeO3/Pt
and Cu0.97Zn0.03GeO3/Pt, respectively. In a-e, the SSE signal is normalized by
the heating power and the detector resistance. f and g, $T$ dependence of
magnetic susceptibility for Cu0.99Zn0.01GeO3 and Cu0.97Zn0.03GeO3,
respectively. h, $T$ dependence of SSE signal in different samples at
$\mu_{0}H=2$ T. Error bars represent standard deviation.  Figure 4: Magnetic
field dependence of SSE voltage.  a, Magnetic field ($\mu_{0}H$) dependence of
the magnitude of $|\tilde{V}_{\mathrm{SSE}}|$ at $2.5$ K. b, Calculation
results of the $\mu_{0}H$ dependence of the triplon spin current (see text and
Supplementary Note G2 for details).

## References

  * [1] Uchida, K. _et al._ Observation of the spin Seebeck effect.  _Nature_ 455, 778–781 (2008).
  * [2] Takahashi, S., Saitoh, E. & Maekawa, S.  Spin current through a normal-metal/insulating-ferromagnet junction.  _Journal of Physics: Conference Series_ 200, 062030 (2010).
  * [3] Saitoh, E., Ueda, M., Miyajima, H. & Tatara, G.  Conversion of spin current into charge current at room temperature: Inverse spin-Hall effect.  _Applied Physics Letters_ 88, 182509 (2006).
  * [4] Valenzuela, S. O. & Tinkham, M.  Direct electronic measurement of the spin Hall effect.  _Nature_ 442, 176–179 (2006).
  * [5] Ito, N. _et al._ Spin Seebeck effect in the layered ferromagnetic insulators CrSiTe3 and CrGeTe3.  _Physical Review B_ 100, 060402 (2019).
  * [6] Cornelissen, L. J., Liu, J., Duine, R. A., Youssef, J. B. & van Wees, B. J.  Long-distance transport of magnon spin information in a magnetic insulator at room temperature.  _Nature Physics_ 11, 1022–1026 (2015).
  * [7] Seki, S. _et al._ Thermal Generation of Spin Current in an Antiferromagnet.  _Phys. Rev. Lett._ 115, 266601 (2015).
  * [8] Wu, S. M. _et al._ Antiferromagnetic Spin Seebeck Effect.  _Phys. Rev. Lett._ 116, 097204 (2016).
  * [9] Wu, S. M., Pearson, J. E. & Bhattacharya, A.  Paramagnetic Spin Seebeck Effect.  _Physical Review Letters_ 114, 186602 (2015).
  * [10] Hirobe, D. _et al._ One-dimensional spinon spin currents.  _Nature Physics_ 13, 30–34 (2017).
  * [11] Hirobe, D., Kawamata, T., Oyanagi, K., Koike, Y. & Saitoh, E.  Generation of spin currents from one-dimensional quantum spin liquid.  _Journal of Applied Physics_ 123, 123903 (2018).
  * [12] Pytte, E.  Peierls instability in Heisenberg chains.  _Physical Review B_ 10, 4637–4642 (1974).
  * [13] Hase, M., Terasaki, I. & Uchinokura, K.  Observation of the spin-Peierls transition in linear Cu2+ ( Spin-$\frac{1}{2}$) chians in an inorganic compound CuGeO3 .  _Physical Review Letters_ 70, 3651–3654 (1993).
  * [14] Uchinokura, K.  Topical review: spin-Peierls transition in CuGeO3 and impurity-induced ordered phases in low-dimensional spin-gap systems.  _Journal of Physics: Condensed Matter_ 14, R195–R237 (2002).
  * [15] Nishi, M., Fujita, O. & Akimitsu, J.  Neutron-scattering study on the spin-Peierls transition in a quasi-one-dimensional magnet CuGeO3.  _Physical Review B_ 50, 6508–6510 (1994).
  * [16] Lorenzo, J. E. _et al._ Soft longitudinal modes in spin-singlet CuGeO3.  _Physical Review B_ 50, 1278–1281 (1994).
  * [17] Takayoshi, S. & Sato, M.  Coefficients of bosonized dimer operators in spin-$\frac{1}{2}$ XXZ chains and their applications.  _Physical Review B_ 82, 214420 (2010).
  * [18] Brill, T. M. _et al._ High-Field Electron Spin Resonance and Magnetization in the Dimerized Phase of CuGeO3.  _Physical Review Letters_ 73, 1545 (1994).
  * [19] Regnault, L. P., Ain, M., Hennion, B., Dhalenne, G. & Revcolevschi, A.  Inelastic-neutron-scattering investigation of the spin-Peierls system CuGeO3.  _Physical Review B_ 53, 5579–5597 (1996).
  * [20] Fujita, O., Akimitsu, J., Nishi, M. & Kakurai, K.  Evidence for a Singlet-Triplet Transition in Spin-Peierls System CuGeO3.  _Physical Review Letters_ 74, 1677 (1995).
  * [21] Bonner, J. C. & Blöte, H. W. J.  Excitation spectra of the linear alternating antiferromagnet.  _Physical Review B_ 25, 6959–6980 (1982).
  * [22] Rønnow, H. M. _et al._ Neutron Scattering Study of the Field-Induced Soliton Lattice in CuGeO3.  _Physical Review Letters_ 84, 4469 (2000).
  * [23] Hase, M. _et al._ Magnetization of pure and Zn-doped spin-Peierls cuprate CuGeO3 in high magnetic field.  _Physica B: Condensed Matter_ 201, 167–170 (1994).
  * [24] Uchinokura, K., Hase, M. & Sasago, Y.  Magnetic phase transitions in CuGeO3 in high magnetic fields.  _Physica B: Condensed Matter_ 211, 175–179 (1995).
  * [25] Maekawa, S., Valenzuela, S., Saitoh, E. & Kimura, T.  _Spin Current_ (OUP Oxford, Oxford, 2012).
  * [26] Kikkawa, T. _et al._ Critical suppression of spin Seebeck effect by magnetic fields.  _Physical Review B_ 92, 064413 (2015).
  * [27] Uchida, K.-i. _et al._ Observation of longitudinal spin-Seebeck effect in magnetic insulators.  _Applied Physics Letters_ 97, 172505 (2010).
  * [28] Hirobe, D. _et al._ Magnon Pairs and Spin-Nematic Correlation in the Spin Seebeck Effect.  _Physical Review Letters_ 123, 117202 (2019).
  * [29] Lussier, J. G., Coad, S. M., McMorrow, D. F. & Paul, D. M.  The temperature dependence of the spin - Peierls energy gap in CuGeO3.  _Journal of Physics: Condensed Matter_ 8, L59–L64 (1996).
  * [30] Sasago, Y. _et al._ New phase diagram of Zn-doped CuGeO3.  _Physical Review B_ 54, R6835–R6837 (1996).
  * [31] Fukuyama, H., Tanimoto, T. & Saito, M.  Antiferromagnetic Long Range Order in Disordered Spin-Peierls Systems.  _Journal of the Physical Society of Japan_ 65, 1182–1185 (1996).
  * [32] Grenier, B. _et al._ Magnetic susceptibility and phase diagram of CuGe1-xSixO3 single crystals.  _Physical Review B_ 57, 3444–3453 (1998).
  * [33] Abrikosov, A.  _Fundamentals of the Theory of Metals_ (Dover Publications, Mineola, New York, 2017).
  * [34] Adachi, H., Ohe, J.-i., Takahashi, S. & Maekawa, S.  Linear-response theory of spin Seebeck effect in ferromagnetic insulators.  _Physical Review B_ 83, 533 (2011).
  * [35] Jauho, A.-P., Wingreen, N. S. & Meir, Y.  Time-dependent transport in interacting and noninteracting resonant-tunneling systems.  _Physical Review B_ 50, 5528 (1994).
  * [36] Oyanagi, K. _et al._ Spin transport in insulators without exchange stiffness.  _Nature Communications_ 10, 4740 (2019).
  * [37] Shiomi, Y. & Saitoh, E.  Paramagnetic Spin Pumping.  _Physical Review Letters_ 113, 266602 (2014).
  * [38] Kamimura, O., Terauchi, M., Tanaka, M., Fujita, O. & Akimitsu, J.  Electron Diffraction Study of an Inorganic Spin-Peierls System CuGeO3.  _Journal of the Physical Society of Japan_ 63, 2467–2471 (1994).
  * [39] Azuma, M., Hiroi, Z., Takano, M., Ishida, K. & Kitaoka, Y.  Observation of a Spin Gap in SrCu2O3 Comprising Spin-½ Quasi-1D Two-Leg Ladders.  _Physical Review Letters_ 73, 3463–3466 (1994).
  * [40] Kageyama, H. _et al._ Exact Dimer Ground State and Quantized Magnetization Plateaus in the Two-Dimensional Spin System SrCu2(BO3)2.  _Physical Review Letters_ 82, 3168–3171 (1999).

See pages - of Supp.pdf


# Revisiting visualization of spiral states in a wide-gap spherical Couette
flow

Isshin Arai  Department of Pure and Applied Physics, Faculty of Engineering
Science, Kansai University, Osaka, 564-8680, Japan  Tomoaki Itano
<EMAIL_ADDRESS>Department of Pure and Applied Physics, Faculty of
Engineering Science, Kansai University, Osaka, 564-8680, Japan  Masako
Sugihara-Seki  Department of Pure and Applied Physics, Faculty of Engineering
Science, Kansai University, Osaka, 564-8680, Japan

###### Abstract

A pioneering study conducted by Egbers and Rath (1995) [Acta Mech. 111 pp.
125–140 (1995)] experimentally captured spiral waves to elucidate the
transition in the wide-gap spherical Couette flow. However, the physical field
quantities of the spiral waves corresponding to light patterns of various
intensities, as obtained in the experiment, remain unclear, and we have yet to
move beyond the understanding that the reflected light from shear-sensitive
flake tracers responds to a flow that appears at the transition. In this
study, the experiment to visualize spiral waves using aluminum flakes, as
performed by Egbers and Rath, was numerically reproduced by solving the
translational and rotational motions of the particles in a spiral wave. First,
the spiral wave in a spherical Couette flow with an aspect ratio $\eta=1/2$
was numerically calculated using the Newton–Raphson method. Subsequently, the
image that was numerically reproduced from the spiral wave was compared with
an experimentally visualized image. The torque acting on the inner sphere and
the phase angular velocity of the spiral waves with various wavenumbers were
provided. Attempts have been made to determine the instantaneous physical
quantity to which the light and dark patterns obtained in the visualization
corresponded, and the orientation motion of the flakes developed in the
advective history of the flow is essential to yield favorable results.
Exploring the correlation between flow visualization results and shear
structures may provide a new avenue for quantitatively estimating spatial
structures and time scales in complex and quickly time-varying flow fields,
such as turbulence.

Spherical Couette flow, spiral wave, flake visualization

††preprint: APS/123-QED

##  I introduction

The fluid phenomena associated with astrophysical bodies, such as those in
their cores, atmospheres, and oceans, involve flows among concentric double
spheres[1, 2]. Although considering the thermohydrodynamic or electromagnetic
factors for the flow in the core of an astrophysical body is important[3, 4,
5, 6, 7, 8], even considering only the angular velocity difference among
spherical boundaries complicates the flow. We considered a system involving
the Newtonian fluid flow confined in the gap between spheres of radii
$\tilde{r}_{\rm in}$ and $\tilde{r}_{\rm out}$, where the outer spherical
boundary is at rest and the inner spherical boundary rotates at a constant
angular velocity $\tilde{\Omega}_{\rm in}$. This system is referred to as the
spherical Couette flow  (SCF)[9, 10, 11, 12, 13, 14, 15].

The SCF is determined based on two dimensionless parameters, that is, the
aspect ratio $\eta=\tilde{r}_{\rm in}/{\tilde{r}}_{\rm out}$ and Reynolds
number $\mbox{{\rm Re}}=\tilde{r}_{\rm in}^{2}\tilde{\Omega}_{\rm
in}/\tilde{\nu}$, where $\tilde{\nu}$ is the kinematic viscosity of the
working fluids. Egbers and Rath[11] obtained a phase diagram in the
$\eta-\mbox{{\rm Re}}$ parameter space based on the assumption that regardless
of the initial conditions, these parameters only determine the final
equilibrium state via nonlinearity and dissipation. Realized below the
critical Reynolds number determined for each value of $\eta$, the axisymmetric
fundamental flow is not considered a laminar flow because the streamlines are
not parallel to the spherical boundaries. Both spherical boundaries must be
non-slipping to generate the axisymmetric basic flow because the rigidly
rotating flow is realized if the slip condition is imposed on either spherical
boundary. For small values of $\eta$ ($\eta\lesssim 3/4$), also referred to as
the “wide gap case”[12], the flow becomes unstable because of the crossing
fundamental flow with inflection points from mid-latitudes to the poles above
the finite critical Reynolds number, leading to a spiral wave, as described in
[16, 17, 18, 19].

The geomagnetic field on a planet is formed by convection currents in the
liquid core. Considering the planets in the Solar System that form the
geomagnetic field, the radius ratio of the inner core to the outer core varies
from 0.35 for Earth to 0.55 for Mercury. In this study, we focus on a wide
spherical SCF with a radius ratio $\eta=1/2$. Belyaef et al. [10, 20]
estimated the first transitional Reynolds number for $\eta=1/2$ using power
spectra obtained from laser Doppler velocimetry. A recent study [21, 22]
showed that the basic flow becomes unstable against infinitesimal sinusoidal
disturbances for an azimuthal wavenumber of 4 at $\mbox{{\rm Re}}=489$.
Immediately after Re exceeds the critical Reynolds number at which a spiral
wave with a wavenumber of 4 is expected to be sustained, a spiral wave with a
lower wavenumber, such as 3 or 2, appears dominant as the Reynolds number
increases, which has been observed in experiments[23, 24]. The classical
visualization technique can visualize the flow state in a spiral wave as a
spiral pattern[11, 25] with $m$ equally spaced arms extending from the poles
to the equatorial zone in each hemisphere. This phenomenon is also observed as
a polygonal pattern in the flow on a rotating planar disk in a stationary
casing[26]. Belyaef et al. [10, 20] reported hysteresis loops related to
spiral waves with various wavenumbers in multiple ranges of Re. The numerical
experiments based on the shooting method indicate that both spiral waves with
wavenumbers 3 and 4 possess independent basin attractions in the phase
space[27]. However, to date, the reasons why experiments predominantly produce
spiral waves with smaller wavenumber remain unclear.

In the 1970s, when SCF research began to flourish, the transitions in SCF were
measured by visualization technique by the addition of aluminum to the
fluid[23]. We revisit the pioneering study by Egbers and Rath (1995), wherein
spiral waves were experimentally captured. They suspended a small amount of
aluminum flakes in a working fluid at a concentration of 0.05% as tracer
particles to visualize the flow structure. In their experiment, the
viscosities of the working fluids with tracer particles were undetectable in
rheometer measurements. They successfully captured the transition of the
visualization image on the path from the basic state to the spiral wave by
setting a light source and video camera through a half mirror in the polar
direction of the SCF experiment. The change in the captured pattern was not
due to a spatial fluctuation in flake density. This is because, as long as the
particles that are neutrally suspended in the fluid are sufficiently small,
the number density of the particles in the flow must be uniform, regardless of
whether they are spherical or flake-shaped. However, in principle, flakes
rotate under the shear of the flow such that the flow structure as a vorticity
is educed due to the orientation.

However, the physical field quantities of the spiral wave to which the bright
and dark patterns correspond remain unclear. When light is passed through the
suspension with flakes, the intensity of the scattered light is very sensitive
to the orientation of the flakes, indicating that even small velocity
gradients can be measured[23]. Thus, flow velocity gradients may be measured
from light intensity using visualization of flake-mixed fluid experiments;
thus, can we conclude that the spatial pattern of light intensity is solely
determined by the instantaneous velocity gradient? In this study, we virtually
distributed infinitesimal flake particles with a uniform density in an exact
field of spiral waves, which were numerically solved using spectral expansion
and the Newton–Raphson method[27]. By solving the translational and rotational
motions of the particles in the spiral wave using the procedure presented by
Goto et al. [28], we virtually reproduced the light and dark patterns of the
flow field photographed in the experiment by Egbers and Rath[11]. Furthermore,
we discerned the physical quantities to which the pattern of brightness
corresponds.

The remainder of this paper is organized as follows: The next section
describes the numerical method for solving spiral waves via the Newton–Raphson
method. To visualize the spiral state, the equations for both translational
and rotational motions of an infinitesimal drifting planar particle are
presented in a flow and numerical procedure. In Section 3, we compare the
images of the spiral state visualized experimentally in [11] and numerically
using our approach. In Section 4, we discuss the physical quantities obtained
from images using the classical visualization technique. Finally, Section 5
summarizes this study.

##  II Numerical Method

Under the nondimensionalization based on the half length of the gap width,
$\Delta\tilde{r}$, and viscosity, we assume that an inner sphere centered at
the origin rotates at a constant angular velocity $\Omega_{\rm
in}=\frac{\tilde{\nu}}{\Delta\tilde{r}^{2}}\tilde{\Omega}_{\rm in}$ about the
z-axis, and we consider the polar coordinates $(r,\theta,\phi)$. The
nondimensionalized momentum equation of an incompressible Newtonian fluid
confined in the gap $2\Delta\tilde{r}$ between the inner and outer spheres is
given by

$\frac{\partial{\bm{u}}}{\partial
t}+{\bm{u}}\cdot{\bm{\nabla}}{\bm{u}}=-{\bm{\nabla}}p+{\bm{\nabla}}^{2}{\bm{u}}$

under nonslip conditions on both the inner and outer spherical boundaries. The
arbitrary incompressible velocity field ${\bm{u}}$ can be expressed as
${\bm{\nabla}}\times\bigl{(}-{\bm{r}}\varPsi+{\bm{\nabla}}\times({\bm{r}}\varPhi)\bigr{)}$
according to two scalar fields, $\varPsi(r,\theta,\phi)$ and
$\varPhi(r,\theta,\phi)$, commonly known as toroidal–poloidal
decomposition[23]. Each scalar field may be described according to the product
of a series of spherical harmonics in $\theta,\phi$ directions and modified
Chebyshev polynomials in the radial direction. If the solution is a rotating
wave solution with constant phase angular velocity $\omega$, the time-
derivative term in the governing equation $\frac{\partial}{\partial t}$ can be
replaced by $-\omega\frac{\partial}{\partial\phi}$. Thus, all the terms in the
governing equations of $\varPsi$ and $\varPhi$ specifying a spiral state can
be numerically evaluated. The condition satisfied by the toroidal and poloidal
components of a spiral wave is equivalent to an algebraic equation with
degrees of freedom representing its flow state. The algebraic equation was
solved using the Newton–-Raphson method with the aid of open numerical
libraries[29, 30, 31].

When the spiral wave (or secondary wave) was first experimentally discovered
in the pioneering works, it was captured as a periodic light reflection and
was perceived as a wave of disturbance with wavenumber $m$ that could be
superimposed. On the other hand, it is a solution with the $m$-fold symmetry
to the Navier–Stokes equation for us, which cannot be superimposed, and its
amplitude is uniquely determined by Re and $\eta$; thus, we refer to it as
spiral state rather than as a spiral wave.

We considered two arbitrary material lines ${\bm{a}}$ and ${\bm{b}}$ with
infinitesimal length floating at a point in the spiral state. In the flow, the
location, length, and direction change with time $t$. The time evolutions of
${\bm{a}}(t)$ and ${\bm{b}}(t)$ are given by
$\dot{{\bm{a}}}={\bm{a}}\cdot{\bm{\nabla}}{\bm{u}}+{\rm O}(|{\bm{a}}|^{2})$,
where ${\bm{u}}$ is the velocity field at ${\bm{r}}(t)$, which is the position
of the material line at time $t$[32]. The positions of the material lines
follow the order of $\dot{{\bm{r}}}={\bm{u}}({\bm{r}}(t),t)$. Owing to the
flow incompressibility, the time variation of the infinitesimal surface
element formed by these material lines ${\bm{S}}(t)={\bm{a}}\times{\bm{b}}$
satisfies the equation $\dot{{\bm{S}}}=-{\bm{\nabla}}{\bm{u}}\cdot{\bm{S}}$.
However, the areas of the small and thin flake-shaped particles distributed in
the flow do not change over time. The unit normal vector of the flake
${\bm{s}}$ is represented by ${\bm{s}}(t)={\bm{S}}/|{\bm{S}}|$ and
$d|{\bm{S}}|^{2}/dt=-2{\bm{S}}\cdot{\bm{\nabla}}{\bm{u}}\cdot{\bm{S}}$.
Therefore, the orientation of the flake at time $t$ is given by the time
integration of the equation:
$\dot{{\bm{s}}}={\bm{s}}\times{\bm{s}}\times{\bm{\nabla}}({\bm{u}}\cdot{\bm{s}})$,
as shown in [28]. We obtained the orientation of the flakes drifting in the
spiral state of the SCF by numerically solving the equation [24].

For illustration, we consider a flake located in a simple shear flow
${\bm{u}}=y{\bm{e}}_{x}$, where the flake moving at a constant speed along the
x-axis changes its orientation with time. The components of ${\bm{s}}$ are
denoted as $(s_{x},s_{y},s_{z})$, and the time evolution of each component is
described according to the ordinary differential equations
$\dot{s_{x}}=s_{x}^{2}s_{y}$, $\dot{s_{y}}=s_{x}(s_{y}^{2}-1)$, and
$\dot{s_{z}}=s_{x}s_{y}s_{z}$. The system possesses two conserved quantities:
$s_{x}^{2}+s_{y}^{2}+s_{z}^{2}$ and $s_{x}/s_{z}$. The steady solution is
$s_{x}=0$. Under arbitrary initial conditions, ${\bm{s}}$ moves away from
$s_{z}=\pm 1$ and asymptotically approaches $s_{y}=\pm 1$ over time. Thus, the
attractors, repellers, and limit cycles may exist in the phase space of the
flake orientations (unit spherical surface) for an arbitrary steady flow. For
general time-dependent flow fields, invariant sets solved in an instantaneous
flow field may vary with time. This example indicates that even if the
orientation of the flakes initially scattered at a point in the flow is
stochastically isotropic, they may align in a specific direction after a
certain time. Conversely, the orientations of the flakes concentrated in a
specific direction may become isotropic after a certain period.

As incident collimated rays from a laser or other light sources outside the
container pass through the fluid in a spherical Couette flow container, some
rays may be reflected by the flakes. When a ray reflected by the flakes enters
the camera, the image is brighter at the position of the flakes. Based on the
assumption of no light attenuation in the fluid, the condition for the flake
to completely reflect a ray from the light source to the camera was evaluated.
The direction of the incident ray on the flake is represented as a unit vector
${\bm{I}}$, which is independent of the position of the flake. The reflected
light is observable if the orientation of flake ${\bm{s}}$ is in the
bisectional direction between $-{\bm{I}}$ and the optical axis of the camera.
Generally, the direction of a ray ${\bm{R}}$ reflected from the flake is
${\bm{R}}={\bm{I}}-2({\bm{I}}\cdot{\bm{s}}){\bm{s}}$. The alignment of the
reflected rays against the optical axis of the camera must be correlated with
the intensity distribution in the recorded images. Thus, the intensity of the
reflected light captured by the camera may be evaluated using the magnitude of
${\bm{C}}\cdot{\bm{R}}$, where ${\bm{C}}$ is a unit vector along the optical
axis of the camera lens with respect to the reflecting flake. If the flakes
are isotropically oriented in a localized region, the light intensities are
scattered throughout the region such that the region has a gritty pattern. On
the other hand, if the flakes are aligned in a direction in the localized
region, the light intensities captured by a camera in the region remain
consistent regardless of the flakes reflecting light such that the region may
appear with a smooth pattern.

##  III Results

Figure 1: The normalized value of the torque acting on the spherical boundary
against the Reynolds number, which shows a bifurcation diagram of the
axisymmetric and 2-, 3-, and 4-fold states.

For $\eta=0.5$, the 4-fold state bifurcates at $\mbox{{\rm Re}}\sim 489$ from
the axisymmetric state, and a 3-fold state occurs at $491$[27], followed by
the 2-fold state bifurcating around $\mbox{{\rm Re}}\sim 600$. All states
originated in the axisymmetric state, which first lost stability at
$\mbox{{\rm Re}}=489$ against a perturbation with a wavenumber of 4. The
torque values exceeded the Stokes flow by a factor of 3 [33] in the range of
$\mbox{{\rm Re}}>500$. The torque values of the $m$-fold states are smaller
than those of the axisymmetric state, and the torque slightly decreases as the
wavenumber $m$ decreases. This may be attributed to the fact that a spiral
state with a smaller wavenumber is predominantly obtained in experiments at a
high Re. This is contrasting to previous experimental results on narrow-gap
SCF[34], where the frictional torque of states successively bifurcated from
the axisymmetric state is larger than that of the predecessor, such as the
axisymmetric state, as measured in a turbulent transition. Intuition indicates
that the transition from laminar to turbulent flow becomes more complex (as
the wave number increases) with increasing Re. Simultaneously, it is
accompanied by a greater dissipation and friction. However, the spiral state
bifurcation in SCF suggests that this intuition is not always correct.

The phase angular velocities of the spiral states calculated using the
Newton–Raphson method have similar values regardless of the wavenumber and are
found to be in the range of approximately $1/6\leq\omega/\Omega_{\rm in}\leq
1/7$. As with the torque graphs, confirming the difference among the branches
with different wavenumbers in the graphs indicating the relationship between
Re and phase speed is difficult. In addition, the value of the dimensionless
frequency proposed by Belyaef[10] may be utilized to better represent the
difference. This quantity is the frequency of the flickering of the light
reflected from the flow, which can be experimentally measured (normalized by
the rotation rate of the inner sphere).

Figure 2:  Two-, three- and four-fold states and branches using dimensionless
frequency  Figure 3:  The isosurface of $u_{r}$ and $\varPhi$ of the 4-fold
state at $\mbox{{\rm Re}}=600$ was captured from $+z$ direction. The blue and
pink isosurfaces correspond to the region with $u_{r}>1$ and $\varPhi>0.1$,
respectively.

Fig.3 illustrates the overall view of the 4-fold spiral state taken from the
$+z$ direction, where the inner spherical boundary rotates anticlockwise. The
blue isosurface corresponds to the region where $u_{r}>1$ and pink isosurface
corresponds to the region where $\varPhi>0.1$. Both images imply the term
“spiral,” which describes the flow structure consisting of equally spaced arms
extending from the poles to the equatorial zone in each hemisphere. Through an
experimental visualization technique using aluminum flakes drifting through a
horizontal plane illuminated by a laser sheet, Yoshikawa et al.[24] recorded a
spiral wave with a wavenumber of 3 on the horizontal cross-section over the
critical Reynolds number. By solving the equation of motion for infinitesimal
planar particles advecting in the flow field of a spiral wave, they also
virtually reproduced a sequence of visual distributions of reflected light.
This is in good agreement with the pictures recorded experimentally. We
applied the same method to reproduce an image of the 4-fold spiral state
(Fig.4) captured by Egbers and Rath in [11], which was the pioneering
observation to capture the transition from the axisymmetric state in a wide-
gap case.

Figure 4:  Fig. 8(b) in [11] reproduced with permission from Springer Nature,
which is an experimental visual observation at $\mbox{{\rm Re}}=1583$ using
both a white light source and video camera via a half mirror positioned at
$+z$. The captured image was identified as “non-axisymmetric secondary wave”
with four spiral arms for $\eta=2/3$. Suspended aluminum flakes advected in
the working fluid are employed as tracer particles

Fig.4 is a copy of Fig.8 in [11], where the 4-fold spiral state, illuminated
by a light source from the direction of a pole, was captured using a camera
through a half mirror placed at the pole. The incident rays from the polar
direction are partly absorbed in the bulk fluid, partly reflected by the
surfaces of the flakes floating in the bulk fluid, and partly transmitted
through the inner fluid layer. By simplifying the aforementioned complex
process followed by the light rays, we consider the image of the 4-fold spiral
state captured by the camera as a superposition of the rays reflected toward
the optical axis of the camera by the flakes existing on various $z$ planes.

Figure 5:  The numerical visualization of the 4-fold state solved at
$\mbox{{\rm Re}}=600$ was virtually reproduced in the case of
${\bm{I}}=-{\bm{e}}_{z}$ and ${\bm{C}}=+{\bm{e}}_{z}$, which shows a camera
and light source positioned at $+z$ direction. Inner products with values
varying from 0 to 1 are illuminated in a black-to-white grayscale. All initial
orientations of flakes are random with respect to space. From top to bottom
$z=2,3,\sqrt{12}$.

The light intensity of the rays reflected by $32000$ flakes distributed at
different values of $z$ was calculated for a 4-spiral state solved at
$(\eta,\mbox{{\rm Re}})=(1/2,600)$ using the computational visualization
method described in the previous section, where we assumed the direction of
the incident ray and optical axis of the camera as ${\bm{I}}=-{\bm{e}}_{z}$
and ${\bm{C}}=+{\bm{e}}_{z}$, respectively. Figs.5(a), (b) and (c) correspond
to the light intensity reflected at $z=2,3$, and $\sqrt{12}$, respectively.
Flakes were isotropically oriented using random numbers at $3T_{\rm in}$
before the time to drift to each position on each $z$, where $T_{\rm
in}=2\pi/\Omega_{\rm in}$. At the poles, there is a flow from the outer sphere
to the inner sphere in the radial direction due to the Ekman pumping. This
flow is expected to incline the normal vector of the flake equatorward, except
near the inner sphere. The flakes move along the flow while satisfying
${\bm{u}}\cdot{\bm{s}}=0$ for relatively long sections, resulting in the
appearance of a small round shadow near the pole in Figs.5(b) and (c).

Fig.5(a) shows that a square streaky shadow surrounding the inner sphere. The
square pattern appears to be related to the cross-sectional image of the shear
arms extending from the poles in the spiral state, as shown in Fig.3. However,
the physical quantity of the state related to the shadow region cannot be
easily explained, and this will be discussed in the next section. The bright
regions in both $z=3$ and $z=\sqrt{12}$, which exist in the foreground from
the camera side, are superimposed on the square shadow in Fig.5(a). The square
pattern is invisible to the camera located at the pole direction. By contrast,
the horn-shaped shadows extending outward from the four vertices of the square
streaky shadow shown in Fig.5(a) do not disappear even if the bright images of
$z=3$ and $z=\sqrt{12}$ are superimposed. The four horn-shaped shadow patterns
shown in Fig.4 in [11] correspond to those in Fig.5(a). This suggests that the
inner sphere in [11] also rotates anticlockwise, as same as in the case
performed in this study. Moreover, the intensity of light at $z=\sqrt{12}$ is
stronger than that at $z=3$. This suggests that the normal vector of the
flakes is aligned poleward in response to the flow, since the radial component
of the flow velocity is smaller than the other components along streamlines
from the periphery to the pole near the stationary outer sphere.

Notably, Fig.5 shows the cross-section of the flow presented in Fig.3, where
feather-like structures can be observed. Fig.3 shows an anticlockwise
extension from the base of the wings, whereas Fig.5 illustrates a clockwise
extension. This also indicates that the visualization results do not simply
correspond to the flow structure.

##  IV Discussion

Flow visualizations using flake particles as tracers have been performed in
fluid dynamics experiments [35, 36]. If the density ratio of the tracer to the
working fluid is in unity and the Stokes number is sufficiently small, the
distribution of tracers in the domain is expected to always be uniform from a
theoretical viewpoint. Thus, the spatial fluctuation of the brightness
observed experimentally was considered to originate not from the distribution
but from the orientation of the particles. We primarily focused on the
features of the flow field that can be read from the spatial pattern of the
brightness in the spiral state of the SCF. In the following subsections, we
discussed this in detail, especially with respect to SCF.

In [28], the authors conducted flow visualization using flake particles for
the flow in a spherical vessel rotating under a small precession, which seems
to be an almost rigidly rotating flow. They reported that the spatial pattern
of brightness visualized in their experiment originated not from the
fluctuation of the number density of flakes but from the fluctuation of the
orientation of flakes. They successfully reproduced the experimentally
observed pattern by numerically simulating the motion and rotation of flakes
distributed in an almost rigid rotating flow.

Notably, their numerical results imply that the bright regions in a visualized
image are the superposition of those with flake orientation is isotropic and
those with flakes tend to align in the direction in which the incident rays
are reflected into the line of sight. If the flakes are anisotropically
oriented in a specific region, the incident light from a specific direction to
that region is reflected only in a specific direction, and it is rare for the
reflecting rays to meet the optical axis of the camera. In some partially
observed regions, the direction of the reflected light matches the optical
axis. These regions are identified by the camera as strongly shrinking areas
but only in limited areas. By contrast, in the region where the flakes are
isotropically oriented, the incident light is equally reflected in all
directions. Then, a part of the reflected light can be captured by a camera
located in any direction. In their experiment[28], most of the volumetric
regions in the flow where flakes were well oriented must have been observed to
be dark, except for a few exceptional areas (subtle bright sector-shaped
pattern) where the direction of reflection matched the optical axis of the
camera. The regions that are bright due to flake isotropic orientation possess
a planar structure in the bulk fluid, whereas the regions that are bright,
owing to an orientation consistent with the optical axis, possess a volumetric
structure.

Suppose a flake advecting through a flow field is traveling at a phase
velocity. We assume that the time variation of the velocity gradient perceived
by the flake is small from a Lagrangian perspective and the magnitude of
vorticity is negligible compared to that of strain. The Poincare–-Benedixon
theorem[37, 38] guarantees the existence of limit cycles in a dynamic system
in a two-dimensional closed domain if no fixed point exists. The normal vector
of the flake moves only in the unit sphere; thus, the theorem may be directly
applied to the present stationary problems. As the eigenvector of
${\bm{\nabla}}{\bm{u}}$ is constant in this case and the equation for
${\bm{s}}$ is essentially linear, we can deduce that the orientation of the
flake approaches the direction of the eigenvector associated with the smallest
strain eigenvalue.

The time constant required for the orientation to align may be estimated as
the inverse of the absolute of the smallest (negative) real component among
the eigenvalues of the velocity gradient tensor, that is , $1/t_{\rm
align}\sim|{\rm min}({\rm E.V.of}{\bm{\nabla u}})|$. If the eigenvalues
include a few complex conjugates, the argument below holds only when the
remaining eigenvalue is negative. In a region with multiple fine-scale eddies,
the orientation of the favorable eigenvector may significantly change within a
short distance for which flakes advect during the interval required for the
orientation of the flakes to align. Suppose that the eigenvector of the
instantaneous velocity gradient favorable by a flake at ${\bm{x}}$ is
represented by ${\bm{e}}_{1}({\bm{x}})$ and that the flake moves from
${\bm{x}}$ to ${\bm{x}}^{\prime}={\bm{x}}+\int_{0}^{t_{\rm
align}}{\bm{u}}({\bm{x}}(t),t)dt$ in $t_{\rm align}$. In the region where the
condition $|{\bm{e}}_{1}({\bm{x}})\cdot{\bm{e}}_{1}({\bm{x}}^{\prime})|<1$ is
realized, the orientations of the flakes are kept insufficiently aligned along
the short path. By contrast, for our spiral state of the SCF, flakes drift
along spatially smooth but undulated streamlines, where
${\bm{e}}({\bm{x}})\cdot{\bm{e}}({\bm{x}}^{\prime})\sim 1$ may be realized in
parts of the entire domain. The 4-fold spiral state rotates at a slow phase
speed around the rotation axis, and the history experienced by the flakes
along the trajectories, which crosses from a vortical structure to its
neighboring structure in the spiral state, allows them to align in the
preferred direction.

According to Fig.5, the present spiral state of the SCF includes some regions
where the brightness contour appears to be gritty. In such regions, the flakes
are not oriented in a specific direction. However, the background gradation,
which is behind the gritty pattern associated with isotropic distribution,
varies smoothly from white to black in space. Such a pattern of light and dark
suggests that the flake particles are well aligned. Furthermore, this implies
that the 4-fold arm-like shadow patterns obtained in the visualization
experiments by Egbers& Rath[11] were not the regions where the flake
orientation remained isotropically distributed, as identified by the three
sharp circles by Goto et al[28].

Figure 6:  Similar to Fig.5, all initial directions of flakes were set
constant ($\theta=0$).

The same procedure shown in Fig.5(a) was carried out with initially
distributing orientations of the flakes as a constant,
${\bm{s}}={\bm{e}}_{z}$, as illustrated in Fig.6. A comparison between
Figs.5(a) and 6 indicates several flow characteristics in the present 4-fold
spiral state. First, the pattern of smoothly varying background brightness is
consistent in both figures. This suggests that a sufficiently long time had
passed for the initial orientation to not largely affect the visualization
results. Both figures show that the orientation of the flakes generating the
background was sufficiently acclimatized to the 4-fold state after the flakes
experienced the shear of the spiral state. Second, a gritty pattern associated
with the isotropic distribution of the flake orientation remains, as shown in
Fig.5. In the isotropic distribution, the orientation of the flakes is not
fully aligned in some regions, although they are subjected to flow shear. If
the integral time is longer, some regions in gritty may disappear and others
may persist. The degree of gritty depends on its position, and the region in
gritty possesses not a volumetric but planar structure, as shown in Fig.6.
Third, gritty patterns indicating isotropic distribution completely disappear
around the axis (origin), thin layers near the outer spherical boundary, and
the large black square.

A common question that would arise is what brightness distribution represents
in the flow field. As the flake particles traverse their orbits in a
Lagrangian manner, they rotate under the effect of the shear from the flow
field, which varies in time. In principle, predicting some characteristics of
an instantaneous velocity field from a bright image is difficult, even if the
flow is almost steady. Considering this, we investigated the correspondence
between the brightness distribution drawn by the flakes and the instantaneous
flow field.

Figure 7:  Velocity field of 4-fold state at $\mbox{{\rm Re}}=600$ taken at
$z=2$. Arrows indicate velocity fields $(u_{x},u_{y})$, and contour indicates
the sign of $u_{z}$ ($u_{z}>0$ (white)).

Fig.7 is the projection of the vector of flow field at a cross-section $z=2$,
where the shading and contours represent $u_{z}$ and the light colored areas
correspond to $u_{z}>0$. A narrow square shadow tangential to the maximum
diameter of the inner sphere is shown in Fig.5(a), exhibiting the brightness
distribution at the cross-section $z=2$. From a comparison with Fig.7, no
common feature is confirmed at the location of this square. If we consider a
windmill-like shape consisting of a field $u_{z}<0$, the slope against the
$x$-axis is different from the square formed by the brightness distribution in
Fig.5(a). Apparently, velocity and brightness are not directly correlated at
each point.

Figure 8:  The Q criterion of the 4-fold state at $\mbox{{\rm Re}}=600$ taken
at $z=2$. Contour indicates the sign of Q ($Q>0$ (white)).

A correlation between the brightness and velocity gradient is more probable.
The vorticity magnitude may be utilized to identify vortex structures;
however, it does not always yield the desired coherent structures. We focused
on the Q criterion[39], which is a commonly used method for vortex
visualization, to identify the vortices in a flow with a positive second
invariant of the velocity gradient tensor. Fig.8 shows the flow field
characterized in terms of the Q criterion, where the light colored area
surrounded by the curves $Q=0$ are expected to possess a vortical structure.
The shearing of the stretch region between neighboring vortical structures may
align the orientation of the flakes. Because our spiral state does not contain
sharp or fine eddy structures as in the turbulent flows but only vague shear
fluctuations extending from the poles to the equatorial plane in a spiral
manner, the characterization of the flow field based on Q values is not
evident. We observed no significant correlation between the luminance
distribution and Q values, but a correlation is evident from the comparison of
Fig.8 to the velocity field presented in Fig.7. The velocity field and
vortical structures estimated using the Q value are correlated.

Neither the velocity field nor the Q-value could reproduce the luminance
distribution produced by the flakes. The present example is a counterexample,
that is, for a general flow field, the brightness distribution of flakes
obtained in the visualization experiments may represent neither the vortex
structure nor the region of stretching at that instance, even if the intensity
of the scattered light reflected from flakes is very sensitive to small
velocity gradients[23]. . Because we could not correlate it with the
instantaneous field, we proposed a physical quantity that considers history.

Figure 9:  Shearing $\log{p}$ of the flow considering the background. Flakes
in the white regions were advected with a relatively strong shear.

The flakes are sufficiently small and stick to the flow to perform
translational and rotational motions, so that the orientation of the flakes is
strongly correlated with the stretching of the fluid elements. Taking this
into account, suppose a flake was present at ${\bm{x}}={\bm{x}}(t)$ at a time
$t$, then the following positive quantity is considered,
$p^{2}({\bm{x}})=\frac{1}{T}\int_{0}^{-T}\sum_{k=1,2,3}\lambda_{k}^{2}({\bm{x}}(t),t)dt$,
where $\lambda_{k}$ is the eigenvalue of strain tensor at ${\bm{x}}(t)$. The
contour map of $\log{p}$ in Fig.9 provides a pattern similar in Fig.5,albeit
with differences in the magnitude of the values. The value of $p$ takes on
large values near the inner and outer spheres, which correlated with a pattern
with the regions of background gradation varying smoothly in space, as seen in
Fig.5.

##  V Summary

The Newton–Raphson method was used to numerically solve the spiral state in a
spherical Couette flow with an aspect ratio of $\eta=1/2$. By simulating the
translational and rotational motions of an infinitesimally small neutrally
floating flaky particle with negligible inertia moving through the obtained
state, an image obtained in the experimental visualization of the spiral state
using thin aluminum flakes according to the experiment performed by Egbers and
Rath [Acta Mech. 111 pp. 125–140 (1995)] was reproduced. The numerically
reproduced image agreed with the experimentally visualized images. Several
instantaneous physical variables calculated from the flow field of the spiral
state could not reproduce the light and dark patterns owing to the flakes
obtained during visualization. Thus, incorporating the effects of the
advection history of the flakes in the spiral state is necessary.

The flicker contrast of the reflected light might increase with increasing
Reynolds number, even for spiral states with the same wavenumber. Also, as the
degree of flake orientation alignment increases, the intensity of the
reflections will change accordingly. This means that the intensity of flow
shear may be quantitatively estimated from the visualization. Furthermore, the
degree of flake alignment may fall rapidly as the flow transitions to
turbulence consisting of small eddies, and the distance and time scale of the
turbulence may be estimated. These points will be pursued in future studies.

###### Acknowledgements.

The authors would like to thank Mr. F. Goto and Mr. T. Inagaki for their pilot
surveys and Mr. Yoshikawa for his valuable comments on the draft. They also
acknowledge ORDIST in Kansai University and the RIMS Joint Research Activities
in Kyoto University for providing a space for their research and
communication. This work was supported in part by a Grant-in-Aid for
Scientific Research(C) and JSPS KAKENHI Grant Nos.20K04294 and 24K07331. This
study also benefited from the interaction within RISE-2018 No.824022 ATM2BT of
the European Union Horizon 2020-MSCA program, which includes Kansai
University.

## References

  * Fowler [2004] C. M. R. Fowler,  _The Solid Earth_ (Cambridge University Press, 2004).
  * D.A.Godfrey [1988] D.A.Godfrey, A hexagonal feature around saturn’s north pole, Icarus 76, 335 (1988).
  * Chandrasekhar [1961] S. Chandrasekhar,  _Hydrodynamic and Hydromagnetic Stability_ (Clerendon Press, Oxford, 1961).
  * Busse [1975] F. H. Busse, Patterns of convection in spherical shells, J. Fluid Mech. 72(1), 67 (1975).
  * Zebib  _et al._ [1983] A. Zebib, G. Schubert, J. Dein, and R. Pariwal, Character and stability of axisymmetric thermal convection in spheres and spherical shells, Geophys. Astrophys. Fluid Dyn. 23, 1 (1983).
  * Feudel  _et al._ [2011] F. Feudel, K. Bergemann, L. Tuckerman, C. Egbers, B. Futterer, M. Gellert, and R. Hollerbach, Convection patterns in a spherical fluid shell, Phys. Rev. E 83, 046304 (2011).
  * Kida  _et al._ [1997] S. Kida, K. Araki, and H. Kitauchi, Periodic reversals of magnetic field generated by thermal convection in a rotating spherical shell, J. Phys. Soc. Jpn 66(7), 2194 (1997).
  * Sakuraba and Kono [1999] A. Sakuraba and M. Kono, Effect of the inner core on the numerical solution of the magnetohydrodynamic dynamo, Physics of the Earth and Planetary Interiors 1111, 105 (1999).
  * Munson and Menguturk [1975] B. R. Munson and M. Menguturk, Viscous incompressible flow between concentric rotating spheres, J. Fluid Mech. 69, 705 (1975).
  * Belyaev  _et al._ [1984] Y. N. Belyaev, A. A. Monakhov, S. A. Scherbakov, and I. M. Yavorskaya, Some routes to turbulence in spherical couette flow, in  _Laminar-Turbulent Transition IUTAM-Symp._ (1984) pp. 669–676.
  * Egbers and Rath [1995] C. Egbers and H. J. Rath, The existence of taylor vortices and wide-gap instabilities in spherical couette flow, Acta Mechanica 111, 125 (1995).
  * Nakabayashi  _et al._ [2002] K. Nakabayashi, Y. Tsuchida, and Z. Zheng, Characteristics of disturbances in the larminar-turbulent transition of spherical couette flow. 1\. spiral taylor-g ortler vortices and traveling waves for narrow gaps, Phys. Fluids 14, 3963 (2002).
  * Wicht [2014] J. Wicht, Flow instabilities in the wide-gap spherical couette system, J. Fluid Mech. 738, 184 (2014).
  * Barik  _et al._ [2018] A. Barik, S. Triana, M. Hoff, and J. Wicht, Triadic resonances in the wide-gap spherical couette system, J. Fluid Mech. 843, 211 (2018).
  * Hoff and Harlander [2019] M. Hoff and U. Harlander, Stewartson-layer instability in a wide-gap spherical couette experiment: Rossby number dependence, J. Fluid Mech. 878, 522 (2019).
  * Araki  _et al._ [1997] K. Araki, J. Mizushima, and S. Yanase, The nonaxisymmetric instability of the wide-gap spherical couette flow, Phys. Fluids 9, 1197 (1997).
  * Hollerbach  _et al._ [2006] R. Hollerbach, M. Junk, and C. Egbers, Non-axisymmetric instabilities in basic state spherical couette flow, Fluid Dyanmics Research 38, 257 (2006).
  * Abbas  _et al._ [2018a] S. Abbas, L. Yuan, and A. Shah, Simulation of spiral instabilities in wide-gap spherical couette flow, Fluid Dyn. Res. 50, 025507 (2018a).
  * Abbas  _et al._ [2018b] S. Abbas, L. Yuan, and A. Shah, Existence regime of symmetric and asymmetric taylor vortices in wide-gap spherical couette flow, Journal of the Brazilian Society of Mechanical Sciences and Engineering 11, 154 (2018b).
  * Belyaef and Yavorskaya [1991] Y. N. Belyaef and I. M. Yavorskaya, Spherical couette flow: transitions and onset of chaos, (translated from) Izvestiya Akademii Nauk SSSR, Mekh. Zhid. i Gaza 1, 10 (1991).
  * Dumas and Leonard [1994] D. Dumas and A. Leonard, A divergence-free spectral expansions method for three dimensional flows in spherical-gap geometries, J. Comput. Phys. 111, 205 (1994).
  * Junk and Egbers [2000] M. Junk and C. Egbers, Isothermal spherical couette flow, in  _Physics of Rotating Fluids_ , Lecture Notes in Physics, Vol. 549, edited by C. Egbers and G. Pfister (Springer, 2000) pp. 215–233.
  * Belyaev  _et al._ [1978] Y. N. Belyaev, A. A. Monakhov, and I. M. Yavorskaya, Stability of spherical couette flow in thick layers when the inner sphere revolves, Izvestiya Akademii Nauk SSSR, Mekchanika Zhidkosti i Gaza 2, 9 (1978).
  * Yoshikawa  _et al._ [2023] K. Yoshikawa, T. Itano, and M. Sugihara-Seki, Numerical reproduction of the spiral wave visualized experimentally in a wide-gap spherical couette flow, Phys. Fluids 35, https://doi.org/10.1063/5.0144365 (2023).
  * Wulf  _et al._ [1999] P. Wulf, C. Egbers, and H. J. Rath, Routes to chaos in wide-gap spherical couette flow, Phys. Fluids 11, 1359 (1999).
  * Mizushima  _et al._ [2009] J. Mizushima, G. Sugihara, and T. Miura, Two modes of oscillatory instability in the flow between a pair of corotating disks, Phys. Fluids 21, 014101 (2009).
  * Goto  _et al._ [2021] F. Goto, T. Itano, M. Sugihara-Seki, and T. Adachi, Bifurcation aspect of polygonal coherence over transitional reynolds numbers in wide-gap spherical couette flow, Phys. Rev. Fluids 6, 113903 (2021).
  * Goto  _et al._ [2011] S. Goto, S. Kida, and S. Fujiwara, Flow visualization using reflective flakes, J. Fluid Mech. 683, 417 (2011).
  * Schaeffer [2013] N. Schaeffer, Efficient spherical harmonic transforms aimed at pseudospectral numerical simulations, Geochemistry, Geophysics, Geosystems 14, 751 (2013).
  * Frigo and Johnson [2005] M. Frigo and S. G. Johnson, The design and implementation of FFTW3, Proceedings of the IEEE 93, 216 (2005), special issue on “Program Generation, Optimization, and Platf orm Adaptation”.
  * Anderson  _et al._ [1999] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen,  _LAPACK Users’ Guide_ , 3rd ed. (Society for Industrial and Applied Mathematics, Philadelphia, PA, 1999).
  * Batchelor [1967] G. K. Batchelor,  _An introduction to fluid dyanmics_ (cambridge university press, Cambridge, UK, 1967).
  * Landau and Lifshitz [1987] L. Landau and E. Lifshitz,  _Fluid Mechanics: Vol.6 (Course of Theoretical Physics) 2nd Edition_ (Butterworth-Heinemann, 1987).
  * Nakabayashi [1978] K. Nakabayashi, Frictional moment of flow between two concentric spheres, one of which rotates, J. Fluids Engineering 100, 97 (1978).
  * Dyke [1982] M. V. Dyke,  _An Album of Fluid Motion_ (Parabolic press, Stanford, California, 1982).
  * Samimy  _et al._ [2003] M. Samimy, K. Breuer, L. Leal, and P. Steen,  _A Gallery of Fluid Motion_ (Cambridge University press, Cambridge, UK, 2003).
  * Guckenheimer and Holmes [1983] J. Guckenheimer and P. Holmes,  _Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, Applied Mathematical Sciences 42_ (Springer-Verlag, New York, 1983).
  * Hale and Koçak [1991] J. K. Hale and H. Koçak,  _Dyanmics and Bifurcations, Texts in Applied Mathematics 3_ (Springer-Verlag, New York, 1991).
  * Hunt  _et al._ [1988] J. C. R. Hunt, A. A. Wray, and P. Moin, Eddies, streams, and convergence zones in turbulent flows, Center for Turbulence Research CTR-S88, 193 (1988).