# Disorder effects in the Kitaev-Heisenberg model

Ayushi Singhania  Institute for Theoretical Solid State Physics, IFW Dresden,
01069 Dresden, Germany  Jeroen van den Brink  Institute for Theoretical Solid
State Physics, IFW Dresden, 01069 Dresden, Germany  Department of Physics,
Technical University Dresden, 01069 Dresden, Germany  Satoshi Nishimoto
Institute for Theoretical Solid State Physics, IFW Dresden, 01069 Dresden,
Germany  Department of Physics, Technical University Dresden, 01069 Dresden,
Germany

###### Abstract

We study the interplay of disorder and Heisenberg interactions in Kitaev model
on honeycomb lattice. The effect of disorder on the transition between Kitaev
spin liquid and magnetic ordered states as well as the stability of magnetic
ordering is investigated. Using Lanczos exact diagonalization we discuss the
consequences of two types of disorder: (i) random-coupling disorder and (ii)
singular-coupling disorder. They exhibit qualitatively similar effects in the
pure Kitaev-Heisenberg model without long-range interactions. The range of
spin liquid phases is reduced and the transition to magnetic ordered phases
becomes more crossover-like. Furthermore, the long-range zigzag and stripy
orderings in the clean system are replaced by their three domains with
different ordering direction. Especially in the crossover range the
coexistence of magnetically ordered and Kitaev spin-liquid domains is
possible. With increasing the disorder strength the area of domains becomes
smaller and the system goes into a spin-glass state. However, the disorder
effect is different in magnetically ordered phases caused by long-range
interactions. The stability of such magnetic ordering is diminished by
singular-coupling disorder, and accordingly, the range of spin-liquid regime
is extended. This mechanism may be relevant to materials like $\alpha$-RuCl3
and H3LiIr2O6 where the zigzag ground state is stabilized by weak long-range
interactions. We also find that the flux gap closes at a critical disorder
strength and vortices appears in the flux arrangement. Interestingly, the
vortices tend to form kinds of commensurate ordering.

##  I Introduction

Quantum spin liquids (QSLs) are entangled states of matter, where quantum
fluctuations prevent the formation of magnetic order down to lowest
temperatures. They have attracted substantial attention in the developments of
topological quantum computing due to its remarkable properties, including long
range entanglement, topological ground state and fractionalized excitations
[1, 2]. One of the most promising theoretical models to realise QSLs, was
proposed by Kitaev in his seminal paper [3]. The model is based on a two
dimensional honeycomb lattice with bond anisotropic interactions among quantum
spin-1/2 particles. Remarkably, such Kitaev-type interactions are predicted to
be realized in transition metal oxides with spin-orbit couplings [4]. Since
then, layered honeycomb materials A2IrO3 (A = Na, Li) [5, 6] and
$\alpha-$RuCl3 [7] have been thoroughly investigated as the candidate Kitaev
materials [8, 9]. Even though thermodynamic and dynamical properties of these
materials match the theoretical prediction of Kitaev model, all the candidate
materials are known to magnetically order at low temperatures [10, 11].
Magnetic orders in such materials are manifested by the interactions beyond
the Kitaev interactions [12, 13, 14].

In addition, presence of disorder like vacancies, impurities, lattice
distortions, and stacking faults, inevitable in real materials, are often
responsible for the instability of QSL [15, 16, 17]. Recently, a new class of
intercalated compounds have been synthesized which are even more susceptible
to disorder H3LiIr2O6, Cu2IrO3, and Ag3LiIr2O6 [18, 19]. Role of disorder in
Kitaev materials has gained special interest since the experimental
observation of low temperature divergence in specific heat of the proximate
Kitaev material, H3LiIr2O6, is understood as the consequence of bond disorder
in the Kitaev model [20, 21]. Theoretically, it is known that bond randomness
in the Kitaev model affects the low-energy dynamical spin-correlators [22].
Thermodynamic properties like specific heat and susceptibility, spin
transport, and spin dynamics at low energy excitations of the Kitaev model
with bond randomness and site dilution have been also studied at finite but
low temperatures [23, 24, 24, 16, 25, 26, 27]. Moreover, the effect of bond
randomness on the zero-temperature topological properties of the Kitaev model
[28, 29] as well as the spin excitations in a diluted Kitaev system containing
spin vacancies and bond randomness [30] have been discussed.

While a lot of studies have been devoted to the thermodynamics (relevant for
experiments) of disordered Kitaev models, very little is understood about the
ground state properties of the Kitaev model in the presence of disorder.
Furthermore, presence of any type of disorder such as impurity or stacking
faults must affect the interactions beyond the Kitaev limit. However, even
less is known about the disorder effects, for example, in the presence of
Heisenberg interaction. Once the Heisenberg interaction is added in the Kitaev
model (referred as Kitaev-Heisenberg model), variety of magnetically ordered
phases are stabilized as a consequence of order by disorder mechanism.
Understanding the effect of disorder in the Kitaev-Heisenberg model is crucial
from the point of view of interplay between frustration and disorder.

Motivated by this, we employ numerical diagonalization of 24-site cluster on a
honeycomb lattice, to study the effects of disorder in the ground state of
Kitaev-Heisenberg model. This serves as a putative minimal model for several
QSL candidates. We emphasise on two different types of bond-disorder
introduced in both Kitaev $and$ Heisenberg interactions. In the pristine
Kitaev-Heisenberg model, Kitaev QSL state survives the small perturbations
like Heisenberg interaction[31]. Here, we show that the QSL remains robust
against weak disorder and Heisenberg interaction, the stability of QSL region
is narrowed with increasing disorder strength. The long-range zigzag and
stripy orderings are replaced by a mixture of their three domains with
different ordering direction. The area of domains is reduced with increasing
disorder strength and the system goes into a spin-glass state. Moreover, we
show that a flux state can be induced by bond-disorder in the ground state.
The transition from flux-free state to random flux state in the presence of
bond disorder has already been predicted just in the Kitaev limit [22, 23].
Interestingly, this random-flux state persists even for a small admixing of
nearest-neighbor Heisenberg interaction. We also find that disorder in
presence of further neighbor interaction can indeed result in a QSL or
disordered state. This is relevant to experimental observations of H3LiIr2O6
[32].

The paper is organized as follows: In Sec. II, we begin by introducing the
type of disorders considered in the Kitaev-Heisenberg model and the details of
numerical methods applied. In Sec. III, we present the results of disorder
effects on ground state phases in K-H parameter space. Moreover, we also
discuss the effect of further neighbour interactions on the ground state of
disordered KH model. Section IV is devoted to the summary.

Figure 1: (a) 24-site periodic cluster used in our exact-diagonalization
calculations. A black rectangle shows the unit cell. Numbers 1-6 correspond to
labeling for the hexagonal plaquette operator [Eq. (3)]. (b) Corresponding
Brillouin zone, where inner hexagon shows first Brillouin zone. Magnetic Bragg
peaks for zigzag and stripy states are also shown.

##  II Model and method

###  II.1 Kitaev-Heisenberg Model

The Kitaev honeycomb model is described by anisotropic Ising-type interactions
among spin-1/2 degrees of freedom on a honeycomb lattice, where the
interaction in spin space ($\gamma=x,y,z$) is dictated by the bond direction.
The bond-directional Ising-type interactions are called Kitaev interactions.
The model is exactly solvable with the ground state being QSL and its
elementary excitations described by Majorana fermions coupled to the $Z_{2}$
gauge field. In general, the systems exhibiting the Kitaev interactions are
referred as Kitaev materials. However, as revealed from extensive research,
real Kitaev materials possess interactions beyond Kitaev, such as diagonal
exchange couplings, i.e., Heisenberg terms, as well as off-diagonal couplings,
i.e., so-called $\Gamma$ terms [12, 13]. Nevertheless, a QSL ground state
state is maintained even in the presence of weak Heisenberg terms [31]. In
order to derive a common understanding on the disorder effects in Kitaev
systems, we here restrict ourselves to models with the Heisenberg and Kitaev
terms, i.e., the Kitaev-Heisenberg model.

####  II.1.1 Clean case

The Hamiltonian of the Kitaev-Heisenberg model reads as

$\mathcal{H}_{\rm KH}=\sum_{\langle
ij\rangle\gamma}2K_{ij}^{\gamma}S_{i}^{\gamma}S_{j}^{\gamma}+\sum_{\langle
i,j\rangle}J_{ij}\mathbf{S}_{i}\cdot\mathbf{S}_{j}$ (1)

where $K_{ij}^{\gamma}$ is the Kitaev interaction of nearest neighbor spins on
three different bonds $\gamma=x,y,z$ and $J_{ij}$ is the Heisenberg
interaction between nearest neighbours. We assume all the Kitaev couplings are
isotropic with $K_{ij}^{x}=K_{ij}^{y}=K_{ij}^{z}=K$ and the nearest-neighbor
Heisenberg coupling is constant with $J_{ij}=J$ in the clean case. For
simplicity, we parameterize as following: $J=A\cos\phi$ and $K=A\sin\phi$,
where $A=\sqrt{J^{2}+K^{2}}$ is used as an energy unit. The ground state phase
diagram is known to host four ordered and two QSL phases as a function of
$\phi$: Néel ($-0.19\pi\lesssim\phi\lesssim 0.49\pi$), AFM QSL
($0.49\pi\lesssim\phi\lesssim 0.51\pi$), zigzag ($0.51\pi\lesssim\phi\lesssim
0.90\pi$), FM ($0.90\pi\lesssim\phi\lesssim 1.40\pi$), FM QSL
($1.40\pi\lesssim\phi\lesssim 1.58\pi$), and stripy
($1.58\pi\lesssim\phi\lesssim 1.89\pi$) [31]. In this paper, we consider two
types of bond disorders explained below, and mainly focus on the effect of
disorders around the two Kitaev QSL phases, i.e., around $\phi=\pi/2,3\pi/2$.

####  II.1.2 Random-coupling disorder

The first type of disorder is a random distribution of bond strength (called
“random-coupling disorder”), which commonly appear due to the structural
disorders in crystals. This can be achieved by changing the parameter
$A\rightarrow(A+\delta A_{ij}$) for all nearest-neighbor bonds, such that
$J\rightarrow(A+\delta A_{ij})\cos\phi$ and $K\rightarrow(A+\delta
A_{ij})\sin\phi$. Namely, the Kitaev and Heisenberg couplings of 36 nearest-
neighbor bonds in our 24-site periodic cluster (see Fig 1(a)) take different
values from each other. Nevertheless, note that the ratio between $K$ and $J$
remains as $K/J=\tan\phi$ for each bond. The random deviation from the
original bond, $\delta A_{ij}$, is defined by a Gaussian distribution
$N(0,\sigma^{2})$ where $\sigma$ is the standard deviation. The disorder
strength is controlled by $\Delta=2\sigma$. The cases of $\Delta=0.125$,
$0.25$, and $0.5$ are studied in this paper. For this kind of bond disorder, a
box probability has been also frequently used in the previous studies but just
for the Kitaev limit, $K/J\to\infty$. Although the effect of disorder with a
box probability may be somewhat a little more drastic than that with a
Gaussian distribution, their results are expected to be qualitatively similar.

