#  Critical Slowing Down of Quadrupole and Hexadecapole Orderings in Iron
Pnictide Superconductor

Ryosuke Kurihara1 E-mail<EMAIL_ADDRESS>Keisuke
Mitsumoto1 Mitsuhiro Akatsu2 Yuichi Nemoto1,5
Terutaka Goto1,5 Corresponding author, Emeritus Prof. of Niigata Univ., E-mail
<EMAIL_ADDRESS>Yoshiaki Kobayashi3,5 and Masatoshi Sato3,4,5
Emeritus Prof. of Nagoya Univ., E-mail<EMAIL_ADDRESS>1 Graduate School
of Science and Technology1 Graduate School of Science and Technology Niigata
University Niigata University Niigata 950-2181 Niigata 950-2181 Japan
2 Faculty of Science Japan
2 Faculty of Science Niigata University Niigata University Niigata 950-2181
Niigata 950-2181 Japan
3 Department of Physics Japan
3 Department of Physics Nagoya University Nagoya University Nagoya 464-8602
Nagoya 464-8602 Japan
4 Research Center for Neutron Science and Technology Japan
4 Research Center for Neutron Science and Technology Comprehensive Research
Organization for Science and Society (CROSS) Comprehensive Research
Organization for Science and Society (CROSS) Ibaraki 300-0811 Ibaraki 300-0811
Japan
5 JST Japan
5 JST Transformative Research-Project on Iron Pnictides (TRIP) Transformative
Research-Project on Iron Pnictides (TRIP) Chiyoda Chiyoda Tokyo 102-0075 Tokyo
102-0075 Japan Japan

###### Abstract

Ultrasonic measurements have been carried out to investigate the critical
dynamics of structural and superconducting transitions due to degenerate
orbital bands in iron pnictide compounds with the formula Ba(Fe1-xCox)2As2.
The attenuation coefficient $\alpha_{\mathrm{L}[110]}$ of the longitudinal
ultrasonic wave for $(C_{11}+C_{12}+2C_{66})/2$ for $x=0.036$ reveals the
critical slowing down of the relaxation time around the structural transition
at $T_{\mathrm{s}}=65$ K, which is caused by ferro-type ordering of the
quadrupole $O_{x^{\prime 2}-y^{\prime 2}}$ coupled to the strain
$\varepsilon_{xy}$. The attenuation coefficient $\alpha_{66}$ of the
transverse ultrasonic wave for $C_{66}$ for $x=0.071$ also exhibits the
critical slowing down around the superconducting transition at
$T_{\mathrm{SC}}=23$ K, which is caused by ferro-type ordering of the
hexadecapole
$H_{z}^{\alpha}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}=O_{x^{\prime}y^{\prime}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}O_{x^{\prime
2}-y^{\prime 2}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}+O_{x^{\prime 2}-y^{\prime
2}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}O_{x^{\prime}y^{\prime}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}$
of the bound two-electron state coupled to the rotation $\omega_{xy}$. It is
proposed that the hexadecapole ordering associated with the superconductivity
brings about spontaneous rotation of the macroscopic superconducting state
with respect to the host tetragonal lattice.

##  1 Introduction

Since the discovery of superconductivity in iron pnictide compounds with the
formula LaFeAs(O1-xFx) with a high transition temperature of
$T_{\mathrm{SC}}=26$ K by Hosono and his coworkers in 2008 [1], many
researchers have been involved in the search for various iron-based composites
showing high superconducting transition temperatures and the inherent
mechanism of the superconductivity. Among the many iron-based compounds, the
family of LaFeAsO with the 1111-type ZrCuSiAs structure shows transition
temperatures as high as $T_{\mathrm{SC}}\sim 56$ K [2]. The family of BaFe2As2
with the 122-type ThCr2Si2 structure exhibits superconductivity with
$T_{\mathrm{SC}}\sim 38$ K upon the substitution of barium for potassium [3].
The other family of LiFeAs with the 111-type PbFCl structure shows
superconductivity with $T_{\mathrm{SC}}\sim 18$ K [4], and FeSe with the
11-type PbO structure shows superconductivity with $T_{\mathrm{SC}}\sim 10$ K
[5]. These compounds have a common lattice structure with a two-dimensional
square of iron layers, which gives the electronic band structure due to 3$d$
orbitals of Fe2+ ions favorable for achieving superconductivity with high
transition temperatures.

The 122-type compounds of BaFe2As2 showing the superconductivity under either
chemical doping or applying pressure have received particular attention [6, 7,
8], because high-quality single crystals with a fair size are available for
experiments. The end material of BaFe2As2 exhibits a structural phase
transition from the tetragonal phase with space group $D_{4h}^{17}$ $(I4/mmm)$
to the orthorhombic phase with space group $D_{2h}^{23}$ $(Fmmm)$
simultaneously accompanied by antiferro-magnetic ordering with a stripe-type
spin structure at transition temperatures of
$T_{\mathrm{s}}=T_{\mathrm{N}}=140$ K [3, 9, 10]. The chemical doping by
substituting Co ions with 3$d^{7}$ electrons for Fe ions with $3d^{6}$
electrons in Ba(Fe1-xCox)2As2 reduces both the structural and antiferro-
magnetic transition temperatures while slightly splitting the two transition
temperatures so that $T_{\mathrm{s}}>T_{\mathrm{N}}$ [11, 12, 13, 14]. The Co-
ion doped compounds with $x>0.03$ reveal the superconductivity below the
successive structural and magnetic transitions. With increasing the Co dopant
concentration $x$ to the quantum critical point (QCP) of
$x_{\mathrm{QCP}}=0.061$, the structural and antiferro-magnetic orderings
disappear and the superconducting phase manifests itself. Upon further doping
over $x_{\mathrm{QCP}}$, the optimized superconducting transition temperature
of as high as $T_{\mathrm{SC}}=23$ K emerges in the compound with $x=0.071$.

In order to clarify the interplay of the structural transition to the
superconductivity in the iron pnictide Ba(Fe1-xCox)2As2, the softening of the
elastic constant $C_{66}$ as a precursor of the structural transition has been
intensively investigated. By using resonance ultrasonic spectroscopy,
Fernandes et al. first showed the softening of $C_{66}$ in Ba(Fe1-xCox)2As2
for the end material and the over-doped compound with $x=0.08$, the latter
exhibiting a superconducting transition at $T_{\mathrm{SC}}=16$ K, and
proposed that the softening of $C_{66}$ is caused by spin nematic fluctuations
of the Fe $3d^{6}$ electrons [15]. By using the ultrasonic pulse-echo method,
Goto et al. showed the softening of $C_{66}$ by $21\%$ for $x=0.071$ from 300
K down to the optimized superconducting transition temperature
$T_{\mathrm{SC}}=23$ K, and deduced that the electric quadrupole $O_{x^{\prime
2}-y^{\prime 2}}$ associated with the degenerate $y^{\prime}z$ and
$zx^{\prime}$ orbitals of the Fe2+ ion coupled to the elastic strain
$\varepsilon_{xy}$ of the transverse ultrasonic wave plays a role in the
appearance of the superconductivity [16]. Here, the $x$ and $y$ ($x^{\prime}$
and $y^{\prime}$) coordinates for neighboring Ba-Ba (Fe-Fe) directions are
adopted. Yoshizawa et al. systematically investigated the elastic constant
$C_{66}$ for compounds with various Co concentrations and clarified the
quantum criticality of the softening of $C_{66}$ around the QCP of
$x_{\mathrm{QCP}}=0.061$ [17].

Angle-resolved photoemission spectroscopy experiments on Ba0.6K0.4Fe2As2 [18],
muon resonance experiments on Ba1-xKxFe2As2 [19], microwave penetration depth
measurements of PrFeAsO1-y [20], and NMR relaxation rate measurements of
LnFe1-yMyAsO1-xFx [21] commonly show mostly isotropic and nodeless
superconducting energy gaps with an s-like shape. The robustness of the
superconducting transition temperature $T_{\mathrm{SC}}$ for nonmagnetic
impurity doping to the system is favorable for the sign-conserving s-wave
state of s++ [21, 22, 23]. Taking the quadrupole fluctuations associated with
the elastic softening of $C_{66}$ into account, the s++ state for the
superconducting energy gap has been theoretically deduced [24, 25, 26]. On the
other hand, neutron scattering experiments on Ba(Fe0.925Co0.075)2As2 have
shown a magnetic excitation peak around $\boldsymbol{q}=(1/2,1/2,1)$,
indicating a role of the antiferro-magnetic fluctuations in the
superconductivity [27]. Accordingly, the sign-reversing s±-wave state has been
presented while emphasizing the spin fluctuation effects due to the antiferro-
magnetic interaction [28, 29, 30]. The clarification of the superconducting
mechanism of the iron pnictides is still an important issue in solid-state
physics.

Since ultrasonic waves with frequencies as high as 100 MHz easily penetrate
into metals, attenuation measurements are useful for examining order parameter
dynamics around phase transitions in metals. There are a few reports on
ultrasonic attenuation measurements around the structural and superconducting
transitions of the iron pnictides [31, 32]. However, the attenuation of
ultrasonic waves relating the elastic constant $C_{66}$ has not been reported
so far to the best of our knowledge. In the present paper, we show the
critical dynamics around the structural and superconducting transitions in the
iron pnictide Ba(Fe1-xCox)2As2 by means of ultrasonic measurements.

This paper is organized as follows. In Sect. 2, the experimental procedures of
sample preparation and ultrasonic measurements are shown. In Sect. 3, we
present experimental results of the critical slowing down of order parameters
around the superconducting transition for $x=0.071$ and the structural
transition for $x=0.036$. The critical temperatures characterizing the phase
diagram of the system are shown. In Sect. 4, we demonstrate the theory
treating the quantum degrees of freedom of the electric quadrupoles
$O_{x^{\prime 2}-y^{\prime 2}}$ and $O_{x^{\prime}y^{\prime}}$ and the angular
momentum $l_{z}$ associated with the degenerate $y^{\prime}z$ and
$zx^{\prime}$ orbitals and also their interactions with the transverse
acoustic phonons accompanying rotation and strain. In addition, we show
plausible models of the quadrupole ordering for the structural transition and
the hexadecapole ordering for the superconductivity in the iron pnictide
Ba(Fe1-xCox)2As2. We show that the superconducting Hamiltonian due to the
quadrupole interaction gives the superconducting ground state bearing the
hexadecapole. In Sect. 5, we present conclusions.

##  2 Experiments

Single crystals of Ba(Fe1-xCox)2As2 with Co concentration $x$ were grown by
the Bridgman method. Energy-dispersive X-ray spectrometry showed that nominal
concentrations of $x_{\mathrm{nominal}}=0.03,0.07$, and $0.10$ corresponded to
actual concentrations of $x=0.017,0.036$, and $0.071$, respectively. A Laue
X-ray camera was used to determine the crystallographic orientation of the
crystals. Samples with sizes of $2.30\times 1.05\times 0.44$ mm3 for the end
material $x=0$, $3.11\times 5.35\times 0.77$ mm3 for $x=0.017$, $2.13\times
2.06\times 0.26$ mm3 for $x=0.036$, and $1.52\times 2.61\times 0.83$ mm3 for
$x=0.071$ were used. The elastic constant $C=\rho_{\mathrm{M}}(x)v^{2}$ was
calculated from the ultrasonic velocity $v$. Mass densities of
$\rho_{\mathrm{M}}(x)=6.50$ g/cm3 for $x=$ 0, 0.017, and 0.036 and
$\rho_{\mathrm{M}}(x=0.071)=6.51$ g/cm3 were estimated using the Vegard’s law
for the lattice constants in Ba(Fe1-xCox)2As2 [33]. Piezoelectric ultrasonic
transducers with a 36∘ Y-cut LiNbO3 plate for the generation and detection of
longitudinal waves and an X-cut plate for transverse waves were glued on
plane-parallel surfaces of the samples. LiNbO3 plates with thicknesses of 200
and 100 $\mu$m were used to generate ultrasonic waves with frequencies as high
as 260 MHz. Because the present crystal shows a large amount of softening with
a considerably small elastic constant $C_{66}$, even a small external stress
$\sigma_{xy}$ might easily deform the sample with the tetragonal lattice to an
orthorhombic deformed lattice with finite $\varepsilon_{xy}$. This unwanted
lattice deformation might easily lift the degenerate $y^{\prime}z$ and
$zx^{\prime}$ orbital states. This lifting due to careless sample treatment
would prevent accurate measurements of the softening of $C_{66}$ and the
critical divergence of the ultrasonic attenuation across the structural and
superconducting transitions. In order to avoid applying undesired stress
$\sigma_{xy}$ to the samples, we carefully kept each sample in a brass disk
during the ultrasonic experiments. A vector-type detector based on the
ultrasonic pulse-echo method was used for simultaneous measurements of the
attenuation coefficient $\alpha$ and velocity $v$. A 3He cryostat (Oxford
Heliox TL) was employed for the ultrasonic measurements down to 250 mK.

##  3 Results

###  3.1 Attenuation around superconducting transition for $x=0.071$

In order to investigate the dynamical features of the superconducting
transition while clearly distinguishing them from the structural and
antiferro-magnetic orderings, we focused on the over-doped compound $x=0.071$,
which exhibits the optimized superconducting transition temperature
$T_{\mathrm{SC}}=23$ K but neither a structural transition nor antiferro-
magnetic ordering. We measured the attenuation coefficient $\alpha_{66}$ using
the transverse ultrasonic waves with the propagation vector
$\boldsymbol{q}//\left[100\right]$ and polarization vector
$\boldsymbol{\xi}//\left[010\right]$ with frequencies of 119, 167, and 215
MHz. The attenuation coefficient $\alpha_{66}$ for $x=0.071$ in Fig. 1(a)
increases with decreasing temperature below 80 K in the normal phase and
reveals critical divergence with approaching the superconducting transition
point of $T_{\mathrm{SC}}=23$ K. With further lowering the temperature below
$T_{\mathrm{SC}}=23$ K, the attenuation coefficient $\alpha_{66}$ rapidly
decreases. In both the normal and superconducting phases, the frequency
dependence of the attenuation coefficient $\alpha_{66}$ obeys the $\omega^{2}$
law, consistent with the low-frequency regime of the Debye formula.

Figure 1:  (Color online) (a) Temperature dependence of the ultrasonic
attenuation coefficient $\alpha_{66}$ measured using transverse ultrasonic
waves with frequencies of 119, 167, and 215 MHz, and the elastic constant
$C_{66}$ in Ba(Fe1-xCox)2As2 with $x=0.071$ having the superconducting
transition temperature $T_{\mathrm{SC}}=23$ K. The solid line is the softening
of $C_{66}$ in the normal phase fit by
$C_{66}=C_{66}^{0}\bigl{(}T-T_{\mathrm{s}}^{0}\bigr{)}/\bigl{(}T-\mathit{\Theta}_{\mathrm{Q}}\bigr{)}$
with $T_{\mathrm{s}}^{0}=-26.5$ K and $\mathit{\Theta}_{\mathrm{Q}}=-47$ K.
(b) Temperature dependence of the relaxation time $\tau$ obtained by the
attenuation coefficient $\alpha_{66}$ for frequencies 119, 167, and 215 MHz.
Solid lines show fits by $\tau=\tau_{0}\varepsilon^{-z\nu}$ for the reduced
temperature $\varepsilon=\left|T-T_{\mathrm{c}}^{0}\right|/T_{\mathrm{c}}^{0}$
with $T_{\mathrm{c}}^{0}=23$ K and $z\nu=1$. The inset in (b) shows the
temperature dependence of the relaxation rate $\tau^{-1}$. Solid lines in both
the normal and superconducting phases show fits by
$\tau^{-1}=\tau_{0}^{-1}\varepsilon^{z\nu}$.

The elastic constant $C_{66}$ in the normal phase for $x=0.071$ exhibits the
softening, which is again shown in Fig. 1(a) [16]. This softening is caused by
the interaction of the quadrupole $O_{x^{\prime 2}-y^{\prime 2}}$ to the
strain $\varepsilon_{xy}$ expressed as

$\displaystyle H_{\mathrm{QS}}=-g_{x^{\prime 2}-y^{\prime 2}}O_{x^{\prime
2}-y^{\prime 2}}\varepsilon_{xy}.$ (1)

Here, $g_{x^{\prime 2}-y^{\prime 2}}$ is the quadrupole-strain coupling
constant. We expect the Curie-type behavior of the quadrupole susceptibility
proportional to the reciprocal temperature as
$\chi_{\mathrm{Q}}=|\langle\psi_{l}|O_{x^{\prime 2}-y^{\prime
2}}|\psi_{l}\rangle|^{2}/T$ for suffixes $l=y^{\prime}z$ and $zx^{\prime}$ of
the degenerate orbitals. The softening of $C_{66}$ in the normal phase is
written by [34]

$\displaystyle C_{66}=C_{66}^{0}-\frac{n_{\mathrm{Q}}g_{x^{\prime 2}-y^{\prime
2}}^{2}\chi_{\mathrm{Q}}}{1-g^{\prime}_{x^{\prime 2}-y^{\prime
2}}\chi_{\mathrm{Q}}}=C_{66}^{0}\left(\frac{T-T_{\mathrm{s}}^{0}}{T-\mathit{\Theta}_{\mathrm{Q}}}\right).$
(2)

Here, $C_{66}^{0}$ is the background elastic constant when softening is absent
and $g^{\prime}_{x^{\prime 2}-y^{\prime 2}}$ is the mutual interaction
coefficient between the quadrupoles $O_{x^{\prime 2}-y^{\prime 2}}$ at
different sites. The solid line in Fig. 1(a) is the softening of $C_{66}$ fit
by Eq. (2) in the normal phase. We obtain negative values of the quadrupole
interaction energy $\mathit{\Theta}_{\mathrm{Q}}=-47.0$ K and the critical
temperature
$T_{\mathrm{s}}^{0}=\mathit{\Theta}_{\mathrm{Q}}+\mathit{\Delta}_{\mathrm{Q}}=-26.5$
K. This fitting gives the quadrupole-strain interaction energy
$\mathit{\Delta}_{\mathrm{Q}}=20.5$ K, which represents the coupling energy
between the quadrupoles at different sites mediated by the strain
$\varepsilon_{xy}$.

The negative quadrupole interaction energy
$\mathit{\Theta}_{\mathrm{Q}}=-47.0$ K indicates the antiferro-type quadrupole
interaction for $x=0.071$. The negative critical temperature of
$T_{\mathrm{s}}^{0}=-26.5$ K corresponds to the fictitious critical
temperature of the structural instability for $C_{66}\rightarrow 0$.
Consequently, the superconducting transition temperature $T_{\mathrm{SC}}=23$
K is definitely distinguished from the fictitious critical temperature
$T_{\mathrm{s}}^{0}=-26.5$ K. We conclude that the marked increase in the
attenuation coefficient $\alpha_{66}$ for $x=0.071$ in Fig. 1(a) is caused by
the critical slowing down of the order parameter around the superconducting
transition, which is strictly distinguished from the quadrupole $O_{x^{\prime
2}-y^{\prime 2}}$ interacting to the strain $\varepsilon_{xy}$ in Eq. (1).

Because two electrons are accommodated in the degenerate $y^{\prime}z$ and
$zx^{\prime}$ orbitals of an Fe2+ ion, we tentatively adopt
$n_{\mathrm{Q}}=2N_{\mathrm{Fe}}=3.92\times 10^{22}$ cm-3 in Eq. (2) as the
maximum number of electrons, that play a role in the softening of $C_{66}$.
Here, $N_{\mathrm{Fe}}$ is the number of Fe2+ ions per unit volume. The solid
line in Fig. 1(a) for the softening of $C_{66}$ fit by Eq. (2) gives the
quadrupole-strain coupling constant $g_{x^{\prime 2}-y^{\prime 2}}=1045$ K per
electron. This coupling constant is comparable with the results for manganese
compounds of $g=1167$ K for La0.88Sr0.12MnO3 and $g=1020$ K for
Pr0.65Ca0.35MnO3 with considerably extended 3$d$ orbitals [35, 36], but is
considerably larger than $g\sim 100$ K for rare-earth compounds of the form
RB6 with well-screened 4$f$ orbitals in the inner shell [37].

