Observation of wall-vortex composite defects in a spinor Bose-Einstein condensate

TL;DR

This paper reports the experimental observation of wall-vortex composite defects in a spinor Bose-Einstein condensate, revealing their formation, dynamics, and nucleation mechanisms during phase transitions.

Contribution

It provides the first direct observation of spin domain walls bounded by half-quantum vortices in a spin-1 BEC, demonstrating their formation and evolution during a controlled phase transition.

Findings

Spin domain walls bounded by half-quantum vortices observed.

Wall-vortex composite defects dynamically split due to snake instability.

Singly charged quantum vortices form during strong quenches.

Abstract

We report the observation of spin domain walls bounded by half-quantum vortices (HQVs) in a spin-1 Bose-Einstein condensate with antiferromagnetic interactions. A spinor condensate is initially prepared in the easy-plane polar phase, and then, suddenly quenched into the easy-axis polar phase. Domain walls are created via the spontaneous $\mathbb{Z}_2$ symmetry breaking in the phase transition and the walls dynamically split into composite defects due to snake instability. The end points of the defects are identified as HQVs for the polar order parameter and the mass supercurrent in their proximity is demonstrated using Bragg scattering. In a strong quench regime, we observe that singly charged quantum vortices are formed with the relaxation of free wall-vortex composite defects. Our results demonstrate a nucleation mechanism for composite defects via phase transition dynamics.

Equations

$$\mathbf{\psi }=\left[\begin{array}{c}{\psi }_{+1}\\ {\psi }_0\\ {\psi }_{-1}\end{array}\right]=\sqrt{n}{e}^{i\varphi }\left[\begin{array}{c}-\frac{d_x-i{d}_y}{\sqrt{2}}\\ d_z\\ \frac{(d_x+i{d}_y)}{\sqrt{2}}\end{array}\right],$$

$$\mathbf{\psi }$$

$$\mathcal{E}(\mathbf{\psi })$$

$$\mathbf{\psi }={\psi }_{\mathrm{T}}\hat{\mathbf{d}}$$

$$\mathcal{E}$$

$${\xi }_q=\frac{\hslash }{\sqrt{2mq}}.$$

$$\mathcal{E}\to \sum _{\nu =x,y,z}\frac{{\hslash }^2}{2m}{\left|{\partial }_{\nu }\mathbf{\psi }\right|}^2-\mu {\left|\mathbf{\psi }\right|}^2+\frac{c_0}{2}{\left|\mathbf{\psi }\right|}^4,$$

$$\mathbf{\psi }(\rho ,\theta )=\sqrt{n(\rho )}{e}^{i\frac{\theta }{2}}\hat{\mathbf{d}}(\theta ),\text{leavevmode}\text{leavevmode}\hat{\mathbf{d}}(\theta )={\left[\cos \frac{\theta }{2},0,\sin \frac{\theta }{2}\right]}^{\mathrm{T}}$$

Full text

Topological defects in a continuous ordered system are a splendid manifestation of symmetry breaking, with their fundamental types, such as walls, strings, and monopoles, inevitably determined by the topology of the order parameter space. However, if there is a hierarchy of energy (length) scales with different symmetries, composite defects, such as domain walls bounded by strings and strings terminated by monopoles may exist in the system Kibble_book . In cosmology, it has been noted that such composite defects can be nucleated through successive phase transitions with different symmetry breaking in grand unification theories; furthermore, composite defect formation has been proposed as a possible mechanism for galaxy formation Kibble_book ; Vilenkin_phyrep85 and baryogenesis Menahem_nuclphys92 in the early Universe.