####  II.1.3 Singular-coupling disorder

The second type of disorder is a random mixture of two kinds of bond strengths
(called “singular-coupling disorder”). This can be achieved by replacing a
part of the original bonds with $r$-times stronger bonds, where the fraction
of replaced bonds to the total bonds is defined by $n$. In other words,
randomly selected fraction $n$ of original bonds are changed as $J\rightarrow
rJ$ and $K\rightarrow rK$. We call $n$ disorder density and $r$ disorder
strength. Since we use a finite-size periodic cluster as described below, it
would be reasonable to keep equal number of replaced $x$, $y$ and $z$ bonds.
In fact, we vary the number of replaced bonds in steps of 3. The case of $n=0$
and/or $r=1$ corresponds to the clean limit. The case of $n=1$ also
corresponds to the clean limit but the energy scales as $rA$ instead of $A$ in
the original clean limit $n=0$. In this paper, the effect of singular-coupling
disorder is widely studied from $n=0$ to $1$ for $r=2$, $4$, and $10$.
Interestingly, the appearance of this type of disorder with $r=4$ was
predicted due to the presence of hydrogen cation vacancies in H3LiIr2O6 [33].

###  II.2 Method

We employ the Lanczos exact diagonalization method of the model described by
Eq. (1) in the presence of disorder. A 24-site cluster with periodic boundary
conditions is used [Fig 1(a)]. For each parameter set, we perform the
simulations for 30-50 random realizations and take an average over them if
needed. This setup is sufficient to investigate the instability of Kitaev QSL
around the Kitaev limit because long-range magnetic orderings are not
considered as discussed below. We calculate the ground-state energy $E_{0}$,
spin-spin correlation functions
$\langle\mathbf{S}_{i}\cdot\mathbf{S}_{j}\rangle$, static spin structure
factor $S({\bf Q})$, and the expectation value of flux operator $\langle
W\rangle$.

Peak positions in the second derivative of ground-state energy with respect to
particular parameters provide valuable information on the position of phase
boundaries. Also, the width of the peaks may reflect the ’sharpness’ of phase
transitions; roughly speaking, a sharp (broad) peak indicates a first-order
(second-order or continuous) transition. To identify the magnetic structure
for each phase, we calculate the static spin structure factor

$S({\bf Q})=\sum_{ij}\langle{\tilde{\bf S}}_{i}\cdot{\tilde{\bf
S}}_{j}\rangle\exp[i{\bf Q}\cdot({\bf r}_{i}-{\bf r}_{j})],$ (2)

where $N$ is the number of sites in the periodic cluster ($N=24$) and ${\bf
r}_{i}$ is the position of site $i$. For ordered phases the magnetic structure
can be determined from the reciprocal-space Bragg-peak positions. The
Brillouin zone and magnetic Bragg peaks for zigzag and stripy states are shown
in Fig 1(b). Note that in our calculations for the clean limit any spatial and
spin symmetries are not broken, so that the structure factor exhibits all
Bragg peaks for three ordering directions. However, as shown below, the
translation and rotation symmetries can be broken when the disorder is
introduced. In QSL phase $S({\bf Q})$ is nearly structureless and the
intensity is relatively small for ${\bf Q}$’s. Furthermore, to see how close
is the system to the Kitaev limit, we calculate the expectation value of flux
operator. The flux operator is defined on a hexagonal plaquette by

$W_{p}=2^{6}S_{1}^{x}S_{2}^{y}S_{3}^{z}S_{4}^{x}S_{5}^{y}S_{6}^{z},$ (3)

where a possible numbering of sites is denoted in Fig 1(a). In the Kitaev
limit, the flux operator is a conserved quantity, commuting with the
Hamiltonian, having values exactly $\pm 1$. Away from the Kitaev QSL regime,
this value drops to zero. Therefore, we can also use this quantity as a
measure of overlap between an observed QSL and ideal Kitaev QSL.

##  III Results

Figure 2: Second derivative of the ground-state energy as a function of $\phi$
around (a) AFM ($\phi=0.5\pi$) and (b) FM ($\phi=1.5\pi$) Kitaev points.
Averaged expectation value of the plaquette operator over all 12 plaquettes
around (c) AFM and (d) FM Kitaev points. For the disordered cases both the
quantities are averaged over 50 disorder realizations.

###  III.1 Random-coupling disorder

We begin by discussing the consequence of the random-coupling disorder. Since
we are particularly interested in the possible emergence of QSL by disorder,
we focus on the parameter regions around the Kitaev points
$\phi=0.5\pi,1.5\pi$; where any magnetic order can be avoided due to the fully
dominant Kitaev interactions over the Heisenberg ones. Fig. 2(a,b) show second
derivative of the ground-state energy with respect to $\phi$, i.e.,
$-\partial^{2}E_{0}/\partial\phi^{2}$, as a function of $\phi$ for several
values of the disorder strength $\Delta$ around the AFM ($\phi=0.5\pi$) and FM
($\phi=1.5\pi$) Kitaev points. Similarly, expectation values of the flux
operator (3), i.e., $\langle W_{p}\rangle$, are plotted in Fig. 2(c,d). In the
presence of disorder, the averaged values of
$-\partial^{2}E_{0}/\partial\phi^{2}$ and $\langle W_{p}\rangle$ are taken
over 50 random realisations.

Around AFM Kitaev point, assuming that a peak position in
$-\partial^{2}E_{0}/\partial\phi^{2}$ vs. $\phi$ indicates a critical point, a
transition from Néel to Kitaev QSL phases occurs at
$\phi=0.486\pi(\equiv\phi_{\rm Neel-QSL})$ and a transition from Kitaev QSL to
zigzag phases at $\phi=0.514\pi(\equiv\phi_{\rm QSL-zigzag})$ in the clean
limit ($\Delta=0$). We find that the phase transitions are hardly affected by
the disorder strength, $\Delta=0.25$. However, with further increasing
$\Delta$ to $0.5$, the deviation from the clean limit becomes prominent. The
critical points are changed to $\phi_{\rm Neel-QSL}=0.488\pi$ and $\phi_{\rm
QSL-zigzag}=0.511\pi$. Thus, the range of Kitaev QSL has shrunk by the
disorder. Even more interestingly, both peaks in
$-\partial^{2}E_{0}/\partial\phi^{2}$ at $\Delta=0.5$ are lower and broader
compared to the sharp peaks indicating first-order transitions in the clean
limit ($\Delta=0$). This means that the disorder makes the phase transitions
more crossover-like. It is also confirmed by a rather blunt behavior of
$\langle W_{p}\rangle$ around the critical points. Accordingly, the system may
be characterized by a short-range Néel phase in a range with peak width around
$\phi=\phi_{\rm Neel-QSL}$ and by a mixture of small zigzag domains and QSL
regions in a range with peak width around $\phi=\phi_{\rm QSL-zigzag}$. The
possibility of the formation of zigzag domains is discussed below. On the
other hand, Néel domain walls would be not created because the random-coupling
disorder does not lift the degeneracy of spin-rotation symmetry unlike that of
spatial-rotation symmetry.

The effect of disorder around FM Kitaev point seems to be even more sensitive
than around AFM Kitaev point. A transition from FM to Kitaev QSL phases occurs
at $\phi=1.40\pi(\equiv\phi_{\rm FM-QSL})$ and a transition from Kitaev QSL to
stripy phases at $\phi=1.58\pi(\equiv\phi_{\rm QSL-stripy})$ in the clean
limit; whereas, they are $\phi_{\rm FM-QSL}=1.42\pi$ and $\phi_{\rm QSL-
stripy}=1.56\pi$ at $\Delta=0.5$. Although the shifts of critical $\phi$
values may look tiny, we can find them somewhat considerable by thinking in
terms of $K/J$ as follows: The condition to achieve a QSL is modified from
$K/J>3.1$ or $K/J<-3.9$ in the clean limit to $K/J>3.9$ or $K/J<-5.2$ at
$\Delta=0.5$. Nevertheless, we can say that the effect of the random-coupling
disorder on magnetic order is relatively weak, which is consistent with the
previous study [34]. Then, as in the case around AFM Kitaev limit, peaks in
$-\partial^{2}E_{0}/\partial\phi^{2}$ are broadened. Thus, the system may be
characterized by a short-range FM phase in a range with peak width around
$\phi=\phi_{\rm FM-QSL}$ and by a mixture of small stripy domains and QSL
regions in a range with peak width around $\phi=\phi_{\rm QSL-stripy}$. The
possibility of the formation of stripy domains is discussed below.

It is also interesting that $|\langle W\rangle|$ keeps nearly $1$ in the both
AFM and FM Kitaev QSL phases despite the presence of sizeable disorder. In the
Kitaev limits $\phi=\pm 0.5\pi$ ($K=\pm 1$, $J=0$), the flux operator is still
commutative with the Hamiltonian even for finite $\Delta$. Thus, as confirmed
in Fig. 2(c,d), the ground state is still characterized by $|\langle
W\rangle|=1$ ($\forall$ hexagonal plaquette). Up to $\Delta=0.5$, we found
$\langle W\rangle$ to be the unity for all studied disorder samples.
Furthermore, as naively expected, we find that the nearest-neighbour spin-spin
correlations are finite and longer-range ones are zero as is the case with the
clean limit. Nevertheless, the value of the nearest-neighbour correlations can
deviate from the clean limit value
$|\langle\mathbf{S}_{i}\cdot\mathbf{S}_{j}\rangle_{\rm NN}|=0.1323$ depending
on the coupling value, namely,
$|\langle\mathbf{S}_{i}\cdot\mathbf{S}_{j}\rangle_{\rm NN}|>0.1323$ for the
bond with $\delta A>0$ and
$|\langle\mathbf{S}_{i}\cdot\mathbf{S}_{j}\rangle_{\rm NN}|<0.1323$ for the
bond with $\delta A<0$.

Figure 3: Static structure factors for (a) zigzag ($\phi=0.52\pi$) and (e)
stripy ($\phi=1.65\pi$) states in the clean limit ($\Delta=0$). Examples of
static spin structure factors with rotation-symmetry breaking for (b-d) zigzag
and (f-h) stripy states at $\Delta=0.5$. (i,j) Energy distributions of 50
disorder samples for $\phi=0.52\pi$ and $\phi=1.65\pi$ at $\Delta=0.5$.

