# Structural evolution of the kagome superconductors $A$V3Sb5 ($A$ = K, Rb,
and Cs) through charge density wave order

Linus Kautzsch  Materials Department, University of California Santa Barbara,
Santa Barbara, CA, 93106, United States  Brenden R. Ortiz  Materials
Department, University of California Santa Barbara, Santa Barbara, CA, 93106,
United States  Krishnanand Mallayya  Department of Physics, Cornell
University, Ithaca, NY, 14853, United States  Jayden Plumb  Materials
Department, University of California Santa Barbara, Santa Barbara, CA, 93106,
United States  Ganesh Pokharel  Materials Department, University of California
Santa Barbara, Santa Barbara, CA, 93106, United States  Jacob P. C. Ruff
CHESS, Cornell University, Ithaca, NY, 14853, United States  Zahirul Islam
Advanced Photon Source, Argonne National Laboratory, Lemont, IL, 60439, United
States  Eun-Ah Kim  Department of Physics, Cornell University, Ithaca, NY,
14853, United States  Ram Seshadri  Materials Department, University of
California Santa Barbara, Santa Barbara, CA, 93106, United States  Stephen D.
Wilson<EMAIL_ADDRESS>Materials Department, University of California
Santa Barbara, Santa Barbara, CA, 93106, United States

###### Abstract

The kagome superconductors KV3Sb5, RbV3Sb5, and CsV3Sb5 are known to display
charge density wave (CDW) order which impacts the topological characteristics
of their electronic structure. Details of their structural ground states and
how they evolve with temperature are revealed here using single crystal X-ray
crystallographic refinements as a function of temperature, carried out with
synchrotron radiation. The compounds KV3Sb5 and RbV3Sb5 present
2$\times$2$\times$2 superstructures in the $Fmmm$ space group with a staggered
tri-hexagonal deformation of vanadium layers. CsV3Sb5 displays more complex
structural evolution, whose details have been unravelled by applying machine
learning methods to the scattering data. Upon cooling through the CDW
transition, CsV3Sb5 displays a staged progression of ordering from a
2$\times$2$\times$1 supercell and a 2$\times$2$\times$2 supercell into a final
2$\times$2$\times$4 supercell that persists to $T$ = 11 K and exhibits an
average structure where vanadium layers display both tri-hexagonal and Star of
David patterns of deformations. Diffraction from CsV3Sb5 under pulsed magnetic
fields up to $\mu_{0}H$ = 28 T suggest the real component of the CDW state is
insensitive to external magnetic fields.

††preprint: APS/123-QED

##  I Introduction

The charge density wave (CDW) instability is central to many of the
unconventional properties reported in the $A$V3Sb5 ($A=$ K, Rb, and Cs) family
of kagome superconductors [1, 2, 3]. Within the kagome lattice of these
materials, the CDW is hypothesized to harbor both real and imaginary
components within the resulting order parameter, with the former corresponding
to a real space charge inhomogeneity and the latter corresponding to a
“chiral” flux phase that breaks time reversal symmetry [4]. The real component
of the CDW results in a superstructure deformation of the crystal lattice.
Precise details of the structural pattern of this deformation has been an
active area of recent investigation.

Understanding the pattern of the real space structural deformation stands to
inform a number of theoretical models of the anomalous CDW transtion in AV3Sb5
compounds. For instance, models of the band folding through the transition and
calculations of the resulting Berry curvature may help inform the origin of
the large anomalous Hall effect reported in the CDW state [5, 6]. At high
temperatures, the structure of all three variants assumes an ideal
(undistorted) kagome network of vanadium ions in the space group $P6/mmm$ as
shown in Fig. 1 [1]. The lattice deformation that results upon cooling through
the CDW transition primarily arises from the displacement of the kagome
network of vanadium atoms [7]. The most energetically favored possibilities
are predicted to be one of two breathing modes: either the Star of David-type
(SoD) deformation or its inverse, a tri-hexagonal (TrH) deformation [8, 9,
10]. Both are triple-q distortion modes, suggestive of the influence of
nesting effects between the saddle-points ($M$-points) close to the Fermi
level in AV3Sb5.

Figure 1: (a) and (b) show the refined room temperature ($T=290$ K) structure
of the $A$V3Sb5 ($A=$ K, Rb, and Cs) family in the hexagonal space group
$P6/mmm$.

An added complexity in solving the pattern of CDW order in AV3Sb5 is that
interlayer interactions are strong enough to promote a fully three-dimensional
CDW. This has been observed in scanning tunneling microscopy (STM)
measurements as a half unit cell phase-shift between layers resolved across
step-edges [11] as well as direct observation of superlattice reflections with
a finite $q_{z}$ propagation wave vector in x-ray diffraction measurements [7,
12]. The out-of-plane $q_{z}$ periodicity varies between AV3Sb5 variants with
the unit cell doubled along the c-axis in A= K, Rb variants [8, 12] and the
unit cell at least partially quadrupled in A= Cs [7, 13, 14]. The increased
interlayer spacing driven by the larger Cs interstitial layer seemingly
modifies the interactions between kagome planes and drives the formation of a
delicate three-dimensional CDW phase whose modulation between kagome planes
potentially hosts phase coexistence and metastable states.

Metastability in the CDW phase of CsV3Sb5 has been primarily observed with
x-ray diffraction measurements. Initial reports found conflicting observations
of $2\times 2\times 2$ [12] and $2\times 2\times 4$ CDW superlattices [7],
suggesting sensitivity to growth conditions. Furthermore, recent x-ray
diffraction studies have reported the coexistence of $2\times 2\times 2$ and
$2\times 2\times 4$ structures with different onset temperatures and thermal
quenching/annealing behaviors [13, 14]. Conclusions regarding which state is
stabilized by annealing also vary, likely reflecting the added variable of
disorder introduced via the growth process in the ground state selection. For
instance, recent studies have shown that small amounts of kagome-plane dopants
can quench the $2\times 2\times 4$ state, leaving only a quasi-two-dimensional
$2\times 2\times 2$ phase remaining [15].

The above complexity of the CDW state in CsV3Sb5 is relevant as it is the most
easily grown and widely studied AV3Sb5 variant. In our prior synchrotron
study, an average structural solution for CsV3Sb5 was presented presuming a
minimal model where the system maintains hexagonal symmetry ($P\bar{3}$). This
average structure suggested a superstructure comprised of modulated TrH and
SoD distortions in a $2\times 2\times 4$ cell [7]. This prior solution could
also be indexed as an orthorhombic cell; however twinning/domain formation was
not directly resolved in the diffraction data. Recently, however, STM [16] and
optical measurements [17] have resolved three sets of orthorhombic twin
domains in the CDW state of all three $A$V3Sb5 compounds. This provides an
excellent basis for developing twinned orthorhombic structural models of these
systems via synchrotron x-ray diffraction measurements.

Here, we report the structural deformations resulting from CDW order in all
three AV3Sb5 compounds determined via single crystal synchrotron x-ray
diffraction measurements. Structural refinements using an orthorhombic, three
domain, twinned model are presented for all three compounds. For KV3Sb5 and
RbV3Sb5, a 2$\times$2$\times$2 superstructure is observed and structural
refinement, combined with recent NMR analysis [18, 19, 20], allows for a
staggered TrH deformation of the vanadium layers to be resolved and the atomic
positions determined. We also revisit the average structure of CsV3Sb5 in an
orthorhombic cell using the same three domain structure, where we identify a
2$\times$2$\times$4 superlattice in the ground state. The average low-
temperature structural solution remains a modulation of SoD and TrH-type
deformations, consistent with prior reports [7, 21, 22]. Curiously, we also
resolve a cross-over regime below $T=94$ K where a 2$\times$2$\times$1 and a
2$\times$2$\times$2 lattice form and then diminish at lower temperatures.
Finally, as a test for potential time reversal symmetry breaking in the charge
density wave order parameter, we explore the field dependence of the
superstructure in CsV3Sb5 under the application of a $\mu_{0}H=$ 28 T pulsed
magnetic field.

##  II Methods

###  II.1 Structural solutions

Synchrotron x-ray diffraction experiments were carried out at the QM2 beam
line at CHESS. The incident x-ray wavelength of $\lambda=0.41328$ $\rm\AA$ was
selected using a double-bounce diamond monochromator. A stream of cold flowing
helium gas was used to cool the sample. The diffraction experiment was
conducted in transmission geometry using a 6-megapixel photon-counting pixel-
array detector with a silicon sensor layer. Data was collected in full 360∘
sample rotations with a step size of 0.1∘. Scattering planes in reciprocal
space were visualized using the NeXpy software package. The diffraction data
was indexed and integrated using the APEX3 software package including
absorption and extinction correction. Crystallographic structural solutions
were determined using the SHELX software package [23]. An unsupervised machine
learning tool using Gaussian Mixture Models called X-ray Temperature
Clustering (X-TEC) [24] was used to identify classes of peaks in reciprocal
space with distinct temperature dependencies.

###  II.2 Pulsed magnetic field experiments

Pulsed magnetic field experiments were carried out at the beam line 6 ID-C of
the Advanced Photon Source at Argonne National Laboratory using a dual
cryostat, single solenoid pulsed-magnet system described in Ref. [25]. X-rays
with a wavelength of $\lambda=0.56356$ $\rm\AA$ were used, and the
crystallographic $c$-axis of CsV3Sb5 was aligned parallel to the direction of
the incident x-ray beam and the direction of the magnetic field. A low-noise
Lambda 750K detector was used with $55\times 55$ $\mu$ m2 CdTe pixels
providing near 100% efficiency. The detector was run in 12-bit mode at a 1 kHz
repetition rate with XRD pattern collected every 1 ms. Multiple frames were
collected over the course of repeated 9 ms pulses with a maximum magnetic
field of $H_{\mathrm{max}}=28$ T.