In our attempt to clarify the order parameter dynamics around the
superconducting transition, we analyzed the frequency dependence of the
attenuation coefficient $\alpha_{66}$ in Fig. 1(a) in terms of the Debye
formula expressed as

$\displaystyle\alpha_{66}=\frac{C_{66}\left(\infty\right)-C_{66}\left(0\right)}{2\rho_{\mathrm{M}}{v_{66}\left(\infty\right)}^{3}}\frac{\omega^{2}\tau}{1+\omega^{2}\tau^{2}}.$
(3)

Here, $C_{66}\left(\infty\right)$ and $C_{66}\left(0\right)$ stand for the
high- and low-frequency limits of the elastic constant $C_{66}$, respectively.
We regard the elastic constant experimentally observed in Fig. 1(a) as the
low-frequency limit $C_{66}\left(0\right)$ in Eq. (3). The background
$C_{66}^{0}$ used in the analysis based on the quadrupole susceptibility of
Eq. (2) corresponds to the high-frequency limit of the elastic constant
$C_{66}\left(\infty\right)$ for the sound velocity
$v_{66}\left(\infty\right)=\sqrt{C_{66}\left(\infty\right)/\rho_{\mathrm{M}}}$
in Eq. (3) [16]. $\omega$ is the angular frequency of the employed ultrasonic
waves, and $\tau$ is the relaxation time of the order parameter fluctuation of
the system. We show the temperature dependence of the relaxation time $\tau$
obtained by Eq. (3) in Fig. 1(b) and that of the relaxation rate $\tau^{-1}$
in the inset of Fig. 1(b). The rapid increase in the relaxation time $\tau$
around $T_{\mathrm{SC}}=23$ K indeed indicates the critical slowing down of
the order parameter accompanying the superconducting transition in the system.

The correlation length $\zeta$ of the order parameter
$Q\left(\boldsymbol{r}\right)$ for an appropriate phase transition is widely
used to describe the correlation function as
$G\left(\boldsymbol{r}\right)=\left\langle
Q\left(\boldsymbol{r}\right)Q\left(0\right)\right\rangle\propto
r^{-\left(d-1\right)}e^{-r/\zeta}$. Here, $|\boldsymbol{r}|=r$ stands for the
spatial distance between the order parameters and $d$ denotes the dimension of
the system. In the vicinity of the critical temperature $T_{\mathrm{c}}^{0}$
of the phase transition, the correlation length $\zeta$ becomes infinite due
to the increase in the local order as $\zeta=\zeta_{0}\varepsilon^{-\nu}$ for
a reduced temperature of
$\varepsilon=\bigl{|}T-T_{\mathrm{c}}^{0}\bigr{|}/T_{\mathrm{c}}^{0}$. The
critical index $\nu=1/2$ is expected from mean field theory [38].

When an ultrasonic pulse wave enters a crystal, the equilibrium state of the
system is instantaneously perturbed to a nonequilibrium state. After the
relaxation time $\tau$, the system returns to the equilibrium state. In the
vicinity of the critical point, however, the relaxation time $\tau$ becomes
infinite due to the divergence of the correlation length as
$\tau\propto\zeta^{z}$ for dynamical critical index $z$. Presumably, the
critical slowing down of the relaxation time $\tau$ is explained by the
critical index $z\nu$ as [39, 40, 41, 42]

$\displaystyle\tau=\tau_{0}\left|\frac{T-T_{\mathrm{c}}^{0}}{T_{\mathrm{c}}^{0}}\right|^{-z\nu}=\tau_{0}\varepsilon^{-z\nu}.$
(4)

In mean field theory, the dynamical critical index $z=2$ is expected.
Actually, the temperature dependence of the relaxation time $\tau$ in the
normal phase of Fig. 1(b) is well reproduced by the solid line for the
critical index $z\nu=1$, the critical temperature $T_{\mathrm{c}}^{0+}=23.0$
K, and the attempt time $\tau_{0}^{+}=6.0\times 10^{-11}$ s. In the
superconducting phase, however, we obtain $z\nu=1/3$,
$T_{\mathrm{c}}^{0-}=23.5$ K, and $\tau_{0}^{-}=5.8\times 10^{-11}$ s. The
distinct deviation of the critical index of $z\nu=1/3$ from $z\nu=1$ of mean
field theory may be caused by the inherent property that the hexadecapole
ordering appears in accordance with the superconductivity in the present iron
pnictide. The analysis of the experimental results in Fig. 1(b) gives a ratio
of the attempt times of $\tau_{0}^{+}/\tau_{0}^{-}=1.03$, which is
distinguished from the ratio of $\tau_{0}^{+}/\tau_{0}^{-}=2$ expected from
mean field theory. Because the ultrasonic echo signal almost disappears in the
vicinity of the superconducting transition point due to the critical slowing
down, the absolute values of the attenuation coefficients inevitably include
experimental errors.

###  3.2 Attenuation around structural transition for $x=0.036$

Next, we examine the critical slowing down around the structural transition
from the tetragonal to orthorhombic phases of the under-doped compound
$x=0.036$. This is worth comparing with the critical slowing down around the
superconducting transition of the over-doped compound $x=0.071$ presented in
Sect. 3.1.

Figure 2:  (Color online) Temperature dependence of elastic constants of
Ba(Fe1-xCox)2As2 with $x=0.036$. The considerable softening of $C_{66}$ in (e)
and $C_{\mathrm{L}[110]}=\left(C_{11}+C_{12}+2C_{66}\right)/2$ in (d) are
found to be a precursor of the structural transition at $T_{\mathrm{s}}=65$ K.
Anomalies associated with the antiferro-magnetic transition at
$T_{\mathrm{N}}=39$ K and the superconducting transition at
$T_{\mathrm{SC}}=16.4$ K are observed. The softening of $C_{66}$ in (e) is
reproduced by a solid line fit by
$C_{66}=C_{66}^{0}\bigl{(}T-T_{\mathrm{s}}^{0}\bigr{)}/\bigl{(}T-\mathit{\Theta}_{\mathrm{Q}}\bigr{)}$
with $T_{\mathrm{s}}^{0}=63$ K and $\mathit{\Theta}_{\mathrm{Q}}=47$ K.

The temperature dependence of the elastic constants for $x=0.036$ is shown in
Fig. 2. We denote the propagation vector as $\boldsymbol{q}$ and the
polarization vector as $\boldsymbol{\xi}$ for the ultrasonic waves measured in
Fig. 2. The elastic constant $C_{66}$ in Fig. 2(e) measured by the transverse
ultrasonic wave reveals considerable softening of 85$\%$ with decreasing
temperature from 300 K down to the structural transition temperature of
$T_{\mathrm{s}}=65$ K. As indicated by arrows in Fig. 2(e), the antiferro-
magnetic ordering at $T_{\mathrm{N}}=39$ K and the superconducting transition
at $T_{\mathrm{SC}}=16.4$ K are also observed as anomalies in $C_{66}$. The
softening of $C_{66}$ in the tetragonal phase is reproduced by the solid line
in Fig. 2(e), which is obtained by Eq. (2) for the quadrupole interaction
energy of $\mathit{\Theta}_{\mathrm{Q}}=47$ K, the quadrupole-strain
interaction energy $\mathit{\Delta}_{\mathrm{Q}}=16$ K, and the critical
temperature of
$T_{\mathrm{s}}^{0}=\mathit{\Theta}_{\mathrm{Q}}+\mathit{\Delta}_{\mathrm{Q}}=63$
K. The background of $C_{66}^{0}=2.87\times 10^{10}$ J/m3 is shown by the
dashed line in Fig. 2(e). The critical temperature $T_{\mathrm{s}}^{0}=63$ K
is in agreement with the experimentally observed structural transition
temperature $T_{\mathrm{s}}=65$ K. We deduce the quadrupole-strain coupling
constant to be $g_{x^{\prime 2}-y^{\prime 2}}=920$ K per electron by adopting
the electron number $n_{\mathrm{Q}}=2N_{\mathrm{Fe}}$ in Eq. (2). The positive
values of the quadrupole interaction energy $\mathit{\Theta}_{\mathrm{Q}}=47$
K and the critical temperature $T_{\mathrm{s}}^{0}=63$ K are consistent with
the ferro-quadrupole ordering $O_{x^{\prime 2}-y^{\prime 2}}$ accompanying the
structural transition from the tetragonal to orthorhombic phase.

The elastic constant
$C_{\mathrm{L}\left[110\right]}=\left(C_{11}+C_{12}+2C_{66}\right)/2$ in Fig.
2(d) measured by a longitudinal ultrasonic wave with the propagation direction
$\boldsymbol{q}//\left[110\right]$ parallel to the polarization direction
$\boldsymbol{\xi}$ also shows marked softening with decreasing temperature
down to the structural transition at $T_{\mathrm{s}}=65$ K. With further
decreasing temperature in the distorted orthorhombic phase, the graph of the
elastic constant $C_{\mathrm{L}\left[110\right]}$ slightly bend at the
antiferro-magnetic transition temperature $T_{\mathrm{N}}=39$ K and distinctly
decreases at the superconducting transition temperature $T_{\mathrm{SC}}=16.4$
K. The longitudinal ultrasonic wave for $C_{\mathrm{L}\left[110\right]}$
induces elastic strain of
$\varepsilon_{\mathrm{L}\left[110\right]}=\varepsilon_{\mathrm{B}}/3-\varepsilon_{u}/\bigl{(}2\sqrt{3}\bigr{)}+\varepsilon_{xy}$.
The coupling of the strain $\varepsilon_{xy}$ to the quadrupole $O_{x^{\prime
2}-y^{\prime 2}}$ induces the softening in $C_{\mathrm{L}\left[110\right]}$,
partly consisting of $C_{66}$. The volume strain
$\varepsilon_{\mathrm{B}}=\varepsilon_{xx}+\varepsilon_{yy}+\varepsilon_{zz}$
and the tetragonal strain
$\varepsilon_{u}=\bigl{(}2\varepsilon_{zz}-\varepsilon_{xx}-\varepsilon_{yy}\bigr{)}/\sqrt{3}$
scarcely affect the softening.

Figure 3:  (Color online) (a) Temperature dependence of the ultrasonic
attenuation coefficient $\alpha_{\mathrm{L}[110]}$ measured using longitudinal
ultrasonic waves with frequencies of 112, 186, and 260 MHz, and the elastic
constant $C_{\mathrm{L}[110]}=\left(C_{11}+C_{12}+2C_{66}\right)/2$ in
Ba(Fe1-xCox)2As2 with $x=0.036$. (b) Temperature dependence of the relaxation
time $\tau$ obtained from the results for 112, 186, and 260 MHz ultrasonic
waves. Solid lines indicate fits by $\tau=\tau_{0}\varepsilon^{-z\nu}$ for the
reduced temperature
$\varepsilon=\left|T-T_{\mathrm{c}}^{0}\right|/T_{\mathrm{c}}^{0}$ with
$T_{\mathrm{s}}^{0}=65$ K and the critical index $z\nu=1$ for both tetragonal
and orthorhombic phases. The inset shows the temperature dependence of the
relaxation rate $\tau^{-1}$ and the solid lines are fits by
$\tau^{-1}=\tau_{0}^{-1}\varepsilon^{z\nu}$.

We initially attempted to measure the attenuation coefficient $\alpha_{66}$ by
using the pure transverse ultrasonic wave. However, the attenuation of
$\alpha_{66}$ near the structural transition was too large to measure. In the
present experiments, therefore, we used the longitudinal ultrasonic wave,
which shows relatively moderate damping, to measure the attenuation
coefficient $\alpha_{\mathrm{L}\left[110\right]}$. In Fig. 3(a), we show the
attenuation coefficient $\alpha_{\mathrm{L}\left[110\right]}$ for $x=0.036$
acquired using longitudinal ultrasonic waves with frequencies of 112, 186, and
260 MHz together with the elastic constant $C_{\mathrm{L}\left[110\right]}$.
The attenuation coefficient $\alpha_{\mathrm{L}\left[110\right]}$ increases
considerably with approaching the structural transition point of
$T_{\mathrm{s}}=65$ K from both sides. With decreasing temperature in the
orthorhombic phase, the attenuation $\alpha_{\mathrm{L}\left[110\right]}$
exhibits a distinct peak at the superconducting transition point of
$T_{\mathrm{SC}}=16.4$ K. The missing ultrasonic echo signal due to
considerable damping of the longitudinal ultrasonic wave prevents us from
acquiring $\alpha_{\mathrm{L}\left[110\right]}$ approaching the structural
transition. Nevertheless, the distinct tendency of the divergence of
$\alpha_{\mathrm{L}\left[110\right]}$ around the structural transition at
$T_{\mathrm{s}}=65$ K indicates the critical slowing down of the relaxation
time $\tau$ for the order parameter fluctuation.

The frequency dependence for $\alpha_{\mathrm{L}\left[110\right]}$ for the
under-doped compound $x=0.036$ in Fig. 3(a) is analyzed in terms of the Debye
formula in Eq. (3), where we read the attenuation as
$\alpha_{\mathrm{L}\left[110\right]}$ and the elastic constant as
$C_{\mathrm{L}\left[110\right]}$. We depict the temperature dependence of the
relaxation time $\tau$ in Fig. 3(b) and the relaxation rate $\tau^{-1}$ in the
inset of Fig. 3(b). The increase in the relaxation time $\tau$ around the
structural transition at $T_{\mathrm{s}}=65$ K is caused by the critical
slowing down of the quadrupole $O_{x^{\prime 2}-y^{\prime 2}}$, which is the
order parameter of the structural transition. The ultrasonic frequencies of
$\omega/2\pi$ up to 260 MHz used in the present experiments are much higher
than the relaxation rate $\tau^{-1}<100$ MHz in the vicinity of the structural
transition point $T_{\mathrm{s}}$ for the narrow reduced temperature region of
$\varepsilon=|T-T_{\mathrm{c}}^{0}|/T_{\mathrm{c}}^{0}<0.038$. As a result,
the ultrasonic echo signal disappears on both sides around the structural
transition point of $T_{\mathrm{s}}=65$ K as depicted in Fig. 3(a). The solid
lines in Fig. 3(b) are fits for the relaxation time $\tau$ in terms of Eq.
(4). Supposing that the critical temperature $T_{\mathrm{c}}^{0}=65$ K
coincides with the structural transition point and the critical indices for
both phases are consistent with $z\nu=1$ of mean field theory, we obtain the
fit shown by the solid line in Fig. 3(b) using Eq. (4). This gives the attempt
time $\tau_{0}^{+}=0.380\times 10^{-9}$ s for the tetragonal phase and
$\tau_{0}^{-}=1.15\times 10^{-9}$ s for the orthorhombic phase.

The quadrupole-strain interaction of Eq. (1) induces the ferro-quadrupole
ordering accompanying the structural transition. It is expected that the long-
range Coulomb force between electrons bearing the quadrupole exhibits critical
phenomena described by mean field theory. Actually, the adoption of the
critical index $z\nu=1$ for both tetragonal and orthorhombic phases consistent
with mean field theory reproduces the experimental results in Fig. 3. The
small-echo signal due to the considerable damping of the longitudinal
ultrasonic waves gives inevitable errors in the absolute-value of the
attenuation coefficient $\alpha_{\mathrm{L}[110]}$ in Fig. 3. The dispersive
scattering of the ultrasonic wave by the domain wall due to the orthorhombic
distortion is also included in the attenuation. These ambiguity might have
resulted in the deviation of the experimentally determined ratio of
$\tau_{0}^{+}/\tau_{0}^{-}=0.33$ from $\tau_{0}^{+}/\tau_{0}^{-}=2$ expected
from mean field theory. The attempt time $\tau_{0}^{+}=0.55\times 10^{-9}$ s
of the structural transition for $x=0.036$ in Fig. 3(b) is one order of
magnitude larger than the attempt time $\tau_{0}^{+}=6.0\times 10^{-11}$ s of
the superconducting transition for $x=0.071$ in Fig. 1(b). This notable result
suggests that the order parameter showing the critical slowing down around the
superconducting transition for $x=0.071$ is distinguished from the ferro-
quadrupole ordering accompanying the structural transition for $x=0.036$.

###  3.3 Critical temperatures of the system

The ultrasonic measurements of the critical divergence in the attenuation
coefficients and the softening in the elastic constants provided us with the
critical temperatures, which characterize the structural and superconducting
transitions of the present iron pnictide Ba(Fe1-xCox)2As2. In Fig. 4, we plot
the critical temperatures obtained by the experimental results for $x=0.036$
and 0.071 together with those for $x=0.017$ and the end material $x=0$ [43].

Figure 4:  (Color online) Critical temperatures relating to structural (ST)
and superconducting (SC) transitions in Ba(Fe1-xCox)2As2 determined by the
present ultrasonic measurements. The quadrupole interaction (Q intr.) energy
$\mathit{\Theta}_{\mathrm{Q}}$, quadrupole-strain interaction (Q-S intr.)
energy $\mathit{\Delta}_{\mathrm{Q}}$, and structural transition temperature
$T_{\mathrm{s}}^{0}=\mathit{\Theta}_{\mathrm{Q}}+\mathit{\Delta}_{\mathrm{Q}}$
associated with the ferro-quadrupole ordering are determined by the elastic
constant $C_{66}\rightarrow 0$. The critical temperature $T_{\mathrm{c}}^{0}$
is determined by the critical slowing down of the relaxation time
$\tau\rightarrow\infty$ measured by the ultrasonic attenuation coefficient
$\alpha$. Structural transition temperatures $T_{\mathrm{s}}$ associated with
ferro-quadrupole (FQ) ordering are indicated by black open circles.
Superconducting transition temperatures $T_{\mathrm{SC}}$ associated with
ferro-hexadecapole (FH) ordering are indicated by black open squares.
Superconducting transition temperatures indicated by black open rhombuses are
results of Ni $et$ $al$. [11].

The softening of the elastic constant $C_{66}$ analyzed by the quadrupole
susceptibility $\chi_{\mathrm{Q}}$ in Eq. (2) for $O_{x^{\prime 2}-y^{\prime
2}}$ gives the critical temperature
$T_{\mathrm{s}}^{0}=\mathit{\Theta}_{\mathrm{Q}}+\mathit{\Delta}_{\mathrm{Q}}$
corresponding to the structural instability point as $C_{66}\rightarrow 0$. As
shown by blue closed circles and the blue solid line in Fig. 4, the critical
temperature of $T_{\mathrm{s}}^{0}=135$ K of the end material $x=0$ decreases
to $T_{\mathrm{s}}^{0}=100$ K for $x=0.017$ and $T_{\mathrm{s}}^{0}=63$ K for
$x=0.036$. Beyond the QCP of $x_{\mathrm{QCP}}=0.061$, the critical
temperature changes its sign to $T_{\mathrm{s}}^{0}=-26.5$ K for $x=0.071$.
The quadrupole interaction energy $\mathit{\Theta}_{\mathrm{Q}}$ shown by red
closed circles and the red solid line decreases from
$\mathit{\Theta}_{\mathrm{Q}}=115$ K for the end material $x=0$ to
$\mathit{\Theta}_{\mathrm{Q}}=80$ K for $x=0.017$,
$\mathit{\Theta}_{\mathrm{Q}}=47$ K for $x=0.036$, and changes to a negative
value of $\mathit{\Theta}_{\mathrm{Q}}=-47$ K for $x=0.071$. On the other
hand, the quadrupole-strain interaction energy $\mathit{\Delta}_{\mathrm{Q}}$
shown by green closed circles and the green solid line remains at the positive
value $\mathit{\Delta}_{\mathrm{Q}}\sim 20$ K. The positive critical
temperatures $T_{\mathrm{s}}^{0}$ of the under-doped compounds are comparable
with the ferro-quadrupole ordering $O_{x^{\prime 2}-y^{\prime 2}}$ associated
with the structural transition from the tetragonal to orthorhombic phase,
while the negative critical temperature $T_{\mathrm{s}}^{0}=-26.5$ K for the
over-doped compound $x=0.071$, suggesting absence of the ferro-quadrupole
ordering, implies the fictitious lattice instability point.