Let us next consider how the ordered phases are affected by the random-
coupling disorder. In the clean limit, the zigzag and stripy ordered states
are associated with spatial rotation symmetry breaking. Yet, the rotation
symmetry is not broken in our calculations with 24-site periodic cluster and
three states with different ordering directions are energetically degenerate.
Fig. 3(a,e), show that the Bragg peaks for all three ordering directions in
zigzag and stripy order appear with the same intensity in the static structure
factor. One of the degenerate three states is then selected as a ground state
when its global magnetic ordering is constructed. However, the rotation
symmetry can be explicitly broken once the random-coupling disorder is
introduced. Examples of static spin structure factors with broken rotation
symmetry for the zigzag state ($\phi=0.52\pi$) and stripy state
($\phi=1.65\pi$) are plotted in Fig. 3(b-d) and (f-h), respectively. In a
disordered system, one of the symmetry-broken state with lowest energy tends
to be locally realized depending on the disorder distribution. This means in
the ordered phase, the three domains with different ordering directions can
coexist as a result of disorder. In Fig. 3(i,j) their energy distributions of
50 disorder samples are shown. For both of $\phi=0.52\pi$ and $1.65\pi$, the
variance of energy distribution is increased with increasing $\Delta$. The
local ordering structure with lower energy is strongly stabilized and that
with higher energy is easily excluded. Thus, the averaged area of zigzag (or
stripy) domains becomes smaller and the system is in more like a spin-glass
state for larger $\Delta$. It is because the short-range correlation is
dominated only by the local disorder structure when the disorder is strong. In
addition, a coexistence of four kinds of domains could be possible in the
crossover regions between QSL and magnetic ordered phases, namely, one QSL
domain and three ordering domains as mentioned above [35].

###  III.2 Singular-coupling disorder

Figure 4: (a-c) Second derivative of the ground-state energy and (d-f)
expectation value of the flux operator 3 as a function of $\phi$ around FM
Kitaev point ($\phi=1.5\pi$). Results are shown for various disorder strength
$r=2$, $4$, and $10$, and disorder density $n=0-1$.

Next, we turn to the case of singular-coupling disorder. For this case we can
study the effect of disorder as functions of the disorder strength $r$ and the
disorder density $n$. As in the case of random-coupling disorder, the impact
of singular-coupling disorder on the two phase transitions around the FM
Kitaev point $\phi=1.5\pi$ is investigated by estimating
$-\partial^{2}E_{0}/\partial\phi^{2}$ with $\phi$. Here, it is convenient to
renormalize the energy unit as $A\rightarrow A/(nr+(1-n))$. In Fig. 4(a-c), we
plot $-\partial^{2}E_{0}/\partial\phi^{2}$ as a function of $\phi$ with fixed
$r=2$, $4$, and $10$. The disorder density $n$ is varied from $0$ to $1$ at
intervals of $3/36$. Results shown in Fig. 4 are averaged values of
$-\partial^{2}E_{0}/\partial\phi^{2}$ and $\langle W_{\rm p}\rangle$ over 30
disorder samples. Interestingly, $-\partial^{2}E_{0}/\partial\phi^{2}$ is
nearly symmetric about $n=18/36$ for any $r$. This means that the qualitative
properties are determined only by the ratio of distinct bonds. At $n=0$ and
$n=1$, where the system is clean, the range of FM QSL is indicated by two
sharp peaks at $\phi=1.40\pi$ and $1.58\pi$. The transitions to FM
($\phi\lesssim 1.40\pi$) and stripy ($\phi\gtrsim 1.58\pi$) phases are of the
first order.

For weak disorder strength ($r=2$), we find that the QSL range is slightly
reduced in the whole density region $0<n<1$ and it is narrowest at $n=0.5$.
This reduction seems to increase continuously with increasing $r$. As a
result, for strong disorder strength ($r=10$), the QSL range almost vanishes
and a nearly direct transition between FM and stripy phases may occur around
$n=0.5$. The remaining QSL seems to be still of the Kitaev-type because of
$\langle W_{p}\rangle\approx 1$. Another effect of singular-coupling disorder
is that the peak of $-\partial^{2}E_{0}/\partial\phi^{2}$ indicating a
transition is significantly broadened even for weak disorder. This implies
that peak position for single disorder sample differs from each other. Thus,
transition in the disordered system becomes continuous or even a crossover-
like one. As in the case of random-coupling disorder, the system may be
characterised by a short-range FM or a short-range stripy ordering in the
crossover region. This might be interpreted as a QSL in the broad sense of the
term. It is rather surprising that the broadening of peak in
$-\partial^{2}E_{0}/\partial\phi^{2}$ seems to be almost independent of $r$.

Figure 5: (a-c) Representative static structure factors of rotation-symmetry
broken stripy ordering at $\phi=1.64\pi$, $n=0.5$ for $r=10$. (d) QSL-like
structure factor appearing at $\phi=1.55\pi$, $n=18/36$ for $r=10$. (e-h)
Energy distributions of 30 samples for various parameters, where the energy
value for the clean limit is set to be zero and the class interval is $\Delta
E_{0}=0.05$.

Furthermore, in order to see how the ordered state is affected by the
singular-coupling disorder, we study the static structure factor in the stripy
phase. Three representative structure factors at $\phi=1.64\pi$, $n=0.5$ for
$r=10$ are shown in Fig. 5(a-c), indicating a stripy state with broken
rotation symmetry. Since most of samples for any $r$ basically exhibit one of
the three structures, the system consists of the three stripy domains with
different ordering directions. As discussed in Sec. III.1, the variation of
energy distribution is related to the averaged area of stripy domains. In Fig.
5(e-g) the energy distributions of 30 disorder samples in the stripy phase
($\phi=1.64\pi$) are shown. As expected, the variance of energy distribution
increases with increasing $r$, and however, it clearly decreases with
increasing $n$. As a result, their energy distributions are completely
different between small and large $n$ regions, although the second derivative
of energy is symmetric about $n=0.5$. This may be explained by considering
random-singlet formations. At lower $n$, the energy can be easily lowered
because the stronger $r$ bonds are mostly isolated and spin-singlets are
formed on them. The direction of stripy ordering is closely related to the
distribution of random-singlets. Therefore, the area of stripy domains tends
to be small at smaller $n$ [Fig. 5(e)], and it would be relatively large at
larger $n$ [Fig. 5(g)]. For example, at $n=30/36$ the variance seems to be
very small and the energies of disorder samples are almost degenerate for any
$n$. This situation is close to the clean limit. Thus, the system consists of
large stripy domains if they exist. On the other hand, at $n=6/36$ and $r=10$
a wide distribution of the variance may indicate the appearance of small
stripy domains. This is because a local stripy configuration can be realised
against a global ordering if it significantly lowers the total energy. In sum,
a spin-glass state due to a mixture of tiny stripy domains may be achieved at
large $r$ and small $n$.

Let us look at the crossover region between stripy and Kitaev QSL phases in a
little more detail. We take the case of $r=10$ and $n=18/36$ for example. A
structureless structure factor as in Fig. 5(d) appears other than the three
symmetry broken ordered ones. Also, the variation of energy distribution is
relatively larger than that in the deep stripy phase [Fig. 5(h)]. Thus, a
coexistence of four kinds of domains may be possible in the crossover region
between Kitaev QSL and stripy phases, namely, one QSL domain and three stripy
ordering domains. At $\phi=1.55\pi$, 10 of 30 samples exhibit such a QSL-like
structure factor, and the remaining 20 samples exhibit either of stripy ones.
Assuming that $\langle W\rangle=0$ for the stripy samples and $\langle
W\rangle=1$ for the QSL-like samples, we may roughly expect $\langle
W\rangle\sim 0.33$. This is comparable to the actual value $\langle
W\rangle\sim 0.4$ for $\phi=1.55\pi$ at $r=10$ and $n=18/36$.

Figure 6: Typical examples of flux configurations for (a) flux-free, (b,c) 1/6
flux density, and (d-f) 1/3 flux density states. Number of samples exhibiting
flux states out of 30 samples, as a function of disorder density $n$ for (g)
$r=2$, (h) $r=4$, and (i) $r=10$

.

#### Appearance of flux state by disorder

We find an intriguing evolution of flux states with the disorder strength $r$
and disorder density $n$ in the Kitaev limit $\phi=3\pi/2$. The introduction
of disorder may mimic thermal ﬂuctuations in real materials. Since the
relations $W_{p}^{2}=1$, $[\mathcal{H},W_{p}]=0$, and
$[W_{p},W_{p^{\prime}}]=\delta_{pp^{\prime}}$ are fulfilled at $\phi=3\pi/2$
even in the presence of disorder, Therefore, the eigenstate for the system (1)
can be still characterized by a set of $\langle W_{p}\rangle$ for all
plaquettes, where $\langle W_{p}\rangle$ has a value of 1 or -1. For each
disorder sample we see the $\langle W_{p}\rangle$ values for 12 plaquettes in
our 24-site periodic cluster.

For the clean case ($r=1$), it is known that the ground state is given by
$\langle W_{p}\rangle=1$ ($\forall p$), the so-called flux-free state [Fig.
6(a)]. Then, the excited flux states can be created by flipping $\langle
W_{p}\rangle$ from $1$ to $-1$ for even number of plaquettes. The flux gap has
been estimated to be $E_{\rm gap}=0.066K$ as a lowest-energy excitation when
the values of $\langle W_{p}\rangle$ for two neighboring plaquettes are
flipped [22]. The ﬂux density is deﬁned by $n_{\rm F}=N_{\rm F}/N_{\rm p}$,
where $N_{\rm F}$ is the numbers of plaquettes exhibiting $\langle
W_{p}\rangle=-1$ and $N_{\rm p}$ is the total number of plaquettes. In
particular, the flux density $n_{\rm F}=1/6$ is achieved by a 24-site periodic
cluster containing 2 plaquettes with $\langle W_{p}\rangle=-1$ and 10
plaquettes with $\langle W_{p}\rangle=1$, i.e., $N_{\rm F}=2$ and $N_{\rm
p}=12$. Similarly, $n_{\rm F}=1/3$ corresponds to $N_{\rm F}=4$ and $N_{\rm
p}=12$ in our calculations. We refer to $n_{\rm F}=1/6$ and $n_{\rm F}=1/3$ as
2-flux and 4-flux state, respectively. Typical examples of flux configurations
for (excited) 2-flux and 4-flux states are illustrated in Fig. 6(b-f). The
appearance of fluxes in the ground state implies the closing of flux gap by
disorder. It is consistent with the previous study [29].