Figure 2: Refined structural models for KV3Sb5 and RbV3Sb5. (a) shows the z =
0.5 layer comprised of vanadium and antimony atoms of the refined room
temperature structure of KV3Sb5 and RbV3Sb5 in the hexagonal space group
$P6/mmm$. (b) and (c) show the suggested TrH model for the 2$\times$2$\times$2
superstructure in the orthogonal space group $Fmmm$. The vanadium planes at
$z=0$ and $z=0.5$ in the $Fmmm$ structure are displayed.

##  III Experimental Results

###  III.1 Structure of KV3Sb5 and RbV3Sb5

We begin by discussing diffraction experiments on KV3Sb5 at $T=10$ K.
Scattering data indicate a 2$\times$2$\times$2 superstructure, and, similar to
prior measurements [7], the area detector resolution does not allow for the
momentum-space resolution of twinning in the distorted state. Recent scanning
optical studies have shown the presence of three structural twins that form
below the CDW transition [17]. These domains have their principle axes rotated
by 120∘ in the ab-plane, consistent with an underlying orthorhombic lattice.
This would imply a pseudohexagonal twinning with three orthorhombic twins that
are rotated 120∘ with respect to each other in the $ab$-plane, which in turn
raises the apparent symmetry of the diffraction pattern to hexagonal.
Theoretical investigations identified the space groups $Cmcm$ [26] and $Fmmm$
[10] as possible starting orthorhombic solutions to the twinned
pseudohexagonal structure.

As reflection intensities indicate a Laue class of $mmm$, model refinements
were performed in both the face-centered orthorhombic space group $Fmmm$ and
the base-centered space group $Cmmm$. Refinement in both space groups yields
identical deformations of the vanadium layers with the distinction being that
$Cmmm$ allows for unique distortions within neighboring vanadium layers.
Refinement in $Cmmm$ also implicitly allows for any intensity at positions
forbidden in $Fmmm$ to be accounted for outside of twinning models. Analysis
in these two possible groups gave comparable or slightly better $R_{1}$ values
for the higher symmetry $Fmmm$ space group, and this group was used for the
final low-temperature structural determination.

KV3Sb5 data refined in $Fmmm$ is well-represented by two structural models
(Tab. 1) consisting of either a staggered TrH ($R_{1}=6.22$, Fig. 2) or a
staggered SoD distortion ($R_{1}=6.31$) of the vanadium layers. We note here
that the assignments of TrH or SoD-type distortions in this paper are based on
the dominant displacement type of V atoms, and that in all solutions there is
detectable variance amongst V-V distances within a single plane (ie. there is
not a single V-V distance in a given distorted kagome plane). This is further
parameterized in Discussion section of this manuscript. The pattern of a given
distortion type is shifted by half of a lattice constant ($0.5a$) with respect
to neighboring layers, and the two solutions differ in the sign of the
displacement vectors of the vanadium atoms in the $ab$-plane.

X-ray refinement alone is unable to discriminate between the distortion types,
which is a result of the previously discussed twinning. However, recent NMR
studies identify a TrH-type distortion and allow us to break the degeneracy
between these solutions [18, 19]. The resulting low-temperature structure
generated from the staggered TrH solution of KV3Sb5 is plotted in Fig. 2 (b).
This is consistent with recent studies [13, 27, 9], and crystallographic
information files (CIF) for both TrH and SoD distortion types can be found in
the supplemental information. After collecting data at 10 K, the sample was
warmed to 290 K and the high temperature structure measured. This undistorted
state was refined in $P6/mmm$ at 290 K with bond lengths illustrated in Fig. 2
(a).

Similar analysis was employed for the x-ray diffraction data of RbV3Sb5, where
RbV3Sb5 also forms a 2$\times$2$\times$2 superstructure in the CDW state. At
$T=10$ K, the data is again best represented when indexed in the space group
$Fmmm$, and our refinements show superior solutions using a staggered TrH
model ($R_{1}=6.40$) versus a staggered SoD model ($R_{1}=7.41$). This is
consistent with NMR results [18], and the TrH solution is selected as the
correct structure for presentation and the resulting bond lengths in the
vanadium network are shown in Fig. 2 (b). For reference, CIF files for both
TrH and SoD distortion types can also be found in the supplemental information
[28]. After measurement at 10 K, the 290 K structure was determined at 290 K
in $P6/mmm$ with bond lengths illustrated in Fig. 2 (a).

###  III.2 Structure of CsV3Sb5

Figure 3: Temperature dependence of the unit cell parameters of CsV3Sb5 upon
traversing the CDW transition. The evolution of the $a$-, $b$-, and $c$-axes
are plotted in panel (a). (b) shows $a$ and $b$ where $a_{\mathrm{hex}}=a/2$
for $T\leq 90$ K and $b_{\mathrm{hex}}=\sqrt{b^{2}+4a_{hex}^{2}}/4$ for $T\leq
90$ K and (c) shows $c$ where $c_{\mathrm{hex}}=c/2$ for $T=90$ K and
$c_{\mathrm{hex}}=c/4$ for $T<90$ K. The subcell parameters undergo no
significant changes. Figure 4: X-TEC analysis of CsV3Sb5 diffraction data. (a)
X-TEC identifies three distinct temperature trajectories (blue, green, and
yellow clusters) whose intensities exhibit a sharp onset near the CDW
transition temperature. The lines denote the mean and shading denotes the
standard deviation of the re-scaled (z-scored) intensities $\tilde{I}(T)$
within each cluster. The grey clusters identify the diffuse background
scattering. (b) The average intensity of the blue, green, and yellow cluster
reveals the different onset temperatures for the $2\times 2\times 4$ CDW
($T_{2\times 2\times 4}\approx 85K$) and $2\times 2\times 2$ CDW ($T_{2\times
2\times 2}\approx 90K$) state and the $2\times 2\times 1$ state emerging
between 95 K and 80 K. $\Delta I(T)$ is the average intensity after
subtracting its minimum value. (c) A section of the (H, 2.5, L) plane, in
reciprocal lattice units (r.l.u.) and indexed in P6/mmm, showing the pixels
colored by the cluster assignments of their intensities. (d,e,f) Distribution
of the out-of-plane momentum $L$ of the peaks in the respective clusters.

Now turning to analysis of the structure of CsV3Sb5, temperature dependent
synchrotron x-ray diffraction experiments were carried out on a single crystal
of CsV3Sb5. The final data set was collected at $T=290$ K where data were
indexed in the hexagonal space group $P6/mmm$. The Bragg peaks of the CsV3Sb5
crystal studied are narrow with minimal stacking disorder, though the
crystallinity is anisotropic. Quantifying this with Gaussian fits to
representative Bragg peaks shows minimum correlation lengths in the $ab$-plane
of $\xi_{min}=890$ $\rm\AA$ or $\xi_{min}=1000$ $\rm\AA$ (depending on the
pixel-limited cut direction) and correlation lengths of $\xi=400$ $\rm\AA$
along the $c$-axis, indicating excellent crystallinity of the sample [28]. The
lower correlation length along the $c$-axis (interlayer direction) is typical
in layered materials and beyond the resolution of experiments using relaxed
collimation.

Figure 5: Refined structural model for CsV3Sb5. (a) shows the kagome network
(z = 0.5) of the 290 K room temperature structure in the hexagonal space group
$P6/mmm$. (b), (c), (d) and (e) represent a single-phase model of the
2$\times$2$\times$4 superstructure in the orthorhombic space group $Cmmm$.
Three kagome layers form with TrH deformations that are capped by a layer with
a predominant SoD deformation.

Temperature-dependent diffraction data were collected on warming after using
an initial cooling rate of approximately 10 K/min to cool to 10 K. Data were
then collected by warming at $\approx 5$ K/min with pauses for measurement
scans (12 min/scan) performed at select temperatures. Collecting a data series
using this cooling/warming profile show that at temperatures below $T=94$ K,
superstructure reflections appear as expected. Just below the CDW transition,
a clear staging behavior is observed. At 90 K, superlattice reflections
consistent with a 2$\times$2$\times$2 supercell are observed due to the
appearance of half-integer $L$ peaks as well as half-integer peaks in the
($H$, $K$)-plane. Below this temperature, at 85 K, weak quarter-integer $L$
peaks also appear, signifying the onset of an enlarged 2$\times$2$\times$4
supercell. The cell parameters through the CDW transition are plotted in Fig.
3

Figure 6: Magnetic field dependence of x-ray diffraction of the
Q=(0.5,0.5,0.25) CDW superlattice reflection in CsV3Sb5 in a pulsed magnetic
field experiment with $H_{\mathrm{max}}=28$ T at $T=11$ K and $T=82.5$ K. (a)
and (b) CDW reflection at $T=82.5$ K and $t=5$ ms from 1 sweep (1 ms exposure
in (a)) and integrated over 61 sweeps (61 ms exposure in (b)). (c) Magnetic
field $H$ as a function of the experiment time in milliseconds. (d) Sum of the
integrated intensity of the (0.5,0.5,0.25) superlattice reflection from a
total of 45 pulsed field sweeps at $T=11$ K and 61 sweeps at $T=82.5$ K.