The infinite divergence of the relaxation time $\tau$ due to the critical
slowing down provides the critical temperature $T_{\mathrm{c}}^{0}$ at which
the long-range ordering appears due to freezing of the order parameter
fluctuation [42, 41]. In the case of the structural transition for $x=0.036$,
the divergence of the relaxation time $\tau$ gives the critical temperature of
$T_{\mathrm{c}}^{0}=65$ K, which coincides with the structural transition
temperature $T_{\mathrm{s}}=65$ K experimentally observed, and is in agreement
with the critical temperature
$T_{\mathrm{s}}^{0}=\mathit{\Theta}_{\mathrm{Q}}+\mathit{\Delta}_{\mathrm{Q}}=63$
K obtained by fitting the elastic softening of $C_{66}$ by Eq. (2). This
result confirms that the critical divergence of the ultrasonic attenuation
coefficient $\alpha_{\mathrm{L}\left[110\right]}$ for $x=0.036$ is caused by
the critical slowing down of the ferro-quadrupole order parameter
$O_{x^{\prime 2}-y^{\prime 2}}$. The critical temperature for
$T_{\mathrm{c}}^{0}=65$ K for $x=0.036$ is shown in Fig. 4 by a purple closed
asterisk together with the critical temperature of $T_{\mathrm{s}}^{0}=63$ K.
Note that the quadrupole-strain interaction energy
$\mathit{\Delta}_{\mathrm{Q}}\sim 20$ K, which is almost independent of the Co
concentration, also promotes the appearance of the structural transition at
the critical temperature
$T_{\mathrm{s}}^{0}=\mathit{\Theta}_{\mathrm{Q}}+\mathit{\Delta}_{\mathrm{Q}}$.
The quadrupole-strain interaction energy $\mathit{\Delta}_{\mathrm{Q}}$ is
shown by green closed circles in Fig. 4.

In the case of the superconducting transition for $x=0.071$, the critical
temperature $T_{\mathrm{c}}^{0}=23$ K determined by the critical slowing down
of the relaxation time $\tau$ is in good agreement with the superconducting
transition temperature $T_{\mathrm{SC}}=23$ K. In Fig. 4, we show the former
critical temperature $T_{\mathrm{c}}^{0}=23$ K for the relaxation time $\tau$
by a purple closed asterisk together with the latter superconducting
transition temperature $T_{\mathrm{SC}}$ = 23 K by a black open square. The
critical temperature $T_{\mathrm{s}}^{0}=-26.5$ K indicated by a blue filled
circle corresponds to the fictitious lattice instability point as
$C_{66}\rightarrow 0$. The fictitious instability point of
$T_{\mathrm{s}}^{0}=-26.5$ K is definitely distinguished from the critical
temperature $T_{\mathrm{c}}^{0}=23$ K of the relaxation time $\tau$ coinciding
with the superconducting transition temperature $T_{\mathrm{SC}}$. This result
evidences that the order parameter exhibiting the critical slowing down around
the superconducting transition does not correspond to the electric quadrupole
$O_{x^{\prime 2}-y^{\prime 2}}$ but to other quantum degrees of freedom, which
interact with the transverse ultrasonic waves used in the experiments. It is
of great importance to identify the order parameter, which brings about the
critical slowing down of the relaxation time $\tau$ around the superconducting
transition. In Sect. 4.6, we give the electric hexadecapole carried by a two-
electron state as a plausible order parameter associated with the
superconducting transition.

##  4 Theory

###  4.1 Transverse acoustic wave

The relaxation time $\tau$ determined by the ultrasonic attenuation reveals
the critical slowing down around the superconducting transition for the
compound $x=0.071$ as well as the structural transition for $x=0.036$. In our
attempt to properly explain the ultrasonic experiments, we treat the couplings
between the ultrasonic waves and electrons accommodated in the degenerate
$y^{\prime}z$ and $zx^{\prime}$ orbitals of Fe ion. We will show plausible
order parameters associated with the superconducting phase and the
structurally distorted phase in the system.

The transverse ultrasonic wave with the propagation vector
$\boldsymbol{q}=\bigl{(}q_{x},0,0\bigr{)}$ and the polarization vector
$\boldsymbol{\xi}=\bigl{(}0,\xi_{y},0\bigr{)}$ employed in the present
experiments is expressed in terms of the plane wave as

$\displaystyle\xi_{y}=\xi_{y}^{0}\exp\left[i\left(q_{x}x-\omega
t\right)\right].$ (5)

Here, $\xi_{y}^{0}$ is the amplitude of the ultrasonic wave, $q_{x}$ is the
wavevector component, and $\omega$ is the angular frequency. The elastic
constant is determined by the sound velocity $v_{x}=\omega/q_{x}=v_{66}$ as
$C_{66}=\rho v_{x}^{2}$. The transverse ultrasonic wave induces the
deformation tensor

$\displaystyle\frac{\partial\xi_{y}}{\partial
x}=iq_{x}\xi_{y}^{0}\exp\left[i\left(q_{x}x-\omega
t\right)\right]=iq_{x}\xi_{y}.$ (6)

The deformation tensor of Eq. (6) associated with the transverse ultrasonic
wave consists of the strain $\varepsilon_{xy}$ expressed as a symmetric tensor
and the rotation $\omega_{xy}$ as an antisymmetric tensor

$\displaystyle\left\\{\begin{array}[]{c}\varepsilon_{xy}\\\ \omega_{xy}\\\
\end{array}\right\\}=\frac{1}{2}\left(\frac{\partial\xi_{y}}{\partial
x}\pm\frac{\partial\xi_{x}}{\partial y}\right).$ (9)

In experiments using the transverse ultrasonic wave with
$\boldsymbol{q}=(q_{x},0,0)$ and $\boldsymbol{\xi}=(0,\xi_{y},0)$, we tune the
perturbation parameter of
$\delta=q_{x}\xi_{y}^{0}=(\omega/v_{x})\xi_{y}^{0}=(2\pi/\lambda_{x})\xi_{y}^{0}$
to an infinitesimal value by controlling the amplitude of the generated
ultrasonic wave $\xi_{y}^{0}$ and the frequency $\omega$ for a given phase
velocity $v_{x}=v_{66}$. The rotation
$\omega_{xy}=q_{x}\xi_{y}^{0}=\sin\theta$ and the strain
$\varepsilon_{xy}=q_{x}\xi_{y}^{0}=\sin\theta$ are caused by slight twisting
of the $x$\- and $y$-axes by an angle of $\pi/2\pm 2\theta$. This will later
be illustrated in Figs. 7(a) and 7(b). The transverse ultrasonic wave
simultaneously induces the rotation $\omega_{xy}$ and strain
$\varepsilon_{xy}$, which are strictly distinguished from each other from a
symmetrical viewpoint. Concerning the space group $D_{4h}^{17}$ of the
tetragonal lattice, the rotation $\omega_{xy}$ is compatible with $A_{2}$
symmetry, while the strain $\varepsilon_{xy}$ is compatible with $B_{2}$ [44]
symmetry. Thus, it is expected that the rotation $\omega_{xy}$ and strain
$\varepsilon_{xy}$ interact with the electronic states of the system in
different manners. When the electron-phonon interaction is absent, both the
rotation $\omega_{xy}$ and strain $\varepsilon_{xy}$ propagate in the lattice
with the same velocity $v_{66}$.

In our attempt to explain the critical divergence of the ultrasonic
attenuation coefficients and the marked elastic softening, we suppose that the
interactions of thermally excited transverse acoustic phonons with electronic
states play a role in the appearance of the structural and superconducting
transitions. In the interactions of the transverse acoustic phonons with the
electronic states of the iron pnictide, we take the Hamiltonian for the
harmonic oscillators with polarization vectors in the $xy$ plane as

$\displaystyle H_{\mathrm{ph}}$
$\displaystyle=\sum_{\boldsymbol{q}}\left[\hbar\omega_{x}\left(\boldsymbol{q}\right)\left(a_{x,\boldsymbol{q}}^{\dagger}a_{x,\boldsymbol{q}}+\frac{1}{2}\right)\right.$
$\displaystyle\left.\qquad\qquad+\hbar\omega_{y}\left(\boldsymbol{q}\right)\left(a_{y,\boldsymbol{q}}^{\dagger}a_{y,\boldsymbol{q}}+\frac{1}{2}\right)\right].$
(10)

Here, $a_{i,\boldsymbol{q}}$ and $a_{i,\boldsymbol{q}}^{\dagger}$ respectively
denote annihilation and creation operators of the transverse acoustic phonon
with the polarization direction $\xi_{i}$ for $i=x$ and $y$ and momentum
$\hbar\boldsymbol{q}$. $\hbar\omega_{i}\left(\boldsymbol{q}\right)$ stands for
the phonon energy. The strain $\varepsilon_{xy}\left(\boldsymbol{r},t\right)$
and rotation $\omega_{xy}\left(\boldsymbol{r},t\right)$ at position
$\boldsymbol{r}$ and time $t$ are expressed in terms of the phonon operators
$a_{i,\boldsymbol{q}}$ and $a_{i,\boldsymbol{q}}^{\dagger}$ as [45]

$\displaystyle\left\\{\begin{array}[]{c}\varepsilon_{xy}\left(\boldsymbol{r},t\right)\\\
\omega_{xy}\left(\boldsymbol{r},t\right)\\\ \end{array}\right\\}$
$\displaystyle=\frac{1}{2}\left[\frac{\partial\xi_{y}\left(\boldsymbol{r},t\right)}{\partial
x}\pm\frac{\partial\xi_{x}\left(\boldsymbol{r},t\right)}{\partial y}\right]$
(13)
$\displaystyle=\frac{i}{2}\sum_{\boldsymbol{q}}\sqrt{\mathstrut\frac{\hbar}{2V\rho_{\mathrm{M}}\omega_{y}\left(\boldsymbol{q}\right)}}q_{x}$
$\displaystyle\quad\times\left[a_{y,\boldsymbol{q}}e^{i\boldsymbol{q}\cdot\boldsymbol{r}}e^{-i\omega_{y}\left(\boldsymbol{q}\right)t}-a_{y,\boldsymbol{q}}^{\dagger}e^{-i\boldsymbol{q}\cdot\boldsymbol{r}}e^{i\omega_{y}\left(\boldsymbol{q}\right)t}\right]$
$\displaystyle\pm\frac{i}{2}\sum_{\boldsymbol{q}}\sqrt{\mathstrut\frac{\hbar}{2V\rho_{\mathrm{M}}\omega_{x}\left(\boldsymbol{q}\right)}}q_{y}$
$\displaystyle\quad\times\left[a_{x,\boldsymbol{q}}e^{i\boldsymbol{q}\cdot\boldsymbol{r}}e^{-i\omega_{x}\left(\boldsymbol{q}\right)t}-a_{x,\boldsymbol{q}}^{\dagger}e^{-i\boldsymbol{q}\cdot\boldsymbol{r}}e^{i\omega_{x}\left(\boldsymbol{q}\right)t}\right].$
(14)

Here, $\rho_{M}$ is the mass density and $V$ is the volume of the system.

###  4.2 Degenerate $\psi_{y^{\prime}z}$ and $\psi_{zx^{\prime}}$ orbitals

Many reports of band-structure calculations on the iron pnictide compounds
have shown inherent 3$d$ electronic structures favorable for the
superconductivity as well as the structural transition [46, 47, 48, 49, 50].
Three sheets of hole surfaces with a cylinder-type structure exist around the
$\Gamma$-point of the zone center and two electron pockets exist around the
X-points of the zone boundary. We focus on the degenerate $y^{\prime}z$ and
$zx^{\prime}$ bands occupied up to half filling, which are lifted by the
transverse ultrasonic waves and acoustic phonons for $C_{66}$ with a small-
wavenumber limit of $|\boldsymbol{q}|\rightarrow 0$.

We show the electronic states of the degenerate $y^{\prime}z$ and
$zx^{\prime}$ orbitals in the crystalline electric field (CEF) Hamiltonian on
an Fe2+ ion with $D_{2d}$ site symmetry, which consists of an electric
monopole, quadrupole, and hexadecapole expressed in terms of spherical
harmonics as

$\displaystyle H_{\mathrm{CEF}}$ $\displaystyle=A_{0}^{0}$
$\displaystyle\quad+A_{2}^{0}\left(\frac{3z^{2}-r^{2}}{2r^{2}}\right)$
$\displaystyle\qquad+A_{4}^{0}\left(\frac{35z^{4}-30z^{2}r^{2}+3r^{4}}{8r^{4}}\right)$
$\displaystyle\quad\qquad+A_{4}^{4}\left(\frac{\sqrt{35}}{8}\frac{x^{\prime
4}-6x^{\prime 2}y^{\prime 2}+y^{\prime 4}}{r^{4}}\right).$ (15)

Here, $A_{l}^{m}$ denotes CEF parameters for the electric multipole
potentials. The 3$d$ orbitals of the Fe2+ ion in the CEF Hamiltonian of Eq.
(4.2) splits into five orbital states with energy $E_{\mathrm{CEF}}^{l}$ as

$\displaystyle E_{\mathrm{CEF}}^{l}=\int
d\boldsymbol{r}\psi_{l}^{\ast}\left(\boldsymbol{r}\right)H_{\mathrm{CEF}}\psi_{l}\left(\boldsymbol{r}\right).$
(16)

Here, $l$ denotes the 3$d$ orbital suffix. The point charge model of the CEF
Hamiltonian gives the low-lying singlet states of $\psi_{3z^{2}-r^{2}}$ with
$A_{1}$ symmetry and $\psi_{x^{\prime 2}-y^{\prime 2}}$ with $B_{1}$ symmetry
and the mid-lying doublet state of $\psi_{y^{\prime}z}$ and
$\psi_{zx^{\prime}}$ with $E$ symmetry and the excited singlet of
$\psi_{x^{\prime}y^{\prime}}$ with $B_{2}$ symmetry [43]. Actually, the band
calculations show that the three hole sheets around the $\Gamma$-point consist
of the doublet state of $\psi_{y^{\prime}z}$ and $\psi_{zx^{\prime}}$ and the
singlet state $\psi_{x^{\prime 2}-y^{\prime 2}}$.

In the present investigation, we suppose that the structural and
superconducting transitions of the iron pnictide manifest themselves as a
result of the spontaneous symmetry breaking associated with the degenerate
$y^{\prime}z$ and $zx^{\prime}$ bands, which are mapped on the special unitary
group SU(2). In order to explore the quadrupole-strain interaction in the iron
pnictide, we specially treat the degenerate $\psi_{y^{\prime}z}$ and
$\psi_{zx^{\prime}}$ orbitals of the $E$ symmetry with half filling as [50,
51]

$\displaystyle\psi_{y^{\prime}z}\left(\boldsymbol{r}\right)$
$\displaystyle=\frac{i}{\sqrt{2}}f_{d}\left(r\right)\left[Y_{2}^{1}\left(\theta,\varphi\right)+Y_{2}^{-1}\left(\theta,\varphi\right)\right]$
$\displaystyle=\sqrt{\frac{15}{4\pi}}f_{d}\left(r\right)\frac{y^{\prime}z}{r^{2}},$
(17) $\displaystyle\psi_{zx^{\prime}}\left(\boldsymbol{r}\right)$
$\displaystyle=\frac{-1}{\sqrt{2}}f_{d}\left(r\right)\left[Y_{2}^{1}\left(\theta,\varphi\right)-Y_{2}^{-1}\left(\theta,\varphi\right)\right]$
$\displaystyle=\sqrt{\frac{15}{4\pi}}f_{d}\left(r\right)\frac{zx^{\prime}}{r^{2}}.$
(18)

Here, the polar coordinate $\left(r,\theta,\varphi\right)$ represents the
position vector $\boldsymbol{r}=\left(x^{\prime},y^{\prime},z\right)$ of an
electron. $f_{d}\left(r\right)$ is the radius function of a 3$d$ electron, and
$Y_{2}^{\pm 1}\left(\theta,\varphi\right)$ is the spherical harmonics with
angular momentum $l=2$ and azimuthal quantum number $m=\pm 1$.

Since the direct product of the $E$ doublet is reduced as $E\otimes
E=A_{1}\oplus A_{2}\oplus B_{1}\oplus B_{2}$, we deduce that the degenerate
$y^{\prime}z$ and $zx^{\prime}$ orbitals carry the electric quadrupoles
$O_{3z^{2}-r^{2}}=\bigl{(}3z^{2}-r^{2}\bigr{)}/\bigl{(}2r^{2}\bigr{)}$ with
$A_{1}$ symmetry,
$O_{x^{\prime}y^{\prime}}=\sqrt{3}x^{\prime}y^{\prime}/r^{2}$ with $B_{1}$
symmetry, and $O_{x^{\prime 2}-y^{\prime 2}}=\sqrt{3}\bigl{(}x^{\prime
2}-y^{\prime 2}\bigr{)}/\bigl{(}2r^{2}\bigr{)}$ with $B_{2}$ symmetry and the
angular momentum $l_{z}=-i\bigl{(}x^{\prime}\partial/\partial
y^{\prime}-y^{\prime}\partial/\partial
x^{\prime}\bigr{)}=-i\partial/\partial\varphi$ with $A_{2}$ symmetry. The
quadrupole $O_{x^{\prime 2}-y^{\prime 2}}$ couples to the strain
$\varepsilon_{xy}$ of the transverse ultrasonic wave of $C_{66}$, while
$O_{x^{\prime}y^{\prime}}$ couples to the strain $\varepsilon_{x^{2}-y^{2}}$
of the transverse ultrasonic wave of $\left(C_{11}-C_{12}\right)/2$. The CEF
Hamiltonian of Eq. (4.2) includes $O_{3z^{2}-r^{2}}$ with full symmetry.

Employing the identity matrix $\tau_{0}$ and the Pauli matrices $\tau_{x}$,
$\tau_{y}$, and $\tau_{z}$ corresponding to the generator elements of the
special unitary group SU(2), we present the quadrupoles $O_{3z^{2}-r^{2}}$,
$O_{x^{\prime 2}-y^{\prime 2}}$, and $O_{x^{\prime}y^{\prime}}$ and the
angular momentum $l_{z}$ in terms of the matrices for the orbital state
$\psi_{y^{\prime}z}$ of Eq. (4.2) and $\psi_{zx^{\prime}}$ of Eq. (4.2) as

$\displaystyle O_{3z^{2}-r^{2}}$
$\displaystyle=\frac{1}{7}\bordermatrix{&\psi_{y^{\prime}z}&\psi_{zx^{\prime}}\cr&1&0\cr&0&1\cr}=\frac{1}{7}\tau_{0},$
(19) $\displaystyle O_{x^{\prime}y^{\prime}}$
$\displaystyle=\frac{\sqrt{3}}{7}\left(\begin{array}[]{cc}0&1\\\ 1&0\\\
\end{array}\right)=\frac{\sqrt{3}}{7}\tau_{x},$ (22) $\displaystyle l_{z}$
$\displaystyle=-\left(\begin{array}[]{cc}0&-i\\\ i&0\\\
\end{array}\right)=-\tau_{y},$ (25) $\displaystyle O_{x^{\prime 2}-y^{\prime
2}}$ $\displaystyle=-\frac{\sqrt{3}}{7}\left(\begin{array}[]{cc}1&0\\\ 0&-1\\\
\end{array}\right)=-\frac{\sqrt{3}}{7}\tau_{z}.$ (28)

The Pauli matrices $\tau_{x}$, $\tau_{y}$, and $\tau_{z}$ obey the commutation
relation

$\displaystyle\left[\tau_{i},\tau_{j}\right]=2i\varepsilon_{ijk}\tau_{k}\ \ \
\left(i,j,k=x,y,z\right).$ (29)