Flux states appears as a ground state on introduction of singular-coupling
disorder. The number of samples exhibiting flux states out of 30 samples is
shown as a function of disorder density $n$ for several disorder strengths in
Fig. 6(g-i). For a weak disorder strength ($r=2$) only $\sim 10-20\%$ of
disorder samples exhibit flux states. Interestingly, all of them are the 4-ﬂux
states shown in Fig. 6(e). Since the effect of disorder is most highlighted
around $n=0.5$, we find the probability of observing a flux state is also
maximum around $n=0.5$. With further increasing $r$ up to $4$ the probability
of finding flux states increases. The 4-flux configuration still accounts for
a large fraction of observed flux states, although a different pattern for
4-flux configurations such as in Fig. 6(f) as well as the 2-flux configuration
like in Fig. 6(b,c) appear as minority. In other words, ordered 4-flux state
would tend to be energetically more favourable than the other flux states.
This may be related to the “Majorana insulating” state [36]. Looking at all of
the flux configurations, the fluxes seem to keep distance from each other.
Thus, a kind of repulsion may exist between fluxes, and the energetically-
favoured ordered-flux state could be a consequence of the repulsive
interaction (also see the Appendix VI.1, VI.2). This also implies no
possibility of a phase-separation-like flux configuration induced by disorder.
Note that, however, the ordered-flux state is not globally stabilized in our
disordered system because various flux patterns locally exist depending on the
disorder realization, for example, as shown in in Fig. 6(b-f).

For a very strong disorder strength ($r=10$) the 4-flux and 2-flux states
appear with similar probability, and the flux-free state still exists.
Furthermore, various flux configurations are mixed. For example, a periodic
configurations like in Fig. 6(f) are also contained. These results indicate a
trend to the random-flux states [22, 37]. This can be naively understood if we
can assume that the effect of disorder is essentially the same as thermal
ﬂuctuations. Our focus here has been only the ground state, we show that
random flux state appears as a consequence of disorder itself, which indicates
that the flux gap closes at some disorder strength, $r$ for a particular
disorder configuration. Eventually, the presence of disorder in a random flux
sector can explain the abundance of low energy states in H3LiIr2O6 compound
[32]. Thus, the low-lying flux and spin excitations [38] of the Kitaev-
Heisenberg model should be studied in future.

###  III.3 Further neighbor interactions

Figure 7: (a) Second derivative of the ground-state energy as a function of
$J_{2,3}$ for the clean and disordered systems. The values for the system is
averaged ones over 30 disorder samples. (b) Second derivative of the ground-
state energy as a function of $J_{2,3}$ and (c) the corresponding static spin
structure factors for two representative disorder samples. (d) Expectation
values of the flux operator, $\langle W_{\rm p}\rangle$, as a function of
$J_{2,3}$ for the clean and disordered systems. For disordered systems, the
values for each disorder sample are also plotted.  Figure 8: (a) Static spin
structure factors for the clean case (upper panel) and the disordered case
(lower panel) for various values of $J_{2,3}$, where $K=-6$ and $J=1$ are
fixed. (b) Possible ground-state phase diagrams for the clean and disordered
systems as a function of $J_{2,3}$.

So far we have considered the effect of disorder on quantum phase transitions
between Kitaev QSL and magnetic ordered states in the Kitaev-Heisenberg model
[Eq. (1)], where the ordered states are caused by the competition of nearest-
neighbor exchange coupling $J$ and Kitaev term $K$. However, in real materials
like Na2IrO3 and $\alpha$-RuCl3, the experimentally observed zigzag ordering
may be possibly stabilised by long-range neighbor couplings $J_{2}$ and
$J_{3}$ [39]. Therefore, it would be informative to study the effect of
disorder on zigzag state stabilised by $J_{2}$ and $J_{3}$. Thus, the
Hamiltonian is modified as

$\mathcal{H}=\mathcal{H}_{\rm
KH}+\sum_{(i,j)}\sum_{\alpha=2,3}J_{\alpha}\mathbf{S}_{i}\cdot\mathbf{S}_{j}$
(4)

where $J_{2}$ and $J_{3}$ are uniform Heisenberg exchange interactions on 2nd
and 3rd nearest neighbours, respectively. We fix $K=-6$, $J=1$ as a typical
parameter set for the Kitaev materials. We omit the $\Gamma$ terms from our
model [40] and set $J_{2}=J_{3}(\equiv J_{2,3})$ for simplicity. We now focus
on the case of $n=1/3$, i.e., 12 bonds out of original 36 bonds are randomly
replaced by singular couplings with disorder strength $r=4$. This value of $r$
was actually predicted as a singular-coupling disorder due to the missing of
hydrogen cations in H3LiIr2O6 [33].

Let us first see the second derivative of ground-state energy as a function of
$J_{2,3}$, i.e., $-\partial^{2}E_{0}/\partial J_{2,3}^{2}$. The results for
the clean and disordered cases are compared in Fig. 7(a). In the clean case,
two phase boundaries are signalled by the sharp peaks. We thus find three
phases depending on the value of $J_{2,3}$: FM ($J_{2,3}\lesssim-0.34$),
Kitaev-type QSL ($-0.34\lesssim J_{2,3}\lesssim 0.29$), and zigzag
($J_{2,3}\gtrsim 0.29$). At $J_{2,3}=0$, although the ground state is a QSL,
the system is in the vicinity of the stripy phase (the phase boundary is
$K=-4$, $J=1$), For the disordered case, the averaged values over 30 disorder
samples are plotted. Compared to the clean case, the heights of both peaks are
much reduced and the widths are significantly broadened. Furthermore, it is
interesting that the phase between FM and zigzag, seems to be extended by the
disorder since the left (right) peak position is shifted from
$J_{2,3}\approx-0.34$ ($J_{2,3}\approx 0.29$) to $J_{2,3}\approx-0.5$
($J_{2,3}\approx 0.5$).

In order to take a more detailed look, let us see the result for each disorder
sample. We find that most of disorder samples exhibit three or four peaks,
i.e., the appearance of additional one or two peaks to the clean case, in
$-\partial^{2}E_{0}/\partial J_{2,3}^{2}$. Typical examples of
$-\partial^{2}E_{0}/\partial J_{2,3}^{2}$ for two disorder samples are plotted
in Fig. 7(b). For both samples there exist five ‘regions’ separated by four
peaks. Fig. 7(c) shows the static spin structure factor at representative
$J_{2,3}$ values for each region: $J_{2,3}=-0.75$, $-0.5$, $0$, $0.5$, and
$0.75$. Incidentally, in the clean case, the system is in a FM phase for
$J_{2,3}=-0.75$ and $-0.5$; in a QSL phase at $J_{2,3}=0$; and in a zigzag
phase at $J_{2,3}=0.5$ and $0.75$. The two sets of $S({\bf Q})$ for disordered
samples show similar qualitative features at the same $J_{2,3}$ value. At
$J_{2,3}=-0.75$ we find a sharp peak at ${\bf Q}={\bf 0}$ indicating FM state,
as in the clean case. At $J_{2,3}=-0.5$ a stripy structure appears in addition
to the FM peak. Therefore, the state may be characterized as a quasi-FM order
with short-range stripy fluctuations. Note that the position of stripy Bragg
peaks reflects one of three possible ordering directions in each disorder
sample. At $J_{2,3}=0$ it is interesting that a QSL-like structure, i.e.,
structureless distribution of the intensity over the first Brillouin zone, and
a stripy Bragg peaks coexist. Since the system is in the proximity to stripy
phase in the clean case, a short-range stripy correlations could be enhanced
by a pinning effect with disorder. Nevertheless, as seen in Fig. 7(d) the
value of $\langle W_{\rm p}\rangle$ is still significantly larger than zero
for most of the disorder samples, so that the system is broadly characterized
as a Kitaev QSL. At $J_{2,3}=0.5$ the intensity is widely distributed and
there is no sharp Bragg peaks indicating a specific ordering. Since the value
of $\langle W\rangle$ is nearly zero for most of the disorder samples, the
system is in a non-Kitaev QSL state. In this QSL phase, unlike a very rapid
decay of the spin-spin correlation functions in the Kitaev QSL phase, the
correlation survives at longer distance (see the Appendix VI.3 for more
details).

At $J_{2,3}=0.75$ the structure factor for each disorder sample exhibits a
symmetry-broken zigzag state with one of three possible ordering directions.
Thus, the system consists of three zigzag domains. As discussed in Sec. III.2
the size of domains depends on the disorder density and strength. Also, it is
worth mentioning that the QSL range is somewhat expanded by disorder. It may
be because the disorder-induced stripy fluctuations push out the phase
boundaries of neighboring ordered phases.

In Fig. 8(a) the evolution of the static structure factors with $J_{2,3}$ is
compared for the clean and disordered case. For the disordered case, we show
the structure factor averaged over ones over 30 disorder samples. Although the
overall explanation has been already given above for two disorder samples,the
effect of disorder can be more obviously seen in the averaged $S({\bf Q})$. Of
particular interest is that the sharp Bragg peak features in the clean system
are totally collapsed in the disordered system at $J_{2,3}=-0.4$, $0.3$, and
$0.5$. This means that even short-range magnetic orderings are destroyed by
disorder. Accordingly, the $J_{2,3}$-range of QSL phase is extended from
$-0.34\lesssim J_{2,3}\lesssim 0.29$ to $-0.45\lesssim J_{2,3}\lesssim 0.55$
by the disorder as shown in prospective ground-state phase diagrams [Fig.
8(b)]. This is in a strong contrast to the results for the case of NN Kitaev-
Heisenberg model, where the QSL phase is always shrank by any types of
disorder. As mentioned above, this value of $r=4$ corresponds to the singular-
coupling disorder in H3LiIr2O6 [33], however, the actual values of $J_{2}$,
$J_{3}$, and $n$ are unknown. Nevertheless, since the qualitative features
would be unaffected by changing $n$, assuming typical values of $J_{2}$ and
$J_{3}$ as several tens of percent of $J$ we can speculate that either spin
glass or QSL is likely as low-temperature state of H3LiIr2O6. Moreover, we
find random-singlet formations between second- and third-neighbor spins in our
system. This is qualitatively consistent with heat capacity measurement under
magnetic field [32]. For more quantitative considerations the $\Gamma$ terms
should be also taken into account because they are expected to be sizeable

##  IV Conclusion

Using Lanczos exact diagonalization technique we have studied the Kitaev-
Heisenberg model on a honeycomb lattice in the presence of bond disorder.
Especially, the effects of disorder on the transition between Kitaev spin
liquid and magnetic ordered states as well as the stability of magnetic
ordering have been investigated. We have considered two types of disorder:
random-coupling disorder and singular-coupling disorder. The former is
achieved by a random distribution of bond strength, and the latter by a random
mixture of two kinds of bond strengths. They exhibit qualitatively similar
effects in the pure Kitaev-Heisenberg model without long-range interactions.
The range of Kitaev spin liquid phase is reduced and the transition to
magnetic ordered phases becomes more crossover-like. Accordingly, although the
regions of zigzag and stripy phases are enlarged by disorder, their long-range
orderings in the clean system are replaced by a mixture of the three domains
with different ordering directions. With increasing the disorder strength the
area of domains becomes smaller and the system goes into a spin-glass state.
Particularly in the crossover regime, the coexistence of magnetic domains and
Kitaev spin-liquid domains may be possible. On the other hand, qualitative
trend is different if magnetic ordered phases caused by long-range
interactions are doped by disorder. The stability of such magnetic ordering is
diminished by singular-coupling disorder, and as a result, the range of spin-
liquid is extended. This mechanism may be relevant to materials like
$\alpha$-RuCl3 and H3LiIr2O6 where the zigzag ground state is stabilized by
small long-range interactions.