Spinful superfluid systems with multiple symmetry breaking provide an experimental platform for studying the physics of composite defects and, thus, to examine the cosmological scenario. In superfluid 3He-B, it has been observed that a spin-mass vortex, on which a planar soliton terminates, can survive after phase transitions by being pinned on the vortex lattice Kondo_prl92 ; Eltsov_prl00 or nafen Makinen_arxiv18 . Composite defects have also been theoretically studied in the atomic Bose-Einstein condensate (BEC) system. Vortex confinement with a domain wall was predicted to occur in a two-component BEC under coherent intercomponent coupling Son_pra02 . In particular, for a spin-1 Bose gas with antiferromagnetic interactions, half-quantum vortices (HQVs) joined by a spin domain wall were anticipated to be responsible for the emergence of an exotic 2D superfluid phase with spin-singlet pair correlations Mukerjee_prl06 ; James_prl11 .

In this Letter, we report the experimental observation of wall-vortex composite defects in a quasi-2D antiferromagnetic spin-1 BEC. The composite defects are nucleated via a two-step instability mechanism in quantum quench dynamics from the easy-plane polar (EPP) phase into the easy-axis polar (EAP) phase. In the first step, spontaneous $\mathbb{Z}_{2}$ symmetry breaking causes domain wall formation, the core of which is occupied by the EPP phase. In the second step, the snake instability splits the domain walls into segments, with each segment forming a composite defect, which is a domain wall terminating on a HQV Seo_prl15 . The mass supercurrent in proximity to the wall end point is demonstrated using Bragg scattering Seo_scirep17 . We also observe that singly charged quantum vortices (QVs) can be formed by the relaxation of free composite defects. Our results directly demonstrate the existence of wall-vortex composite defects and their nucleation mechanism via phase transition dynamics in a spinful superfluid system.

The experiment is performed with a BEC of 23Na atoms in the $F=1$ hyperfine state having an antiferromagnetic spin interaction coefficient $c_{2}>0$ Stenger98 . The ground state of a spin-1 antiferromagnetic BEC is a polar state with $\langle\mathbf{F}\rangle=0$ Ho98 ; Ohmi98 , where $\mathbf{F}=(F_{x},F_{y},F_{z})$ is the spin operator of the particle. The order parameter of the BEC is parametrized with the superfluid phase $\phi$ and a real unit vector $\hat{\bm{d}}=(d_{x},d_{y},d_{z})$ for the spin director, and is expressed as

[TABLE]

where $\psi_{m_{z}=0,\pm 1}$ is the condensate wave function of the $|m_{z}\rangle$ Zeeman component and $n$ is the particle density. In the presence of an external magnetic field, e.g., along the $z$ axis, uniaxial spin anisotropy is imposed by the quadratic Zeeman energy $E_{z}=q(1-d_{z}^{2})$ and the ground state of the system is the EAP state with $\hat{\bm{d}}=\pm\hat{\bm{z}}$ for $q>0$ and the EPP state with $\hat{\bm{d}}\perp\hat{\bm{z}}$ for $q<0$.

As a means of creating wall-vortex composite defects, we employ the quantum quench dynamics from the EPP phase to the EAP phase via a sudden change of spin anisotropy, which can be implemented by dynamically controlling the $q$ value Gerbier_pra06 ; Zhao_pra14 . Because of the positional difference between the two phases in the order parameter space, the quench dynamics involves spontaneous $\mathbb{Z}_{2}$ symmetry breaking, as \hat{\bm{d}}$$\perp$$\hat{\bm{z}}\rightarrow\hat{\bm{d}}=$$\pm\hat{\bm{z}}. For q$$>0, the initial EPP state is dynamically unstable so that spin fluctuations will be exponentially amplified via the spin exchange process of |$$+1\rangle|$$-$$1\rangle\rightarrow|0\rangle|0\rangle Kawaguchi_phyrep12 . The microscopic origin of the $\mathbb{Z}_{2}$ symmetry breaking arises from the two equivalent choices for the phase of the $|0\rangle$ component. A rapid quench can give rise to a complex network of domain walls in a uniform system according to the Kibble-Zurek mechanism Kibble80 ; Zurek85 .