Here, $\varepsilon_{ijk}$ is the Levi$-$Civita symbol. As will be shown in
Sects. 4.6 and 4.10, the commutation relation among the quadrupoles
$O_{x^{\prime 2}-y^{\prime 2}}$ and $O_{x^{\prime}y^{\prime}}$ and the angular
momentum $l_{z}$ of Eq. (29) brings about quantum fluctuations, which play a
significant role in the manifestation of the specific superconductivity
accompanying the hexadecapole ordering in the vicinity of the QCP.

Figure 5:  (Color online) Schematic view of interactions between electric
quadrupole $O_{x^{\prime}y^{\prime}}$ and strain $\varepsilon_{x^{2}-y^{2}}$
with $B_{1}$ symmetry in (a), $O_{x^{\prime 2}-y^{\prime 2}}$ and
$\varepsilon_{xy}$ with $B_{2}$ symmetry in (b), and electric hexadecapole
$H_{z}^{\alpha}\propto\partial l_{z}/\partial t$ and rotation $\omega_{xy}$
with $A_{2}$ symmetry in (c) of the space group $D_{4h}^{17}$ in the present
iron pnictide. The coordinates of $x$ and $y$ ($x^{\prime}$ and $y^{\prime}$)
are adopted for the neighboring Ba-Ba (Fe-Fe) directions.

The elastic strains $\varepsilon_{\mathrm{\Gamma}_{\gamma}}$ induced by the
thermally excited transverse acoustic phonons or the experimentally generated
transverse ultrasonic waves perturb the CEF Hamiltonian as

$\displaystyle
H_{\mathrm{CEF}}\left(\boldsymbol{r},\varepsilon_{\mathrm{\Gamma}_{\gamma}}\right)=H_{\mathrm{CEF}}\left(\boldsymbol{r}\right)+\sum_{\mathrm{\Gamma}_{\gamma}}\frac{\partial
H_{\mathrm{CEF}}\left(\boldsymbol{r}\right)}{\partial\varepsilon_{\mathrm{\Gamma}_{\gamma}}}\varepsilon_{{\mathrm{\Gamma}}_{\gamma}}.$
(30)

Because electrons accommodated in the degenerate $y^{\prime}z$ and
$zx^{\prime}$ orbitals bear the electric quadrupoles, the second term of Eq.
(30) proportional to the strain $\varepsilon_{\mathrm{\Gamma}_{\gamma}}$ is
expressed in terms of the quadrupole-strain interaction as [52, 53, 54]

$\displaystyle H_{\mathrm{QS}}=-g_{x^{\prime 2}-y^{\prime 2}}O_{x^{\prime
2}-y^{\prime
2}}\varepsilon_{xy}-g_{x^{\prime}y^{\prime}}O_{x^{\prime}y^{\prime}}\varepsilon_{x^{2}-y^{2}}.$
(31)

As the present ultrasonic experiments give a large coupling constant of
$g_{x^{\prime 2}-y^{\prime 2}}\sim 1\times 10^{3}$ K per electron, the first
term $-g_{x^{\prime 2}-y^{\prime 2}}O_{x^{\prime 2}-y^{\prime
2}}\varepsilon_{xy}$ in Eq. (31) is important for the critical divergence of
the attenuation coefficient $\alpha_{66}$ and the softening of the elastic
constant $C_{66}$ around the structural transition. The little softening of
$\left(C_{11}-C_{12}\right)/2$ in the present experiments implies that the
second term
$-g_{x^{\prime}y^{\prime}}O_{x^{\prime}y^{\prime}}\varepsilon_{x^{2}-y^{2}}$
in Eq. (31) scarcely affects the system. We schematically show the quadrupole
$O_{x^{\prime}y^{\prime}}$ conjugate with the strain
$\varepsilon_{x^{2}-y^{2}}$ in Fig. 5(a) and the quadrupole $O_{x^{\prime
2}-y^{\prime 2}}$ conjugate with $\varepsilon_{xy}$ in Fig. 5(b).

Next, we consider the interaction of the rotation $\omega_{xy}$ associated
with the transverse ultrasonic wave with one-electron states having the
angular momentum $l_{z}$ in the CEF Hamiltonian. The rotation operator of
$\exp\bigl{[}-il_{z}\omega_{xy}\bigr{]}$ twists the $\psi_{y^{\prime}z}$ and
$\psi_{zx^{\prime}}$ orbital states in the CEF Hamiltonian given by Eq. (4.2)
as [55]

$\displaystyle\left\langle
l\left|H_{\mathrm{CEF}}\bigl{(}\omega_{xy}\bigr{)}\right|l^{\prime}\right\rangle=\int
d\boldsymbol{r}\psi_{l}^{\ast}\left(\boldsymbol{r}\right)e^{il_{z}\omega_{xy}}H_{\mathrm{CEF}}e^{-il_{z}\omega_{xy}}\psi_{l^{\prime}}\left(\boldsymbol{r}\right).$
(32)

Here, ket $|l\rangle=\psi_{l}\left(\boldsymbol{r}\right)$ stands for the
orbital states $l=y^{\prime}z$ of Eq. (4.2) and $zx^{\prime}$ of Eq. (4.2).
Supposing infinitesimal rotation of $\omega_{xy}\ll 1$ associated with the
ultrasound in experiments, the rotation-operated Hamiltonian of Eq. (32) is
expanded in terms of power series up to the second order of $\omega_{xy}$ as

$\displaystyle
e^{il_{z}\omega_{xy}}H_{\mathrm{CEF}}e^{-il_{z}\omega_{xy}}\approx$
$\displaystyle
H_{\mathrm{CEF}}+i\left[l_{z},H_{\mathrm{CEF}}\right]\omega_{xy}$
$\displaystyle-\frac{1}{2}\left[l_{z},\bigl{[}l_{z},H_{\mathrm{CEF}}\bigr{]}\right]\bigl{(}\omega_{xy}\bigr{)}^{2}.$
(33)

Here, $\left[A,B\right]=AB-BA$ denotes the Poisson bracket. Applying the
differential operation of $l_{z}=-i\left(x^{\prime}\partial/\partial
y^{\prime}-y^{\prime}\partial/\partial
x^{\prime}\right)=-i\partial/\partial\varphi$ to the CEF Hamiltonian of Eq.
(4.2), we obtain the perturbation Hamiltonian as

$\displaystyle
H_{\mathrm{rot}}^{\mathrm{CEF}}\bigl{(}\omega_{xy}\bigr{)}=4A_{4}^{4}H_{z}^{\alpha}\omega_{xy}-8A_{4}^{4}H_{4}^{4}\bigl{(}\omega_{xy}\bigr{)}^{2}.$
(34)

The first term of Eq. (34) represents the linear coupling of the rotation
$\omega_{xy}$ to the electric hexadecapole
$H_{z}^{\alpha}=\sqrt{35}x^{\prime}y^{\prime}\bigl{(}x^{\prime 2}-y^{\prime
2}\bigr{)}/\bigl{(}2r^{4}\bigr{)}$ with $A_{2}$ symmetry of the $D_{2d}$ point
group [56]. In Fig. 5(c), we schematically show the hexadecapole
$H_{z}^{\alpha}$ carried by the one-electron state and the rotation
$\omega_{xy}$ in the interaction Hamiltonian
$H_{\mathrm{rot}}^{\mathrm{CEF}}\bigl{(}\omega_{xy}\bigr{)}$ of Eq. (34). The
second term in Eq. (34) is the coupling of the square of the rotation
$\bigl{(}\omega_{xy}\bigr{)}^{2}$ to the electric hexadecapole
$H_{4}^{4}=\sqrt{35}\bigl{(}x^{\prime 4}-6x^{\prime 2}y^{\prime 2}+y^{\prime
4}\bigr{)}/\bigl{(}8r^{4}\bigr{)}$ with full symmetry. The hexadecapole
$H_{4}^{4}$ has already appeared in the CEF Hamiltonian of Eq. (4.2).

The coefficient proportional to the rotation $\omega_{xy}$ in Eq. (4.2) is
identified with the Heisenberg equation of the angular momentum $l_{z}$ for
the CEF Hamiltonian of Eq. (4.2). This gives the time derivative of the
angular momentum $l_{z}$, which is identical to the torque $\tau_{xy}$ for a
quantum system, as

$\displaystyle\ i\hbar\frac{\partial l_{z}}{\partial
t}=\left[l_{z},H_{\mathrm{CEF}}\right]=-ig_{\mathrm{H}}H_{z}^{\alpha}=i\tau_{xy}.$
(35)

Note that the coupling parameter of $g_{\mathrm{H}}=4A_{4}^{4}$ in Eq. (35)
has already been given in the CEF Hamiltonian of Eq. (4.2). Since the angular
momentum $l_{z}$ is constant for the motion of the one-electron states in the
CEF Hamiltonian of Eq. (4.2), $l_{z}$ commutes with $H_{\mathrm{CEF}}$ as
$\left[l_{z},H_{\mathrm{CEF}}\right]=0$. Consequently, the hexadecapole
$H_{z}^{\alpha}$ proportional to the time derivative of the angular momentum
$\partial l_{z}/\partial t$ of Eq. (35) is expected to vanish. This is well
known as the rotational invariance for the electronic states in a central
force potential such as the CEF potential. The rotation invariance can be
naturally understood by the fact that the hexadecapole interaction of Eq.
(4.2) does not change the energy of an electron moving along the contour line
of the CEF potential.

The rotation $\omega_{xy}$ due to a transverse ultrasonic wave should also
affect the quadrupole-strain interaction of Eq. (31) within the bilinear
coupling between the strain and rotation as

$\displaystyle H_{\mathrm{rot}}^{\mathrm{QS}}\bigl{(}\omega_{xy}\bigr{)}$
$\displaystyle=i\left[l_{z},H_{\mathrm{QS}}\right]\omega_{xy}$
$\displaystyle=2g_{x^{\prime 2}-y^{\prime
2}}O_{x^{\prime}y^{\prime}}\varepsilon_{xy}\omega_{xy}$
$\displaystyle\qquad-2g_{x^{\prime}y^{\prime}}O_{x^{\prime 2}-y^{\prime
2}}\varepsilon_{x^{2}-y^{2}}\omega_{xy}.$ (36)

Here, the commutation relation among the Pauli matrices of Eq. (29) is
employed. The perturbation energies for the orbital states $l=y^{\prime}z$ and
$zx^{\prime}$ with energy $E_{l}^{0}$ due to the perturbation Hamiltonian of
Eq. (4.2) are written as

$\displaystyle E_{l}\bigl{(}\omega_{xy}\bigr{)}$
$\displaystyle=E_{l}^{0}+\left\langle
l\left|H_{\mathrm{rot}}^{\mathrm{QS}}\bigl{(}\omega_{xy}\bigr{)}\right|l\right\rangle$
$\displaystyle=E_{l}^{0}+2g_{x^{\prime 2}-y^{\prime 2}}\left\langle
l\left|O_{x^{\prime}y^{\prime}}\right|l\right\rangle\varepsilon_{xy}\omega_{xy}$
$\displaystyle\qquad\qquad-2g_{x^{\prime}y^{\prime}}\left\langle
l\left|O_{x^{\prime 2}-y^{\prime
2}}\right|l\right\rangle\varepsilon_{x^{2}-y^{2}}\omega_{xy}.$ (37)

The null of the diagonal elements of the Pauli matrix of
$\boldsymbol{O}_{x^{\prime}y^{\prime}}$ of Eq. (22) diminishes the second term
$2g_{x^{\prime 2}-y^{\prime 2}}\left\langle
l\left|O_{x^{\prime}y^{\prime}}\right|l\right\rangle$ in the second line in
Eq. (4.2). Although the matrix of $\boldsymbol{O}_{x^{\prime 2}-y^{\prime 2}}$
of Eq. (28) possesses diagonal elements, we may ignore the third term
$-2g_{x^{\prime}y^{\prime}}\left\langle l\left|O_{x^{\prime 2}-y^{\prime
2}}\right|l\right\rangle$ in the second line in Eq. (4.2) because the little
softening of $(C_{11}-C_{12})/2$ indicates a small quadrupole-strain coupling
constant of $g_{x^{\prime}y^{\prime}}<<1$. Presumably, the strain-rotation
bilinear coupling terms of Eq. (4.2) have little effect on the present system
with degenerate $y^{\prime}z$ and $zx^{\prime}$ orbitals.

The applied magnetic field, however, breaks the rotational invariance for the
one-electron states in the CEF potential. This effect of rotation on the
elastic properties in applied magnetic fields has been theoretically
investigated [52, 53] and experimentally verified by ultrasonic measurements
on the antiferro-magnetic compound MnF2 [57] and the rare-earth compounds TmSb
[58], CeAl2 [59], CeB6 [60], Ce0.5La0.5B6 [61], and HoVO4 [55]. Furthermore,
the effect of rotation on the quantum oscillation of elastic constants of
transverse ultrasonic waves for Fermi surfaces with an ellipsoidal shape has
been discussed [62]. In Sect. 4.6, we show the coupling of the rotation
associated with the transverse acoustic phonons to the hexadecapole carried by
two-electron states.

###  4.3 Quadrupole-strain interaction

As shown in Eq. (20), the quadrupoles carried by electrons in the CEF
Hamiltonian interact to the strains of transverse ultrasonic waves. This
quadrupole-strain interaction brings about the sizable softening of $C_{66}$
and the divergence of the attenuation coefficients $\alpha_{66}$ and
$\alpha_{\mathrm{L}[110]}$ around the structural transition. In order to
properly describe the structural transition of the iron pnictide, we will
consider the interaction of the quadrupole $O_{x^{\prime 2}-y^{\prime
2}}\left(\boldsymbol{r}_{i}\right)$ of band electrons to the transverse
acoustic phonons carrying the strain
$\varepsilon_{xy}\left(\boldsymbol{r}_{i}\right)$ for the position vector
$\boldsymbol{r}_{i}$ per electron as

$\displaystyle H_{\mathrm{QS}}=-g_{x^{\prime 2}-y^{\prime
2}}\sum_{l_{1},l_{2}=y^{\prime}z,zx^{\prime}}\sum_{\sigma=\uparrow,\downarrow}$
$\displaystyle\int
d\boldsymbol{r}_{i}\mathit{\Psi}_{l_{1},\sigma}^{\ast}\left(\boldsymbol{r}_{i}\right)O_{x^{\prime
2}-y^{\prime 2}}\left(\boldsymbol{r}_{i}\right)$
$\displaystyle\times\mathit{\Psi}_{l_{2},\sigma}\left(\boldsymbol{r}_{i}\right)\varepsilon_{x^{\prime
2}-y^{\prime 2}}\left(\boldsymbol{r}_{i}\right).$ (38)

The electron field operators of
$\mathit{\Psi}_{y^{\prime}z,\sigma}\left(\boldsymbol{r}_{i}\right)$ and
$\mathit{\Psi}_{zx^{\prime},\sigma}\left(\boldsymbol{r}_{i}\right)$ at
position $\boldsymbol{r}_{i}$ acting on the degenerate $y^{\prime}z$ and
$zx^{\prime}$ bands with spin orientation $\sigma$ in Eq. (4.3) are written as
[K.Yoshida]

$\displaystyle\mathit{\Psi}_{l,\sigma}\left(\boldsymbol{r}_{i}\right)=d_{i,l,\sigma}\psi_{l}\left(\boldsymbol{r}_{i}\right)v_{\sigma}\left(\boldsymbol{r}_{i}\right),$
(39)
$\displaystyle\mathit{\Psi}_{l,\sigma}^{\ast}\left(\boldsymbol{r}_{i}\right)=d_{i,l,\sigma}^{\dagger}\psi_{l}^{\ast}\left(\boldsymbol{r}_{i}\right)v_{\sigma}^{\ast}\left(\boldsymbol{r}_{i}\right).$
(40)

Here, $l$ denotes the band suffix of $y^{\prime}z$ and $zx^{\prime}$, and
$v_{\sigma}\left(\boldsymbol{r}_{i}\right)$ is the spin function of
$\alpha\left(\boldsymbol{r}_{i}\right)$ for up-spin $\sigma=\uparrow$ or
$\beta\left(\boldsymbol{r}_{i}\right)$ for down-spin $\sigma=\downarrow$. The
annihilation operator $d_{i,l,\sigma}$ and creation operator
$d_{i,l,\sigma}^{\dagger}$ are expressed by Fourier transforms as below

$\displaystyle d_{i,l,\sigma}$
$\displaystyle=\sum_{\boldsymbol{k}}d_{\boldsymbol{k},l,\sigma}e^{i\boldsymbol{k}\cdot\boldsymbol{r}_{i}},$
(41) $\displaystyle d_{i,l,\sigma}^{\dagger}$
$\displaystyle=\sum_{\boldsymbol{k}}d_{\boldsymbol{k},l,\sigma}^{\dagger}e^{-i\boldsymbol{k}\cdot\boldsymbol{r}_{i}}.$
(42)

Here, $\boldsymbol{k}$ is the wavevector of an electron. We express the
Hamiltonian $H_{\mathrm{K}}$ for the electrons accommodated in the degenerate
$y^{\prime}z$ and $zx^{\prime}$ bands, which play a significant role in the
appearance of the quadrupole ordering in the system as

$\displaystyle H_{\mathrm{K}}$
$\displaystyle=\sum_{\boldsymbol{k}}\sum_{\sigma}\left[\varepsilon_{y^{\prime}z}\left(\boldsymbol{k}\right)d_{\boldsymbol{k},y^{\prime}z,\sigma}^{\dagger}d_{\boldsymbol{k},y^{\prime}z,\sigma}\right.$
$\displaystyle\left.\qquad\qquad\qquad+\varepsilon_{zx^{\prime}}\left(\boldsymbol{k}\right)d_{\boldsymbol{k},zx^{\prime},\sigma}^{\dagger}d_{\boldsymbol{k},zx^{\prime},\sigma}\right].$
(43)

Here, $\varepsilon_{l}\left(\boldsymbol{k}\right)$ for suffix $l=y^{\prime}z$
or $zx^{\prime}$ is the energy from the Fermi level.

The quadrupole-strain interaction of Eq. (4.3) is rewritten in terms of
Fourier transforms as

$\displaystyle H_{\mathrm{QS}}=-G_{x^{\prime 2}-y^{\prime
2}}\sum_{\boldsymbol{k},\boldsymbol{q}}O_{x^{\prime 2}-y^{\prime
2},\boldsymbol{k},\boldsymbol{q}}\varepsilon_{xy}\left(\boldsymbol{q}\right).$
(44)

Here, we use the coupling constant $G_{x^{\prime 2}-y^{\prime
2}}=\sqrt{3}g_{x^{\prime 2}-y^{\prime 2}}/7$. The interaction Hamiltonian of
Eq. (44) means that the electrons with wavevector $\boldsymbol{k}$ bearing the
quadrupole $O_{x^{\prime 2}-y^{\prime 2},\boldsymbol{k},\boldsymbol{q}}$ are
scattered by the strain $\varepsilon_{xy}\left(\boldsymbol{q}\right)$ of the
transverse acoustic phonons with wavevector $\boldsymbol{q}$. The quadrupole
$O_{x^{\prime 2}-y^{\prime 2},\boldsymbol{k},\boldsymbol{q}}$ in Eq. (44) is
described in terms of the annihilation and creation operators of Eqs. (41) and
(42) as

$\displaystyle O_{x^{\prime 2}-y^{\prime
2},\boldsymbol{k},\boldsymbol{q}}=\sum_{\sigma}\left(-d_{\boldsymbol{k}+\boldsymbol{q},y^{\prime}z,\sigma}^{\dagger}d_{\boldsymbol{k},y^{\prime}z,\sigma}+d_{\boldsymbol{k}+\boldsymbol{q},zx^{\prime},\sigma}^{\dagger}d_{\boldsymbol{k},zx^{\prime},\sigma}\right).$
(45)

Note that the quadrupole $O_{x^{\prime 2}-y^{\prime
2},\boldsymbol{k},\boldsymbol{q}}$ in Eq. (45) is simply described by the
difference between the occupation numbers for the degenerate $y^{\prime}z$ and
$zx^{\prime}$ bands at the long-wavelength limit of
$\left|\boldsymbol{q}\right|=2\pi/\lambda\rightarrow 0$ of the transverse
acoustic phonons carrying the strain
$\varepsilon_{xy}\left(\boldsymbol{q}\right)$ as