We have also found that the flux states are induced by a certain amount of
disorder strength. Interestingly, the fluxes tend to form kinds of
commensurate ordering at intermediate disorder strength. This may be related
to the Majorana insulating” state. This may indicate an existence of repulsive
interaction between fluxes. With further increasing disorder strength the
system goes into the random-flux states. At any disorder strength we have seen
no signature for a phase-separation-like flux configuration.

##  V Acknowledgments

We acknowledge Joji Nasu, Rajyavardhan Ray, and Pranay Patil for fruitful
discussions. We thank Ulrike Nitzsche for technical assistance. This project
is funded by IFW Excellence Programme 2020 and the German Research Foundation
(DFG) via the projects A05 of the Collaborative Research Center SFB 1143
(project-id 247310070) and through the Würzburg-Dresden Cluster of Excellence
on Complexity and Topology in Quantum Matter ct.qmat (EXC 2147, Project No.
39085490).

##  VI Appendix

###  VI.1 Excitation of free-flux state by applying local flux field

Figure 9: Critical flux field for various clusters. The used clusters are
shown in Fig. 10(a).

It is known that the ground state of the Kitaev model on a honeycomb lattice
is characterized by $\langle W_{p}\rangle=1$ ($\forall p$), the so-called
flux-free state. The flux operator is defined as
$W_{p}=2^{6}S_{1}^{x}S_{2}^{y}S_{3}^{z}S_{4}^{x}S_{5}^{y}S_{6}^{z}$ (also see
the main text). An excited flux state can be created by flipping $\langle
W_{p}\rangle$ from $1$ to $-1$ for even number of plaquettes. We here consider
the Kitaev model with $K=-2$. To study flux configuration in the excited flux
state, we apply flux field $h_{\rm f}W_{p}$ to a single plaquette on finite-
size clusters shown in Fig. 10(a). This is a kind of the pinning-field method
which is usually used to study magnetically ordered state. With increasing
$h_{\rm f}$, the value of $\langle W_{p}\rangle$ for the field-applied
plaquette is flipped from $1$ to $-1$ at some $h_{\rm f}$. We call it critical
flux field. In Fig. 9 the critical flux field is plotted as a function of
$1/N$. Above the critical flux field, the clusters 24a and 8b contain four
fluxes, and the other clusters contain only two fluxes. For the clusters with
two fluxes, the position of fluxes may be not uniquely fixed due to the
symmetry of cluster. For example, let us see the cluster 24b in Fig. 10(a). A
red plaquette corresponds to $\langle W_{p}\rangle=-1$ where the flux field is
applied. The remaining $\langle W_{p}\rangle=-1$ is put into either of two
thin red plaquettes.

Interestingly, the flux configuration for 24a and 18a exhibits a staggered
order. (Note that the cluster denoted by 18a is not perfectly consistent with
a staggered order because three fluxes are needed to complete it.
Nevertheless, the possible positions of fluxes are of the staggered order [see
Fig. 10(a)].) As seen in Fig. 9, the critical flux fields for 24a and 18a are
lower than those for the other clusters with the same system size. Especially,
the critical flux field for 24a is appreciably lower than that for24b and 24c.
This may mean that the staggered ordered state is more stable than the other
flux state. This is consistent with the fact that the staggered ordered state
has significantly lower energy than the random flux state at the flux density
$n_{\rm F}=1/3$ [36]. The clusters having relatively large critical flux
field, i.e., 12b, 20b, and 24c, might be related to a random flux state.

Furthermore, some clusters having non-flux-free ground state are shown in Fig.
10(b) as additional information. For the clusters 12b, 16b, 16c, and 20c, the
flux-free state and full-flux state are (nearly) degenerate in the ground
state. For the cluster 14b, the ground state contains five $\langle
W_{p}\rangle=1$ and two $\langle W_{p}\rangle=-1$. For the cluster 8a, the
ground state is full-flux state. These results might be useful to consider the
finite-size effects in the Majorana basis for future studies.

Figure 10: (a) Clusters used for estimating the critical flux field. The red
plaquette corresponds to $\langle W_{p}\rangle=-1$, and either of the thin red
plaquettes is occupied by $\langle W_{p}\rangle=-1$ in the excited flux state
for each cluster with two fluxes. (b) Clusters having non-flux-free ground
state (see text).

###  VI.2 Repulsive interaction between fluxes

Table 1: Critical flux field $h_{\rm f,c}$(NN) and $h_{\rm f,c}$(NNN) for the clusters denoted by 24b, 24c, and 18a. The effective repulsive interaction is estimated as $V_{\rm eff}=h_{\rm f,c}$(NN)-$h_{\rm f,c}$(NNN).  cluster |  $h_{\rm f,c}$(NN) |  $h_{\rm f,c}$(NNN) |  $V_{\rm eff}$
---|---|---|---
24b |  0.0404969256 |  0.0337186000 |  0.006778326
24c |  0.0471150167 |  0.0401504917 |  0.006964525
18a |  0.0338108475 |  0.0271080335 |  0.006702814

If a flux-ordered state such as staggered order is more stable than random
flux state, there must be a kind of repulsive interaction between fluxes. We
now estimate the effective value of the repulsive interaction by applying flux
field $h_{\rm f}W_{p}$ to two plaquettes on the clusters exhibiting excited
flux states with two nonneighboring fluxes. Such clusters are 24b, 24c, and
18a in Fig. 10(a). First, we obtain two critical flux fields: (i) $h_{\rm
f,c}$(NN) with flux field applied to two neighboring plaquettes and (ii)
$h_{\rm f,c}$(NNN) with flux field applied to two next-nearest-neighbor
plaquettes. Since there is a repulsion between fluxes, $h_{\rm f}$(NN) is
larger than $h_{\rm f}$(NNN). The obtained values are shown in Table 1. Thus,
we may estimate the effective value of the repulsive interaction $V_{\rm eff}$
from the difference between $h_{\rm f,c}$(NN) and $h_{\rm f,c}$(NNN). The
obtained value of $V_{\rm eff}$ are very close for three clusters, i.e.,
$V_{\rm eff}\approx 0.007$ in unit of $|2K|$, although a finite-size scaling
analysis is required to further confirm this value.

###  VI.3 spin-spin correlations in the non-Kitaev QSL phase

As shown in FIG. 8, the region of $-0.45\lesssim J_{2,3}\lesssim 0.55$ for the
disordered system is characterized as a QSL based on the structureless static
structure factor. Part of the QSL phase, i.e., at $0.35\lesssim
J_{2,3}\lesssim 0.55$, is recognized as non-Kitaev one due to the reduced
$\langle W_{p}\rangle$ value. Here, we look at the spin-spin correlation
functions $\langle{\bf S}_{i}\cdot{\bf S}_{j}\rangle$ to further support this
discrimination of Kitaev and non-Kitaev QSLs. In a pure Kitaev model, only the
NN correlations are finite and longer-range ones are zero; that is faithfully
reproduced by the 24-site calculations. Even away from the pure Kitaev limit,
the spin-spin correlation decays rapidly with distance if the system is in the
Kitaev QSL state.

In Fig. 11 the absolute values of spin-spin correlation functions for nearest-
neighbor, 2nd-neighbor, 3rd-neighbor, and 4th-neighbor bonds are plotted as a
function of $J_{2,3}$ for typical three disorder samples. Since the system is
disordered, the correlations are averaged over all possible bonds. In the
Kitaev QSL region, we can see a rapid decay of $\overline{|\langle{\bf
S}_{i}\cdot{\bf S}_{j}\rangle|}$. Whereas in the non-Kitaev QSL region, we
find an enhancement of the 4th-neighbor correlation, and it becomes larger
than that for the 3rd-neighbor bond. This is clearly in contrast to the Kitaev
QSL feature. Interestingly, this reverse of 3rd-neighbor and 4th-neighbor
correlations in the non-Kitaev QSL region is confirmed in all of our 30
disorder samples.

Figure 11:  (a,c,e) Expectation value of flux operator and (b,d,f) spin-spin
correlation functions for nearest-neighbor, 2nd-neighbor, 3rd-neighbor, and
4th-neighbor bonds as a function of $J_{2,3}$ for three disorder samples. The
inset in (b) shows the spin pairs considered and the yellow highlighted area
shows the non-Kitaev spin liquid.