The spatial structure of a domain wall is described in Fig. 1(a), which is formed at the interface between two domains with opposite spin directions. Here $\hat{\bm{d}}$ is denoted by a double-head arrow and the superfluid phase is indicated by the color of each arrow head, reflecting the discrete symmetry of the order parameter under the operation of $(\phi,\hat{\bm{d}})\rightarrow(\phi+\pi,-\hat{\bm{d}})$ Zhou_prl01 . In the wall region, $\hat{\bm{d}}$ continuously flips to the opposite direction and the |$$\pm 1$$\rangle components are present, sandwiched by the $|0\rangle$ component. The wall thickness is determined by the competition between the quadratic Zeeman energy and the gradient energy associated with the vector field $\hat{\bm{d}}(\bm{r})$, giving a characteristic length scale of $\xi_{q}=\hbar/\sqrt{2mq}$ for $q\ll\mu$ with the particle mass $m$ and the chemical potential $\mu$ supple . An interesting observation is that the two domains separated by the wall comprise only the $|0\rangle$ component, which means that they can be continuously connected to each other by varying $\phi$ without flipping $\hat{\bm{d}}$, thus, allowing spatial termination of the domain wall as shown in Fig. 1(b). In this case, the wall end point exhibits a superfluid phase winding of $\pi$, forming a HQV Seo_prl15 . This is the wall-vortex composite defect expected in the EAP phase. The spatial structure of the composite defect is analogous to that of the spin-mass vortex, also referred to as a $\theta$ soliton, in superfluid 3He Kondo_prl92 ; Makinen_arxiv18 ; supple .

We prepare a condensate containing $N_{c}\approx 8.0\times 10^{6}$ atoms in the |F$$=$$1,m_{F}$$=$$0\rangle hyperfine spin state in an optical dipole trap with trapping frequencies of (\omega_{x},\omega_{y},\omega_{z})=2\pi\times(3.8,5.5,402)\leavevmode\nobreak\Hz. The Thomas-Fermi radii for the trapped condensate are $(R_{x},R_{y},R_{z})\approx(230,160,2.2)\leavevmode\nobreak\ \mu$m. The external magnetic field is $B_{z}=33$ mG, giving $q/h=0.3$ Hz, and the field gradient is controlled to be less than 0.1 mG/cm Kang_pra17 . The EPP-to-EAP quench dynamics is initiated by rotating $\hat{\bm{d}}$ from $\hat{\bm{z}}$ to the $xy$ plane by applying a short rf pulse and then suddenly changing the $q$ value to a target value $q_{f}>0$ using a microwave dressing technique Seo_prl15 .

The postquench evolution of the BEC is examined by measuring the spatial density distributions of the three spin components at a variable hold time $t$ with taking an absorption image after Stern-Gerlach (SG) spin separation for 24 ms time-of-flight Kang_pra17 . The spin healing length is $\xi_{s}=\hbar/\sqrt{2mc_{2}n_{0}}\approx 4.0\leavevmode\nobreak\ \mu$m for the peak atom density $n_{0}$, and our highly oblate sample with $R_{z}<\xi_{s}$ constitutes a quasi-2D system for spin dynamics. In our experiment, $q_{f}$, which represents the initial excitation energy per particle with respect to the ground state, is much smaller than \mu\approx h\times 880\leavevmode\nobreak\Hz, so incurrence of density perturbations is energetically improbable. Note that $q_{f}\ll\mu$ sets a clear hierarchy of energy scales in the system supple .

Figure 2 shows several images of the quenched condensate after various hold times for $q_{f}/h=1.0$ Hz. In the early stage of the quench dynamics, a few line defects are clearly observed to appear across the condensate [Fig. 2(a)]. The $|0\rangle$ component shows density-depleted trenches and the trench regions are filled by both of the |$$\pm$$1\rangle components, consistent with the spin distribution for the domain wall described in Fig. 1(a). The high visibility of the trench in the $|0\rangle$ component after such a long time-of-flight reflects the nature of a topological soliton. At later $t>0.2$ s, the line defects are typically observed to end in the middle of the condensate [Figs. 2(b) and 2(c)]. The profile of the total condensate density was confirmed to remain unperturbed by imaging without SG spin separation. Such a smooth wall termination is the key characteristic of the wall-vortex composite defects. As $t$ increases, the end point seems to recede toward the boundary of the condensate with the decrease of the domain wall length [Fig. 2(d)], which we attribute to the wall tension due to the quadratic Zeeman energy. It is noteworthy that ring-shaped domain walls were sporadically observed supple , which is reminiscent of a 2D skyrmion that is also a topologically allowed defect for the polar phase Choi_prl12 .