$\displaystyle O_{x^{\prime 2}-y^{\prime 2},\boldsymbol{k},\boldsymbol{q}=0}$
$\displaystyle=\sum_{\sigma}\left(-d_{\boldsymbol{k},y^{\prime}z,\sigma}^{\dagger}d_{\boldsymbol{k},y^{\prime}z,\sigma}+d_{\boldsymbol{k},zx^{\prime},\sigma}^{\dagger}d_{\boldsymbol{k},zx^{\prime},\sigma}\right)$
$\displaystyle=\sum_{\sigma}\left(-n_{\boldsymbol{k},y^{\prime}z,\sigma}+n_{\boldsymbol{k},zx^{\prime},\sigma}\right).$
(46)

In Sect. 4.6, we will use the quadrupole
$O_{x^{\prime}y^{\prime},\boldsymbol{k},\boldsymbol{q}}$ expressed as

$\displaystyle
O_{x^{\prime}y^{\prime},\boldsymbol{k},\boldsymbol{q}}=\sum_{\sigma}\left(d_{\boldsymbol{k}+\boldsymbol{q},y^{\prime}z,\sigma}^{\dagger}d_{\boldsymbol{k},zx^{\prime},\sigma}+d_{\boldsymbol{k}+\boldsymbol{q},zx^{\prime},\sigma}^{\dagger}d_{\boldsymbol{k},y^{\prime}z,\sigma}\right).$
(47)

The strain $\varepsilon_{xy}\left(\boldsymbol{q}\right)$ in Eq. (44) is
expressed by the annihilation and creation operators of phonons defined in Eq.
(4.1) as

$\displaystyle\varepsilon_{xy}\left(\boldsymbol{q}\right)=$
$\displaystyle\frac{i}{2}\sqrt{\mathstrut\frac{\hbar}{2V\rho_{\mathrm{M}}\omega_{y}\left(\boldsymbol{q}\right)}}q_{x}\left(a_{y,\boldsymbol{q}}-a_{y,-\boldsymbol{q}}^{\dagger}\right)$
$\displaystyle+\frac{i}{2}\sqrt{\mathstrut\frac{\hbar}{2V\rho_{\mathrm{M}}\omega_{x}\left(\boldsymbol{q}\right)}}q_{y}\left(a_{x,\boldsymbol{q}}-a_{x,-\boldsymbol{q}}^{\dagger}\right).$
(48)

The strain $\varepsilon_{xy}$ induced by a transverse ultrasonic wave with
frequency as low as 100 MHz is identified with
$\varepsilon_{xy}\left(\boldsymbol{q}\right)$ of Eq. (4.3) in the long-
wavelength limit of $\left|\boldsymbol{q}\right|=2\pi/\lambda\rightarrow 0$ as
$\varepsilon_{xy}=\varepsilon_{xy}\left(\boldsymbol{q}=0\right)$. In the case
of the quadrupole-strain interaction of Eq. (44) mediated by thermally excited
transverse acoustic phonons, the emission of phonons carrying the strain
$\varepsilon_{xy}\left(\boldsymbol{q}\right)$ scatters the electron state of
momentum $\hbar\boldsymbol{k}$ to the state of
$\hbar\boldsymbol{k}+\hbar\boldsymbol{q}$ and the absorption of phonons
scatters the electron state of $\hbar\boldsymbol{k}$ to the state of
$\hbar\boldsymbol{k}-\hbar\boldsymbol{q}$.

The canonical transformation for the electron-phonon scattering processes by
the quadrupole-strain interaction $H_{\mathrm{QS}}$ of Eq. (44) provides us
with the electron interaction mediated by the transverse acoustic phonons as
[45]

$\displaystyle H_{\mathrm{ind}}^{\mathrm{QQ}}=$
$\displaystyle-\sum_{\boldsymbol{k},\boldsymbol{k}^{\prime},\boldsymbol{q}}\sum_{\sigma,\sigma^{\prime}}D_{y^{\prime}z}^{\mathrm{QQ}}\left(\boldsymbol{k},\boldsymbol{q}\right)d_{\boldsymbol{k}-\boldsymbol{q},y^{\prime}z,\sigma}^{\dagger}d_{\boldsymbol{k},y^{\prime}z,\sigma}d_{\boldsymbol{k}^{\prime}+\boldsymbol{q},y^{\prime}z,\sigma^{\prime}}^{\dagger}d_{\boldsymbol{k}^{\prime},y^{\prime}z,\sigma^{\prime}}$
$\displaystyle+\sum_{\boldsymbol{k},\boldsymbol{k}^{\prime},\boldsymbol{q}}\sum_{\sigma,\sigma^{\prime}}D_{y^{\prime}z}^{\mathrm{QQ}}\left(\boldsymbol{k},\boldsymbol{q}\right)d_{\boldsymbol{k}-\boldsymbol{q},y^{\prime}z,\sigma}^{\dagger}d_{\boldsymbol{k},y^{\prime}z,\sigma}d_{\boldsymbol{k}^{\prime}+\boldsymbol{q},zx^{\prime},\sigma^{\prime}}^{\dagger}d_{\boldsymbol{k}^{\prime},zx^{\prime},\sigma^{\prime}}$
$\displaystyle+\sum_{\boldsymbol{k},\boldsymbol{k}^{\prime},\boldsymbol{q}}\sum_{\sigma,\sigma^{\prime}}D_{zx^{\prime}}^{\mathrm{QQ}}\left(\boldsymbol{k},\boldsymbol{q}\right)d_{\boldsymbol{k}-\boldsymbol{q},zx^{\prime},\sigma}^{\dagger}d_{\boldsymbol{k},zx^{\prime},\sigma}d_{\boldsymbol{k}^{\prime}+\boldsymbol{q},y^{\prime}z,\sigma^{\prime}}^{\dagger}d_{\boldsymbol{k}^{\prime},y^{\prime}z,\sigma^{\prime}}$
$\displaystyle-\sum_{\boldsymbol{k},\boldsymbol{k}^{\prime},\boldsymbol{q}}\sum_{\sigma,\sigma^{\prime}}D_{zx^{\prime}}^{\mathrm{QQ}}\left(\boldsymbol{k},\boldsymbol{q}\right)d_{\boldsymbol{k}-\boldsymbol{q},zx^{\prime},\sigma}^{\dagger}d_{\boldsymbol{k},zx^{\prime},\sigma}d_{\boldsymbol{k}^{\prime}+\boldsymbol{q},zx^{\prime},\sigma^{\prime}}^{\dagger}d_{\boldsymbol{k}^{\prime},zx^{\prime},\sigma^{\prime}}.$
(49)

Here, the coupling coefficient
$D_{l}^{\mathrm{QQ}}\left(\boldsymbol{k},\boldsymbol{q}\right)$ in Eq. (4.3)
for $l=y^{\prime}z$ and $zx^{\prime}$ is given by

$\displaystyle
D_{l}^{\mathrm{QQ}}\left(\boldsymbol{k},\boldsymbol{q}\right)=-\frac{1}{2}G_{x^{\prime
2}-y^{\prime 2}}^{2}$
$\displaystyle\times\left\\{\frac{\hbar}{2V\rho_{\mathrm{M}}\omega_{y}\left(\boldsymbol{q}\right)}q_{x}^{2}\frac{\hbar\omega_{y}\left(\boldsymbol{q}\right)}{\left[\varepsilon_{l}\left(\boldsymbol{k}\right)-\varepsilon_{l}\left(\boldsymbol{k}-\boldsymbol{q}\right)\right]^{2}-\hbar^{2}\omega_{y}\left(\boldsymbol{q}\right)^{2}}\right.$
$\displaystyle\left.+\frac{\hbar}{2V\rho_{\mathrm{M}}\omega_{x}\left(\boldsymbol{q}\right)}q_{y}^{2}\times\frac{\hbar\omega_{x}\left(\boldsymbol{q}\right)}{\left[\varepsilon_{l}\left(\boldsymbol{k}\right)-\varepsilon_{l}\left(\boldsymbol{k}-\boldsymbol{q}\right)\right]^{2}-\hbar^{2}\omega_{x}\left(\boldsymbol{q}\right)^{2}}\right\\}.$
(50)

The four independent scattering processes involving virtually excited one-
phonon states due to the strain $\varepsilon_{xy}\left(\boldsymbol{q}\right)$
in Eq. (44) are schematically pictured in Fig. 6. The processes in Fig. 6(a)
indicate the scattering of the electrons accommodated in the same band
$y^{\prime}z$ and those in Fig. 6(d) indicate the scattering in the same band
$zx^{\prime}$. The processes in Figs. 6(b) and 6(c) indicate the scattering
between two electrons in the different bands $y^{\prime}z$ and $zx^{\prime}$.

Figure 6:  The indirect quadrupole interaction
$H_{\mathrm{ind}}^{\mathrm{QQ}}$ of Eq. (4.3) gives four scattering processes
between the one-electron states through the strain
$\varepsilon_{xy}\left(\boldsymbol{q}\right)$ of a transverse acoustic phonon,
while indirect hexadecapole interaction $H_{\mathrm{ind}}^{\mathrm{HH}}$ of
Eq. (4.7) gives four scattering processes between the two-electron states
through the rotation $\omega_{xy}\left(\boldsymbol{q}\right)$ of a transverse
acoustic phonon.

Supposing the identity of the energy
$\varepsilon_{y^{\prime}z}\left(\boldsymbol{k}\right)=\varepsilon_{zx^{\prime}}\left(\boldsymbol{k}\right)$
for electrons with a small wavevector $\boldsymbol{k}$ located near the Fermi
level, we take the indirect quadrupole interaction coefficient
$D_{y^{\prime}z}^{\mathrm{QQ}}\left(\boldsymbol{k},\boldsymbol{q}\right)=D_{zx^{\prime}}^{\mathrm{QQ}}\left(\boldsymbol{k},\boldsymbol{q}\right)=D^{\mathrm{QQ}}\left(\boldsymbol{k},\boldsymbol{q}\right)$
in the tetragonal phase. Using the quadrupole $O_{x^{\prime 2}-y^{\prime
2},\boldsymbol{k},\boldsymbol{q}}$ expressed in Eq. (45), we reduce the
quadrupole interaction Hamiltonian of Eq. (4.3) as

$\displaystyle
H_{\mathrm{ind}}^{\mathrm{QQ}}=-\sum_{\boldsymbol{k},\boldsymbol{k}^{\prime},\boldsymbol{q}}D^{\mathrm{QQ}}\left(\boldsymbol{k},\boldsymbol{q}\right)O_{x^{\prime
2}-y^{\prime 2},\boldsymbol{k},-\boldsymbol{q}}O_{x^{\prime 2}-y^{\prime
2},\boldsymbol{k}^{\prime},\boldsymbol{q}}.$ (51)

In order to understand the frequency dependence of the ultrasonic attenuation
of $\alpha_{\mathrm{L}[110]}$ in Fig. 3(a), we will examine the indirect
quadrupole interaction coefficient
$D^{\mathrm{QQ}}\left(\boldsymbol{k},\boldsymbol{q}\right)$ in Eq. (51), which
strongly depends on the magnitude of wavevectors $\boldsymbol{k}$ and
$\boldsymbol{k}^{\prime}$ for the interacting electrons and $\boldsymbol{q}$
for the transverse acoustic phonons participating in the scattering. Taking
the electron energies of
$\varepsilon_{y^{\prime}z}\left(\boldsymbol{k}\right)=\varepsilon_{zx^{\prime}}\left(\boldsymbol{k}\right)=\varepsilon\left(\boldsymbol{k}\right)=\hbar^{2}\left|\boldsymbol{k}\right|^{2}/2m_{e}$
with effective electron mass $m_{e}$ and transverse acoustic phonon energies
of
$\hbar\omega_{x}\left(\boldsymbol{q}\right)=\hbar\omega_{y}\left(\boldsymbol{q}\right)=\hbar\omega\left(\boldsymbol{q}\right)=\hbar
v_{66}|\boldsymbol{q}|$ with ultrasonic velocity $v_{66}$, we obtain the
indirect quadrupole interaction coefficient
$D^{\mathrm{QQ}}\left(\boldsymbol{k},\boldsymbol{q}\right)$ of Eq. (51) as

$\displaystyle D^{\mathrm{QQ}}\left(\boldsymbol{k},\boldsymbol{q}\right)$
$\displaystyle=-G_{x^{\prime 2}-y^{\prime
2}}^{2}\frac{\hbar}{2V\rho_{\mathrm{M}}\omega\left(\boldsymbol{q}\right)}q^{2}$
$\displaystyle\qquad\quad\times\frac{\hbar\omega\left(\boldsymbol{q}\right)}{\left[\varepsilon\left(\boldsymbol{k}\right)-\varepsilon\left(\boldsymbol{k}-\boldsymbol{q}\right)\right]^{2}-\hbar^{2}\omega\left(\boldsymbol{q}\right)^{2}}.$
(52)

In the small-wavenumber limit of $|\boldsymbol{q}|\rightarrow 0$ for the
transverse acoustic phonons, the indirect quadrupole interaction coefficient
of Eq. (4.3) is reduced to the simple formula

$\displaystyle
D^{\mathrm{QQ}}\left(\boldsymbol{k},\boldsymbol{q}=0\right)=-\frac{G_{x^{\prime
2}-y^{\prime
2}}^{2}}{2V\rho_{\mathrm{M}}}\frac{1}{\left(\displaystyle\frac{\hbar}{m_{e}}\right)^{2}\left(k^{2}-\displaystyle\frac{m_{e}^{2}v_{66}^{2}}{\hbar^{2}}\right)}.$
(53)

The indirect quadrupole interaction coefficient of Eq. (53) for electrons with
a small wavenumber of
$-\big{|}\boldsymbol{k}_{\mathrm{b}}^{\mathrm{Q}}\bigr{|}<k<\big{|}\boldsymbol{k}_{\mathrm{b}}^{\mathrm{Q}}\bigr{|}=m_{e}v_{66}/\hbar$
has a positive sign,
$D^{\mathrm{QQ}}\left(\boldsymbol{k},\boldsymbol{q}=0\right)>0$, showing the
ferro-type quadrupole interaction, while the interaction coefficient for a
relatively large wavenumber
$k>\big{|}\boldsymbol{k}_{\mathrm{b}}^{\mathrm{Q}}\bigr{|}=m_{e}v_{66}/\hbar$
possesses a negative sign,
$D^{\mathrm{QQ}}\left(\boldsymbol{k},\boldsymbol{q}=0\right)<0$, indicating
the antiferro-type quadrupole interaction. The ferro-type interaction of
$D^{\mathrm{QQ}}\left(\boldsymbol{k},\boldsymbol{q}=0\right)>0$ brings about
the ferro-quadrupole ordering accompanying the structural transition, as
actually observed in the under-doped compound $x=0.036$.

It is worth estimating the boundary wavenumber
$k_{\mathrm{b}}^{\mathrm{Q}}=\left|\boldsymbol{k}_{\mathrm{b}}^{\mathrm{Q}}\right|$
and the corresponding energy for the present iron pnictide. In the compound
with $x=0.036$, the ultrasonic velocity $v_{66}=772$ m/s at 65 K in the
vicinity of the structural transition (1780 m/s at 100 K far above the
transition) gives the boundary wavenumber
$k_{\mathrm{b}}^{\mathrm{Q}}=6.7\times 10^{6}$ m-1 ($15\times 10^{6}$ m-1),
the boundary wavelength $\lambda_{\mathrm{b}}^{\mathrm{Q}}=0.94\times 10^{-6}$
m ($0.41\times 10^{-6}$ m), and the corresponding frequency
$f_{\mathrm{b}}^{\mathrm{Q}}=\varepsilon\left(k_{\mathrm{b}}^{\mathrm{Q}}\right)/h=0.41$
GHz (2.2 GHz). Here, we take the electron rest mass as $m_{e}$. In the low-
energy regime below 2 GHz $\approx$ 100 mK, therefore, the electrons bearing
the quadrupoles $O_{x^{\prime 2}-y^{\prime 2},\boldsymbol{k},-\boldsymbol{q}}$
and $O_{x^{\prime 2}-y^{\prime 2},\boldsymbol{k}^{\prime},\boldsymbol{q}}$ in
Eq. (51) interact with each other through the strain
$\varepsilon_{xy}\left(\boldsymbol{q}\right)$ in the manner of the ferro-type
quadrupole interaction. Actually, the relaxation rate $\tau^{-1}$ observed
around the structural transition for $x=0.036$ in Fig. 3(b) is always smaller
than the estimated crossover frequency of about 2 GHz. This implies that the
transverse acoustic waves with frequencies below 2 GHz exhibit the full extent
of elastic softening and considerable damping in the vicinity of the
structural transition point, while the acoustic waves with frequencies higher
than 2 GHz exhibit less softening and damping. As shown in Fig. 2(e), the
ultrasonic measurements with frequencies less than 260 MHz actually exhibited
the full amount of softening in $C_{66}$, while the Raman scattering
measurements of the transverse acoustic phonons with frequencies as high as 20
GHz in the end material $x=0$ exhibited only 25% of the full amount of
softening [64].

###  4.4 Quadrupole interaction

The electrons accommodated in the degenerate $y^{\prime}z$ and $zx^{\prime}$
bands carry the quadrupoles. In addition to the indirect quadrupole
interaction $H_{\mathrm{ind}}^{\mathrm{QQ}}$ mediated by the strain
$\varepsilon_{xy}$ presented in Sect. 4.3, we expect direct quadrupole
interactions through the long-range Coulomb potential between electrons at
positions $\boldsymbol{r}_{i}$ and $\boldsymbol{r}_{j}$, expressed as

$\displaystyle\left\langle l_{1}\ \sigma,l_{2}\
\sigma^{\prime}\left|H_{\mathrm{Coulomb}}\right|l_{3}\ \sigma^{\prime},l_{4}\
\sigma\right\rangle$ $\displaystyle\ \ \ \ \ \ \ \ \ \ =\sum_{i\leq j}\iint
d\boldsymbol{r}_{i}d\boldsymbol{r}_{j}\psi_{l_{1}}^{\ast}\bigl{(}\boldsymbol{r}_{i}\bigr{)}v_{\sigma}^{\ast}\bigl{(}\boldsymbol{r}_{i}\bigr{)}\psi_{l_{2}}^{\ast}\bigl{(}\boldsymbol{r}_{j}\bigr{)}v_{\sigma^{\prime}}^{\ast}\bigl{(}\boldsymbol{r}_{j}\bigr{)}$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\times\frac{e^{2}}{\left|\boldsymbol{r}_{i}-\boldsymbol{r}_{j}\right|}\psi_{l_{3}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}v_{\sigma^{\prime}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}\psi_{l_{4}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}v_{\sigma}\bigl{(}\boldsymbol{r}_{i}\bigr{)}.$
(54)

Here, $l_{n}$ ($n=1,2,3,4$) denote suffixes of $y^{\prime}z$ and $zx^{\prime}$
orbital bands and $v_{\sigma}$ is the spin function in Eqs. (39) and (40). The
electric multipole expansion for the Coulomb interaction in Eq. (4.4) for an
isotropic space with $r_{i}<r_{j}$ provides terms consisting of the monopole,
dipole, and quadrupole as [56, 65]