## References

  * Balents [2010] L. Balents, Spin liquids in frustrated magnets, Nature 464, 199 (2010).
  * Takagi  _et al._ [2019] H. Takagi, T. Takayama, G. Jackeli, G. Khaliullin, and S. E. Nagler, Concept and realization of kitaev quantum spin liquids, Nature Reviews Physics 1, 264 (2019).
  * Kitaev [2006] A. Kitaev, Anyons in an exactly solved model and beyond, Annals of Physics 321, 2 (2006).
  * Jackeli and Khaliullin [2009] G. Jackeli and G. Khaliullin, Mott insulators in the strong spin-orbit coupling limit: from heisenberg to a quantum compass and kitaev models, Physical review letters 102, 017205 (2009).
  * Singh  _et al._ [2012] Y. Singh, S. Manni, J. Reuther, T. Berlijn, R. Thomale, W. Ku, S. Trebst, and P. Gegenwart, Relevance of the heisenberg-kitaev model for the honeycomb lattice iridates a 2 iro 3, Physical review letters 108, 127203 (2012).
  * Chaloupka  _et al._ [2010] J. Chaloupka, G. Jackeli, and G. Khaliullin, Kitaev-heisenberg model on a honeycomb lattice: possible exotic phases in iridium oxides a 2 iro 3, Physical review letters 105, 027204 (2010).
  * Koitzsch  _et al._ [2016] A. Koitzsch, C. Habenicht, E. Müller, M. Knupfer, B. Büchner, H. C. Kandpal, J. Van Den Brink, D. Nowak, A. Isaeva, and T. Doert, J eff description of the honeycomb mott insulator $\alpha$\- rucl 3, Physical Review Letters 117, 126403 (2016).
  * Winter  _et al._ [2017] S. M. Winter, A. A. Tsirlin, M. Daghofer, J. van den Brink, Y. Singh, P. Gegenwart, and R. Valentí, Models and materials for generalized kitaev magnetism, Journal of Physics: Condensed Matter 29, 493002 (2017).
  * Trebst [2017] S. Trebst, Kitaev materials, arXiv preprint arXiv:1701.07056 (2017).
  * Choi  _et al._ [2012] S. Choi, R. Coldea, A. Kolmogorov, T. Lancaster, I. Mazin, S. Blundell, P. Radaelli, Y. Singh, P. Gegenwart, K. Choi, _et al._ , Spin waves and revised crystal structure of honeycomb iridate na 2 iro 3, Physical review letters 108, 127204 (2012).
  * Williams  _et al._ [2016] S. Williams, R. Johnson, F. Freund, S. Choi, A. Jesche, I. Kimchi, S. Manni, A. Bombardi, P. Manuel, P. Gegenwart, _et al._ , Incommensurate counterrotating magnetic order stabilized by kitaev interactions in the layered honeycomb $\alpha$\- li 2 iro 3, Physical Review B 93, 195158 (2016).
  * Rau  _et al._ [2014] J. G. Rau, E. K.-H. Lee, and H.-Y. Kee, Generic spin model for the honeycomb iridates beyond the kitaev limit, Physical review letters 112, 077204 (2014).
  * Hwan Chun  _et al._ [2015] S. Hwan Chun, J.-W. Kim, J. Kim, H. Zheng, C. C. Stoumpos, C. Malliakas, J. Mitchell, K. Mehlawat, Y. Singh, Y. Choi, _et al._ , Direct evidence for dominant bond-directional interactions in a honeycomb lattice iridate na2iro3, Nature Physics 11, 462 (2015).
  * Sears  _et al._ [2015] J. A. Sears, M. Songvilay, K. Plumb, J. Clancy, Y. Qiu, Y. Zhao, D. Parshall, and Y.-J. Kim, Magnetic order in $\alpha$\- rucl 3: A honeycomb-lattice quantum magnet with strong spin-orbit coupling, Physical Review B 91, 144420 (2015).
  * Freitas and Pereira [2022] L. R. Freitas and R. G. Pereira, Gapless excitations in non-abelian kitaev spin liquids with line defects, Physical Review B 105, L041104 (2022).
  * Andrade  _et al._ [2020] E. C. Andrade, L. Janssen, and M. Vojta, Susceptibility anisotropy and its disorder evolution in models for kitaev materials, Physical Review B 102, 115160 (2020).
  * Kao and Perkins [2021] W.-H. Kao and N. B. Perkins, Disorder upon disorder: Localization effects in the kitaev spin liquid, Annals of Physics 435, 168506 (2021).
  * Bahrami  _et al._ [2021] F. Bahrami, E. M. Kenney, C. Wang, A. Berlie, O. I. Lebedev, M. J. Graf, and F. Tafti, Effect of structural disorder on the kitaev magnet ag 3 liir 2 o 6, Physical Review B 103, 094427 (2021).
  * Bahrami  _et al._ [2022] F. Bahrami, M. Abramchuk, O. Lebedev, and F. Tafti, Metastable kitaev magnets, Molecules 27, 871 (2022).
  * Kitagawa  _et al._ [2018] K. Kitagawa, T. Takayama, Y. Matsumoto, A. Kato, R. Takano, Y. Kishimoto, S. Bette, R. Dinnebier, G. Jackeli, and H. Takagi, A spin–orbital-entangled quantum liquid on a honeycomb lattice, Nature 554, 341 (2018).
  * Knolle  _et al._ [2019] J. Knolle, R. Moessner, and N. B. Perkins, Bond-disordered spin liquid and the honeycomb iridate h 3 liir 2 o 6: Abundant low-energy density of states from random majorana hopping, Physical review letters 122, 047202 (2019).
  * Zschocke and Vojta [2015] F. Zschocke and M. Vojta, Physical states and finite-size effects in kitaev’s honeycomb model: Bond disorder, spin excitations, and nmr line shape, Physical Review B 92, 014403 (2015).
  * Nasu and Motome [2021] J. Nasu and Y. Motome, Spin dynamics in the kitaev model with disorder: Quantum monte carlo study of dynamical spin structure factor, magnetic susceptibility, and nmr relaxation rate, Physical Review B 104, 035116 (2021).
  * Nasu and Motome [2020] J. Nasu and Y. Motome, Thermodynamic and transport properties in disordered kitaev models, Physical Review B 102, 054437 (2020).
  * Willans  _et al._ [2010] A. Willans, J. Chalker, and R. Moessner, Disorder in a quantum spin liquid: Flux binding and local moment formation, Physical review letters 104, 237203 (2010).
  * Willans  _et al._ [2011] A. Willans, J. Chalker, and R. Moessner, Site dilution in the kitaev honeycomb model, Physical Review B 84, 115146 (2011).
  * Kao  _et al._ [2021] W.-H. Kao, J. Knolle, G. B. Halász, R. Moessner, and N. B. Perkins, Vacancy-induced low-energy density of states in the kitaev spin liquid, Physical Review X 11, 011034 (2021).
  * Yamada [2020] M. G. Yamada, Anderson–kitaev spin liquid, npj Quantum Materials 5, 1 (2020).
  * Dantas and Andrade [2022] V. Dantas and E. C. Andrade, Disorder, low-energy excitations, and topology in the kitaev spin liquid, Physical Review Letters 129, 037204 (2022).
  * Do  _et al._ [2020] S.-H. Do, C. Lee, T. Kihara, Y. Choi, S. Yoon, K. Kim, H. Cheong, W.-T. Chen, F. Chou, H. Nojiri, _et al._ , Randomly hopping majorana fermions in the diluted kitaev system $\alpha$-ru 0.8 ir 0.2 cl 3, Physical review letters 124, 047204 (2020).
  * Chaloupka  _et al._ [2013] J. Chaloupka, G. Jackeli, and G. Khaliullin, Zigzag magnetic order in the iridium oxide na 2 iro 3, Physical review letters 110, 097204 (2013).
  * Kimchi  _et al._ [2018] I. Kimchi, J. P. Sheckelton, T. M. McQueen, and P. A. Lee, Scaling and data collapse from local moments in frustrated disordered quantum spin systems, Nature communications 9, 1 (2018).
  * Yadav  _et al._ [2018] R. Yadav, R. Ray, M. S. Eldeeb, S. Nishimoto, L. Hozoi, and J. van den Brink, Strong effect of hydrogen order on magnetic kitaev interactions in ${\mathrm{h}}_{3}{\mathrm{liir}}_{2}{\mathrm{o}}_{6}$, Phys. Rev. Lett. 121, 197203 (2018).
  * Laflorencie  _et al._ [2006] N. Laflorencie, S. Wessel, A. Läuchli, and H. Rieger, Random-exchange quantum heisenberg antiferromagnets on a square lattice, Phys. Rev. B 73, 060403 (2006).
  * Andrade and Vojta [2014] E. C. Andrade and M. Vojta, Magnetism in spin models for depleted honeycomb-lattice iridates: Spin-glass order towards percolation, Physical Review B 90, 205112 (2014).
  * Koga  _et al._ [2021] A. Koga, Y. Murakami, and J. Nasu, Majorana correlations in the kitaev model with ordered-flux structures, Phys. Rev. B 103, 214421 (2021).
  * Kao and Perkins [2022] W.-H. Kao and N. B. Perkins, Spin- and flux-gap renormalization in the random kitaev spin ladder, Phys. Rev. B 106, L100402 (2022).
  * Song  _et al._ [2016] X.-Y. Song, Y.-Z. You, and L. Balents, Low-energy spin dynamics of the honeycomb spin liquid beyond the kitaev limit, Phys. Rev. Lett. 117, 037209 (2016).
  * Yadav  _et al._ [2016] R. Yadav, N. A. Bogdanov, V. M. Katukuri, S. Nishimoto, J. Van Den Brink, and L. Hozoi, Kitaev exchange and field-induced quantum spin-liquid states in honeycomb $\alpha$-rucl3, Scientific reports 6, 1 (2016).
  * [40] For example, the Kitaev and Heisenberg interactions were estimated as $K=-6.15$ meV and $J=1.65$ meV ($|K/J|\sim 3.7$) for H3LiIr2O6 [33]. They lead to a stripy phase if the other interactions are ignored. However, the system is shifted to a QSL state by considering the $\Gamma$ terms. Thus, the assumption of a QSL state in the vicinity of stripy phase would be a good starting point at $J_{2}=J_{3}=0$.


# Finite amplitude wave propagation through bubbly fluids

Siew-Wan Ohl  Department of Soft Matter, Institute of Physics & Faculty of
Natural Sciences, University of Magdeburg, Universitätplatz 2, 39108
Magdeburg, Germany.  Research Campus STIMULATE, University of Magdeburg, Otto-
Hahn-Straße 2, 39106 Magdeburg, Germany.  Juan Manuel Rosselló  Department of
Soft Matter, Institute of Physics & Faculty of Natural Sciences, University of
Magdeburg, Universitätplatz 2, 39108 Magdeburg, Germany.  Daniel Fuster
Institut Jean Le Rond D’Alembert UPMC, Tour 55-65, 75005 Paris, France.
Claus-Dieter Ohl  Department of Soft Matter, Institute of Physics & Faculty of
Natural Sciences, University of Magdeburg, Universitätplatz 2, 39108
Magdeburg, Germany.

###### Abstract

A few bubbles in a liquid could drastically reduce the acoustic wave speed.
Wood’s equation models the linear sound speed, while the speed of an ideal
shock waves is derived as a function of the pressure ratio across the shock.
The common finite amplitude waves lie, however, in between these limits. We
show that in a bubbly medium, the high frequency components of finite
amplitude waves dissipate quickly but a low-frequency part remains. We
demonstrate that the  the transition between linear and the shock wave regimes
can be smoothly connected through a Mach number based on the collapse velocity
of the bubbles.

Since the observation of damping of sound in frothy liquid in 1910 [1] the
physics of sound propagation in bubbly liquid has attracted many research
works. The simplest models have the assumptions that the bubbly liquid is a
continuum, the sound wave has a small amplitude, the bubbles are uniformly
distributed [2, 3, 4, 5]. Wijngaarden [5] could show that a void fraction of
0.3% reduces the speed of sound to about 200 m/s from 1500 m/s in the pure
liquid.

Wave velocity and attenuation depend on frequency in a complex way. Effective
medium approaches estimate the dispersion relation of pure gas bubbles [3, 6,
7] and gas/vapor bubbles [8] in the linear regime. At low frequencies,
dispersion effects are negligible, and the effective speed approaches an
asymptotic value known as Wood’s limit [9]:

$\frac{1}{c^{2}}=\frac{(1-\alpha)^{2}}{c_{l}^{2}}+\frac{\alpha^{2}}{c_{g}^{2}}+\frac{\rho_{l}\alpha(1-\alpha)}{\kappa\,p}\quad,$
(1)

where $c_{l}$ and $c_{g}$ are the speed of sound in the liquid and gas
respectively, $\alpha$ is the void fraction, and $\kappa$ the ratio of
specific heats [10]. Silberman [11], among others, have measured the sound
velocity in a standing wave tube filled with a mixture of bubbles and water,
and verified Wood’s derivation.

For the strongly nonlinear regime where the small amplitude assumption is not
valid, Campbell and Pitcher [12] connected Eq. (1) with the propagation speed
of a shock wave, $U$. They used the shock relations for continuity, momentum,
and energy in a bubbly liquid and obtained