It is clear that the nature of the charge density wave and the potential
coexistence/competition between different states involves subtle interactions
and changes in the diffraction patterns. In addition, the large temperature-
dependent and Q-dependent data sets contain substantial quantities of
information that can be overlooked during standard data reduction and
integration processes. Recently, an unsupervised machine learning tool called
X-ray Temperature Clustering (X-TEC) was developed to enable analysis of
parametric, complex data sets. X-TEC uses a Gaussian Mixture Model to identify
peaks in reciprocal space with distinct temperature dependencies. In order to
distinguish the functional form of the intensity-temperature trajectory rather
than their magnitudes, the intensity of each momenta $\vec{q}$ is re-scaled
(z-scored) as
$\tilde{I}_{\vec{q}}(T)=\left(I_{\vec{q}}(T)-\mu_{\vec{q}}\right)/\sigma_{\vec{q}}$
where $\mu_{\vec{q}}$ is the mean over temperature $T$, and $\sigma_{\vec{q}}$
is the standard deviation in $T$ [24].

By applying X-TEC on the CsV3Sb5 data, $\sim 15,000$ peaks were identified
with three distinct temperature trajectories color-coded as yellow, green and
blue clusters. These clusters exhibit an onset behavior around the CDW
transition temperature as shown in Figs. 4(a) and (b). By looking at the
locations of these clusters in reciprocal space (Fig. 4(c)), we find that all
three clusters have the same in-plane momenta at $(H,K)\equiv(0.5,0.5)$, but
differ in their out-of-plane momentum $L$. While the green clusters are
located at $L=1/4$ (Fig. 4(e)) corresponding to the $2\times 2\times 4$
structure, the blue clusters derive from the $L=1/2$ peaks (Fig. 4(f)) of the
$2\times 2\times 2$ structure, and the yellow clusters to an additional class
of $L=0$ peaks (Fig. 4 (d)) representing a $2\times 2\times 1$ structure. From
their intensity-temperature trajectories (Fig. 4(b)), we find that the onset
of CDW proceeds in two stages. First, the $2\times 2\times 1$ and $2\times
2\times 2$ cells appear below temperature $\approx 95$ K, followed by the
onset of a $2\times 2\times 4$ state at $\approx 85$ K. Remarkably, with the
onset of the $2\times 2\times 4$ peaks, the $2\times 2\times 1$ peaks are
strongly suppressed and vanish below 80 K, and, additionally, the $2\times
2\times 2$ peaks initially decrease with the onset of $2\times 2\times 4$
order and then recover their intensity at lower temperatures. This suggests an
exchange of scattering weight between order parameters.

Table 1: Parameters of structural models for the room temperature and the low-temperature structures of KV3Sb5, RbV3Sb5, and CsV3Sb5. All .cif files from the refinements can be found in the supplemental information [28]. | KV3Sb5 | RbV3Sb5 | CsV3Sb5
---|---|---|---
Wavelength ($\rm\AA$) | 0.41328
Temperature (K) | 290 | 10 | 290 | 10 | 290 | 90 | 11
Space group | $P6/mmm$ | $Fmmm$ | $P6/mmm$ | $Fmmm$ | $P6/mmm$ | $Fmmm$ | $Cmmm$
Type |  | $2\times 2\times 2$ |  | $2\times 2\times 2$ |  | $2\times 2\times 2$ | $2\times 2\times 4$
Distortion |  | TrH | SoD |  | TrH | SoD |  | TrH | SoD | TrH/SoD
$N_{\mathrm{measured}}$ | 2751 | 17014 | 2239 | 16168 | 2888 | 6788 | 24048
$N_{\mathrm{independent}}$ | 232 | 2401 | 201 | 2444 | 238 | 2159 | 5002
$N_{\mathrm{sig.}}$ [$I>2\sigma(I)$] | 231 | 1617 | 201 | 1906 | 238 | 641 | 3098
$a$ ($\rm\AA$) | 5.4260(9) | 10.9576(9) | 5.4941(6) | 11.0065(9) | 5.5236(6) | 11.0275(9) | 11.005(2)
$b$ ($\rm\AA$) | 5.4260(9) | 18.9813(16) | 5.4941(6) | 19.0575(16) | 5.5236(6) | 19.1082(16) | 19.046(4)
$c$ ($\rm\AA$) | 8.845(2) | 17.8552(15) | 9.1071(14) | 18.1010(15) | 9.3623(15) | 18.6213(15) | 37.176(7)
$\rho$ (g/cm3) | 5.90 | 5.73 | 5.91 | 5.93 | 6.00 | 6.06 | 7.22
$R_{1}>4\sigma$ (%) | 8.62 | 6.22 | 6.31 | 10.42 | 6.40 | 7.41 | 5.91 | 7.71 | 7.73 | 10.0
$R_{1}$ all (%) | 8.63 | 6.44 | 6.56 | 10.42 | 6.74 | 7.82 | 5.91 | 9.61 | 9.62 | 11.55
Twin fraction 1 | - | 0.35 | 0.36 | - | 0.32 | 0.32 | - | 0.31 | 0.30 | 0.31
Twin fraction 2 | - | 0.33 | 0.33 | - | 0.33 | 0.33 | - | 0.36 | 0.36 | 0.34
Twin fraction 3 | - | 0.32 | 0.31 | - | 0.34 | 0.34 | - | 0.33 | 0.34 | 0.35

The 2$\times$2$\times$4 superstructure persists down to temperatures of 11 K.
We note here that the in-plane and out-of-plane correlation lengths of the CDW
superlattice peaks here are identical to the primary structural peaks. This
indicates that the 2$\times$2$\times$4 CDW state is long-range with
correlation lengths constrained only by the sample’s underlying crystallinity.
Notably, our diffraction data using this cooling/warming profile do not
resolve anomalies at $T=60$ K as seen in magnetotransport experiments [29, 30]
or a crossover to a 2$\times$2$\times$2 state similar to Ref. [13]. This
difference potentially arises from subtle effects attributed to the different
cooling profiles during measurement and potentially from sample dependent
effects imparted during the crystal growth process.

We now discuss the structural refinement in the 2$\times$2$\times$2 regime at
$T=90$ K and the 2$\times$2$\times$4 superstructure at lower temperatures.
Using the same twinning approach described earlier, the 2$\times$2$\times$2
state at 90 K is best fit in the $Fmmm$ group rendering either a staggered SoD
or TrH distortion (CIFs provided in the supplemental material [28]). This is
consistent with the low temperature state recently reported by Stahl et al.
[13] where only a 2$\times$2$\times$2 superstructure was observed. Upon
cooling into the 2$\times$2$\times$4 state, the data at 10 K are instead best
refined in the lower symmetry, base-centered $Cmmm$ group. Fig. 5 shows the
resulting model of the superstructure consisting of vanadium layers with both
SoD and TrH deformations. The three TrH layers are staggered along the
$c$-axis by $0.5a$ (Fig. 5 (d),(e),(f)) and are then capped by a SoD layer
(Fig. 5 (c)). The model assumes a uniform phase at this temperature and
represents the measured diffraction data well, especially considering the
large number of significant reflections (3098 shown in Tab. 1). This mixed
distortion-type result is consistent with previous refinements in a
$P\overline{3}$ average structure containing both SoD and TrH layers [3] as
well as recent angle-resolved photoemission data [21, 22]. Mixed reports from
NMR measurements [19, 20] provide a picture interpreted as either SoD or TrH
distortion types, and it is unclear whether a mixed distortion-type structure
can also be modeled via the data. Moreover, our X-TEC analysis indicates that
a 2$\times$2$\times$2 phase competes with and potentially coexists with the
2$\times$2$\times$4 phase at lowest temperatures and that the stacking
sequence of the TrH- and SoD-like layers may be influenced by this
coexistence.

###  III.3 Pulsed magnetic field experiment

Scanning tunneling microscopy measurements have reported that a magnetic field
couples to the apparent Fourier weight attributed to CDW order at the surface
in $A$V3Sb5 compounds [8], suggesting the presence of a chiral CDW state that
breaks time reversal symmetry. As means of testing whether a magnetic field
couples to the real component of the CDW order parameter in the bulk, the
magnetic field dependence of the superlattice reflections was investigated. A
pulsed magnetic field of $H_{\mathrm{pulse}}=28$ T was used, and the geometry
of the inner bore of the magnet coil allows for the measurement of the
Q$=(0.5,0.5,0.25)$ superlattice reflection in transmission geometry with the
field applied parallel to the $c$-axis. Scattering data were collected at
$T=11$ K and at $T=82.5$ K after cooling the sample into the CDW state under
zero field conditions. Magnetic field pulses were applied that last for
approximately 10 ms (half-sine, start to finish), and where the maximum
magnetic field of 28 T is reached. Multiple pulses were used to improve
measurement statistics.

The Q$=(0.5,0.5,0.25)$ peak was measured at $T=11$ K for a total of 45 pulses
and at $T=82.5$ K for a total of 61 pulses. For illustration, Fig. 6 (a) shows
the result of 1 ms of data collection during a single field pulse on the
Q$=(0.5,0.5,0.25)$ superlattice peak on a pixel array. Fig. 6 (b) then shows
the same peak after a cumulative sum of 61 pulses. Each time bin corresponds
to a field value shown in Fig. 6 (c).

The integrated intensities of the Q$=(0.5,0.5,0.25)$ peak as a function of
time were summed over the number of individual pulses to increase the signal-
to-noise ratio with the results plotted in Fig. 6 (d). No significant change
of the integrated intensity outside of the experimental error could be
observed at either of the two temperatures explored. The pre-pulse
fluctuations and sharp post-pulse decrease in the integrated intensity are
attributed to mechanical vibrations. Within our experimental resolution, we do
not observe any evidence of magnetic-field dependence of the superlattice peak
intensity, consistent with a recent report [13] and signifying that the real
component of the CDW order parameter is unchanged.

##  IV Discussion and Conclusions