$\displaystyle\frac{e^{2}}{\left|\boldsymbol{r}_{i}-\boldsymbol{r}_{j}\right|}$
$\displaystyle=\frac{e^{2}}{r_{j}}\sum_{k=0}^{\infty}\frac{4\pi}{2k+1}\left(\frac{r_{i}}{r_{j}}\right)^{k}\sum_{m=-k}^{k}Y_{k}^{m}\bigl{(}\theta_{i},\varphi_{i}\bigr{)}Y_{k}^{m\ast}\bigl{(}\theta_{j},\varphi_{j}\bigr{)}$
$\displaystyle=\frac{e^{2}}{r_{j}}+e^{2}\frac{r_{i}}{r_{j}^{2}}\left[\frac{x_{i}}{r_{i}}\frac{x_{j}}{r_{j}}+\frac{y_{i}}{r_{i}}\frac{y_{j}}{r_{j}}+\frac{z_{i}}{r_{i}}\frac{z_{j}}{r_{j}}\right]$
$\displaystyle\quad+e^{2}\frac{r_{i}^{2}}{r_{j}^{3}}\left[O_{3z^{2}-r^{2}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}O_{3z^{2}-r^{2}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}+O_{x^{\prime
2}-y^{\prime 2}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}O_{x^{\prime 2}-y^{\prime
2}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}\right.$
$\displaystyle\qquad\qquad+O_{y^{\prime}z}\bigl{(}\boldsymbol{r}_{i}\bigr{)}O_{y^{\prime}z}\bigl{(}\boldsymbol{r}_{j}\bigr{)}+O_{zx^{\prime}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}O_{zx^{\prime}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}$
$\displaystyle\qquad\qquad\quad\left.+O_{x^{\prime}y^{\prime}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}O_{x^{\prime}y^{\prime}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}\right].$
(55)

According to the multipole expansion of Eq. (4.4), we identify the electric
quadrupole with the order parameter of the structural transition in the
present iron pnictide. In our attempt to clarify the spontaneous symmetry
breaking of the degenerate $y^{\prime}z$ and $zx^{\prime}$ orbital bands in
the tetragonal lattice, we express the direct quadrupole interaction
$H_{\mathrm{C}}^{\mathrm{QQ}}$ in terms of the quadrupole $O_{x^{\prime
2}-y^{\prime 2}}$ with the $B_{2}$ symmetry and $O_{x^{\prime}y^{\prime}}$
with $B_{1}$ symmetry as

$\displaystyle H_{\mathrm{C}}^{\mathrm{QQ}}=$ $\displaystyle-\sum_{i\leq
j}J_{\mathrm{C}}\bigl{(}\boldsymbol{r}_{i}-\boldsymbol{r}_{j}\bigr{)}$
$\displaystyle\times\left[O_{x^{\prime 2}-y^{\prime
2}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}O_{x^{\prime 2}-y^{\prime
2}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}+O_{x^{\prime}y^{\prime}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}O_{x^{\prime}y^{\prime}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}\right].$
(56)

Here, we adopt the Coulomb interaction coefficient
$J_{\mathrm{C}}\bigl{(}\boldsymbol{r}_{i}-\boldsymbol{r}_{j}\bigr{)}$, which
depends on the relative position $\boldsymbol{r}_{i}-\boldsymbol{r}_{j}$. The
quadrupole interaction of $O_{3z^{2}-r^{2}}$ with $A_{1}$ symmetry is excluded
because symmetry breaking is not expected. The case of quadrupoles
$O_{y^{\prime}z}$ and $O_{zx^{\prime}}$ being absent in the $y^{\prime}z$ and
$zx^{\prime}$ bands is beyond the scope of this study. Note that the monopole
interaction among the electrons in Eq. (4.4) is effectively included in the
band model concerned and that the dipole interaction is irrelevant for the
present lattice bearing the inversion symmetry in the $xy$ plane.

Using the expressions $O_{x^{\prime}y^{\prime}}=\sqrt{3}\tau_{x}/7$ in Eq.
(22) and $O_{x^{\prime 2}-y^{\prime 2}}=\sqrt{3}\tau_{z}/7$ in Eq. (28), we
map the quadrupole interaction of Eq. (4.4) on the ideal $xz$ model. In the
present real system, however, the deviation of the anisotropic quadrupole
interaction from the ideal $xz$ model is expected. From the viewpoint of the
symmetry of the tetragonal lattice, the mutual quadrupole interactions of
$O_{x^{\prime 2}-y^{\prime 2}}$ with $B_{2}$ symmetry and
$O_{x^{\prime}y^{\prime}}$ with $B_{1}$ symmetry may be different. Actually,
it has been reported that the mixing between the 4$p$ states of As ions and
the 3$d$ states of Fe ions brings about an anisotropic quadrupole interaction
[66]. The quadrupoles $O_{x^{\prime 2}-y^{\prime 2}}$ at positions
$\boldsymbol{r}_{i}$ and $\boldsymbol{r}_{j}$ interact with each other through
the strain $\varepsilon_{xy}$ of the transverse acoustic phonons for $C_{66}$,
which exhibits the large amount of softening. The quadrupoles
$O_{x^{\prime}y^{\prime}}$ may also interact with each other mediated by the
strain $\varepsilon_{x^{2}-y^{2}}$ of the transverse acoustic phonons for
$(C_{11}-C_{12})/2$, which monotonically increases. The former quadrupole
interaction of $O_{x^{\prime 2}-y^{\prime 2}}$ plays a significant role in the
structural transition, while the latter quadrupole interaction of
$O_{x^{\prime}y^{\prime}}$ has a minor effect on the phase transition.
Consequently, the indirect quadrupole interaction mediated by the transverse
acoustic phonons also indicates the anisotropic nature of the quadrupole
interaction.

Combining the direct quadrupole interaction $H_{\mathrm{C}}^{\mathrm{QQ}}$ of
Eq. (4.4) followed by the Coulomb potential and the indirect quadrupole
interaction $H_{\mathrm{ind}}^{\mathrm{QQ}}$ of Eq. (51) mediated by the
strain $\varepsilon_{xy}$ of the transverse acoustic phonons, we obtain the
anisotropic quadrupole interaction $H_{\mathrm{QQ}}$ specified by parameter
$\gamma$ as

$\displaystyle H_{\mathrm{QQ}}\bigl{(}\gamma\bigr{)}$
$\displaystyle\quad=-\sum_{i\leq
j}J_{\mathrm{Q}}\bigl{(}\boldsymbol{r}_{i}-\boldsymbol{r}_{j}\bigr{)}$
$\displaystyle\qquad\times\left[O_{x^{\prime 2}-y^{\prime
2}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}O_{x^{\prime 2}-y^{\prime
2}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}+\gamma
O_{x^{\prime}y^{\prime}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}O_{x^{\prime}y^{\prime}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}\right].$
(57)

The quadrupole interaction coefficient
$J_{\mathrm{Q}}\bigl{(}\boldsymbol{r}_{i}-\boldsymbol{r}_{j}\bigr{)}$ in Eq.
(4.4) is written as the sum of the Coulomb interaction coefficient
$J_{\mathrm{C}}\bigl{(}\boldsymbol{r}_{i}-\boldsymbol{r}_{j}\bigr{)}$ in the
direct interaction of Eq. (4.4) and the indirect quadrupole interaction
coefficient
$D^{\mathrm{QQ}}\bigl{(}\boldsymbol{r}_{i}-\boldsymbol{r}_{j}\bigr{)}$
mediated by the strain $\varepsilon_{xy}$ of the transverse acoustic phonons
in Eq. (4.3).

$\displaystyle
J_{\mathrm{Q}}\bigl{(}\boldsymbol{r}_{i}-\boldsymbol{r}_{j}\bigr{)}=J_{\mathrm{C}}\bigl{(}\boldsymbol{r}_{i}-\boldsymbol{r}_{j}\bigr{)}+D^{\mathrm{QQ}}\bigl{(}\boldsymbol{r}_{i}-\boldsymbol{r}_{j}\bigr{)}.$
(58)

The quadrupole interaction coefficient
$J_{\mathrm{Q}}\bigl{(}\boldsymbol{r}_{i}-\boldsymbol{r}_{j}\bigr{)}$ and the
anisotropic parameter $\gamma$ in Eq. (4.4) play a substantial role in the
appearance of the ferro-quadrupole ordering associated with the structural
transition and the hexadecapole ordering associated with the superconducting
transition.

###  4.5 Quadrupole ordering

In order to explain the elastic softening of $C_{66}$ and critical slowing
down of the relaxation time $\tau$ around the structural transition in the
iron pnictide, we consider the quadrupole interaction Hamiltonian
$H_{\mathrm{QQ}}(\gamma)$ of Eq. (4.4) for the case of $\gamma=0$ mapped on
the Ising model. A transverse ultrasonic wave with a small-wavenumber limit of
$|\boldsymbol{q}|\rightarrow 0$ induces the strain
$\varepsilon_{xy}=\varepsilon_{xy}\left(\boldsymbol{q}=0\right)$, which is
treated as a classical quantity. According to the quadrupole-strain
interaction of Eq. (44), the strain $\varepsilon_{xy}$ generated by the
ultrasonic wave lifts the degenerate $y^{\prime}z$ and $zx^{\prime}$ bands, as
schematically shown in Fig. 7(a). This ultrasonic perturbation gives a
quadrupole susceptibility proportional to the reciprocal temperature as

$\displaystyle\chi_{\mathrm{Q}}=\frac{N_{\mathrm{Q}}G_{x^{\prime 2}-y^{\prime
2}}^{2}/C_{66}^{0}}{T}=\frac{\mathit{\Delta}_{\mathrm{Q}}}{T}.$ (59)

Here, $N_{\mathrm{Q}}$ is the number of electrons carrying the ferro-type
quadrupole interaction in the small-wavenumber regime of
$-k_{\mathrm{b}}^{\mathrm{Q}}<k<k_{\mathrm{b}}^{\mathrm{Q}}$ for the boundary
wavenumber $k_{\mathrm{b}}^{\mathrm{Q}}$ as discussed for Eq. (53). The
quadrupole-strain interaction energy of
$\mathit{\Delta}_{\mathrm{Q}}=N_{\mathrm{Q}}G_{x^{\prime 2}-y^{\prime
2}}^{2}/C_{66}^{0}$ in Eq. (59) is determined in terms of
$D^{\mathrm{QQ}}\left(\boldsymbol{k},\boldsymbol{q}\right)$ of Eq. (53) for
the small-wavenumber regime of $|\boldsymbol{q}|\rightarrow 0$ as

$\displaystyle\frac{\mathit{\Delta}_{\mathrm{Q}}}{2VN_{\mathrm{Q}}}\frac{C_{66}^{0}}{\rho_{\mathrm{M}}v_{66}^{2}}$
$\displaystyle=\frac{G_{x^{\prime 2}-y^{\prime
2}}^{2}}{2V\rho_{\mathrm{M}}v_{66}^{2}}$
$\displaystyle=\frac{1}{N_{\mathrm{Q}}}\sum_{|\boldsymbol{k}|<k_{\mathrm{b}}^{\mathrm{Q}}}D^{\mathrm{QQ}}\left(\boldsymbol{k},\boldsymbol{q}=0\right)$
$\displaystyle=\widetilde{D}^{\mathrm{QQ}}.$ (60)

Here, $\sum_{|\boldsymbol{k}|<k_{\mathrm{b}}^{\mathrm{Q}}}$ means the sum over
the electron states with wavenumber
$|\boldsymbol{k}|<k_{\mathrm{b}}^{\mathrm{Q}}$, that participate in the ferro-
type quadrupole interaction. $\widetilde{D}^{\mathrm{QQ}}$ in Eq. (4.5) stands
for the effective indirect quadrupole interaction of Eq. (4.3). As shown in
Fig. 4, the ultrasonic experiments give a quadrupole-strain interaction energy
of $\mathit{\Delta}_{\mathrm{Q}}\sim 20$ K that is almost independent of the
Co concentration $x$. The indirect quadrupole interaction energy
$D^{\mathrm{QQ}}$ in Eq. (4.5) is enhanced by the softening of
$\rho_{\mathrm{M}}v_{66}^{2}=C_{66}$ due to the quadrupole-strain interaction
of Eq. (44).

Figure 7:  (Color online) (a) Energy splitting of the degenerate $y^{\prime}z$
and $zx^{\prime}$ bands by the quadrupole-strain interaction $H_{\mathrm{QS}}$
of Eq. (44). (b) Energy splitting of the two-electron states
$\psi_{\pm}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ of Eq. (4.6)
by the hexadecapole-rotation interaction
$H_{\mathrm{rot}}\bigl{(}\omega_{xy}\bigr{)}$ of Eq. (92). The transverse
ultrasonic wave with amplitude $\xi_{y}^{0}$ and frequency $\omega$ induces
infinitesimal strain $\varepsilon_{xy}=\sin\theta$ and rotation
$\omega_{xy}=\sin\theta$ with
$\sin\theta=q_{x}\xi_{y}^{0}=\left(\omega/v_{x}\right)\xi_{y}^{0}=\left(2\pi/\lambda_{x}\right)\xi_{y}^{0}$
for a given phase velocity $v_{x}=v_{66}$.

The softening of $C_{66}$ as a precursor of the structural transition is
expressed as

$\displaystyle C_{66}=C_{66}^{0}\left(1-\widetilde{\chi}_{\mathrm{Q}}\right)$
$\displaystyle=C_{66}^{0}\left(\frac{T-T_{\mathrm{s}}^{0}}{T-\mathit{\Theta}_{\mathrm{Q}}}\right)$
$\displaystyle=C_{66}^{0}\left(1-\frac{\mathit{\Delta}_{\mathrm{Q}}}{T-\mathit{\Theta}_{\mathrm{Q}}}\right).$
(61)

Here, we adopt the renormalized susceptibility $\widetilde{\chi}_{\mathrm{Q}}$
expressed by the quadrupole interaction energy $\mathit{\Theta}_{\mathrm{Q}}$
and the quadrupole-strain interaction energy $\mathit{\Delta}_{\mathrm{Q}}$ as

$\displaystyle\widetilde{\chi}_{\mathrm{Q}}=\frac{\mathit{\Delta}_{\mathrm{Q}}}{T-\mathit{\Theta}_{\mathrm{Q}}}.$
(62)

The elastic instability point due to the full softening $C_{66}\rightarrow 0$
gives the structural transition temperature as
$T_{\mathrm{s}}^{0}=\mathit{\Theta}_{\mathrm{Q}}+\mathit{\Delta}_{\mathrm{Q}}$.
This has already been used in Eq. (2) for the analysis of the experimentally
observed softening of $C_{66}$.

By using the mean field approximation, we obtain the following self-consistent
equation giving the temperature dependence of the quadrupole order parameter
$O_{x^{\prime 2}-y^{\prime 2}}$:

$\displaystyle\left\langle O_{x^{\prime 2}-y^{\prime
2}}\right\rangle=\tanh\left[\widetilde{J}_{\mathrm{Q}}\left\langle
O_{x^{\prime 2}-y^{\prime 2}}\right\rangle/k_{\mathrm{B}}T\right].$ (63)

Here, $\langle A\rangle$ stands for the thermal average. We use the effective
coupling constant $\widetilde{J}_{\mathrm{Q}}$, which is the average of the
Fourier transformed quadrupole interaction of
$J_{\mathrm{Q}}\bigl{(}\boldsymbol{r}_{i}-\boldsymbol{r}_{j}\bigr{)}=J_{\mathrm{C}}\bigl{(}\boldsymbol{r}_{i}-\boldsymbol{r}_{j}\bigr{)}+D^{\mathrm{QQ}}\bigl{(}\boldsymbol{r}_{i}-\boldsymbol{r}_{j}\bigr{)}$
in Eq. (58) over the electron states with wavenumber
$|\boldsymbol{k}|<k_{\mathrm{b}}^{\mathrm{Q}}$. $\widetilde{J}_{\mathrm{Q}}$
is expressed as

$\displaystyle\widetilde{J}_{\mathrm{Q}}$
$\displaystyle=\widetilde{J}_{\mathrm{C}}+\widetilde{D}^{\mathrm{QQ}}$
$\displaystyle=\frac{1}{N_{\mathrm{Q}}}\sum_{|\boldsymbol{k}|<k_{\mathrm{b}}^{\mathrm{Q}}}J_{\mathrm{C}}\bigl{(}\boldsymbol{k},\boldsymbol{q}=0\bigr{)}+\frac{1}{N_{\mathrm{Q}}}\sum_{|\boldsymbol{k}|<k_{\mathrm{b}}^{\mathrm{Q}}}D^{\mathrm{QQ}}\bigl{(}\boldsymbol{k},\boldsymbol{q}=0\bigr{)}.$
(64)

Here, $\widetilde{J}_{\mathrm{C}}$ is the effective energy of the direct
Coulomb interaction of Eq. (4.4). $\widetilde{D}^{\mathrm{QQ}}$ is the
effective indirect quadrupole interaction coefficient in Eq. (4.5). We suppose
that the ferro-quadrupole ordering is caused by the quadrupole interaction
with the positive coupling constant $\widetilde{J}_{\mathrm{Q}}>0$. We expect
there to be a non vanishing solution $\langle O_{x^{\prime 2}-y^{\prime
2}}\rangle\neq 0$ at low temperatures below the ferro-quadrupole ordering
point of $T_{\mathrm{c}}=\widetilde{J}_{\mathrm{Q}}/k_{\mathrm{B}}$.

The elastic instability due to the full softening $C_{66}\rightarrow 0$ brings
about spontaneous lattice distortion accompanying the ferro-quadrupole
ordering as

$\displaystyle\left\langle\varepsilon_{xy}\right\rangle=\frac{N_{\mathrm{Q}}G_{x^{\prime
2}-y^{\prime 2}}\left\langle O_{x^{\prime 2}-y^{\prime
2}}\right\rangle}{C_{66}}=\sqrt{\frac{N_{\mathrm{Q}}\mathit{\Delta}_{\mathrm{Q}}}{C_{66}}}\left\langle
O_{x^{\prime 2}-y^{\prime 2}}\right\rangle.$ (65)

Here, we use the relation between the coupling constant $G_{x^{\prime
2}-y^{\prime 2}}$ and the quadrupole-strain interaction energy
$\mathit{\Delta}_{\mathrm{Q}}$ of $G_{x^{\prime 2}-y^{\prime
2}}=\sqrt{\mathit{\Delta}_{\mathrm{Q}}C_{66}^{0}/N_{\mathrm{Q}}}$. The ferro-
quadrupole ordering of $\bigl{\langle}O_{x^{\prime 2}-y^{\prime
2}}\bigr{\rangle}\neq 0$ associated with the spontaneous strain of
$\bigl{\langle}\varepsilon_{xy}\bigr{\rangle}\neq 0$ in Eq. (65) is caused by
the structural transition from the tetragonal phase with the high-symmetry
space group $D_{4h}^{17}$ to the orthorhombic phase with the low-symmetry
space group $D_{2h}^{23}$. The generators of the symmorphic space group
$D_{4h}^{17}$ consist of the rotation, reflection, and inversion operations
but do not involve the screw and glide operations. The quenching of the
$B_{2g}$ symmetry in the mother phase of $D_{4h}^{17}$ associated with the
ferro-quadrupole ordering of $\langle O_{x^{\prime 2}-y^{\prime 2}}\rangle\neq
0$ and the spontaneous strain $\langle\varepsilon_{xy}\rangle\neq 0$ across
the structural transition loses the symmetry operations, which consist of
rotation through $\pm\pi/2$ about the vertical axis of $C_{4}$, rotations
through $\pi$ about the horizontal $x$\- and $y$-axes of $2C_{2}^{\prime}$,
rotation through $\pm\pi/2$ about the vertical axis followed by inversion of
$2IC_{4}$, and mirror reflection in the vertical plane of
$2\sigma_{\mathrm{v}}$ [51, 67]. Note that the ferro-type quadrupole ordering
denoted in the present paper has been frequently referred to as cooperative
Jahn-Teller effects [68], orbital orderings, or electronic nematic orders [69,
70, 71, 72, 73, 74, 75, 76, 77, 78, 79].