To corroborate the existence of HQVs at the wall end points, we measure the mass superflow distribution using a spatially resolved Bragg scattering method Seo_scirep17 . Before applying a pulse of magnetic field gradient for the SG spin separation, we irradiate two pairs of counterpropagating Bragg laser beams onto the sample in the $xy$ plane for $0.8$ ms. The frequencies of the laser beams are set to be resonant to atoms with a velocity of 0.3 mm/s $\approx 0.4\frac{\hbar}{m\xi_{s}}$ so that the atoms that have such high velocities near the HQV cores may be scattered out of the condensate. The HQV core size is characterized by $\xi_{s}$ Seo_prl15 ; Anglin03 ; Lovegrove12 . Then, the mass circulation around the HQVs can be identified by examining the spatial distribution of the scattered atoms with respect to the wall end points supple ; footnote:Condition .

Two examples of data of the Bragg scattering measurement are provided in Fig. 3. In the case of Fig. 3(a), a single domain wall terminates in the center region and a strong scattering signal is detected at the front side of its end point, consistent with the mass flow expected from a HQV with counterclockwise circulation at the end point [Fig. 3(e)]. Figure 3(b) presents another case in which two end points are close to each other. The Bragg signal shows that a superflow passes through the gap between the two walls, indicating that the two HQVs at their end points have opposite circulations [Fig. 3(f)]. The spatial configuration of the two walls, together with the superflow pattern, conjures up the possibility that they might be formed by breaking a single domain wall that initially traverses the condensate [Fig. 2(a)].

The domain wall can be viewed as a three-component soliton Nistazaki_pra08 ; Xiong_pra10 ; Bersano_prl18 , where a dark soliton of the $|0\rangle$ component with a $\pi$ phase step coexists with the bright solitons of the |$$\pm$$1\rangle components. A dark soliton in a scalar BEC is dynamically unstable due to snake instability, which causes a local Josephson current by breaking the soliton Feder_pra00 ; Anderson_prl01 ; Huang_pra03 . If a similar mechanism is active for the domain wall, the wall can break into many free composite defects, i.e., domain walls bounded by two HQVs at both ends. In the experiment with $q_{f}/h=1.0$ Hz, however, we rarely observed free composite defects detached from the condensate boundary, which means that the domain wall does not suffer much from the snake instability. It is the |$$\pm 1$$\rangle components at the wall core that suppress the snake instability by providing an effective pinning potential, which, thus, stabilizes the domain wall. In other words, in the quench dynamics for large $q_{f}$, the snake instability can be enhanced with the domain wall becoming thinner, thus leading to proliferation of free composite defects.

We perform the same quench experiment with a higher $q_{f}/h=10.6$ Hz. The Bogoliubov analysis of the dynamic instability of the initial EPP state shows that the characteristic length and time scales for the quench dynamics are proportional to $\sqrt{q_{f}}$ and $1/\sqrt{q_{f}}$, respectively, for $q_{f}\ll c_{2}n_{0}$ Kawaguchi_phyrep12 , and it is expected that domain walls will be nucleated in a denser pattern within a faster time scale. Indeed, we observe that a characteristically dense network of thinner domain walls develops within 60 ms and, also, that the domain walls subsequently break into many composite defects [Fig. 4(a)], demonstrating the enhancement of the snake instability.