Tab. 1 summarizes the parameters of the crystallographic refinements carried
out on KV3Sb5, RbV3Sb5, and CsV3Sb5 crystals at key temperatures discussed in
the preceding sections. All refined models represent the collected data well
and the $R$-values are low considering the high number of significant
reflections collected in synchrotron experiments. The CDW distortion results
in an average V-V bond length that is nearly the same in all three AV3Sb5
variants with (V-V)avg=2.7022$\pm 0.022$ in A = K, (V-V)avg=2.6992$\pm 0.028$
in A = Rb, and (V-V)avg=2.6979$\pm 0.036$ in A = Cs. The standard deviations
in these V-V distances highlight the fact that the assignments of TrH or SoD
distortions are only effective classifications and that the V-V bond lengths
have substantial variance within a given layer.

The increase in the standard deviation of distorted bond lengths about a fixed
average value likely indicates a progressively increasing instability of the
TrH state as the A-site changes from K to Rb to Cs. Notably, the Cs system has
the largest room temperature V-V distance, rendering its low-temperature
reconstruction of the kagome network the most dramatic of the series. This is
likely reflected in the larger heat capacity anomaly in the Cs system relative
to the other two as well as its accommodation of a mixed TrH and SoD
distortion.

The mixed SoD and TrH state refined in the CDW state of CsV3Sb5 should also be
viewed through the lens of the experimental parameters used in this study.
Slower cooling or faster quenching may modify the resulting structure in the
CDW state given the seeming competition with the $2\times 2\times 2$ TrH state
and a newly observed $2\times 2\times 1$ state that both appear just below
$T_{CDW}$. It is possible that materials with different levels of strain
frozen during the growth process or dilute impurity content may modify this
behavior as well. Subtle disorder along the interlayer direction in the form
of variable crystallinity between samples does not seem to play a dominant
role as we have screened many crystals with variable quality and they all
exhibit the same $2\times 2\times 4$ low-temperature superlattice using this
same cooling profile.

In summary, we studied the structural deformations related to the CDW
transitions in the systems CsV3Sb5, KV3Sb5, and RbV3Sb5 using synchrotron
x-ray diffraction. Our diffraction data allow for twinned models of
orthorhombic superstructures in the CDW state. We find 2$\times$2$\times$2
superstructures in KV3Sb5 and RbV3Sb5 and our crystallographic refinements
allow us to parameterize the distortion of vanadium-based kagome layers in a
staggered TrH state. Temperature dependent diffraction data from the CsV3Sb5
system indicate a complex evolution of the structural superlattice. Initially,
a 2$\times$2$\times$1 state appears along with a 2$\times$2$\times$2 state at
temperatures below 90 K, and this is followed by the onset of a
2$\times$2$\times$4 superlattice below 85 K. Structural refinement of this
larger superstructure reveals a three-layer staggered TrH deformation capped
by a layer with a SoD type distortion. Scattering data from CsV3Sb5 measured
under a 28 T pulsed magnetic field fails to resolve a change in the CDW
superlattice intensity, setting bounds on the response of the real component
of the CDW state to the application of a magnetic field.

##  V Acknowledgments

SDW acknowledges helpful discussions with Turan Birol, Rafael Fernandes, and
Binhai Yan. This work was supported by the National Science Foundation (NSF)
through Enabling Quantum Leap: Convergent Accelerated Discovery Foundries for
Quantum Materials Science, Engineering and Information (Q-AMASE-i): Quantum
Foundry at UC Santa Barbara (DMR-1906325). The research made use of the shared
facilities of the NSF Materials Research Science and Engineering Center at UC
Santa Barbara (DMR- 1720256). The UC Santa Barbara MRSEC is a member of the
Materials Research Facilities Network. (www.mrfn.org). This work is based upon
research conducted at the Center for High Energy X-ray Sciences (CHEXS) which
is supported by the National Science Foundation under award DMR-1829070. The
high-field pulsed magnet and a choke coil were installed at the Advanced
Photon Source through a partnership with International Collaboration Center at
the Institute for Materials Research (ICC-IMR) and Global Institute for
Materials Research Tohoku (GIMRT) at Tohoku University. The development and
use of X-TEC by K.M. and E.-A.K. was supported by the U.S. Department of
Energy, Office of Basic Energy Sciences, Division of Materials Sciences and
Engineering. KM and E-AK are supported in part by the Gordon and Betty Moore
Foundation’s EPiQS Initiative, Grant GBMF10436.

## References

  * Ortiz  _et al._ [2019] B. R. Ortiz, L. C. Gomes, J. R. Morey, M. Winiarski, M. Bordelon, J. S. Mangum, I. W. Oswald, J. A. Rodriguez-Rivera, J. R. Neilson, S. D. Wilson, _et al._ , New kagome prototype materials: discovery of kv3sb5, rbv3sb5, and csv3sb5, Physical Review Materials 3, 094407 (2019).
  * Ortiz  _et al._ [2020] B. R. Ortiz, S. M. Teicher, Y. Hu, J. L. Zuo, P. M. Sarte, E. C. Schueller, A. M. Abeykoon, M. J. Krogstad, S. Rosenkranz, R. Osborn, _et al._ , Csv3sb5: A $z$2 topological kagome metal with a superconducting ground state, Physical Review Letters 125, 247002 (2020).
  * Ortiz  _et al._ [2021a] B. R. Ortiz, P. M. Sarte, E. M. Kenney, M. J. Graf, S. M. Teicher, R. Seshadri, and S. D. Wilson, Superconductivity in the z 2 kagome metal kv 3 sb 5, Physical Review Materials 5, 034801 (2021a).
  * Park  _et al._ [2021] T. Park, M. Ye, and L. Balents, Electronic instabilities of kagome metals: Saddle points and landau theory, Phys. Rev. B 104, 035142 (2021).
  * Yang  _et al._ [2020] S.-Y. Yang, Y. Wang, B. R. Ortiz, D. Liu, J. Gayles, E. Derunova, R. Gonzalez-Hernandez, L. Šmejkal, Y. Chen, S. S. Parkin, _et al._ , Giant, unconventional anomalous hall effect in the metallic frustrated magnet candidate, kv3sb5, Science advances 6, eabb6003 (2020).
  * Yu  _et al._ [2021] F. Yu, T. Wu, Z. Wang, B. Lei, W. Zhuo, J. Ying, and X. Chen, Concurrence of anomalous hall effect and charge density wave in a superconducting topological kagome metal, Physical Review B 104, L041103 (2021).
  * Ortiz  _et al._ [2021b] B. R. Ortiz, S. M. Teicher, L. Kautzsch, P. M. Sarte, N. Ratcliff, J. Harter, J. P. Ruff, R. Seshadri, and S. D. Wilson, Fermi surface mapping and the nature of charge-density-wave order in the kagome superconductor CsV3Sb5, Phys. Rev. X 11, 041030 (2021b).
  * Jiang  _et al._ [2021] Y.-X. Jiang, J.-X. Yin, M. M. Denner, N. Shumiya, B. R. Ortiz, G. Xu, Z. Guguchia, J. He, M. S. Hossain, X. Liu, _et al._ , Unconventional chiral charge order in kagome superconductor kv3sb5, Nature Materials 20, 1353 (2021).
  * Tan  _et al._ [2021] H. Tan, Y. Liu, Z. Wang, and B. Yan, Charge density waves and electronic properties of superconducting kagome metals, Physical review letters 127, 046401 (2021).
  * Christensen  _et al._ [2021] M. H. Christensen, T. Birol, B. M. Andersen, and R. M. Fernandes, Theory of the charge density wave in a v 3 sb 5 kagome metals, Physical Review B 104, 214513 (2021).
  * Liang  _et al._ [2021] Z. Liang, X. Hou, F. Zhang, W. Ma, P. Wu, Z. Zhang, F. Yu, J.-J. Ying, K. Jiang, L. Shan, Z. Wang, and X.-H. Chen, Three-dimensional charge density wave and surface-dependent vortex-core states in a kagome superconductor ${\mathrm{csv}}_{3}{\mathrm{sb}}_{5}$, Phys. Rev. X 11, 031026 (2021).
  * Li  _et al._ [2021] H. Li, T. T. Zhang, T. Yilmaz, Y. Y. Pai, C. E. Marvinney, A. Said, Q. W. Yin, C. S. Gong, Z. J. Tu, E. Vescovo, C. S. Nelson, R. G. Moore, S. Murakami, H. C. Lei, H. N. Lee, B. J. Lawrie, and H. Miao, Observation of unconventional charge density wave without acoustic phonon anomaly in kagome superconductors ${A\mathrm{V}}_{3}{\mathrm{sb}}_{5}$ ($a=\mathrm{Rb}$, cs), Phys. Rev. X 11, 031050 (2021).
  * Stahl  _et al._ [2022] Q. Stahl, D. Chen, T. Ritschel, C. Shekhar, E. Sadrollahi, M. Rahn, O. Ivashko, M. v. Zimmermann, C. Felser, and J. Geck, Temperature-driven reorganization of electronic order in csv3sb5, Physical Review B 105, 195136 (2022).
  * Xiao  _et al._ [2022] Q. Xiao, Y. Lin, Q. Li, X. Zheng, W. Xia, S. Zhang, Y. Guo, J. Feng, and Y. Peng, Coexistence of multiple stacking charge density waves in kagome superconductor ${\mathrm{csv}}_{3}{\mathrm{sb}}_{5}$, arXiv preprint arXiv:2201.05211 (2022).
  * Kautzsch  _et al._ [2022] L. Kautzsch, Y. M. Oey, H. Li, Z. Ren, B. R. Ortiz, R. Seshadri, J. Ruff, Z. Wang, I. Zeljkovic, and S. D. Wilson, Incommensurate charge-stripe correlations in the kagome superconductor csv3sb5-xsnx, arXiv preprint arXiv:2207.10608 (2022).
  * Li  _et al._ [2022] H. Li, H. Zhao, B. R. Ortiz, T. Park, M. Ye, L. Balents, Z. Wang, S. D. Wilson, and I. Zeljkovic, Rotation symmetry breaking in the normal state of a kagome superconductor kv3sb5, Nature Physics 18, 265 (2022).
  * Xu  _et al._ [2022] Y. Xu, Z. Ni, Y. Liu, B. R. Ortiz, Q. Deng, S. D. Wilson, B. Yan, L. Balents, and L. Wu, Three-state nematicity and magneto-optical kerr effect in the charge density waves in kagome superconductors, Nature Physics 10.1038/s41567-022-01805-7 (2022).
  * Frassineti  _et al._ [2023] J. Frassineti, P. Bonfà, G. Allodi, E. Garcia, R. Cong, B. R. Ortiz, S. D. Wilson, R. De Renzi, V. F. Mitrović, and S. Sanna, Microscopic nature of the charge-density wave in the kagome superconductor rbv3sb5, Physical Review Research 5, L012017 (2023).
  * Wang  _et al._ [2022] Y. Wang, T. Wu, Z. Li, K. Jiang, and J. Hu, The structure of kagome superconductors av3sb5 in the charge density wave states, arXiv preprint arXiv:2210.07585 (2022).
  * Luo  _et al._ [2022] J. Luo, Z. Zhao, Y. Z. Zhou, J. Yang, A. F. Fang, H. T. Yang, H. J. Gao, R. Zhou, and G.-q. Zheng, Possible star-of-david pattern charge density wave with additional modulation in the kagome superconductor csv3sb5, npj Quantum Materials 7, 30 (2022).
  * Kang  _et al._ [2022] M. Kang, S. Fang, J. Yoo, B. R. Ortiz, Y. Oey, S. H. Ryu, J. Kim, C. Jozwiak, A. Bostwick, E. Rotenberg, E. Kaxiras, J. Checkelsky, S. D. Wilson, J.-H. Park, and R. Comin, Microscopic structure of three-dimensional charge order in kagome superconductor av3sb5 and its tunability, arXiv preprint arXiv:2202.01902 (2022).
  * Hu  _et al._ [2022] Y. Hu, X. Wu, B. R. Ortiz, X. Han, N. C. Plumb, S. D. Wilson, A. P. Schnyder, M. Shi, _et al._ , Coexistence of trihexagonal and star-of-david pattern in the charge density wave of the kagome superconductor a v 3 sb 5, Physical Review B 106, L241106 (2022).
  * Sheldrick [2008] G. M. Sheldrick, A short history of shelx, Acta Crystallographica Section A: Foundations of Crystallography 64, 112 (2008).
  * Venderley  _et al._ [2022] J. Venderley, K. Mallayya, M. Matty, M. Krogstad, J. Ruff, G. Pleiss, V. Kishore, D. Mandrus, D. Phelan, L. Poudel, A. G. Wilson, K. Weinberger, P. Upreti, M. Norman, S. Rosenkranz, R. Osborn, and E.-A. Kim, Harnessing interpretable and unsupervised machine learning to address big data from modern X-ray diffraction, Proceedings of the National Academy of Sciences 119, e2109665119 (2022).
  * Islam  _et al._ [2012] Z. Islam, D. Capatina, J. P. Ruff, R. K. Das, E. Trakhtenberg, H. Nojiri, Y. Narumi, U. Welp, and P. C. Canfield, A single-solenoid pulsed-magnet system for single-crystal scattering studies, Review of Scientific Instruments 83, 035101 (2012).
  * Subedi [2022] A. Subedi, Hexagonal-to-base-centered-orthorhombic 4 q charge density wave order in kagome metals kv 3 sb 5, rbv 3 sb 5, and csv 3 sb 5, Physical Review Materials 6, 015001 (2022).
  * Ratcliff  _et al._ [2021] N. Ratcliff, L. Hallett, B. R. Ortiz, S. D. Wilson, and J. W. Harter, Coherent phonon spectroscopy and interlayer modulation of charge density wave order in the kagome metal csv 3 sb 5, Physical Review Materials 5, L111801 (2021).
  * [28] See Supplemental Information for further details.
  * Chen  _et al._ [2021] H. Chen, H. Yang, B. Hu, Z. Zhao, J. Yuan, Y. Xing, G. Qian, Z. Huang, G. Li, Y. Ye, _et al._ , Roton pair density wave in a strong-coupling kagome superconductor, Nature 599, 222 (2021).
  * Xiang  _et al._ [2021] Y. Xiang, Q. Li, Y. Li, W. Xie, H. Yang, Z. Wang, Y. Yao, and H.-H. Wen, Twofold symmetry of c-axis resistivity in topological kagome superconductor csv3sb5 with in-plane rotating magnetic field, Nature communications 12, 1 (2021).