The neutron scattering experiments on the end material BaFe2As2 have shown a
spontaneous strain $\bigl{\langle}\varepsilon_{xy}\bigr{\rangle}$ of
$5.08\times 10^{-3}$ at 5 K in the distorted orthorhombic phase [10], where
the order parameter is fully polarized as $\bigl{\langle}O_{x^{\prime
2}-y^{\prime 2}}\bigr{\rangle}=1$. Adopting the quadrupole-strain interaction
energy $\mathit{\Delta}_{\mathrm{Q}}=21$ K and the elastic constant
$C_{66}^{0}=3.19\times 10^{10}$ J/m3 of BaFe2As2 [43], we estimate the number
of 3$d$ electrons $N_{\mathrm{Q}}$ in Eq. (65) to be $2.84\times 10^{21}$
cm-3, which is approximately one-fourteenth of the total number of electrons
of $2N_{\mathrm{Fe}}=3.92\times 10^{22}$ cm-3. Furthermore, we deduce the
coupling energy $G_{x^{\prime 2}-y^{\prime 2}}$ to be $\sim 4\times 10^{3}$ K
per electron, which is four times larger than the tentatively estimated
$G_{x^{\prime 2}-y^{\prime 2}}$ of $\sim 1\times 10^{3}$ K for the localized
electron picture in Sect. 3.1. The large quadrupole-strain coupling energy of
$G_{x^{\prime 2}-y^{\prime 2}}\sim 4\times 10^{3}$ K is caused by the enhanced
quadrupole due to the extended orbital radius compatible with the itinerant
feature of the 3$d$ electron bands. Taking the number of 3$d$ electrons
$N_{\mathrm{Q}}=2.84\times 10^{21}$ cm-3 into account, we may roughly estimate
the indirect quadrupole interaction energy in BaFe2As2 as $D^{\mathrm{QQ}}\sim
220$ K at $T=T_{\mathrm{s}}=135$ K for $\mathit{\Delta}_{\mathrm{Q}}=21$ K,
$C_{66}^{0}=3.19\times 10^{10}$ J/m3, and
$C_{66}=\rho_{\mathrm{M}}v_{66}^{2}=0.281\times 10^{10}$ J/m3. This value is
comparable with the structural transition temperature $T_{\mathrm{s}}$ of the
under-doped compounds. It is worth referring to the extremely enhanced
quadrupole-strain interaction energy of $2.8\times 10^{5}$ K for the vacancy
orbital of a silicon wafer with a large orbital radius, which was verified by
means of bulk ultrasonic waves and surface acoustic waves [80, 81].

The divergence of the relaxation time $\tau$ observed via the ultrasonic
attenuation coefficient $\alpha_{66}$ for $x=0.036$ in Fig. 3(b) is expressed
in terms of the critical slowing down due to freezing of the ferro-quadrupole
order parameter $O_{x^{\prime 2}-y^{\prime 2}}$ at the structural transition
temperature $T_{\mathrm{s}}^{0}\cong T_{\mathrm{c}}^{0}=65$ K. The relaxation
time $\tau_{\mathrm{Q}}$ is written as

$\displaystyle\tau_{\mathrm{Q}}=\tau_{0}\left|\frac{T-T_{\mathrm{c}}^{0}}{T_{\mathrm{c}}^{0}}\right|^{-1}\propto\frac{1}{1-\widetilde{\chi}_{\mathrm{Q}}}.$
(66)

This expression for $\tau_{\mathrm{Q}}$ based on the degenerate $y^{\prime}z$
and $zx^{\prime}$ bands well reproduces the experimental results analyzed in
terms of Eq. (4) for the critical index $z\nu=$ 1 for both tetragonal and
orthorhombic phases. The critical slowing down of the quadrupole $O_{x^{\prime
2}-y^{\prime 2}}$ in the present iron pnictide agrees with the critical
dynamics of the kinetic Ising model [39]. According to the dissipation
fluctuation theorem, the critical slowing down of the relaxation time
$\tau_{\mathrm{Q}}$ is expressed by the divergence of the susceptibility for
the order parameter of the quadrupole $O_{x^{\prime 2}-y^{\prime 2}}$ due to
the infinite increase in the correlation length $\zeta$ in the vicinity of the
structural transition. By analogy with the ferro-magnetic spin system [40,
42], the development of the correlation function of $\langle O_{x^{\prime
2}-y^{\prime 2}}(\boldsymbol{k},t)O_{x^{\prime 2}-y^{\prime
2}}(-\boldsymbol{k},0)\rangle$ obeying the time decay factor of
$\exp[-t/\tau_{\mathrm{Q}}(\boldsymbol{k})]$ should cause the relaxation time
$\tau_{\mathrm{Q}}(\boldsymbol{k})$ to diverge at the ferro-type ordering
point. Knowledge of the diffusion coefficient is necessary for the numerical
estimation of $\tau_{0}$ in Eq. (66).

###  4.6 Hexadecapole-rotation interaction of two-electron states

In the over-doped compound $x=0.071$ exhibiting superconductivity, the
antiferro-type quadrupole interaction with the negative quadrupole interaction
energy $\mathit{\Theta}_{\mathrm{Q}}=-47$ K is verified by analysis of the
softening of $C_{66}$ in the normal phase. The softening of $C_{66}$ tends to
zero at the fictitious structural instability point at the negative
temperature $T_{\mathrm{s}}^{0}=-26.5$ K. This means that the critical slowing
down of the relaxation time $\tau$ around the superconducting transition at
$T_{\mathrm{c}}^{0}=T_{\mathrm{SC}}=23$ K in Fig. 1(b) is caused by an
appropriate order parameter, which is strictly distinguished from the
quadrupole. In the present section, we introduce two-electron states bound by
the quadrupole interaction and consider the coupling of the hexadecapole
carried by two-electron states to the rotation $\omega_{xy}$ of transverse
ultrasonic waves.

Taking the Pauli exclusion principle into account, we express the energy of
the anisotropic quadrupole interaction of Eq. (4.4) by using the two-electron
state
$\psi_{\Gamma_{\gamma}}^{S,S_{z}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$
of the Slater determinant with the irreducible representation
$\Gamma_{\gamma}$ of the orbital part, and the spin state denoted by the total
spin $S$ and $z$ component $S_{z}$ as

$\displaystyle E_{\mathrm{QQ}}^{\Gamma_{\gamma},S,S_{z}}$
$\displaystyle=\bigl{\langle}\Gamma_{\gamma}\ S\
S_{z}|H_{\mathrm{QQ}}\bigl{(}\gamma\bigr{)}|\Gamma_{\gamma}\ S\
S_{z}\bigr{\rangle}$ $\displaystyle=\iint
d\boldsymbol{r}_{i}d\boldsymbol{r}_{j}\psi_{\Gamma_{\gamma}}^{S,S_{z}\ast}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}H_{\mathrm{QQ}}\bigl{(}\gamma\bigr{)}\psi_{\Gamma_{\gamma}}^{S,S_{z}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}.$
(67)