After the fast wall splitting process, short free defects are clearly identified in the center region of the condensate at $t>0.2$ s. The free defects have a spatial size $\lesssim 5\xi_{s}$ with various shapes. Some of the defects appear very round [Fig. 4(a)A], while others show dumbbell shapes, implying splitting [Fig. 4(a)B]. In the subsequent relaxation evolution, the defect number $N_{d}$ decays with a $1/e$ lifetime of $\approx 5$ s, where it is observed that most long-lived free defects exhibit round shapes [Fig. 4(b)]. At $t>0.8$ s, the fractional population of the round defects increases to over 80%. Here we count a defect as a round one if its aspect ratio ($\geq 1$) is smaller than $1.2$. From in situ magnetization imaging Seo_prl15 , we find that the long-lived defects can have nonzero axial magnetization, i.e., contain unequal |$$\pm$$1\rangle spin populations at their cores [Fig. 4(e)], which indicates that the spin current dynamics is intricately involved in the defect formation process.

The free wall-vortex composite defects are classified into two types according to the net circulation around them. One type has a circulation of $h/m$ with two HQVs having the same circulation, which is topologically identical to a singly charged QV in a coarse-grained view [Fig. 4(c)], and the other involves with two HQVs having opposite circulations, which might be described as a magnetic bubble having linear momentum [Fig. 4(d)]. We refer to them as vortex-vortex (VV) and vortex-antivortex (VAV) types, respectively. Immediately after domain wall splitting, the system can contain statistically equal numbers of the defects for the two types. However, when the free defects become short, comparable to $\xi_{s}$ as observed, the defect dynamics of each type will be different because of different HQV interactions Seo_prl16 . We may expect that small VV-type defects will survive longer with the topological character of singly charged QVs, whereas those of the VAV type will decay faster due to their linear motion in the trapped sample and possible vortex pair annihilation Kwon14 . In the experiment, we identify the long-lived round defects as the VV type by confirming the superflow circulation around them with the Bragg scattering measurement supple ; Lovegrove14 .

Finally, we remark on the peculiarity of defect formation in the EPP-to-EAP quantum phase transition. Since the U(1) symmetry is simply broken in the final EAP ground state in a similar manner to the case of scalar BECs, one may expect a conventional Kibble-Zurek scenario for vortex nucleation in our system. However, we observe a two-step defect formation process: first, domain wall creation via the $\mathbb{Z}_{2}$ symmetry breaking, and then, production of composite defects by a splitting of the domain walls. The two-step scenario is also confirmed in our numerical simulation for a uniform system, based on the Gross-Pitaevskii equation for spin-1 BECs supple ; Takeuchi . Our observations show that defect formation in phase transition dynamics critically depends on the symmetry breaking sequence of the system Vilenkin_phyrep85 .

In conclusion, we observed the creation of wall-vortex composite defects in the EPP-to-EAP quantum quench dynamics of an antiferromagnetic BEC and demonstrated the unconventional mechanism of defect formation in the phase transition dynamics. Our findings provide a different framework for the nucleation of composite defects via the Kibble-Zurek mechanism Kibble_book ; Eltsov_prl00 . Additionally, the observation of the free composite defects encourages the efforts to search for the exotic superfluid phase in 2D antiferromagnetic spinor gases James_prl11 ; Chandrasekharan_prl06 ; Podolsky_prb09 .

A1 Domain wall in the polar phase

To describe the theme of spontaneous symmetry breaking in our system, we introduce a complex vector field in the Cartesian representation for spin-1 spinor BECs as

[TABLE]

In this representation, the energy density ${\cal E}({\bm{\psi}})$ is written as

[TABLE]

In antiferromagnetic BECs with $c_{2}>0$, the ground state is the polar state, where we have the spin density vector ${\bm{S}}\equiv i{\bm{\psi}}\times{\bm{\psi}}^{*}=0$ and the order parameter of the polar phase represented by

[TABLE]