11institutetext:  †Physics, Faculty of Science, University of East Anglia, NR4
7TJ, United Kingdom
⋆School of Mathematics, University of East Anglia, NR4 7TJ, United Kingdom

Vortices and turbulence Multicomponent condensates; spinor condensates
Tunneling, Josephson effect, Bose-Einstein condensates in periodic potentials,
solitons, vortices, and topological excitations

# Relaxation Dynamics of Half-Quantum Vortices in a Two-Dimensional Two-
Component Bose-Einstein Condensate

M. T. Wheeler† H. Salman⋆ M. O. Borgh†

###### Abstract

We study the relaxation dynamics of quantum turbulence in a two-component
Bose-Einstein condensate containing half-quantum vortices. We find a temporal
scaling regime for the number of vortices and the correlation lengths that at
early times is strongly dependent on the relative strength of the inter-
species interaction. At later times we find that the scaling becomes
universal, independent of the inter-species interaction, and approaches that
numerically observed in a scalar Bose-Einstein condensate.

###### pacs:

67.25.dk

###### pacs:

67.85.Fg

###### pacs:

03.75.Lm

##  1 Introduction

Since the realization of superfluidity, quantum turbulence (QT) has been
studied in systems ranging from superfluid liquid Helium [1, 2] to quasi-
particle condensates in solid-state systems [3]. Due to their unprecedented
experimental accessibility, QT in Bose-Einstein condensates (BECs) in dilute,
ultracold atomic gases have attracted considerable theoretical [4, 9, 5, 6, 7,
8] and experimental [10, 11, 12, 13, 14, 15] interest in both 2D and 3D
configurations. In a scalar BEC, the QT state is made up of a large number of
vortices with quantised circulation. The collective behaviour of the vortices
plays a key role in the hydrodynamics, recovering features of classical
turbulence that can exhibit the characteristic Kolmogorov power-law spectrum
[16].

In contrast to the scalar superfluids, multicomponent and spinor BECs are
described by multicomponent order parameters and allow for a wider range of
topological defects, which give rise to novel dynamics [17, 18, 19, 20].
Consequently, there has been increasing interest in the properties of QT and
non-equilibrium dynamics in such systems [23, 21, 24, 22, 25]. The simplest
non-scalar topological excitation appears in a two-component BEC, described by
two complex fields, as the appearance of a phase singularity in only one
component. When the atomic mass and mean density of the components are equal,
such vortices are often referred to as half-quantum vortices (HQVs), due to
their similarities with vortices carrying half a quantum of superfluid
circulation in superfluid 3He [26, 27] and spin-1 BECs [28, 29]. The study of
QT in BECs can be separated into two distinct categories: 1) forced turbulence
where a statistically stationary state is established; 2) decaying turbulence
where a non-equilibrium initial condition, typically involving vortices,
relaxes towards equilibrium. Here, we numerically investigate the spatial and
temporal properties of the relaxation dynamics of a non-equilibrium initial
state in a two-dimensional two-component system containing HQVs. Using a
pseudospin interpretation, we compute the temporal scaling of the correlation
functions associated with the spin- and mass-superfluid ordering. We relate
these to the vortex decay rate and analyse how this depends on the intra-
component interaction strength. We contrast our observations for this system
with similar simulations that have been performed for scalar BECs and reported
in [30, 31, 32].

##  2 The two-component BEC

We consider an untrapped two-component BEC described by the Gross-Pitaevskii
(GP) mean-field theory subject to periodic boundary conditions. The dynamics
of the condensate is described by the two coupled GP equations

$\displaystyle\begin{aligned} i\hbar\frac{\partial\psi_{1,2}}{\partial
t}&=\left(-\frac{\hbar^{2}}{2m_{1,2}}\nabla^{2}+g_{1,2}|\psi_{1,2}|^{2}+g_{12}|\psi_{2,1}|^{2}\right)\psi_{1,2}\end{aligned}$
(1)

where $\psi_{j}$ is the condensate wavefunction and $m_{j}$ $(j=1,2)$ is the
atomic mass for the $j$th component. The strength of inter- and intra-
component interactions are described by $g_{j}$ and $g_{12}$, respectively. We
consider a condensate where $m_{1}=m_{2}=m$, as is the case, e.g., when the
two components are different hyperfine states of the same atomic species, and
also assume $g_{1}=g_{2}=g$. The key parameter is then the ratio of intra- to
inter-species interactions

$\displaystyle\gamma=\frac{g_{12}}{g},$ (2)

which in experiment could be tuned using magnetic [33] or microwave-induced
[34] Feshbach resonances. Here we consider $0<\gamma<1$, such that all
interactions are repulsive, while keeping the condensate stable against
separation of the components.