$U^{2}=c^{2}\left(\frac{p_{1}}{p_{0}}-1\right)\left(1-\left(\frac{p_{0}}{p_{1}}\right)^{1/\kappa}\right)^{-1}\quad,$
(2)

with $p_{0}$ and $p_{1}$ being the pressures in the bubbly liquid ahead and
behind the shock wave, respectively [10]. This expression for $\kappa=1$, i.e.
$U^{2}=(p_{1}/p_{0})c^{2}$, was confirmed by Campbell and Pitcher [12] with a
vertical shock tube experiment where small gas bubbles and void fractions
between $\alpha=5\%$ and $\alpha=30\%$ were used.  It is noted that Eq. (2) is
only valid for sufficiently high void fractions and thus $c\ll c_{l}$, for
which the effective sound speed of Eq. (1) no longer depends on $c_{l}$

$c\approx\sqrt{\frac{\kappa p_{0}}{\rho_{l}(1-\alpha)}}$

Interestingly, this theory predicts that the effective speed of sound of the
medium is proportional to $\sqrt{p_{1}/p_{0}}$.

Between the linear and the strongly nonlinear regimes, models for dispersion
and relaxation of weak and moderate-strength shock waves have been established
and tested experimentally. Kameda et al. [13, 14] confirmed Wijngaarden’s [5]
bubble-liquid mixture model with $\Delta p=p_{1}-p_{0}$ up to $1\,$bar. Ando
_et al._ [15] modeled bubble poly-dispersity’s effect on shock front shape,
and found that the broad bubble size distributions smooth out the profile.
Therefore, the linear theory is found not to be not suitable for predicting
strong shock propagation speed.

Here, we study how shock waves propagate through bubbles in gelatin. Our
research is inspired by Dear and Field’s classic work on shock waves
interaction with gas bubbles in hydrogels [16, 17]. Firstly, we use Direct
Numerical Simulation to understand wave transmission across all acoustic
regimes. After that, the experiments on shock wave propagation through bubbles
in gelatin gels are presented. We found that the experimental results are in
agreement with DNS simulations on high-pressure finite amplitude wave regimes.

Direct Numerical Simulation — The simulations are conducted with a multiphase
compressible solver presented in Fuster and Popinet [18]. We model an
idealized system where a pressure pulse propagates through a mono-dispersed
line of bubbles with size $R_{0}=100\,\mu$m in a viscous fluid (Figure 1(a)).
We consider a periodic axisymmetric configuration where the bubbles are
uniformly distributed along the axis of symmetry at a constant inter-bubble
distance $d$. On the top domain boundary, slip boundary conditions are
applied. The simulation domain extends from $z\in[-d/2:52.5d]$ and from
$r\in[0:d]$ in the radial direction. With this configuration, the effective
void fraction is, therefore, $\alpha=4/3(R_{0}/d)^{3}$.

Figure 1: (a) Domain and boundary conditions used for DNS simulations. (b) A
series of snapshots of the pressure field and interface contour at an initial
void fraction, $\beta=9.2\times 10^{-4}$ and $\Delta p/p_{0}=100$ with
indicated timing. (c) The effect of $\Delta p/p_{0}$ at constant void
fraction. Red line shows that for low Ma, the velocities are predicted by the
Wood’s theory. But for Ma ¿ 0.1, the effective velocity is proportional to Ma
at constant void fraction. Blue dots show simulation results for high
amplitude driving of $p_{1}/p_{0}=100$. (d) Temporal evolution of the non-
dimensionalized peak pressure in the domain
$p_{max}(t)^{\prime}=\frac{p_{max}(t)-p_{0}}{\Delta p}$(solid lines) and
evolution of a measure of the bubble collapse strength
$V^{\prime}=~{}\frac{V_{0}-V_{min}(t)}{V_{0}-min(V_{min}(t))}$ (dashed lines)
for a linear case (red line) and a strongly nonlinear case (blue line).
$V_{min}(t)$ represents the minimum bubble volume among all bubbles at a given
instant.

A Gaussian pulse, denoted by equations $p(z,t=0)=p_{0}+\Delta
pe^{-z^{2}/\lambda^{2}}$ and $u_{z}(z,t=0)=\frac{p(z,0)-p_{0}}{\rho_{l}c_{l}}$
is propagating in the simulation domain from left to right. The pulse shape is
designed to mimic the shock wave profile in the experiment. The initial
amplitude of the pressure pulse is left as a variable parameter, and the
parameter $\lambda$ is set to $\lambda=500~{}\mu$m  imposing a range of
forcing frequencies below $\frac{\omega_{pulse}}{\omega_{0}}\leq\frac{2\pi
c_{0}}{\sqrt{3\kappa p_{0}/\rho_{l}}}\frac{R_{0}}{\lambda}\approx 100$. The
gas inside the bubble is modeled as an ideal gas. The viscous fluid is assumed
to be water-like but with a higher viscosity ($\nu_{l}=2.3\times
10^{-4}m^{2}/s$)(see supplementary material for value derivation). This solver
for the Navier–Stokes equations in a Newtonian fluid has previously been
applied for sub-harmonic emissions of bubbles oscillating in a tube [19].

Figure 1(b) depicts various snapshots of the pressure field and the bubble
interface contour during the propagation of a wave with amplitude $\Delta
p/p_{0}=100$. The precursor wave, which is the high-frequency component,
propagates across the bubbly medium at approximately the speed of sound in the
pure liquid, inducing an initial compression of the bubbles. The term “high
frequency” here is referring to the signal frequencies which are higher than
the bubble resonance frequency, $\omega\gg\omega_{0}$. There is an increasing
delay between the instant at which the precursor wave reaches a given bubble
location and the instant at which the bubble volume becomes minimum.

Figure 1(c) shows the effect of $\Delta p/p_{0}$ at constant void fraction.
The velocities obtained for low amplitudes of the forcing pulse (low Ma) are
independent of the Mach number and are predicted by the Wood theory. For
$Ma>0.1$, the effective velocity is then proportional to the Mach number for a
constant void fraction. For reference, simulations performed at constant
pressure $p_{1}/p_{0}=100$ and different concentrations of bubbles are also
included. As expected, the resulting velocities do not scale with the Wood’s
velocity, obtaining a dependence of $U/c_{wood}$ on the value of the pressure
amplitude used.

Figure 1(d) reveals that the high-frequency content of the signal propagates
at liquid velocity. As the attenuation and dispersion effects become dominant,
the wave quickly attenuates. The bubbles act as an acoustic low-pass filter,
such after 0.01 ms only low-frequency components of the signal remains. In the
nonlinear regime, the low-frequency component remains approximately constant
for all bubbles (dashed blue line in Figure 1(d)).

Experiments — A gelatin gel is used to trap the bubbles that have been
generated by shock waves from a piezoelectric-driven lithotripter. We vary
consistently the void fraction in the gelatin gel by allowing the bubbles to
grow in size over time. The pressure is measured with a needle hydrophone
(diameter 1.1 mm, model Müller-Platte Nadelsonde) inserted into the gelatin
cuvette about $2\,$mm sideways from the focus. The first bubbles are nucleated
in the tensile region with approximately $-10\,$bar amplitude. The bubble
dynamics is captured with a high-speed camera (HPV2, Shimadzu Inc.) together
with a generic camera flashlight at a framing rate of 1 million frames per
second. The timing of the devices is controlled with a digital delay generator
(BNC 525, Berkeley Nucleonics).

Frames from the high-speed recording are analyzed with home-built algorithms
[20] to detect the bubbles, register their locations and shapes, and link them
across frames to follow their dynamics. The shock wave speed was then obtained
by performing a numerical fit on the position where the bubbles reach a
minimum size over time using a linear model. Initially, there is no bubble in
the gelatin gel  (Figure 2(i)), and thus to get the shock speed for the first
shock, we track the bubble nucleation as the tensile phase of the shock wave
passes through. Details about the image processing, as well as gelatin
preparation and elasticity measurement, are provided in the supplementary
material.

Figure 2:  Bubble nucleation and compression by a traveling shock in 4%
gelatin. The frame height is $4.974\,$mm and the time is indicated in
microseconds. (i)$1^{\mathrm{st}}$ shock wave, (ii) $2^{\mathrm{nd}}$ shock
wave, and (iii) $8^{\mathrm{th}}$ shock wave.

Figure 2(i) depicts a high-speed recording of the first pressure pulse
propagating in a virgin gelatin block (4% by weight) from the bottom to the
top in positive $z$-direction. The time between frames is $1\,\mu$s. Within
the first 4 frames, the bubbles are nucleated over a distance of about
$5\,$mm. The approximate location where the rarefaction wave starts is
indicated with a solid white line. In the later frames, the bubbles grow to a
larger size and collapse approximately $150\,\mu$s later (not shown here). In
contrast to the bubbles observed in a liquid, these gas bubbles in gelatin do
not quickly dissolve or rise out of the field of view. Once formed, they
remain at their location. This is utilized by allowing the bubbles to grow in
between shocks for about 10 minutes through the  diffusion of gasses from the
gelatin [21], thus to increase the void fraction. This void fraction growth is
repeated 10 times. After about 100 minutes, the experiment for one gelatin
sample is concluded.

Figure 2(ii) shows the shock wave that passes through the gel that contains
gas bubbles from the first shock. Over time, we see an upwards propagating
region of compressed bubbles, i.e. over $3\,\mu$s in the frames from $1\,\mu$s
$\leq t\leq 4\,\mu$s. The approximate location of the minimum bubble volume is
indicated with the dash horizontal line.  Figure 2(iii) shows the bubble
dynamics induced by the 8th shock wave. We find a considerable delay in the
arrival of the shock wave, which enters the field of view from below only at
$t=5\,\mu$s. This later arrival time is confirmed by the considerably slower
propagation of the shock front. Now the compressed bubbles are seen at the top
of the frame in the right-most frame, i.e. $t=14\,\mu$s. The average speed by
which the compression region propagates upwards can be very roughly estimated
from the distance of $4.5\,$mm over $9\,\,\mu$s. This wave speed is only
$\approx 500\,$m/s! If we compare this value with the speed of the tensile
wave propagating through the virgin gelatin in  Figure 2(i), the wave travels
within $3\,\mu$s from the bottom to the top of the frame, a distance of
$5\,$mm, i.e. at a speed of about $1600\,$m/s. The high-speed recordings are
provided as supplementary materials.

Discussion — The determination of the wave speed in both experiment and
simulation is based on the tracking of the minimum volume bubbles. Figure 3(a)
shows the volume evolution of a number of bubbles in the numerical system. The
location of the minimum bubble radius provides the value of the effective
speed of the wave (the yellow line). For comparison, the wave speed in pure
liquid is also drawn (red line). A similar experimental observation is plotted
in Figure 3(b). The analysis is from the 7th shock through the 4% gelatin
sample. The volume of selected bubbles with increasing distances $z$ are
plotted as a function of time. The minimum volume positions of each bubble are
indicated with an asterisk. The yellow line fitted through these locations
reveals the propagation speed of the wave compressing the bubbles. For this
example, the wave speed is $U=385\pm 15$. The speed of sound in the pure
gelatin sample ($1605\,$m/s) is shown with a red line for comparison, see Ref.
[22].