with a complex scalar field $\psi_{\rm T}=\sqrt{n}e^{i\phi}$ and a real unit vector field $\hat{\bm{d}}=[d_{x},d_{y},d_{z}]^{\rm T}$.

To capture the feature of domain walls, we consider a stationary solution for a planar domain wall perpendicular to the $y$ axis in a uniform system; $\hat{\bm{d}}\to\pm\hat{\bm{z}}$ for $y\to\pm\infty$ [Fig. 1(a)]. The vector $\hat{\bm{d}}$ varies continuously across the wall, with $\hat{\bm{d}}\bot\hat{\bm{z}}$ at the core ($y=0$). We assume $\phi=0$ and $\hat{\bm{d}}\bot\hat{\bm{y}}$ for $\psi_{+1}=-\psi_{-1}$ without a loss of generality. Then, the planar wall is evaluated by the reduced energy density for the real functions $\psi_{x}$ and $\psi_{z}$,

[TABLE]

For the case of $\mu\gg q$, the total density $n$ ($=\psi_{x}^{2}+\psi_{z}^{2}$) is approximately constant; the wall is thus characterized by the parameter $q$. Comparing the gradient term and the quadratic Zeeman term, the thickness of a domain wall is characterized by

[TABLE]

Figure S1 shows numerical solutions for the domain wall for different values of $q$. According to the expression for $\xi_{q}$, the domain wall becomes thicker as $q$ decreases. This behavior is qualitatively consistent with our experimental observations.

A2 Structure of the wall-vortex composite defect

The structure of a wall-vortex composite defect is understood in a similar manner to that reported in the literature for superfluid 3He-B Kondo_prl92 ; Eltsov_prl00 . In the 3He-B system, the composite defect is formed due to the existence of two different length scales; the coherence length and the dipolar healing length. The length of the former is defined by the superfluid condensation energy and is much shorter than that of the latter, which is related to the much weaker spin-orbit interaction. Correspondingly, our system of anti-ferromagnetic spinor BECs has also two different length scales; $\xi=\frac{\hbar}{\sqrt{2m\mu}}$ and $\xi_{q}$. The length $\xi$ is defined by the term of $c_{0}$ (or $\mu$) in Eq. (S1), which is associated with the condensation energy. The large length $\xi_{q}$ is related to the much weaker quadratic Zeeman energy with $q\ll\mu$; $\xi\ll\xi_{q}$. The length (energy) hierarchy supports the coexistence of different kinds of topological defects as composite defects. This argument is connected with the condition $q\ll\mu$ to stabilize the domain wall as a part of the composite object.

The two length scales $\xi$ and $\xi_{q}$ are associated with a vortex and a domain wall, respectively. If the order parameter varies spatially on length scales much shorter than $\xi_{q}$, the quadratic Zeeman term is unimportant compared with the gradient term. Then, the energy density (S1) can be reduced to the simplest form

[TABLE]

where again we assume ${\bm{\psi}}\times{\bm{\psi}}^{*}=0$ by considering the polar phase. This form of energy density supports the existence of a HQV. A linear HQV is realized by considering a combination of two transformations, $\phi\to\phi+\pi$ and $\hat{\bm{d}}\to-\hat{\bm{d}}$, e.g.,

[TABLE]

with the cylindrical coordinate ${\bm{r}}=(\rho,\theta,z)$ [Fig. S2(a)]. Here, the mismatch between $(\phi=0,\hat{\bm{d}}=\hat{\bm{z}})$ and $(\phi=\pi,\hat{\bm{d}}=-\hat{\bm{z}})$ along the $x$ axis for $x>0$ is avoided and the order parameter varies continuous there because of the discrete symmetry under the operation of $(\phi,\hat{\bm{d}})\rightarrow(\phi+\pi,-\hat{\bm{d}})$. Note that $\hat{\bm{d}}$ is ill-defined at the origin because of the discrepancy for its orientation. In the HQV core, the broken-axisymmetry (BA) phase is energetically preferred causing a magnetized vortex core with finite spin density of ${\bm{S}}\bot\hat{\bm{z}}$. The size of the magnetized core for the BA phase is at most on the order of the spin healing length $\xi_{s}=\xi\sqrt{\frac{c_{0}}{c_{2}}}$.