The vortex states of the two-component BEC may be understood as follows: We
write the two-component wavefunction as the vector $(\psi_{1},\psi_{2})^{T}$.
Taking $\theta_{j}=\mathrm{Arg}(\psi_{j})$, this may be decomposed as

$\displaystyle\begin{pmatrix}\psi_{1}\\\
\psi_{2}\end{pmatrix}=\begin{pmatrix}|\psi_{1}|e^{i\theta_{1}}\\\
|\psi_{2}|e^{i\theta_{2}}\end{pmatrix}=e^{i\Theta}\begin{pmatrix}|\psi_{1}|e^{i\Phi}\\\
|\psi_{2}|e^{-i\Phi}\end{pmatrix},$ (3)

where

$\displaystyle\Theta=(\theta_{1}+\theta_{2})/2,\qquad\Phi=(\theta_{1}-\theta_{2})/2.$
(4)

Gradients in $\Phi$ can then be interpreted in terms of pseudospin currents,
while gradients in $\Theta$ may be associated with a total, superfluid mass
current.

Now consider a vortex state consisting of a phase singularity in $\psi_{1}$,
around which $\theta_{1}$ winds by $2\pi$ while $\theta_{2}$ remains
unchanged, such that

$\displaystyle\begin{pmatrix}\psi_{1}\\\
\psi_{2}\end{pmatrix}=\begin{pmatrix}|\psi_{1}|e^{i\phi}\\\
|\psi_{2}|\end{pmatrix}=e^{i\phi/2}\begin{pmatrix}|\psi_{1}|e^{i\phi/2}\\\
|\psi_{2}|e^{-i\phi/2}\end{pmatrix},$ (5)

where $\phi$ is the azimuthal angle around the vortex. The vortex is thus
equivalently described by a $\pi$ change in $\Theta$ (and a simultaneous $\pi$
change in $\Phi$) along a closed path encircling the vortex. Since $\Theta$
can be associated with a total mass current in the two components together,
these vortex states are often referred to as HQVs and we adopt this language
from here on. However, the two-component vortices are topologically distinct
from HQVs in the $A$ and polar phases of superfluid 3He [26, 27] and in the
uniaxial nematic phase of spin-1 BECs [28, 29].

A pseudospin picture also allows us to understand the size of HQV cores in
terms of an energetic hierarchy of length scales arising from the inter- and
intra-component interactions. These length scales are associated,
respectively, with variations of the total superfluid density and of the
density difference between the components. We thus define the density and spin
healing lengths as [35]

$\displaystyle\xi_{d}=\frac{\hbar}{\sqrt{2mgn_{0}}},\qquad\xi_{s}=\xi_{d}\left(\frac{1+\gamma}{1-\gamma}\right)^{1/2}\,,$
(6)

where $n_{0}$ is the number density of each component in a uniform system.
Since a HQV consists of a phase singularity in only one condensate component,
the remaining component is free to fill the vortex core. This can be
interpreted as a variation of the pseuodspin $z$-component, whose size is
determined by the spin healing length. When $\xi_{s}\gtrsim\xi_{d}$, the
vortex core can thus expand, lowering the total energy. Therefore, $\gamma$
directly determines the sizes of the vortex cores in the system. A similar
energetic hierarchy of length scales leads to dramatic defect-core
deformations in spinor BECs [36], including splitting of singly quantised
vortices into HQVs [37, 29].

##  3 Numerical method

To study the dynamics of vortices in a turbulent regime we numerically evolve
the time-dependent two-component Gross-Pitaevskii equations using a split-step
algorithm [38]. We write eq. (1) in terms of the dimensionless variables:
$\tilde{\vect{r}}=\vect{r}/a_{s}$, $\tilde{t}=t/\tau$,
$\tilde{g}=2mg/\hbar^{2}$ and $\tilde{\psi_{j}}=a_{s}\psi_{j}$, where $a_{s}$
is the lattice spacing and $\tau=2ma_{s}^{2}/\hbar$ is the lattice time. The
resulting equations then become

$\displaystyle i\frac{\partial\psi_{1,2}}{\partial
t}=\left(-\nabla^{2}+g|\psi_{1,2}|^{2}+\gamma
g|\psi_{2,1}|^{2}\right)\psi_{1,2},$ (7)

where we have dropped the tildes for notational convenience. Our simulations
were performed on a periodic domain of non-dimensional area $L^{2}$ with side
length $L=N_{s}$ where $N_{s}^{2}$ is the number of grid points. We solve eq.
(7) on a grid of $1024^{2}$ points with $a_{s}=1$. Motivated by similar work
in a scalar BEC [32], we take $N=3.2\times 10^{9}$ atoms per component and
dimensionless $g=L^{2}/4N$. The non-dimensional density healing length is thus
fixed at $N_{s}/(gN)^{1/2}=2$. We now explore the role of the inter-component
interaction by varying $\gamma$ within the range $0<\gamma<1$.

The initial condition for the GP evolution is constructed as a grid of vortex
positions containing $48^{2}$ vortices in each component with the grids of
each component offset in both $x$ and $y$ to avoid overlapping positions. We
then add a small, random displacement to each position to create an irregular
distribution of vortices. This facilitates the development of an initially
chaotic and subsequently turbulent vortex evolution during the relaxation
dynamics. The phase of each component is subsequently constructed as an
alternating $2\pi$ winding around each vortex position using the method
described in ref. [5] that also accounts for the periodic boundary conditions.
An initial short period of imaginary-time propagation, keeping the phase
profile fixed, allows the vortex cores to form. The resulting HQVs consist of
a density depletion in one component at the position of the phase singularity,
which is then filled with atoms of the other component, as illustrated in fig.
1. From this initial state, the system is evolved according to eq. (7). HQVs
with opposite circulation but with the phase singularity in the same component
may annihilate which leads to a decay of the total vortex number.

Figure 1: Density $|\psi_{1}|^{2}$ in a $100\xi_{d}\times 100\xi_{d}$
subdomain of the initial state after a short imaginary-time evolution. We can
identify the vortices in this component by the density depletion (blue).
Density peaks (red) form in $\psi_{1}$ at the positions of vortices in the
$\psi_{2}$ component.

##  4 Results

Figure 2: Density (a)–(c) and pseudo-vorticity (d)–(f) of the $\psi_{1}$
component in a $256\xi_{d}\times 256\xi_{d}$ subregion at time $t=2.5\times
10^{4}\xi_{d}^{2}$, for $\gamma=0.1$ (left), $\gamma=0.6$ (middle) and
$\gamma=0.8$ (right). Vortices in $\psi_{1}$ appear as a density depletion.
For $\gamma\geq 0.6$, bright density peaks show where $\psi_{1}$ atoms fill
the cores of HQVs with the phase singularity in $\psi_{2}$. Vortices with
positive (blue) and negative (red) circulation are identifiable in the pseudo-
vorticity field.

We first investigate the effect of $\gamma$ on the relaxation dynamics of
HQVs. Fig. 2(a)–(c) shows the density of the $\psi_{1}$ component for
$\gamma=0.1,0.6,0.8$. HQVs with a phase singularity in this component are
readily apparent by the corresponding density depletion, and have a core size
that grows with increasing $\gamma$. For $\gamma\geq 0.6$, high density peaks
also become noticeable and correspond to the positions of HQVs with phase
singularity in $\psi_{2}$. This can be understood from the healing lengths,
eq. (6). For small $\gamma$, $\xi_{s}\sim\xi_{d}$. As $\gamma$ increases, the
spin healing length also increases. Consequently, the cores of the HQVs fill
with atoms from the other component as the resulting lowering of the kinetic
energy offsets the cost in interaction energy. This causes the vortex cores to
expand to a size similar to the spin healing length, as borne out by our
simulations.

To track the vortex positions, we evaluate the pseudo-vorticity [39, 40]

$\displaystyle\vect{\omega}_{\mathrm{p_{j}}}=\frac{1}{2}\nabla\times(n\vect{v})_{j}\,,$
(8)

where

$\displaystyle(n\vect{v})_{j}=\frac{1}{i}\left[\psi_{j}^{*}(\nabla\psi_{j})-(\nabla\psi_{j}^{*})\psi_{j}\right],$
(9)

is the mass current of component $j=1,2$. The pseudo-vorticity remains regular
and non-zero within the core of each vortex, and relaxes to zero away from the
vortex singularity (at length scales exceeding the spin-healing length,
$\xi_{s}$), as shown in fig. 2(d)–(f). The sign of the pseudo-vorticity also
determines the charge of the vortex. The pseudo-vorticity shows the vortex
positions particularly sharply for small $\gamma$, where the vortex cores are
small.

We now investigate the spatial properties of our turbulent system. We split
the kinetic energy, $E_{\mathrm{kin}}=E^{v}+E^{q}$, into classical ($E^{v}$),
and quantum-pressure ($E^{q}$) contributions. These are given by

$\displaystyle E^{v}$
$\displaystyle=\frac{1}{4}\int\upd^{2}\mathbf{x}\,\left(|\sqrt{n}_{1}\vect{v}_{1}|^{2}+|\sqrt{n}_{2}\vect{v}_{2}|^{2}\right)\,,$
(10) $\displaystyle E^{q}$
$\displaystyle=\int\upd^{2}\mathbf{x}\,\left(|\nabla\sqrt{n}_{1}|^{2}+|\nabla\sqrt{n}_{2}|^{2}\right)\,,$
(11)

where $n_{j}=|\psi_{j}|^{2}$ for $j=1,2$.

The energy spectra for these contributions involve the Fourier transforms of
the generalised velocities for the incompressible ($i$), compressible ($c$),
and quantum pressure ($q$) parts [21], defined as