From the viewpoint of symmetry for the orbital state, the direct product of
the $E$-doublet in the point group symmetry $D_{2d}$ reduces as $E\otimes
E=A_{1}\oplus A_{2}\oplus B_{1}\oplus B_{2}$. Thus, the symmetric $A_{1}$,
$B_{1}$, and $B_{2}$ orbital states exist upon exchanging $\boldsymbol{r}_{i}$
for $\boldsymbol{r}_{j}$ as well as the antisymmetric $A_{2}$ state.
Furthermore, there are two kinds of spin states consisting of a spin singlet
of total spin $S=0$ that is symmetric upon exchanging $\boldsymbol{r}_{i}$ for
$\boldsymbol{r}_{j}$ and the antisymmetric spin triplet with $S=1$.
Consequently, we deduce the two-electron state of
$\psi_{\Gamma_{\gamma}}^{S,S_{z}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$
in Eq. (4.6), which is denoted by the orbital state of the symmetry
$\Gamma_{\gamma}=A_{1}$, $B_{2}$, and $B_{1}$ with the spin singlet of $S=0$
and the orbital state of the $A_{2}$ symmetry with the spin triplet of $S=1$
for $S_{z}=1$, $0$, and $-1$. As a result, we obtain the Slater determinant in
terms of the one-electron orbital states of
$\psi_{y^{\prime}z}\bigl{(}\boldsymbol{r}_{i}\bigr{)}$ and
$\psi_{zx^{\prime}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}$ and the spin states of
$\alpha(\boldsymbol{r}_{i})$ and $\beta(\boldsymbol{r}_{i})$ as

$\displaystyle\psi_{\Gamma_{\gamma}=A_{1}}^{S=0,S_{z}=0}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$
$\displaystyle=\frac{1}{\sqrt{2}}\left[\psi_{y^{\prime}z}\bigl{(}\boldsymbol{r}_{i}\bigr{)}\psi_{y^{\prime}z}\bigl{(}\boldsymbol{r}_{j}\bigr{)}+\psi_{zx^{\prime}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}\psi_{zx^{\prime}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}\right]$
$\displaystyle\quad\times\frac{1}{\sqrt{2}}\left[\alpha\bigl{(}\boldsymbol{r}_{i}\bigr{)}\beta\bigl{(}\boldsymbol{r}_{j}\bigr{)}-\beta\bigl{(}\boldsymbol{r}_{i}\bigr{)}\alpha\bigl{(}\boldsymbol{r}_{j}\bigr{)}\right],$
(68)
$\displaystyle\psi_{B_{2}}^{0,0}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$
$\displaystyle=\frac{1}{\sqrt{2}}\left[\psi_{y^{\prime}z}\bigl{(}\boldsymbol{r}_{i}\bigr{)}\psi_{y^{\prime}z}\bigl{(}\boldsymbol{r}_{j}\bigr{)}-\psi_{zx^{\prime}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}\psi_{zx^{\prime}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}\right]$
$\displaystyle\quad\times\frac{1}{\sqrt{2}}\left[\alpha\bigl{(}\boldsymbol{r}_{i}\bigr{)}\beta\bigl{(}\boldsymbol{r}_{j}\bigr{)}-\beta\bigl{(}\boldsymbol{r}_{i}\bigr{)}\alpha\bigl{(}\boldsymbol{r}_{j}\bigr{)}\right],$
(69)
$\displaystyle\psi_{B_{1}}^{0,0}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$
$\displaystyle=\frac{1}{\sqrt{2}}\left[\psi_{y^{\prime}z}\bigl{(}\boldsymbol{r}_{i}\bigr{)}\psi_{zx^{\prime}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}+\psi_{zx^{\prime}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}\psi_{y^{\prime}z}\bigl{(}\boldsymbol{r}_{j}\bigr{)}\right]$
$\displaystyle\quad\times\frac{1}{\sqrt{2}}\left[\alpha\bigl{(}\boldsymbol{r}_{i}\bigr{)}\beta\bigl{(}\boldsymbol{r}_{j}\bigr{)}-\beta\bigl{(}\boldsymbol{r}_{i}\bigr{)}\alpha\bigl{(}\boldsymbol{r}_{j}\bigr{)}\right],$
(70)
$\displaystyle\psi_{A_{2}}^{1,1}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$
$\displaystyle=\frac{1}{\sqrt{2}}\left[\psi_{y^{\prime}z}\bigl{(}\boldsymbol{r}_{i}\bigr{)}\psi_{zx^{\prime}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}-\psi_{zx^{\prime}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}\psi_{y^{\prime}z}\bigl{(}\boldsymbol{r}_{j}\bigr{)}\right]$
$\displaystyle\quad\times\alpha\bigl{(}\boldsymbol{r}_{i}\bigr{)}\alpha\bigl{(}\boldsymbol{r}_{j}\bigr{)},$
(71)
$\displaystyle\psi_{A_{2}}^{1,0}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$
$\displaystyle=\frac{1}{\sqrt{2}}\left[\psi_{y^{\prime}z}\bigl{(}\boldsymbol{r}_{i}\bigr{)}\psi_{zx^{\prime}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}-\psi_{zx^{\prime}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}\psi_{y^{\prime}z}\bigl{(}\boldsymbol{r}_{j}\bigr{)}\right]$
$\displaystyle\quad\times\frac{1}{\sqrt{2}}\left[\alpha\bigl{(}\boldsymbol{r}_{i}\bigr{)}\beta\bigl{(}\boldsymbol{r}_{j}\bigr{)}+\beta\bigl{(}\boldsymbol{r}_{i}\bigr{)}\alpha\bigl{(}\boldsymbol{r}_{j}\bigr{)}\right],$
(72)
$\displaystyle\psi_{A_{2}}^{1,-1}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$
$\displaystyle=\frac{1}{\sqrt{2}}\left[\psi_{y^{\prime}z}\bigl{(}\boldsymbol{r}_{i}\bigr{)}\psi_{zx^{\prime}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}-\psi_{zx^{\prime}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}\psi_{y^{\prime}z}\bigl{(}\boldsymbol{r}_{j}\bigr{)}\right]$
$\displaystyle\quad\times\beta\bigl{(}\boldsymbol{r}_{i}\bigr{)}\beta\bigl{(}\boldsymbol{r}_{j}\bigr{)}.$
(73)

In order to properly describe the energy of the anisotropic quadrupole
interaction of Eq. (4.6), we calculate the quadrupole interaction energy of
the Hamiltonian $H_{\mathrm{QQ}}\bigl{(}\gamma\bigr{)}$ of Eq. (4.4) in terms
of the matrix representation for the two-electron states
$\psi_{\Gamma_{\gamma}}^{S,S_{z}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$
of Eqs. (4.6)-(4.6).

$\displaystyle\boldsymbol{H}_{\mathrm{QQ}}\bigl{(}\gamma\bigr{)}=-\sum_{i\leq
j}J_{\mathrm{Q}}^{ij}$ $\displaystyle\times$
$\displaystyle\bordermatrix{&\psi_{A_{1}}^{0,0}&\psi_{B_{2}}^{0,0}&\psi_{B_{1}}^{0,0}&\psi_{A_{2}}^{1,1}&\psi_{A_{2}}^{1,0}&\psi_{A_{2}}^{1,-1}\cr&1+\gamma&0&0&0&0&0\cr&0&1-\gamma&0&0&0&0\cr&0&0&-1+\gamma&0&0&0\cr&0&0&0&-1-\gamma&0&0\cr&0&0&0&0&-1-\gamma&0\cr&0&0&0&0&0&-1-\gamma}.$
(74)

Here, the quadrupole interaction energy of $J_{\mathrm{Q}}^{ij}$ between
electrons at positions $\boldsymbol{r}_{i}$ and $\boldsymbol{r}_{j}$ is
calculated in terms of the radius function $f_{d}\bigl{(}r\bigr{)}$ for the
3$d$ electron in Eqs. (4.2) and (4.2) as

$\displaystyle J_{\mathrm{Q}}^{ij}$
$\displaystyle=\left(\frac{\sqrt{3}}{7}\right)^{2}\iint
dr_{i}dr_{j}J_{\mathrm{Q}}\bigl{(}\boldsymbol{r}_{i}-\boldsymbol{r}_{j}\bigr{)}r_{i}^{2}f_{d}\bigl{(}r_{i}\bigr{)}^{2}r_{j}^{2}f_{d}\bigl{(}r_{j}\bigr{)}^{2}$
$\displaystyle=\frac{3}{49}\iint
drdr^{\prime}J_{\mathrm{Q}}\bigl{(}\boldsymbol{r}-\boldsymbol{R}_{i}-\boldsymbol{r}^{\prime}+\boldsymbol{R}_{j}\bigr{)}|\boldsymbol{r}-\boldsymbol{R}_{i}|^{2}f_{d}\bigl{(}|\boldsymbol{r}-\boldsymbol{R}_{i}|\bigr{)}^{2}$
$\displaystyle\qquad\qquad\qquad\times|\boldsymbol{r}^{\prime}-\boldsymbol{R}_{j}|^{2}f_{d}\bigl{(}|\boldsymbol{r}^{\prime}-\boldsymbol{R}_{j}|\bigr{)}^{2}.$
(75)

In the integral in the lower line in Eq. (4.6), we take coordinates
$\boldsymbol{r}$ and $\boldsymbol{r}^{\prime}$ instead of
$\boldsymbol{r}_{i}=\boldsymbol{r}-\boldsymbol{R}_{i}$ and
$\boldsymbol{r}_{j}=\boldsymbol{r}^{\prime}-\boldsymbol{R}_{j}$ for the
positions $\boldsymbol{R}_{i}$ and $\boldsymbol{R}_{j}$ of Fe2+ ion,
respectively. The intra-atomic quadrupole interaction of $J_{\mathrm{Q}}^{ij}$
for the electrons accommodated in the orbital at the same Fe2+ ion, $i$. $e$.,
$\boldsymbol{R}_{i}=\boldsymbol{R}_{j}$, is much stronger than the inter-
atomic quadrupole interaction for electrons located at different Fe2+ ions,
$i$. $e$., $\boldsymbol{R}_{i}\neq\boldsymbol{R}_{j}$. It is expected that the
intra-atomic quadrupole interaction $J_{\mathrm{Q}}^{i=j}$ for
$\boldsymbol{R}_{i}=\boldsymbol{R}_{j}$, which is mostly independent of the
lattice structure, will show isotropic azimuth angle dependence. This almost
isotropic feature of the quadrupole interaction of $J_{\mathrm{Q}}^{ij}$
dominated by the intra-atomic coupling favors the s-like superconducting
energy gap of the present iron pnictide. This will be discussed in Sect. 4.9.

Ultrasonic measurements under pulsed magnetic fields of up to 62 T reveal
increases of 7.4% for the softened $C_{66}$ at 84 K in the tetragonal phase of
$x=0.036$ and 3.3% for the softened $C_{66}$ at 30 K in the normal phase of
$x=0.071$ [82]. Taking into account the robustness of the softening of
$C_{66}$ against magnetic fields in the normal phase above the structural and
superconducting transition temperatures, we adopt the two-electron states of
the Slater determinant with the spin-singlet states of Eqs. (4.6)-(4.6), but
disregard those with the spin-triplet states of Eqs. (4.6)-(4.6). We
abbreviate the two-electron state with the spin singlet of
$\psi_{\Gamma_{\gamma}}^{0,0}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$
in Eqs. (4.6)-(4.6) to
$\psi_{\Gamma_{\gamma}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$
in the following discussion.

In the case of the highly anisotropic quadrupole interaction with $\gamma=0$,
the quadrupole interaction energies of the
$\psi_{A_{1}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ and
$\psi_{B_{2}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ states in
Eq. (4.6) causes the degeneration of each other and the energy of the
$\psi_{B_{1}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ state
deviates from them. This highly anisotropic case of $\gamma\approx 0$ in Eq.
(4.4) leads to the ferro-quadrupole ordering of $O_{x^{\prime 2}-y^{\prime
2}}$ accompanying the structural transition. With increasing the Co
concentration $x$ to the QCP of $x_{\mathrm{QCP}}=0.061$, the quadrupole
interaction of Eq. (4.4) develops an almost isotropic feature with
$\gamma\stackrel{{\scriptstyle<}}{{{}_{\sim}}}1$, where the energies of the
$\psi_{B_{2}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ and
$\psi_{B_{1}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ states get
closer to each other and the energy of the
$\psi_{A_{1}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ state is
different. Since the quadrupoles $O_{x^{\prime 2}-y^{\prime 2}}$ of Eq. (28)
and $O_{x^{\prime}y^{\prime}}$ of Eq. (22) and the angular momentum $l_{z}$ of
Eq. (25) obey the commutation relation of Eq. (29) among the Pauli matrices,
we expect quantum fluctuation between the quadrupoles $O_{x^{\prime
2}-y^{\prime 2}}$ and $O_{x^{\prime}y^{\prime}}$ in the vicinity of the QCP,
which is particularly important for explaining the superconductivity
accompanying the hexadecapole ordering. The special case of $\gamma=1$ mapped
on the ideal $xz$ model leads to the fully degenerate
$\psi_{B_{2}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ and
$\psi_{B_{1}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ states,
which causes the disappearance of the hexadecapole. This is inconsistent with
the experimental observation of the critical slowing down around the
superconducting transition temperature.

When the rotation $\omega_{xy}$ is induced by either the thermally excited
transverse acoustic phonons or the experimentally generated transverse
ultrasonic waves, the phases of the two-electron state of
$\psi_{{\Gamma}_{\gamma}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$
at positions $\boldsymbol{r}_{i}$ and $\boldsymbol{r}_{j}$ in Eq. (4.6) bound
by the quadrupole interaction of Eq. (4.4) are simultaneously changed by the
rotation operator as follows:

$\displaystyle\psi_{\Gamma_{\gamma}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$
$\displaystyle\rightarrow\exp\left[-il_{z}\bigl{(}\boldsymbol{r}_{i}\bigr{)}\omega_{xy}\right]\exp\left[-il_{z}\bigl{(}\boldsymbol{r}_{j}\bigr{)}\omega_{xy}\right]$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\times\psi_{\Gamma_{\gamma}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}.$
(76)

Thus, an infinitesimal amount of rotation $\omega_{xy}$ perturbs the
quadrupole interaction Hamiltonian of Eq. (4.4) as

$\displaystyle\left\langle\Gamma_{\gamma}\left|H_{\mathrm{QQ}}\bigl{(}\omega_{xy}\bigr{)}\right|\Gamma_{\gamma^{\prime}}\right\rangle$
$\displaystyle=\iint
d\boldsymbol{r}_{i}d\boldsymbol{r}_{j}\psi_{\Gamma_{\gamma}}^{\ast}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}\exp\left\\{i\left[l_{z}\bigl{(}\boldsymbol{r}_{i}\bigr{)}+l_{z}\bigl{(}\boldsymbol{r}_{j}\bigr{)}\right]\omega_{xy}\right\\}$
$\displaystyle\qquad\qquad\times
H_{\mathrm{QQ}}\bigl{(}\gamma\bigr{)}\exp\left\\{-i\left[l_{z}\bigl{(}\boldsymbol{r}_{i}\bigr{)}+l_{z}\bigl{(}\boldsymbol{r}_{j}\bigr{)}\right]\omega_{xy}\right\\}\psi_{\Gamma_{\gamma^{\prime}}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$
$\displaystyle\approx\iint
d\boldsymbol{r}_{i}d\boldsymbol{r}_{j}\psi_{\Gamma_{\gamma}}^{\ast}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}H_{\mathrm{QQ}}\bigl{(}\gamma\bigr{)}\psi_{\Gamma_{\gamma^{\prime}}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$
$\displaystyle+\iint
d\boldsymbol{r}_{i}d\boldsymbol{r}_{j}\psi_{\Gamma_{\gamma}}^{\ast}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}i\left[l_{z}\bigl{(}\boldsymbol{r}_{i}\bigr{)}+l_{z}\bigl{(}\boldsymbol{r}_{j}\bigr{)},H_{\mathrm{QQ}}\bigl{(}\gamma\bigr{)}\right]$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\times\psi_{\Gamma_{\gamma^{\prime}}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}\omega_{xy}$
$\displaystyle+\iint
d\boldsymbol{r}_{i}d\boldsymbol{r}_{j}\psi_{\Gamma_{\gamma}}^{\ast}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$
$\displaystyle\qquad\qquad\times\left(-\frac{1}{2}\right)\left[l_{z}\bigl{(}\boldsymbol{r}_{i}\bigr{)}+l_{z}\bigl{(}\boldsymbol{r}_{j}\bigr{)},\left[l_{z}\bigl{(}\boldsymbol{r}_{i}\bigr{)}+l_{z}\bigl{(}\boldsymbol{r}_{j}\bigr{)},H_{\mathrm{QQ}}\bigl{(}\gamma\bigr{)}\right]\right]$
$\displaystyle\qquad\qquad\qquad\qquad\times\psi_{\Gamma_{\gamma^{\prime}}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}\bigl{(}\omega_{xy}\bigr{)}^{2}.$
(77)

Using the commutation relation of Eq. (29) for the Pauli matrices, we reduce
the perturbation Hamiltonian $H_{\mathrm{rot}}\bigl{(}\omega_{xy}\bigr{)}$
which depends on the rotation $\omega_{xy}$ as follows:

$\displaystyle H_{\mathrm{rot}}\bigl{(}\omega_{xy}\bigr{)}$
$\displaystyle=2\bigl{(}1-\gamma\bigr{)}\sum_{i\leq
j}J_{\mathrm{Q}}\bigl{(}\boldsymbol{r}_{i}-\boldsymbol{r}_{j}\bigr{)}$
$\displaystyle\qquad\times\left[O_{x^{\prime}y^{\prime}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}O_{x^{\prime
2}-y^{\prime 2}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}+O_{x^{\prime 2}-y^{\prime
2}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}O_{x^{\prime}y^{\prime}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}\right]\omega_{xy}$
$\displaystyle+4\bigl{(}1-\gamma\bigr{)}\sum_{i\leq
j}J_{\mathrm{Q}}\bigl{(}\boldsymbol{r}_{i}-\boldsymbol{r}_{j}\bigr{)}$
$\displaystyle\qquad\times\left[O_{x^{\prime 2}-y^{\prime
2}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}O_{x^{\prime 2}-y^{\prime
2}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}-O_{x^{\prime}y^{\prime}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}O_{x^{\prime}y^{\prime}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}\right]\bigl{(}\omega_{xy}\bigr{)}^{2}.$
(78)

From the first term in Eq. (4.6), which linearly depends on the rotation
$\omega_{xy}$, we identify the electric hexadecapole as

$\displaystyle
H_{z}^{\alpha}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}=O_{x^{\prime}y^{\prime}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}O_{x^{\prime
2}-y^{\prime 2}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}+O_{x^{\prime 2}-y^{\prime
2}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}O_{x^{\prime}y^{\prime}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}.$
(79)

The hexadecapole
$H_{z}^{\alpha}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ of Eq.
(79) is invariant upon the exchange of positions $\boldsymbol{r}_{i}$ and
$\boldsymbol{r}_{j}$ but is antisymmetric upon the exchange of coordinates
$x^{\prime}$ and $y^{\prime}$. In Fig. 5(c), we schematically show the
interaction between the hexadecapole
$H_{z}^{\alpha}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ of Eq.
(79) carried by the two-electron states and the rotation $\omega_{xy}$ of the
transverse acoustic phonons. The hexadecapole
$H_{z}^{\alpha}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ operates
on the two-electron state that is bound by the quadrupole interaction of Eq.
(4.4), while the hexadecapole
$H_{z}^{\alpha}\bigl{(}\boldsymbol{r}^{\prime}\bigr{)}$ in Eq. (34) acts on
the one-electron state that is trapped in the CEF Hamiltonian of the central
force. Note that the hexadecapole
$H_{z}^{\alpha}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ and
rotation $\omega_{xy}$ commonly belong to the $A_{2}$ symmetry of the point
group symmetry $D_{2d}$. The coefficient of the second term proportional to
$\bigl{(}\omega_{xy}\bigr{)}^{2}$ in Eq. (4.6) represents the energy
modulation of the anisotropic quadrupole interaction.

Comparing the term proportional to the rotation $\omega_{xy}$ in the
expressions of Eq. (4.6) with that in Eq. (4.6), we write the Heisenberg
equation of the time derivative for the total angular momentum
$l_{z}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}=l_{z}\bigl{(}\boldsymbol{r}_{i}\bigr{)}+l_{z}\bigl{(}\boldsymbol{r}_{j}\bigr{)}$
which is proportional to the hexadecapole
$H_{z}^{\alpha}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ as
follows:

$\displaystyle i\hbar\frac{\partial}{\partial
t}l_{z}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$
$\displaystyle=\left[l_{z}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)},H_{\mathrm{QQ}}\bigl{(}\gamma\bigr{)}\right]$
$\displaystyle=-2i\bigl{(}1-\gamma\bigr{)}J_{\mathrm{Q}}\bigl{(}\boldsymbol{r}_{i}-\boldsymbol{r}_{j}\bigr{)}H_{z}^{\alpha}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}.$
(80)

Furthermore, the time derivative of the angular momentum in Eq. (4.6) is
identified with the torque
$\tau_{xy}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ for the two-
electron state bound by the anisotropic quadrupole interaction of Eq. (4.4) as

$\displaystyle\hbar\frac{\partial}{\partial
t}l_{z}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}=\tau_{xy}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}.$
(81)

Note that the torque
$\tau_{xy}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ in Eq. (81)
vanishes for the special case of the ideal isotropic quadrupole interaction
with $\gamma=1$ of Eq. (4.4) identified with the ideal $xz$ model. The
hexadecapole
$H_{z}^{\alpha}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$
proportional to the time derivative of the angular momentum in Eq. (4.6)
conserves the time-reversal symmetry, while the angular momentum
$l_{z}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ itself breaks the
time-reversal symmetry.

The hexadecapole-rotation interaction linearly coupled to the rotation
$\omega_{xy}$ in Eq. (4.6) has the following matrix elements for the two-
electron states of
$\psi_{\Gamma_{\gamma}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$:

$\displaystyle\left\langle\Gamma_{\gamma}\left|H_{\mathrm{rot}}^{0}\bigl{(}\omega_{xy}\bigr{)}\right|\Gamma_{\gamma^{\prime}}\right\rangle$
$\displaystyle=2(1-\gamma)\sum_{i\leq j}\iint
d\boldsymbol{r}_{i}d\boldsymbol{r}_{j}\psi_{\Gamma_{\gamma}}^{\ast}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}J_{\mathrm{Q}}\bigl{(}\boldsymbol{r}_{i}-\boldsymbol{r}_{j}\bigr{)}H_{z}^{\alpha}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$
$\displaystyle\qquad\qquad\times\psi_{\Gamma_{\gamma^{\prime}}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}\omega_{xy}$
$\displaystyle=2(1-\gamma)\sum_{i\leq j}\iint
d\boldsymbol{r}d\boldsymbol{r}^{\prime}\psi_{\Gamma_{\gamma}}^{\ast}\bigl{(}\boldsymbol{r}-\boldsymbol{R}_{i},\boldsymbol{r}^{\prime}-\boldsymbol{R}_{j}\bigr{)}$
$\displaystyle\qquad\qquad\times
J_{\mathrm{Q}}\bigl{(}\boldsymbol{r}_{i}-\boldsymbol{R}_{i}-\boldsymbol{r}_{j}+\boldsymbol{R}_{j}\bigr{)}H_{z}^{\alpha}\bigl{(}\boldsymbol{r}-\boldsymbol{R}_{i},\boldsymbol{r}^{\prime}-\boldsymbol{R}_{j}\bigr{)}$
$\displaystyle\qquad\qquad\qquad\qquad\times\psi_{\Gamma_{\gamma^{\prime}}}\bigl{(}\boldsymbol{r}-\boldsymbol{R}_{i},\boldsymbol{r}^{\prime}-\boldsymbol{R}_{j}\bigr{)}\omega_{xy}.$
(82)

In the integral of Eq. (4.6), we respectively take coordinates
$\boldsymbol{r}$ and $\boldsymbol{r}^{\prime}$ instead of
$\boldsymbol{r}_{i}=\boldsymbol{r}-\boldsymbol{R}_{i}$ and
$\boldsymbol{r}_{j}=\boldsymbol{r}^{\prime}-\boldsymbol{R}_{j}$ for the
positions $\boldsymbol{R}_{i}$ and $\boldsymbol{R}_{j}$ of Fe2+ ions. As
mentioned for Eq. (4.6), the intra-atomic quadrupole interaction between
electrons accommodated in the same Fe2+ ion with
$\boldsymbol{R}_{i}=\boldsymbol{R}_{j}$ is dominant over the inter-atomic
quadrupole interaction with $\boldsymbol{R}_{i}\neq\boldsymbol{R}_{j}$. The
hexadecapole formed by the intra-atomic quadrupole interaction plays an
important role in the appearance of the hexadecapole ordering in the
superconducting state.

Among the two-electron states listed in Eqs. (4.6)-(4.6), we take the spin-
singlet states consisting of the two-electron states
$\psi_{\Gamma_{\gamma}}^{0,0}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}=\psi_{\Gamma_{\gamma}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$
for the orbital symmetries of $\Gamma_{\gamma}=A_{1}$, $B_{2}$, and $B_{1}$
compatible with the magnetic robustness of the present system, while we
disregard the spin-triplet state of $\psi_{A_{2}}^{S=1,\
S_{z}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ with $S_{z}=1$,
$0$, and $-1$. We deduce the hexadecapole-rotation interaction of Eq. (71) for
the spin-singlet states as

$\displaystyle\boldsymbol{H}_{\mathrm{rot}}^{0}\bigl{(}\omega_{xy}\bigr{)}=-4\bigl{(}1-\gamma\bigr{)}\sum_{i\leq
j}J_{\mathrm{Q}}^{ij}\bordermatrix{&\psi_{A_{1}}&\psi_{B_{2}}&\psi_{B_{1}}\cr&0&0&0\cr&0&0&1\cr&0&1&0}\omega_{xy}.$
(83)

Here, the matrix elements
$\bigl{\langle}\Gamma_{\gamma}|H_{\mathrm{rot}}^{0}\bigl{(}\omega_{xy}\bigr{)}|\Gamma_{\gamma^{\prime}}\bigr{\rangle}$
in Eq. (83) for the two-electron state
$\psi_{\Gamma_{\gamma}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$
with the symmetries of $\Gamma_{\gamma}=A_{1}$, $B_{2}$, and $B_{1}$ and the
spin singlet of S = 0 of Eqs. (4.6)-(4.6) are calculated. The hexadecapole
$H_{z}^{\alpha}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ with the
$A_{2}$ symmetry possesses the off-diagonal elements between the two-electron
states $\psi_{B_{2}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ and
$\psi_{B_{1}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$, but these
matrix elements vanish for the state
$\psi_{A_{1}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$. This is
confirmed by the symmetry property of $B_{2}\otimes B_{1}=A_{2}$.

We have particular interest in the interplay of the rotation $\omega_{xy}$
with the appearance of the hexadecapole ordering and the superconductivity.
Employing the unitary transformation for Eq. (83), we obtain the diagonal
representation of the hexadecapole-rotation interaction as

$\displaystyle\boldsymbol{H}_{\mathrm{rot}}\bigl{(}\omega_{xy}\bigr{)}=-4\left(1-\gamma\right)\sum_{i\leq
j}J_{\mathrm{Q}}^{ij}\bordermatrix{&\psi_{+}&\psi_{-}\cr&1&0\cr&0&-1}\omega_{xy}.$
(84)

Here, we adopt the wave function
$\psi_{h}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ with $h=\pm$
for the two-electron states as

$\displaystyle\psi_{\pm}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$
$\displaystyle=\frac{1}{\sqrt{2}}\left[\psi_{B_{2}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}\pm\psi_{B_{1}}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}\right]$
$\displaystyle=\frac{1}{\sqrt{2}}\left[\frac{\psi_{y^{\prime}z}\bigl{(}\boldsymbol{r}_{i}\bigr{)}\pm\psi_{zx^{\prime}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}}{\sqrt{2}}\psi_{y^{\prime}z}\bigl{(}\boldsymbol{r}_{j}\bigr{)}\right.$
$\displaystyle\left.\qquad\qquad\pm\frac{\psi_{y^{\prime}z}\bigl{(}\boldsymbol{r}_{i}\bigr{)}\mp\psi_{zx^{\prime}}\bigl{(}\boldsymbol{r}_{i}\bigr{)}}{\sqrt{2}}\psi_{zx^{\prime}}\bigl{(}\boldsymbol{r}_{j}\bigr{)}\right]$
$\displaystyle\qquad\times\frac{1}{\sqrt{2}}\left[\alpha\bigl{(}\boldsymbol{r}_{i}\bigr{)}\beta\bigl{(}\boldsymbol{r}_{j}\bigr{)}-\beta\bigl{(}\boldsymbol{r}_{i}\bigr{)}\alpha\bigl{(}\boldsymbol{r}_{j}\bigr{)}\right].$
(85)

By analogy with the Zeeman term for the magnetic dipole moments in a magnetic
field, we identify the wave function
$\psi_{+}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ with the
eigenstate for the hexadecapole corresponding to the right-hand rotation of
$h=+$ and the wave function
$\psi_{-}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ with
eigenstate for the hexadecapole corresponding to the left-hand rotation of
$h=-$. The hexadecapole-rotation interaction of Eq. (84) indicates that the
rotation $\omega_{xy}$ acts as a symmetry-breaking field on the hexadecapole
carried by the two-electron states
$\psi_{+}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ and
$\psi_{-}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ of Eq. (4.6).

In order to properly describe the hexadecapole ordering, it is convenient to
use the one-electron states $\lambda_{\pm}\left(\boldsymbol{r}_{i}\right)$
described in terms of the spherical harmonics of
$Y_{l}^{m}\left(\theta_{i},\varphi_{i}\right)$ with the orbital quantum number
$l=2$ and the azimuthal quantum number $m=\pm 1$ and the radial part of
$f_{d}\left(r_{i}\right)$ as

$\displaystyle\lambda_{\pm}\left(\boldsymbol{r}_{i}\right)$
$\displaystyle=f_{d}\left(r_{i}\right)Y_{2}^{\pm
1}\left(\theta_{i},\varphi_{i}\right)$
$\displaystyle=-\frac{1}{\sqrt{2}}\left[i\psi_{y^{\prime}z}\left(\boldsymbol{r}_{i}\right)\pm\psi_{zx^{\prime}}\left(\boldsymbol{r}_{i}\right)\right].$
(86)

We rewrite the two-electron states
$\psi_{\pm}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ of Eq. (4.6)
in terms of the one-electron wave function
$\lambda_{\pm}\left(\boldsymbol{r}_{i}\right)$ of Eq. (4.6) as

$\displaystyle\psi_{\pm}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$
$\displaystyle=\frac{1}{\sqrt{2}}\left[e^{\mp
i\frac{3\pi}{4}}\lambda_{+1}\bigl{(}\boldsymbol{r}_{i}\bigr{)}\lambda_{+1}\bigl{(}\boldsymbol{r}_{j}\bigr{)}+e^{\pm
i\frac{3\pi}{4}}\lambda_{-1}\bigl{(}\boldsymbol{r}_{i}\bigr{)}\lambda_{-1}\bigl{(}\boldsymbol{r}_{j}\bigr{)}\right]$
$\displaystyle\qquad\qquad\times\frac{1}{\sqrt{2}}\left[\alpha\bigl{(}\boldsymbol{r}_{i}\bigr{)}\beta\bigl{(}\boldsymbol{r}_{j}\bigr{)}-\beta\bigl{(}\boldsymbol{r}_{i}\bigr{)}\alpha\bigl{(}\boldsymbol{r}_{j}\bigr{)}\right].$
(87)

Note that the eigenstates of
$\psi_{\pm}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ of Eq. (4.6)
consist of the two-electron states of
$\lambda_{+1}\bigl{(}\boldsymbol{r}_{i}\bigr{)}\lambda_{+1}\bigl{(}\boldsymbol{r}_{j}\bigr{)}$
and
$\lambda_{-1}\bigl{(}\boldsymbol{r}_{i}\bigr{)}\lambda_{-1}\bigl{(}\boldsymbol{r}_{j}\bigr{)}$,
which are superposed on each other while maintaining an orthogonal relation
with a phase difference of $\pm 3\pi/2$.

Taking the spin-singlet state in Eq. (4.6) into account, we introduce
annihilation operators of $B_{i,j,\pm,\sigma,\overline{\sigma}}$ and creation
operators of $B_{i,j,\pm,\sigma,\overline{\sigma}}^{\dagger}$ for the two-
electron eigenstates
$\psi_{\pm}\bigl{(}\boldsymbol{r}_{i},\boldsymbol{r}_{j}\bigr{)}$ of the
hexadecapole as

$\displaystyle
B_{i,j,\pm,\sigma,\overline{\sigma}}=\frac{1}{\sqrt{2}}\left(e^{\mp
i\frac{3\pi}{4}}l_{j,+1,\overline{\sigma}}l_{i,+1,\sigma}+e^{\pm
i\frac{3\pi}{4}}l_{j,-1,\overline{\sigma}}l_{i,-1,\sigma}\right),$ (88)
$\displaystyle
B_{i,j,\pm,\sigma,\overline{\sigma}}^{\dagger}=\frac{1}{\sqrt{2}}\left(e^{\pm
i\frac{3\pi}{4}}l_{i,+1,\sigma}^{\dagger}l_{j,+1,\overline{\sigma}}^{\dagger}+e^{\mp