To understand the structure of a wall-vortex composite defect, it is instructive to show a way to build a composite defect by putting a HQV in the EAP state. A HQV is a high-energy object and should be deformed due to the quadratic Zeeman energy with $q>0$. The order parameter field in Eq. (S4) is deformed to decrease the area of the region for $d_{x}\neq 0$ since the energy density there is higher than that in the $\hat{\bm{d}}=\pm\hat{\bm{z}}$ region. The phase $\phi$ is not affected directly by the quadratic Zeeman energy. After the deformation, the $\hat{\bm{d}}\bot\hat{\bm{z}}$ region condenses along the $x$ axis for $x>0$ to form a domain wall terminating at the HQV at the origin [Fig. S2(b)]. Inside the wall core, the vector field $\hat{\bm{d}}$ flips in a continuous manner from $\hat{\bm{d}}=\hat{\bm{z}}$ to $\hat{\bm{d}}=-\hat{\bm{z}}$ with the phase $\phi$ fixed approximately. The phase $\phi$ rotates from [math] to $\pi$ about the HQV core, while the vector field $\hat{\bm{d}}$ is fixed to point along the $z$-direction outside the wall core.

A3 Numerical time evolution of the quench dynamics in a uniform system

Such composite defects can be nucleated in the course of our non-equilibrium quench dynamics. Domain walls nucleated via the Kibble-Zurek mechanism for the $\mathbb{Z}_{2}$ symmetry breaking can transform into composite defects owing to the snake instability, which causes a local Josephson current and nucleates a pair of HQVs by breaking a domain wall into two segments.

To demonstrate this scenario more clearly in a uniform system, we numerically solve the Gross-Pitaevskii equation for spin-1 BECs. The numerical simulations were done in a method similar to those of the phase transition dynamics described by the Gross-Pitaevskii equation for multi-component systems Takeuchi . Figure S3 shows a two-dimensional simulation for $q_{f}/h=2.6$ Hz from the initial EPP state of $(\psi_{0},\psi_{\pm 1})=(0,\pm\sqrt{\frac{n_{0}}{2}})$ under the Neumann boundary condition at $x=\pm\frac{L}{2}$ and $y=\pm\frac{L}{2}$ with the system size $L=724\xi$. The time evolution is quite consistent with the two-step scenario of the defect formation and explains well the experimental observations.

A4 Superflow detection with Bragg scattering

The superflow distribution in the BECs containing composite defects is investigated by employing a spatially-resolved Bragg scattering method Seo_scirep17 . When atoms are irradiated by a pair of counterpropagating Bragg laser beams along the $x^{\prime}$-direction, a two-photon process may resonantly occur by imparting momentum $\vec{p}_{0}=2\hbar k_{L}\hat{x}^{\prime}$ and energy ${\varepsilon}=\hbar\delta$ to the atoms, where $k_{L}$ is the wavenumber of the two Bragg beams and $\delta$ is their frequency difference. The resonance condition is determined by $\varepsilon=\frac{p_{0}^{2}}{2m}+\vec{p}_{0}\cdot\vec{v}$ for the atom with velocity $\vec{v}$, giving $\delta_{\textrm{d}}\equiv\delta-\delta_{0}=2k_{L}v_{x^{\prime}}$ with $v_{x^{\prime}}=\vec{v}\cdot\hat{x}^{\prime}$ and $\hbar\delta_{0}=\frac{p_{0}^{2}}{2m}$. Because of the velocity dependence of $\delta$, the spatial distribution of the Bragg scattering response for a BEC directly shows the corresponding velocity region in the BEC.