$\displaystyle\begin{gathered}\vect{w}^{i,c}=\sqrt{n_{1}}\vect{v}_{1}^{i,c}+\sqrt{n_{2}}\vect{v}_{2}^{i,c}\,,\\\
\vect{w}^{q}=2\left(\nabla\sqrt{n_{1}}+\nabla\sqrt{n_{2}}\right)\,.\end{gathered}$
(14)

The incompressible and compressible components of the velocity field are
recovered from a Helmholtz decomposition into a divergence free,
incompressible part $\nabla\cdot\vect{v}^{i}=0$, and an irrotational,
compressible part $\nabla\times\vect{v}^{c}=\vect{0}$. The kinetic energy
spectrum can then be calculated by integrating the corresponding Fourier
transforms over the full $k$-space angle

$\displaystyle E^{\delta}(k)$
$\displaystyle=\frac{1}{4}\int_{0}^{2\pi}\upd\Omega_{\vect{k}}\,|\tilde{\vect{w}}^{\delta}(\vect{k})|^{2},\enskip(\delta=i,c,q),$
(15)

for wave number $k=|\vect{k}|$. The total kinetic energy is given by
integrating over all $k$ and summing over the different contributions:
$E_{\textrm{kin}}=\sum_{\delta}\int\upd kE^{\delta}(k)$ for $\delta=(i,c,q)$.
The occupation numbers corresponding to the different energy contributions are
then

$\displaystyle n^{\delta}(k)=k^{-2}E^{\delta}(k),\enskip(\delta=i,c,q).$ (16)

Fig. 3 shows the occupation number for each energy contribution along with the
total occupation number $n(k)$ for the case of $\gamma=0.6$ at a time
$t=2\times 10^{5}\tau$. The total single-particle spectrum obeys the predicted
scaling $n(k)\sim k^{-4}$ in the infrared (IR) region and $n(k)\sim k^{-2}$ in
the ultraviolet (UV) seen for some turbulent, 2D, scalar BEC systems [31, 32,
42]. Decomposing the kinetic energy into its constituent parts, we see that
the incompressible contribution dominates in the IR and is responsible for the
change in scaling to $k^{-4}$ in this region. This incompressible contribution
is associated with the vortices in the system [24]. At large $k$, the spectrum
is dominated by the compressible and quantum pressure contributions exhibiting
the weak-wave-turbulence scaling $k^{-2}$. This scaling of the energy is
qualitatively insensitive to variations in $\gamma$.

Figure 3: Occupation numbers for different fractions of the system with
$\gamma=0.6$ at $t=5\times 10^{4}\xi_{d}^{2}$: single particle spectrum for
quantum pressure (purple diamonds), incompressible (red diamonds) and
compressible (blue diamonds) contributions. The total occupation number (black
diamonds) for the single particle spectrum is obtained by summing the
corresponding fractions from each condensate component. The single particle
spectrum obeys a $k^{-2}$ scaling (dotted line) in the ultra-violet and a
$k^{-4}$ scaling (dashed line) in the infrared regions.

Next, we consider the time-dependent properties of the turbulent dynamics. For
this purpose, we will concentrate on the correlation functions for the spin
and mass parts of the pseudospinor order parameter. For a homogeneous
turbulent system these are defined, respectively, as [43]

$\displaystyle G_{\Phi}(\vect{r},t)$
$\displaystyle=\frac{2}{n^{2}}\mathrm{Tr}\left[\langle\mathsf{Q}(\vect{0})\mathsf{Q}(\vect{r})\rangle\right]\,,$
(17) $\displaystyle G_{\Theta}(\vect{r},t)$
$\displaystyle=\frac{1}{n^{2}}\langle\alpha^{*}(\vect{0})\alpha(\vect{r})\rangle\,,$
(18)

where $\langle\cdot\rangle$ denotes ensemble averaging. Here, the matrix

$\displaystyle\mathsf{Q}=\begin{pmatrix}Q_{xx}&Q_{xy}\\\
Q_{xy}&-Q_{xx}\end{pmatrix},$ (19)

where $Q_{xx}=\mathrm{Re}\\{\psi_{1}^{*}\psi_{2}\\}$ and
$Q_{xy}=\mathrm{Im}\\{\psi_{1}^{*}\psi_{2}\\}$, is associated with spin
ordering in the system, while $\alpha=-2\psi_{1}\psi_{2}$ is an alignment
parameter. Exploiting the fact that our turbulent system is homogeneous, we
can replace ensemble averages with spatial averages. The spin correlation
function is then equivalently defined as [43]

$\displaystyle
G_{\Phi}(r,t)=\int\upd\Omega_{r}\int\frac{\upd^{2}\vect{x}^{\prime}}{L^{2}}\frac{2\mathrm{Tr}\left[\mathsf{Q}(\vect{\vect{x}^{\prime}})\mathsf{Q}(\vect{x}^{\prime}+\vect{r})\right]}{n^{2}},$
(20)

where $\int\upd\Omega_{r}$ denotes angular integration. We perform the same
averaging for the superfluid correlation function.

Figure 4: (a): Spin-correlation function as a function of time for
$\gamma=0.6$. The spin order decays more slowly as time increases, indicating
domain growth within the system. Inset: collapse of the spin correlation
function when scaled by the spin correlation length. (b): Correlation lengths
corresponding to the spin and superfluid correlation functions as a function
of time for $\gamma=0.3,0.6,0.8$. Larger $\gamma$ give a faster initial
growth, with a universal scaling appearing for $t\gtrsim 2.5\times
10^{3}\xi_{d}^{2}$. The $t^{1/5}$ scaling predicted from the scalar BEC is
indicated for comparison. Figure 5: Mean vortex distance in a scalar BEC for
three different initial conditions using the same parameters as in ref. [32].
The $t^{1/5}$ early-time as well as the $t^{1/2}$ late-time scaling regimes
are recovered.

In fig. 4(a) we plot the results for the spin correlation function at
different times for $\gamma=0.6$. As the time increases, the correlation
function decays over a larger distance, indicating the emergence of long-range
order within the system. We verify the same behaviour for the mass-correlation
function. From these correlation functions we obtain the correlation length,
$L_{\delta}(t)$ for $\delta=\\{\Phi,\Theta\\}$, which we take as the distance
at which the corresponding correlation function decays to a quarter of its
value at $r=0$: $G_{\delta}(L_{\delta},t)=\frac{1}{4}G_{\delta}(0,t)$. The
correlation functions are said to exhibit dynamical scaling when their form at
different times remains self similar. This means that they collapse to a
universal, time-independent function when scaled by the correlation lengths,
i.e. $H_{\delta}(r)=G_{\delta}(r/L_{\delta}(t),t)$. The inset in fig. 4 shows
this collapse of the spin correlation function in our system, indicating that
$G_{\Phi}(r,t)$ does indeed exhibit dynamical scaling. We again verify the
same behaviour for $G_{\Theta}(r,t)$.

Fig. 4(b) shows both correlation lengths $L_{\Phi,\Theta}(t)$ as a function of
time for $\gamma=0.3$, $\gamma=0.6$ and $\gamma=0.8$. After the initial
evolution the temporal scaling of the correlation lengths becomes universal
for all values of $\gamma$. However, the effect of $\gamma$ is apparent in the
early time evolution where a larger $\gamma$ leads to a faster growth of the
correlation lengths. This is indicative of a difference in the decay rate of
the vortices in the early-time dynamics.

Figure 6: Total vortex number in both components (red) as a function of time
for $\gamma=0.3,0.7,0.9$. Larger $\gamma$ leads to a steeper decay due to the
rapid annihilation of opposite-signed vortices in the same component. Overlaid
for comparison is twice the vortex number (black) from a corresponding scalar-
BEC simulation (equivalent to $\gamma=0$) with the same initial vortex
distribution, atom number, and interaction strength $g$ as $\psi_{1}$. Figure
7: Exponent $z$ as a function of $\gamma$ in the interval $2.5\times
10^{2}\xi_{d}^{2}<t<2.5\times 10^{3}\xi_{d}^{2}$. A rapid decrease of the
exponent arises for $\gamma\gtrsim 0.6$.

We can investigate this behaviour by considering the total number of vortices
in the system as a function of time. We extract the mean distance between
vortices as $\ell_{d}=1/\sqrt{N_{\mathrm{vort}}}$, where $N_{\mathrm{vort}}$
is the total number of vortices in the system. As a point of reference, in a
scalar BEC initially containing a large number of vortices, $\ell_{d}\sim
t^{\beta}$ [32], where $\beta$ characterises the vortex annihilation rate. In
particular, after some (possibly short) period of evolution, a $\beta=1/5$
scaling is observed. For comparison, we have indicated this theoretically
expected scaling in fig. 4(b) for our two-component BEC. At late times, a
$\beta=1/2$ scaling appears in the scalar BEC, whose onset is delayed if the
initial vortex distribution is highly clustered [32]. In fig. 5 we reproduce
this late-time scaling using the parameters of ref. [32] for an initial grid
of elementary vortices analogous to our two-component initial state, as well
as for a random vortex distribution with and without noise added to the energy
spectrum. In all cases we recover both the $t^{1/5}$ scaling after initial
evolution and the $t^{1/2}$ late-time scaling, indicating that this behavior
is robust and qualitatively insensitive to details of the initial condition.