Figure 3: (a) DNS simulated spatio-temporal diagram of the evolution of the
bubble radius of the bubbles in the system for $\beta=9.2\times 10^{-4}$ and
$\Delta p/p_{0}=100$. The yellow fitted curve represents the minimum radius
instances of bubbles. (b) Experimental volume dynamics of bubbles at different
heights $z$. The yellow line connects the minimum bubble volumes (red
crosses). Its slope is $385\,$m/s, the red line indicates a slope of
$1605\,$m/s, i.e. shock speed in pure gelatin. Figure 4:  Collected shock wave
velocities for the 4%, 6%, and 8% gelatin as a function of the void fraction.
The graph is compared to the acoustic speed of sound, $c$, from Eq. (1), the
shock velocity, $U$, from Eq. (2), and the results from DNS simulations. The
upper horizontal scale is the Mach number, $Ma=c^{-1}(\Delta p/\rho)^{1/2}$.

Experimental results from 61 individual measurements of the compression wave
speed in 6 gelatin samples with concentrations of $4\,$%, $6\,$%, and $8\,$%
by weight are shown in Figure 4. The lowest velocity measured is $298\,$m/s at
a void fraction of $1.5\,$%. The wave speed is a function of the void fraction
and the pressure amplitude. Thus, the measured pressure is fitted to a
polynomial function of order 2 to obtain an approximated function of the
pressure amplitude as a function of the void fraction. The dashed line in
Figure 4 is obtained from Eq. (2) using this fitted function for the pressure.
Additionally, we plot in Figure 4 Wood’s speed of sound Eq. (1) as a dashed-
dotted line. The latter predicts a much faster reduction of the propagation
speed than measured. Clearly, important characteristics of the finite
amplitude waves need to be accounted for. The expression of Campbell & Pitcher
Eq. (2), however, is in reasonable agreement with the experimental results and
is within the measurement errors for sufficiently high void fractions starting
from about $0.3\,$%. For smaller void fractions, the measured pressure
amplitude is above $70\,$bar. The combination of small void fraction and high-
pressure ratio $p_{1}/p_{0}$ limits the applicability of Eq. (2). Figure 4
also compares the measured wave speed with the low-frequency speed of bubble
compression from the numerical simulations (DNS) (solid line).

We use a Mach number to characterize the different wave regimes. By
considering the importance of nonlinear wave propagation phenomena, the Mach
number is quantified by the ratio between the characteristic flow velocity
$U_{c}$ and the wave propagation velocity in the linear regime of the flow
$c$, that is, the Mach number $Ma=U_{c}/c$. The characteristic flow velocity
$U_{c}$ is derived from the average velocity for a bubble to collapse, which
is $U_{c}=R_{max}/T_{C}$ where $R_{max}$ is the maximum bubble size and
$T_{C}$ the collapse time. Using the approximate relationship of the Rayleigh
(inertial) collapse time $T_{C}\approx R_{max}\sqrt{\rho/\Delta p}$, we obtain
for the average speed of collapse $U_{c}=\sqrt{\Delta p/\rho}$.

It is noted that the response of our system transits between the two limiting
regimes, the linear regime $Ma_{LF}\ll 1$, and the fully supersonic regime
$Ma_{HF}\geq 1$, where Campbell and Pitcher’s theory certainly applies. This
theory is indeed shown to provide a good first-order approximation to the
influence of the bubble concentration and pressure amplitude on the evolution
of the finite amplitude waves in the gelatin gel. This can be justified if we
assume that the energy exchanges between the long wavelengths and the short
wavelengths (or low and high frequencies) can be neglected at the leading
order. These effects only appear as a second-order correction when
$Ma_{LF}\,\geq\,1$. Overall, the DNS provides a better agreement with the
experimental data considering the simplifications involved. This agreement
supports the explanation that with increasing void fraction, a low-frequency
wave with an approximately constant speed compresses the bubbles.

Conclusion — We presented a Direct Numerical Simulation of a pressure pulse
through a line of mono-dispersed bubbles. The bubbles cause the dispersion of
the initial pulse into a high-frequency component (the precursor wave), and a
low-frequency component which is propagated by the collapsing bubbles. The
precursor wave, which travels at sound speed in pure liquid, vanishes quickly.
Essentially, the collapsing bubbles act as a low-band filter that mainly
propagates the low-frequency component of the initial pulse. The increase in
bubble size, and thus the void fraction, reduce both the amplitude and speed
of the propagated wave.

Experimentally, we used high-speed photography to observe the slow-down of the
compressive wave by following the bubble collapses. The gas bubbles are
generated in the rarefaction part of the first lithotripter shock administered
to the gelatin gel. As the bubbles grow by diffusion and thereby increase the
void fraction, both amplitude and velocity of the compressive wave are
reduced. Delays in the arrival of subsequent shock waves are observed by high-
speed photography.

The wave propagation in a bubbly gelatin gel can be characterized with a Mach
number, $\mathrm{Ma}=\frac{1}{c}\sqrt{\Delta p/\rho_{l}}$, which compares the
characteristic speed of the bubble collapses with the linear speed of sound in
the bubbly liquid. At low void fraction and small $\mathrm{Ma}$, the wave
propagation follows Wood [9]’s prediction for small amplitude waves. At the
high void fraction regime, the linear acoustic velocity approach the speed of
collapse, i.e. $\mathrm{Ma}$ closes to 1. In our experimental domain where a
finite amplitude wave is used, the experimental results span a Ma number
between about 0.2 to slightly above 1. The experimental data agree with Wood’s
and Campbell’s models only in the asymptotic limit of small and large void
fractions; the DNS simulation gives reasonable agreement over the whole
measurement range.

An important implication of this study is the effect of the void fraction on
focusing compression waves in tissue should be accounted for in medical
treatments such as in histotripsy [23], and high-intensity focused ultrasound.
While experimentally challenging, we expect that experiments in a liquid such
as water should provide similar results, which is predicted by our
simulations.

###### Acknowledgements.

J. M. R acknowledges support by the Alexander von Humboldt Foundation
(Germany) through a Georg Forster Research Fellowship. We acknowledge
constructive and insightful scientific exchanges with Dr. Evert Klaseboer and
Dr. Sun Qiang.

## References

  * Mallock [1910] H. R. A. Mallock, The damping of sound by frothy liquids, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 84, 391 (1910).
  * Minnaert [1933] M. Minnaert, Xvi. on musical air-bubbles and the sounds of running water, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 16, 235 (1933).
  * Carstensen and Foldy [1947] E. Carstensen and L. Foldy, Propagation of sound through a liquid containing bubbles, The Journal of the Acoustical Society of America 19, 481 (1947).
  * Meyer and Skudrzyk [1953] E. Meyer and E. Skudrzyk, Über die akustischen eigenschaften von gasblasen-schleiern in wasser, Acta Acustica united with Acustica 3, 434 (1953).
  * Wijngaarden [1968] L. V. Wijngaarden, On the equations of motion for mixtures of liquid and gas bubbles, Journal of Fluid Mechanics 33, 465–474 (1968).
  * Wijngaarden [1972] L. v. Wijngaarden, One-dimensional flow of liquids containing small gas bubbles, Annual review of fluid Mechanics 4, 369 (1972).
  * Commander and Prosperetti [1989] K. W. Commander and A. Prosperetti, Linear pressure waves in bubbly liquids: Comparison between theory and experiments, The Journal of the Acoustical Society of America 85, 732 (1989).
  * Fuster and Montel [2015] D. Fuster and F. Montel, Mass transfer effects on linear wave propagation in diluted bubbly liquids, Journal of Fluid Mechanics 779, 598 (2015).
  * Wood [1941] A. Wood,  _A textbook of sound: Bell_ (London, 1941).
  * Wijngaarden [2007] L. v. Wijngaarden, Shock waves in bubbly liquids, in  _Shock Wave Science and Technology Reference Library: Multiphase Flows I_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 2007) pp. 3–33.
  * Silberman [1957] E. Silberman, Sound velocity and attenuation in bubbly mixtures measured in standing wave tubes, The journal of the Acoustical Society of America 29, 925 (1957).
  * Campbell and Pitcher [1958] I. J. Campbell and A. Pitcher, Shock waves in a liquid containing gas bubbles, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 243, 534 (1958).
  * Kameda and Matsumoto [1996] M. Kameda and Y. Matsumoto, Shock waves in a liquid containing small gas bubbles, Physics of Fluids 8, 322 (1996).
  * Kameda  _et al._ [1998] M. Kameda, N. Shimaura, F. Higashino, and Y. Matsumoto, Shock waves in a uniform bubbly flow, Physics of Fluids 10, 2661 (1998).
  * Ando  _et al._ [2011] K. Ando, T. Colonius, and C. E. Brennen, Numerical simulation of shock propagation in a polydisperse bubbly liquid, International Journal of Multiphase Flow 37, 596 (2011).
  * Dear  _et al._ [1988] J. Dear, J. Field, and A. J. Walton, Gas compression and jet formation in cavities collapsed by a shock wave, Nature 332, 505 (1988).
  * Dear and Field [1988] J. Dear and J. Field, A study of the collapse of arrays of cavities, Journal of Fluid Mechanics 190, 409 (1988).
  * Fuster and Popinet [2018] D. Fuster and S. Popinet, An all-mach method for the simulation of bubble dynamics problems in the presence of surface tension, Journal of Computational Physics 374, 752 (2018).
  * Fan  _et al._ [2020] Y. Fan, H. Li, and D. Fuster, Optimal subharmonic emission of stable bubble oscillations in a tube, Physical Review E 102, 013105 (2020).
  * Rosselló and Ohl [2023] J. M. Rosselló and C.-D. Ohl, Clean production and characterization of nanobubbles using laser energy deposition, Ultrasonics Sonochemistry 94, 106321 (2023).
  * Ando and Shirota [2019] K. Ando and E. Shirota, Quasistatic growth of bubbles in a gelatin gel under dissolved-gas supersaturation, Physics of Fluids 31, 111701 (2019), https://doi.org/10.1063/1.5128887 .
  * Gojani  _et al._ [2016] A. Gojani, K. Ohtani, K. Takayama, and S. Hosseini, Shock hugoniot and equations of states of water, castor oil, and aqueous solutions of sodium chloride, sucrose and gelatin, Shock Waves 26, 63 (2016).
  * Maxwell  _et al._ [2011] A. D. Maxwell, T.-Y. Wang, C. A. Cain, J. B. Fowlkes, O. A. Sapozhnikov, M. R. Bailey, and Z. Xu, Cavitation clouds created by shock scattering from bubbles during histotripsy, The Journal of the Acoustical Society of America 130, 1888 (2011), https://doi.org/10.1121/1.3625239 .