Figure S5(a) illustrates the experimental sequence for our Bragg scattering measurement. We first turn off the optical dipole trap (ODT). After a short $200\leavevmode\nobreak\ \mu$s time-of-flight (TOF), we apply two pairs of counterpropagating Bragg laser beams to the sample for $800\leavevmode\nobreak\ \mu$s [Fig. S5(b)]. Then, we apply an external magnetic field gradient for 3 ms to spatially separate the |m_{z}=$$\pm$$1\rangle components from the $|0\rangle$ component. After a subsequent TOF, we take an absorption image of the sample including the scattered atom clouds [Fig. S5(c)]. The Bragg signal $S_{B}(x^{\prime},y^{\prime})$ is constructed as $S_{B}=n_{-}-n_{+}$ with $n_{\pm}(x^{\prime},y^{\prime})$ being the density distribution of the atoms scattered out to the $\pm x^{\prime}$ direction, translated back to the condensate frame [Figs. 3(c), 3(d), and S6(d)] Seo_scirep17 . The original density distribution of the $|0\rangle$ component can be obtained by combining that of the unscattered, remaining condensate with $n_{+}$ and $n_{-}$, which is used for locating composite defects in the BEC [Fig. 3(a) and 3(b)].

In our Bragg scattering experiments, the sample condition is slightly different from that in the experiments described in the main text. The condensate atom number is $N_{c}\approx 5.8\leavevmode\nobreak\ \times 10^{6}$ and the Thomas-Fermi radii of the condensate are $(R_{x},R_{y},R_{z})\approx(200,143,2.0)\leavevmode\nobreak\ \mu$m for the ODT with trapping frequencies of (\omega_{x},\omega_{y},\omega_{z})=2\pi\times(4.5,6.3,460)\leavevmode\nobreak\Hz. The chemical potential and the peak spin interaction energy are \mu=h\times 927\leavevmode\nobreak\Hz and c_{2}n_{0}=h\times 14.7\leavevmode\nobreak\Hz, respectively. The Bragg beams are red-detuned by 1.7\leavevmode\nobreak\GHz from the $|F=1\rangle$ to $|F^{\prime}=2\rangle$ transition and each beam has an intensity of $0.3\leavevmode\nobreak\ \text{mW/cm}^{2}$ with a $1/e^{2}$ width of $1.8\leavevmode\nobreak\ \text{mm}\gg R_{x,y}$. The TOF duration before taking an absorption image was chosen from the range 6 to 11 ms for different measurements.

Figure S6(c) presents one of the rare occasions in the Bragg scattering measurements for low $q_{f}/h=2.6$ Hz, where the quenched BEC contains free composite defects. Two free composite defects are shown to line up, detached from the boundary of the condensate. The circulation directions of the HQVs at their end points are assigned based on the measured Bragg signal distribution around them [Fig. S6(d)] despite the signal near the condensate boundary being faint due to low atom density. The free composite defect on the left side shows a strong negative Bragg signal along its wall, indicating its linear motion in the $+x^{\prime}$-direction, with its moving direction consistent with the circulation directions of its HQVs.

In Fig. S7, we display four example data of the Bragg scattering measurements at long hold times of $t>1\leavevmode\nobreak\ s$ for high $q_{f}/h=10.0$ Hz, where several free composite defects of small size are present in the quenched BEC. Taking into account the short distances between the free composite defects, the frequency difference of the Bragg beams was set to be \delta_{\text{d}}/2\pi=-2.0\leavevmode\nobreak\kHz, corresponding to a higher velocity of $v_{x^{\prime}}=0.6$ mm/s $\approx 0.8\frac{\hbar}{m\xi_{s}}$, to probe a region closer to the HQVs. The Bragg signal around the long-lived, round-shaped defect shows opposite signs at the different lateral sides with respect to the Bragg beam axis line passing through the defect core. This demonstrates the nonzero superflow circulation around the defect, i.e., that the round-shaped defect is of the VV type. We also present a couple of image data for dumbbell-shaped defects in Figs. S7(c) and (d). The Bragg signal indicates that the defect is of a VAV-type having two HQVs with opposite circulations.