Motivated by this previous work, we perform a similar analysis to establish
how these results extend to a two-component system with HQVs and how the
vortex annihilation rate depends on $\gamma$. We focus on the early vortex
evolution, where fig. 4(b) suggests that the $\gamma$-dependence is
significant. Fig. 6 shows $N_{\mathrm{vort}}$ as a function of time for three
different values of $\gamma$. For $\gamma=0.7$ and $0.9$, a new scaling regime
emerges at early times ($2.5\times 10^{2}\xi_{d}^{2}\lesssim t\lesssim
2.5\times 10^{3}\xi_{d}^{2}$), where $N_{\mathrm{vort}}(t)$ decays as $t^{-1}$
($\gamma=0.7$) and $t^{-1.5}$ ($\gamma=0.9$). For $t\gtrsim 2.5\times
10^{3}\xi_{d}^{2}$, $N_{\mathrm{vort}}(t)$ approaches a universal $t^{-2/5}$
scaling corresponding to $\ell_{d}\sim t^{1/5}$, similar to the scalar BEC
also shown. These results imply a better agreement with the theoretical
$t^{1/5}$ scaling than indicated from the correlation lengths [fig. 4(b)].
This suggests that although their growth is driven by vortex annihilation, the
length scales $L_{\Phi,\Theta}(t)$ are not fully equivalent to $\ell_{d}(t)$.
The region of interest in fig. 6 only extends up to $t=5\times
10^{4}\xi_{d}^{2}$ and we therefore expect a universal transition to
$\ell_{d}\sim t^{1/2}$ at times extending beyond the time interval of our
simulations.

Previous work has demonstrated that, for a sufficiently high $\gamma\gtrsim
0.6$, a dipole consisting of HQVs with opposite phase winding in the same
component will shrink in size as the vortices move toward one another and
annihilate [17]. We therefore attribute the different scaling regime at early
times, when the mean inter-vortex separation is small, to this behaviour. This
is further supported by the fact that we do not see such scaling for
$\gamma\lesssim 0.6$, where such rapid annihilation rate is not prevalent.
Within that range of values for $\gamma$, the vortex dynamics begins to
recover the behavior observed in a scalar BEC.

We can model the vortex decay rate by a kinetic-like equation of the form

$\displaystyle\partial_{t}N_{\mathrm{vort}}\sim N_{\mathrm{vort}}^{\eta},$
(21)

where $\eta>1$. The dependence of $N_{\mathrm{vort}}$ on the right-hand side
of the equation indicates that the decay rate is a function of the number of
vortices that are involved in facilitating the annihilation. Using this simple
model, we can derive temporal scaling of the total vortex number as [44]

$\displaystyle N_{\mathrm{vort}}\sim t^{-2/z},$ (22)

where $z=-2(1-\eta)$. We note that an exponent of $z=2$ corresponds to a two-
body collision process whereas $z=5$ corresponds to three-body collisions
[32]. In fig. 7, we quantify the $\gamma$ dependence of the early-time scaling
by considering the exponent $z$ in the region $2.5\times
10^{2}\xi_{d}^{2}<t<2.5\times 10^{3}\xi_{d}^{2}$. We see a rapid decrease of
the exponent after $\gamma>0.6$, when the more rapid annihilation becomes
prevalent. The observed decrease in the value of $z$ with $\gamma$ in our
simulations signals an additional interaction effect not present in the scalar
system.

##  5 Conclusions

We have investigated spatial and temporal aspects of a decaying turbulent two-
component BEC containing HQVs. The occupation-number spectrum is found to show
a scaling behaviour consistent with similar results for a scalar BEC across a
wide range of values of the inter-component interaction strength.

However, we find that a new interaction-dependent scaling regime appears in
the temporal properties of the mass- and spin-correlation functions, as well
as the mean inter-vortex separation. For large values of the relative inter-
component interaction strength, $\gamma\gtrsim 0.6$, these exhibit a
$\gamma$-dependent scaling that is markedly different from the universal
behavior, which conforms to that of a scalar BEC at a similar stage of time
evolution. Modelling the total vortex number using a simple kinetic equation,
we have found that this early-time decay rate for high $\gamma$ cannot be
explained by simple two- or three-body collisions. The observed enhanced
vortex decay rate at early times for large $\gamma$ may be due to the role
played by an additional inter-vortex force that arises between vortices in the
same component. The results suggest that this force is short-range and appears
in addition to the well-known $1/R$ inter-vortex force. This latter force
appears to dominate once the vortex density drops significantly following the
rapid vortex annihilations occurring at early times.

###### Acknowledgements.

The results presented were obtained using the High Performance Computing
Cluster supported by the Research and Specialist Computing Support service at
the University of East Anglia.

## References

  * [1] Barenghi C.F., Skrbek L. Sreenivasan K.R. PNAS11120144647-4652
  * [2] Walmsley P., Zmeev D., Pakpour F. Golov A. PNAS11120144691-4698
  * [3] Kreil A.J.E., Bozhko D.A., Musiienko-Shmarova H.Yu., Vasyuchka V.I., L’vov V.S., Pomyalov A., Hillebrands B. Serga A.A. Phys. Rev. Lett.1212018077203.
  * [4] Kobayashi M. Tsubota M. Phys. Rev. A762007045603.
  * [5] Reeves M., Billam T., Anderson B. Bradley A. Phys. Rev. Lett.1102013104501.
  * [6] Billam T., Reeves M., Anderson B. Bradley A. Phys. Rev. Lett.1122014145301.
  * [7] Simula T., Davis M.J. Helmerson K. Phys. Rev. Lett.1132014165302
  * [8] Baggaley A.W. Barenghi C.F. Phys. Rev. A972018033601
  * [9] Numasato R., Tsubota, M. L’vov, V.S. Phys. Rev. A812010063630
  * [10] Henn E., Seman J., Roati G., Magalhães K. Bagnato V. Phys. Rev. Lett.1032009045301.
  * [11] Kwon W.J., Moon G., Choi J., Seo S. and Shin Y. Phys. Rev. A902014063627
  * [12] Seo S.W., Ko B., Kim J.H. Shin Y. Sci. Rep.720174587
  * [13] Navon N., Eigen C., Zhang J., Lopes R., Gaunt A.L., Fujimoto K., Tsubota M., Smith R.P. Hadzibabic Z. Science3662019382
  * [14] Guillaume G., Reeves M.T., Yu X., Bradley A.S., Baker M.A., Bell T.A., Rubinsztein-Dunlop H., Davis M.J. Neely T.W. Science36420191264.
  * [15] Johnstone S.P., Groszek A.J., Starkey P.T., Billington C.J., Simula T.P. Helmerson, K. Science36420191267.
  * [16] Kobayashi M. Tsubota M. J. Phys. Soc. Jpn.7420053248.
  * [17] Kasamatsu K., Eto M. Nitta M. Phys. Rev. A932016013615.
  * [18] Weiss L.S., Borgh M.O., Blinova A, Ollikainen T, Möttönen M, Ruostekoski J, Hall D.S. Nat. Commun.1020194772.
  * [19] Kobayashi M, Kawaguchi Y, Nitta M Ueda M. Phys. Rev. Lett.1032009115301.
  * [20] Kasamatsu K., Tsubota M. Ueda M. Int. J. Mod. Phys.1920051835.
  * [21] Schmied C.-M, Gasenzer T. Blakie P. Phys. Rev. A1002019033603.
  * [22] Prüfer M., Kunkel P., Strobel H., Lannig S., Linnemann D., Schmied C.-M., Berges J., Gasenzer T. Oberthaler M.K. Nature5632018217.
  * [23] Salman H. Berloff N.G. Physica D23820091482.
  * [24] Karl M., Nowak B. Gasenzer T. Phys. Rev. A882013063615.
  * [25] Hofmann J., Natu S. S. Sarma S. D. Phys. Rev. Lett.1132014095702.
  * [26] Volovik G.E. Mineev V.P. JETP Lett.241976561.
  * [27] Autti S, Dmitriev V.V., Mäkinen J.T., Soldatov A.A., Volovik G.E., Yudin A.N., Zavjalov, V.V. Eltsov V.B. Phys. Rev. Lett.1172016255301.
  * [28] Leonhardt U. Volovik G.E. JETP Lett.72200046.
  * [29] Seo S.W., Kang S., Kwon W.J., Shin Y.-i. Phys. Rev. Lett.1152015015301.
  * [30] Schole J., Nowak B. Gasenzer T. Phys. Rev. A862012013624
  * [31] Nowak B., Schole J., Sexty D. Gasenzer T. Phys. Rev. A852012043627.
  * [32] Karl M. Gasenzer T. New J. Phys.192017093014.
  * [33] Inouye S., Andrews M.R., Stenger J., Miesner H.-J., Stamper-Kurn D.M. Ketterle W. Nature3921998151.
  * [34] Papoular D.J., Shlyapnikov G.V. Dalibard J. Phys. Rev. A812010041603.
  * [35] Eto M., Kasamatsu K., Nitta M., Takeuchi H. Tsubota M. Phys. Rev. A832011063603.
  * [36] Ruostekoski J. Anglin, J.R. Phys. Rev. Lett.912003190402.
  * [37] Lovegrove J., Borgh M.O. Ruostekoski J. Phys. Rev. A862012013613.
  * [38] Javanainen J. Ruostekoski J. J. Phys. A: Math. Gen.392006L179.
  * [39] Villois A., Krstulovic G., Proment D. Salman H. J. Phys. A: Math. Theor.492016415502.
  * [40] Rorai C., Skipper J., Kerr R.M. Sreenivasan K.R. J. Fluid Mech.8082016641.
  * [41] Kawaguchi Y. Ueda M. Phys. Rep.5202012253.
  * [42] Salman H. Maestrini D. Phys. Rev. A942016043642.
  * [43] Symes L. Blakie P. Phys. Rev. A962017013602.
  * [44] Bray A.J. Adv. Phys.431994357-459.