{'fraction_non_alphanumeric': 0.04082747978718274, 'fraction_numerical': 0.020614662530339316, 'mean_word_length': 4.217940115742683, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Vortices in multicomponent Bose-Einstein condensates

We review the topic of quantized vortices in multicomponent Bose-Einstein condensates of dilute atomic gases, with an emphasis on that in two-component condensates. First, we review the fundamental structure, stability and dynamics of a single vortex state in a slowly rotating two-component condensates. To understand recent experimental results, we use the coupled Gross-Pitaevskii equations and the generalized nonlinear sigma model. An axisymmetric vortex state, which was observed by the JILA group, can be regarded as a topologically trivial skyrmion in the pseudospin representation. The internal, coherent coupling between the two components breaks the axisymmetry of the vortex state, resulting in a stable vortex molecule (a meron pair). We also mention unconventional vortex states and monopole excitations in a spin-1 Bose-Einstein condensate. Next, we discuss a rich variety of vortex states realized in rapidly rotating two-component Bose-Einstein condensates. We introduce a phase diagram with axes of rotation frequency and the intercomponent coupling strength. This phase diagram reveals unconventional vortex states such as a square lattice, a double-core lattice, vortex stripes and vortex sheets, all of which are in an experimentally accessible parameter regime. The coherent coupling leads to an effective attractive interaction between two components, providing not only a promising candidate to tune the intercomponent interaction to study the rich vortex phases but also a new regime to explore vortex states consisting of vortex molecules characterized by anisotropic vorticity. A recent experiment by the JILA group vindicated the formation of a square vortex lattice in this system.


Introduction
Due to their being a common feature of fluid motion, vortices appear in diverse fields of science. In a superfluid, that is, a fluid without viscosity, the flow has quantized circulation which originates in quantized vortices 1 . A quantized vortex is commonly viewed as a topological defect of the order parameter in Bose-Einstein condensation. Since the gradient of the phase of the order parameter is proportional to the superfluid velocity, quantized vortices play a crucial role in understanding superfluidity. Moreover, the structure of the order parameters in multicomponent Bose-Einstein condensates, hereafter BECs, are analogous to other studies including superfluid 3 He 2 , unconventional superconductors 3 , quantum Hall systems 4 , nonlinear optics 5 , nuclear physics 6 , and cosmology 7 . Therefore, the investigation of unconventional topological defects in such multicomponent BECs is of great importance and of broad interest.
While quantized vortices have been extensively studied in the field of superfluid helium 1 , they have seen a remarkable resurgence of interest since the realization of Bose-Einstein condensation in dilute alkali-atomic gases 8,9 . This Bose-condensed system has some distinct advantages for the study of quantized vortices. First, the BECs of alkali-atomic gases are so dilute that the interatomic interaction can be accurately parametrized in terms of the s-wave scattering length. Experimentally, the diluteness leads to a relatively large healing length that characterizes the vortex core size, thus enabling us to directly visualize the quantized vortices by optical techniques. Also, because the gas is dilute, the Gross-Pitaevskii equation gives an accurate description of static and dynamic properties of the atomic condensate within a mean-field approximation. This situation differs from superfluid 4 He in which the high density and the strong interaction complicate the theoretical treatments. Second, rapidly developing techniques of atomic, molecular and optical physics allow us to carry out quantum-state engineering of the condensates in a well-controlled manner, thereby opening a new direction for the study of coherent quantum phenomena. For example, since alkali atoms have internal degrees of freedom attributed to the hyperfine spin, an external field can be used to induce conversions between spin states and thus it can be used to control the spatial variation of the wave functions. Moreover, a field-induced Feshbach resonance allows one to tune of the atom-atom interaction and create of a quasi-bound molecular state 10 .
The first experimental detection of a vortex in a dilute alkali-atomic gas BEC was made by Matthews et al. in 1999 11 . Their study involves condensates of 87 Rb atoms residing in two hyperfine states. The generation of a vortex was achieved by a phase imprinting technique proposed by Williams and Holland 12 , in which the interconversion between two components was controlled spatially and temporally by an external coupling field. After turning off the coupling, the vortex state is formed by one circulating component surrounding a nonrotating core of the other component. Leanhardt et al. 13 used "topological phase imprinting" 14 to create a vortex in a spin-1 BEC. These two methods use internal degrees of freedom of the condensate to create vortices, but a simple way to create vortices is to stir a condensate mechanically with a rotating "bucket". Following the latter method, the ENS group observed the formation of a single vortex and multiple vortices in a single-component elongated condensate 15 . The condensate was trapped in a static axisymmetric magnetic trap which was deformed by a nonaxisymmetric attractive dipole potential created with stirring laser beams. This combined potential produces a cigar-shaped harmonic trap with a slightly anisotropic transverse profile. By rotating the orientation of the transverse anisotropy at frequency Ω, they observed the dynamic formation of a vortex above a certain critical value of Ω. A lattice involving more vortices appeared when Ω was further increased. In contrast to indirect visualization methods used in superfluid helium systems, the quantized vortices were directly visualized as hollows in the transverse density profile of the time-of-flight image. Later, groups at MIT 16 , JILA 17 , and Oxford 18 observed vortex lattices with the first two groups succeeding in creating 100 or more vortices.
Theoretical studies mainly used the mean-field Gross-Pitaevskii theory to describe the main features of the vortex states (see Ref. 19 for a review), and several predictions have been shown to agree with experiments. Some of the important problems are concerned with the equilibrium properties of a single vortex, includ-ing its structures and dynamics 20−24 , the critical frequency and the surface-mode instability for vortex nucleation 25−33 , and the nonlinear dynamics of vortex lattice formation 34 Another important issue in vortex physics is to elucidate the vortex phase in multicomponent BECs. Multicomponent order parameters allow the formation of various unconventional topological defects with complex properties that arise from interactions between different order-parameter components 68 . Since it is possible to load and cool atoms in more than one hyperfine spin state or more than one atomic element in the same trap, multicomponent condensates can be realized experimentally 69−77 . Such systems offer an ideal testing ground for the study of such topological defects; moreover the microscopic characterization of these systems should in principle be possible because this system is free from impurities and well controlled by optical techniques.
In this paper we review recent studies on quantized vortices in multicomponent BECs, especially focusing on those in two-component BECs. To the best of our knowledge, there are no reviews of two-component BECs despite the large number of studies of this system 78−96 . Also, while a number of interesting phenomena have been reported on vortices of a spin-1 BEC, we mention this system only briefly and would like to refer the reader to references 97-117 for more details.
This paper is organized as follows. After a brief introduction of quantized vortices in a single-component BEC in Sec. 2, we formulate in Sec. 3 a mean-field theory that describes two-component atomic-gas BECs and give a brief review on the ground state structure without vortices. In Sec. 4, after describing experiments on the creation and observation of vortices in two-component BECs 11 , we discuss the static structure and dynamic properties of the single-vortex state in two-component BECs in a slowly rotating regime. By introducing the "pseudospin" that describes the two-component BECs, we discuss the properties of topologically "trivial" and "nontrivial" spin textures. In particular, we predict a vortex molecule in a twocomponent BEC with an internal coherent coupling. We also mention briefly unconventional vortex states in a spin-1 BEC. In Sec. 5, we consider the vortex states of two-component BECs under a rapid rotation. Combination of an analytical approach and numerical simulations reveals a rich variety of vortex states, which is not found in a single-component BEC. We also discuss the effect of the internal coherent coupling on the structure of vortex states. Finally, we conclude this review article in Sec. 6.

What is a quantized vortex?
Below a critical temperature in an ideal Bose gas, a finite fraction of the particles occupies the same single-particle ground state and forms a BEC. When the particles interact, the single-particle density matrix n (1) (r, r ′ ) = Ψ † (r)Ψ(r) is still defined, whereΨ † (Ψ) is the boson creation (anihilation) operator. When a single-particle state is occupied by a macroscopic number of particles, n (1) (r, r ′ ) has a macroscopic eigenvalue. We may define a condensate wave function Ψ(r, t) as the corresponding eigenstates 118 . In the absence of interactions this wave function reduces to that of the macroscopically occupied single particle state. When a BEC is held in an external potential V (r), the dynamics of Ψ(r, t) is described by the Gross-Pitaevskii (GP) equation 8,9 i ∂Ψ(r, t) ∂t = − 2 2m ∇ 2 + V (r) + g|Ψ(r, t)| 2 Ψ(r, t), where g = 4π 2 a/m represents the strength of interaction characterized by the s-wave scattering length a and particle mass m. The condensate wave function is normalized by the total particle number N as d 3 r|Ψ| 2 = N . Writing Ψ(r, t) = ρ(r, t) exp(iθ(r, t)), the absolute squared amplitude ρ(r, t) = |Ψ(r, t)| 2 gives the condensate density and the gradient of the phase θ(r, t) gives the superfluid velocity v s (r, t) = ( /m)∇θ(r, t). As a result, the circulation of v s around a closed path C in the fluid is quantized as C ds · v s = m C ds · ∇θ = nκ (n = 0, ±1, ±2, . . .) (2) with the quantum of circulation κ = h/m. Such a vortex with a quantized circulation is called a "quantized vortex". Let us review the fundamental properties of quantized vortices. Consider a single straight vortex line along the symmetry axis of a cylindrical vessel that can rotate at frequency Ω = Ωẑ. In this case, the velocity field of the vortex line is where r is the radial distance from the line and θ is the azimuthal angle. The corresponding "vorticity" is There is a singularity on the vortex line (r = 0), where v s diverges. Near r = 0 the superfluid density is significantly supressed, forming the vortex core. In the case of superfluid helum, the size of the vortex core r c is of the order of the interatomic distance, while in a dilute gas this size equals the healing length r c ≃ ξ, where The angular momentum L z of the fluid due to the vortex line is where L is the height of the vessel, R is its radius, andρ is the value of the superfluid density far from the core. In deriving the last equality we assumed r c ≪ R. The energy associated with the vortex line is dominated by the kinetic energy term. A straight calculation yields Since E v ∝ n 2 , vortices with n > 1 are energetically unfavorable. That is, the energy cost to create one n = 2 vortex is higher than that to create two n = 1 vortices, even though they give the same angular momentum for R ≫ 1 (L z ∝ n). Therefore, a stable quantized vortex always has n = 1; the disintegration of a multiply quantized vortex has been intensely studied 119−122 and it was observed experimentally by Shin et al. in an atomic-gas BEC 123 . Using Eqs. (6) and (7) one can easily find the critical frequency of rotation Ω c for the existence of an energetically-stable vortex line. In a frame rotating with Ω, the free energy is given by F = E − ΩL z and the presence of a vortex increases E by E v . For a vortex with n = 1, the critical frequency is obtained by imposing that the change in the energy ΩL z due to the creation of the vortex be equal to E v as For Ω > Ω c the free energy of the system with one vortex is lower than that without a vortex. When the rotation frequency is significantly higher than Ω c , more vortices will appear in the form of a triangular lattice of vortices. For very large Ω the rotation of the superfluid will then look similar to the one of a rigid body characterized by ∇ × v s = 2Ω. Using Eq. (4), one finds that the average vorticity per unit area is given by ∇ × v s = κn vẑ , wheren v is the number of vortices per unit area, so that the density of vortices is related to the rotation frequency Ω as This relation, often referred to as Feynman's rule 124 , can be used to determine the maximum possible number of vortices in a given area as a function of Ω. The above formalism is extended to quantized vortices in a trapped BEC by including the explicit form of the trapping potential V (r). This potential leads to an inhomogeneous condensate density, which plays an important role in the stability and dynamical behavior of the vortices. As an example, consider the structure of a single vortex in a condensate trapped by the axisymmetric harmonic potential V (r, z) = mω 2 ⊥ (r 2 + λ 2 z 2 )/2 with the aspect ratio λ = ω z /ω ⊥ . The equilibrium state is determined by the time-independent GP equation where the time dependence of Ψ 0 is fixed by the chemical potential µ as Ψ(r, t) = Ψ 0 (r)e −iµt/ . The value of µ is determined by the normalization condition |Ψ 0 (r)| 2 = N . The condensate wave function with a singly-quantized vortex line located along the z-axis takes the form Ψ 0 (r) = f (r, z)e iθ with the cylindrical coordinate (r, θ, z), where f is a real function that is related to the condensate density as ρ 0 = f 2 . Then, Eq. (10) becomes Here, the centrifugal term 2 /2mr 2 arises from the azimuthal motion of the condensate which gives rise to a kinetic energy density ρ 0 mv 2 s /2 = ρ 0 2 /2mr 2 . The trapping potential defines the characteristic length scale b ho = /mω ⊥ . Without an interparticle interaction (g = 0), the scale b ho determines the spatial extent of the wave function. An increase in the repulsive atomic interaction expands the condensate, so that the kinetic energy associated with the density variation can be made negligible compared with the trap energy and the interaction energy. The ratio of the interaction energy to the harmonic-oscillator energy gρ 0 / ω ⊥ ∼ N a/b ho serves to quantify the effect of the interaction. In the limit N a/b ho ≫ 1, which is applicable to current experiments on trapped BECs, one can omit the terms involving derivatives with respect to r and z in Eq. (11). This approximation is called the Thomas-Fermi approximation 125 , and the density profile for the vortex state is then given by where Θ(x) is a step function and ρ 0 (0) is the density at the center of the vortexfree Thomas-Fermi profile. We define the Thomas-Fermi radius R ⊥ = 2µ/mω 2 ⊥ and R z = 2µ/mω 2 z and the healing length ξ(0) = 2 /2mgρ 0 (0). Due to the centrifugal term ξ(0) 2 /r 2 , the condensate density vanishes at the center out to a distance of order ξ, whereas the density in the outer region has the form of an upward-oriented parabola. The Thomas-Fermi approximation is valid in the region ξ ≪ r < R ⊥ . The asymptotic form of ρ 0 (r) for r ≪ ξ is 126 ρ 0 (r) = ρ 0 (0)(r/ξ(0)) 2 whose interpolation to Eq. (12) gives This form agrees fairly well with the profile of the single vortex state obtained by the numerical calculation of the GP equation 126 . For a detailed account of the single vortex state, we refer to a review article by Fetter and Svidzinsky 19 (see also Chap.9 in Ref. 8 and Chap.14 in Ref.9).

Coupled Gross-Pitaevskii equations
We now consider two-component BECs as described by the condensate wave functions for the two hyperfine spin states Ψ 1 and Ψ 2 11,69,70 . These wave functions can also describe a system of two atomic species 77,127,128 or two isotopes of an atomic element. The theory for the latter case is described in Refs. 129 and 130, and experiments are described in Ref. 131. In addition, the condensates are confined by trapping potentials V i (r) (i = 1, 2), which is assumed to rotate at a rotation frequency Ω about the z axis as Ω = Ωẑ. Viewed from the rotating frame, the GP energy functional of the two-component BEC is 8,9 E[Ψ 1 , Here, m i is the mass of the i-th atom and ΩL z = i Ω(y∂ x − x∂ y ) is a centrifugal term. The coefficients g i (i = 1, 2) and g 12 represent the atom-atom interaction between atoms iin the same hyperfine state, and atoms in different hyperfine states, respectively. They are expressed in terms of the corresponding s-wave scattering lengths a 1 , a 2 , and a 12 as where m −1 12 = m −1 1 + m −1 2 is the reduced mass. When the trapping potential V i (r) is produced by a magnetic field, each component has its own potential energy because the potential depends on a magnetic moment of the atom. For example, the trapping potential for atoms with the hyperfine states |F i , M F i is expressed as a function of the magnetic field |B(r)| as 127 where µ B is the Bohr magneton, G i is the g-factor of the i-th atom, V 0 and ω ik (k = x, y, z) is the trapping frequency satisfying the relation The potential minima V 0 i for the BECs with different values of G i M F i do not coincide. This positional shift is also influenced by a gravitational effect, by the difference of nuclear magnetic moments, and by nonlinearity in the Zeeman shift 69,70 .
If the potential is created by a pure optical laser 71,73,74,75 , all atoms share the same harmonic potential.
The time-dependent coupled GP equations for two-component BECs can be obtained by using a variational procedure i ∂Ψ i = δE/δΨ * i as where we have introduced the Lagrange multiplier µ i which represents the chemical potential and is determined so as to satisfy the conservation of particle number In the two-component BECs, each component interacts with the other through the intercomponent mean-field coupling ∝ g 12 , which yields structures and dynamics not found in a single-component BEC.

Ground state structure
It is instructive to start with the ground state of two-component BECs without rotation (Ω = 0). The equilibrium solutions of the coupled GP equations (21) and (22) show a rich variety of stable structures 78,130,127,132−142 , depending on the various parameters of the system. In particular, the intercomponent interaction g 12 plays an important role in determining the structure of the ground state; when the intercomponent interaction is strongly repulsive, the two components phase separate 136,137 . The qualitative features of the ground state for a trapped two-component BEC can be understood by using the Thomas-Fermi approximation 78,127,140 in which one neglects the quantum pressure terms ( 2 /2m i √ ρ i )∇ 2 √ ρ i with the condensate density ρ i = |Ψ i | 2 in Eqs. (21) and (22). First consider the simple case in the spatial region where the two components coexist (ρ where, we defined the following dimensionless parameters For the two components to coexist, the right-hand sides of Eqs. (23) and (24) must be positive. Then, the parameter Γ 12 must be positive, that is, The coexisting region is given by the region delimited by the surface defined as the solution to the equation µ i − g 12 µ j /g j − Γ j V i (r) = 0 127 . Otherwise, that is when there is no coexisting region, which implies that the two components phase-separate. For Γ 12 > 0, there are two possible density profiles that depend on the signs of the parameters Γ 1 and Γ 2 . These signs determine the curvature of the density profile at the center of the trapping potential because Γ 1 and Γ 2 are the coefficients of the x 2 , y 2 and z 2 terms in Eqs. (23) and (24). For Γ 1 > 0 and Γ 2 > 0, both density profiles are written by upward-convex parabolas in the coexisting region. For Γ 1 > 0 and Γ 2 < 0 (Γ 1 < 0 and Γ 2 > 0), one density profile is an upward-oriented parabola, but the other is a usual parabola. In the peripheral region in which the density of one component vanishes, the other component follows the single-component Thomas-Fermi profile ρ (0) i = (µ i − V i (r))/g i . This profile is connected to the density profile of Eqs. (23) and (24) by adjusting the chemical potential µ i (i = 1 or 2) subject to the normalization condition d 3 rρ 78 . For Γ 12 < 0 the two components phase-separate as the coexisting regions do not exist. In the region in which only one density is nonvanishing (ρ (0) i = 0 and ρ (0) k = 0, for i = k), the density profile is given by the single-component Thomas-Fermi profile ρ (0) i = (µ i − V i (r))/g i . Since the Thomas-Fermi approximation fails to describe the domain boundary region, where the quantum pressure term cannot be neglected, numerical analysis of Eqs. (21) and (22) is necessary to determine the explicit density profile in this regime 135 . The numerical calculation shows that there are two characteristic configurations of the phase-separated two-component BEC in the trapping potential. One configuration consists of a core of one component at the center of the trap and a shell of the other around the one component, a structure that preserves the spatial symmetry of the trapping potential 78,135 . The other configuration breaks the spatial symmetry by displacing the center of each condensate from the center of the trapping potential 134,138 . Whether the ground state takes the former or the latter configuration depends not only on the intercomponent interaction but also on the particle number, intracomponent interaction, and the shape of the trapping potential 139,141,142 .
Compared to single-component BECs, two-component BECs exhibit a rich variety of ground state structures. These structures and their stability depend upon various parameters, particularly on the condition of the phase separation given by 2vorrev Vortices in multicomponent BECs 11 Eq. (29). Two-component mixtures of BECs with different spin states have been created in the laboratory for the systems of 87 Rb (Refs. 69, 70) and 23 Na (Ref. 72). For the case of 87 Rb, the nearly-equal scattering lengths (a 1 = 5.67 nm, a 2 = 5.34 nm, and a 12 = 5.50 nm) and the small relative displacement of the trapping minima lead to a weak separation between the two components. In contrast, for the 23 Na-mixture, the value a 1 = 2.75 nm, a 2 = 2.65 nm, and a 12 = 2.75 nm causes a spin domain to form 71,72 . Current experiments are aiming to create two-component BECs with two atomic species 128,143 ; a mixture of the 41 K and 87 Rb BEC has been observed 77 . In this case, the intercomponent scattering length a 12 may take values different from the intracomponent ones a 1 or a 2 . Furthermore, we should note that an additional possibility to change the value of a 12 by the use of the Feshbach resonance 129,144−147 .

Creation of a quantized vortex
Although most theoretical studies on quantized vortices have focused on singlecondensate systems, the first observation of a quantized vortex in a gaseous BEC was achieved in a two-component system 11,12 . In the JILA experiment 11 , the vortex was created in a two-component BEC consisting of 87 Rb atoms with the hyperfine spin states |F = 1, M F = −1 ≡ |1 and |F = 2, M F = 1 ≡ |2 . Since the scattering lengths between atoms of |1 and |1 , |2 and |2 , and |1 and |2 are all different, the two states are not equivalent, and the two-component condensate is characterized by the two-component order parameter.
The advantage of using the |1 and |2 states of 87 Rb atoms to perform quantum state engineering is that, first, their magnetic moments are nearly equal, so that they can be confined simultaneously in almost identical magnetic potentials. This is advantageous to investigate the intrinsic interaction phenomena between two condensates. Second, nearly equal s-wave scattering lengths for two 87 Rb atoms in |1 and |2 result in the fortuitously small inelastic spin-exchange rate 148 , which provides a stable binary system with a long lifetime. In contrast, the inelastic rate is very large in 23 Na atoms. Third, we can use a microwave rf-pulse to change the populations of atoms from the |1 state to the |2 state 11,149,150 . If a microwave is applied continuously, this system is referred to as an "internal" Josephson effect because the coexisting components are coupled through their internal degrees of freedom 151−154 . Thus, this system offers new possibilities not only for quantummechanical state engineering, but also for exploring new vortex structures in coherently coupled multicomponent BECs as described in Sec. 4.3. In experiments at JILA experiments 70 , atoms are initially trapped in one state, say, the |1 state, and then cooled below the Bose-Einstein transition temperature. When the atoms in the |1 state forms a condensate, a two-photon microwave field is applied, inducing transitions between the |1 state and the |2 state. As a result, the atoms undergo coherent oscillations between the two hyperfine levels with an effective Rabi frequency ω eff R = ω 2 R + ∆ 2 , where the Rabi frequency ω R gives the rate at which population would oscillate between the two states when ∆ = 0 and depends on the microwave power and the atom-field coupling constant. The detuning parameter ∆ denotes the difference between the frequency of the coupling rf field and the frequency difference between the two internal atomic states.
The effect of internal coupling can be formulated by introducing the energy density which allows atoms to change their internal state coherently. Since the microwave coupling fields are time dependent, it is useful to use the frame in which the field is constant (a frame rotating with the laser field). Then, the coupled time-dependent GP equations read The last terms in Eqs. (31) and (32) cause interconversions between the components. For a homogenous system without V i but having uniform phases of the condensate wave functions for both components, the interconversion takes place at the same rate everywhere. However, the time variation of the spatially inhomogeous potential V i changes the way the interconversion is made, which is a key ingredient of the phase imprinting method for vortex creation 12 .
To proceed, let us assume that the external potential consists of a radially isotropic harmonic potential and a time-dependent perturbation V 1 (r, t) = mω 2 ⊥ r 2 /2 +Ĥ 1 and V 2 (r, t) = mω 2 ⊥ r 2 /2 −Ĥ 1 witĥ where κ characterizes the strength of the perturbation and f (r) and g(r) describe the spatial dependence of the perturbation. The form ofĤ 1 can be determined from the symmetry of the quantum state to be prepared. To create a vortex state with one unit of angular momentum, one can take κ = mω 2 ⊥ r 0 , f (r) = x, and g(r) = y in Cartesian coordinates 12 . This form of perturbation effectively confines the two hyperfine states in separate axisymmetric harmonic potentials with the same trap frequency ω ⊥ . The trap centers are offset spatially in the x-y plane by a distance r 0 from the center and rotate about the symmetry axis at angular frequency Ω ′ .
The underlying physics for creating a vortex can be understood by considering the co-rotating frame with the trap centers at the angular frequency Ω ′ of the perturbationĤ 1 . In this frame, the energy of a vortex with one unit of angular momentum is shifted by Ω ′ relative to its value in the laboratory frame. When this energy shift is compensated for by the sum of the detuning energy ∆ and the small chemical potential difference between vortex and non-vortex states, a resonant 2vorrev Vortices in multicomponent BECs 13 transfer of population can occur. From an initial nonrotating component |1 , the |2 component is transferred to the state with unit angular momentum by controlling precisely the time of turning off of the coupling drive; this scheme was confirmed to work by numerical simulation of Eqs. (31) and (32) 12 . A typical structure shows that the |2 component has a vortex with a single quantized circulation at the center, whereas the nonrotating |1 component is located at the center and provides a "pinning" potential that stabilizes the vortex core [see Figs. 3 and 4].
Following the above method, Matthews et al. created a vortex for the first time in trapped BECs 11 . To realize the above configuration, they used a microwave and a movable laser beam with a spatially inhomogeneous profile on the condensate 11 . The laser beam was rotated around the initial condensate as in Fig. 1(a), which reproduces the perturbation given by Eq. (33). The resonant transfer of population from the nonrotating condensate into the vortex state was accomplished by detuning the microwave frequency and rotating the laser beam at an appropriate frequency for the resonant coupling. As shown in Fig. 1(b) for large detunings, the energy resonance condition means that atoms can change only the internal state through the coupling of the time-varying perturbation, and therefore must follow any selection rule that the spatial symmetry of that perturbation might impose. It should be emphasized that it is not simply a mechanical force of the optical field that excites the vortex. If one changes the sign of detuning ∆ at a fixed trap rotation frequency, a vortex with opposite circulation will be created. Vortices with oppo-site circulations experience opposite energy shifts in transforming to the rotating frame, and therefore they must be detuned in an opposite sense to obtain resonant coupling.
Experimentally, it is possible to put the initial condensate into either the |1 or |2 state, and then make a vortex in the |2 or |1 state. In the case of 87 Rb atoms, the intracomponent and intercomponent scattering lengths are in the proportion a 1 : a 12 : a 2 = 1.03 : 1 : 0.97 with the average of the three being 55Å. For the case without vorticity, the atoms with the larger scattering length a 1 in the state |1 formed a lower-density shell around the atoms with the smaller scattering length  The starting point for the calculation is the coupled GP equations (21) and (22). Let us assume that the condensate wave functions have one quantized vortex at the center as Ψ i (r) = f i (r, z)e iγiθ ; then Eqs. (21) and (22) become where ρ i = f 2 i is the condensate density and the trapping potential is given by V i = m i (ω 2 i⊥ r 2 + ω 2 iz z 2 )/2. The axisymmetric vortex state is obtained by solving Eqs. (34) and (35) with (γ 1 , γ 2 ) = (1, 0) or (0, 1).
Under the Thomas-Fermi approximation, the quantum pressure terms (derivatives with respective to the coordinates r and z) in Eqs. (34) and (35) are neglected. Then, the density profiles of the condensates can be calculated in a manner similar to that for the non-vortex case shown in Sec. 3. In the region where only one wave 2vorrev Vortices in multicomponent BECs 15 function is nonvanishing (ρ j = 0 and ρ k = 0 for j = k), Eqs. (34) and (35) are decoupled and the solution is On the other hand, in the region where the wave functions overlap (f i > 0), one obtains Here, the chemical potential has been rewritten as µ i − γ i Ω → µ i . The term 2 γ i /(2m i r 2 ) represents the density singularity of the vortex core. To obtain the density profile, we have to calculate the chemical potentials µ i from the normalization condition d 3 rρ i = N i . This leads to a problem involving two coupled integral equations that can be solved numerically. Choosing the parameters corresponding to those of the JILA experiment 11 , namely, a 1 = 5.67 nm, a 2 = 5.34 nm, a 12 = 5.50 nm, ω ⊥ = 7.8Hz, ω z /ω ⊥ = 1/ √ 8, and N 1 = N 2 = 5 × 10 5 , we calculate the density profile for (γ 1 , γ 2 ) = (0, 0) as well as (1, 0) and (0, 1) by both using the Thomas-Fermi approximation and numerical simulation, and the results are shown in Fig. 3. When neither component has nonzero vorticity [ Fig. 3(a)], the components phase separate due to the intercomponent repulsion. Specifically, to decrease the intracomponent mean field energy, the Ψ 1 component with the larger intracomponent scattering length forms a shell outside the Ψ 2 component in the central region. The density contours in the (r, z) plane show that |Ψ 1 | 2 extends more along the z axis than the r axis. When the Ψ 1 component has a vortex, a density hole associated with the vortex core appears in the density of Ψ 1 and penetrates through the z-axis as shown in Fig. 3(b). Then, the Ψ 2 component fills the vortex core because of the intercomponent repulsion, forming a "cigar-shaped" profile. As a result, the characteristic scale of the core size becomes larger than the healing length ξ = /(2mg 1 ρ 1 ), which characterizes the core size of a single-component BEC. In this case, the centrifugal force associated with the vortex makes the Ψ 1 component expand radially, which allows a decrease in the intracomponent mean-field energy of the Ψ 1 component rather than that of the nonvortex state. On the other hand, when the Ψ 2 component with the smaller intracomponent scattering length has a vortex, the structure is different from the (γ 1 , γ 2 ) = (1, 0) solution as shown in Fig.  3(c). First, the core size is smaller than that of the (1, 0) solution. Second, while a fraction of the Ψ 1 component fills the vortex core, the rotating Ψ 2 component is surrounded by a shell of the excessive Ψ 1 component, despite the centrifugal force. As discussed below, the latter configuration is dynamically unstable 82,83,84 . The above consideration shows that the Thomas-Fermi approximation captures some qualitative features of the solutions 78,80,81 . However, one cannot neglect the quantum pressure term for the description of the vortex states in this system; in particular, the quantum pressure is necessary to determine the boundary profile between segregated components 136,137 .

Pseudospin representation: Nonlinear sigma model
To analyze the vortex structure, we present a 2-D (two-dimensional) model 92 that is based on the concept of "pseudospin" 79,91,110,150 . The spinor order parameter of twocomponent BECs allows us to analyze this system as a spin-1/2 BEC 79,93,94,110,150 .
To simplify the problem, we assume that all atoms have the same mass (m 1 = m 2 = m) and the trapping potentials are equal (V 1 = V 2 = V ). We map the Hilbert space of the system into pseudospin space by introducing the normalized complex-valued spinor Here, ρ T = ρ 1 + ρ 2 is the total density, and the spinor thus satisfies In terms of the spinor, the spin density is defined as S =χ(r)σχ(r), where σ is the Pauli matrix. Explicit expressions of S = (S x , S y , S z ) read where the modulus of the total spin is |S| = 1 everywhere. We assume that Ψ 1 and Ψ 2 represent the up and down components of the spin-1/2 spinor, respectively. The nonzero spin projection on the x-y plane implies a relative phase coherence between the up and down components. The transformation of Eq. (14) to the pseudospin representation is straightforward 92 . The result is where we define new coupling constants as and an effective velocity field as which depends on the gradient of the total phase Θ = θ 1 + θ 2 and that of the pseudospin. The form of Eq. (43) is analogous to the classical nonlinear sigma model (NLσM) for Heisenberg ferromagnets in which only the (∇S) 2 term appears. A NLσM allows one to study a various type of "topological excitations" such as monopoles, domain walls, and skyrmions 155 . Because we have changed merely the representation of the Gross-Pitaevskii energy functional, Eq. (43) is still an exact description of the two-component BEC. The four degrees of freedom of the original condensate wave functions Ψ 1 and Ψ 2 (their amplitudes and phases) have now been expressed in terms of the total density ρ T , the total phase Θ, and two of the spin density components (S x , S y , S z ) (one of which is determined besides its sign by the other two because of |S| = 1).
Unlike the classical NLσM case, the system here has several unique features that are revealed in Eq. (43): (i) the total density ρ T , which is a prefactor of the (∇S) 2 term and effectively 'stiffens' the pseudospin, is spatially inhomogeneous because of the trapping potential; (ii) there is a hydrodynamic kinetic-energy term mρ T (v eff − Ω × r) 2 /2, which is associated with the topological excitations; (iii) if g 1 = g 2 = g 12 , the anisotropic terms involving the coefficients c 1 and c 2 effectively break the "SU(2)-invariance" 79,92,93,156 because then the rotation of S removes the energy degeneracy. The coefficient c 1 can be interpreted as a longitudinal Vortices in multicomponent BECs 19 (pseudo)magnetic field that tends to align the spin along the z-axis. The term with the coefficient c 2 may determine the spin-spin interaction associated with S z , which is antiferromagnetic for c 2 > 0 and ferromagnetic for c 2 < 0 88 . If m 1 = m 2 or V 1 = V 2 , other anisotropic terms appear in Eq. (43).
In the following, we confine ourselves to the 2-D problem for simplicity. Twodimensional calculations of the GP equations (21) and (22) are carried out by separating the degrees of freedom of the original wave function as Ψ i (r) = ψ i (x, y)φ(z). This separation yields the dimensionless 2-D coupling constants u i = 4πa i ηN and where the length and time are measured in units of b ho = /mω ⊥ and ω −1 ⊥ , which are the characteristic scales of the radial trapping potential V = mω 2 ⊥ r 2 /2. Also, the NLσM (43) becomes with c 0 = (u 1 + u 2 + 2u 12 )/4, c 1 = (u 1 − u 2 )/2, and c 2 = (u 1 + u 2 − 2u 12 )/4.

A skyrmion: Numerical analysis
First, we consider the axisymmetric vortex state for the case g 1 = g 2 = g 12 (c 1 = c 2 = 0), in which case the system possesses SU(2) symmetry. In the pseudospin picture, the axisymmetric vortex observed in Ref. 11 can be interpreted as a type of spin texture called a "skyrmion" 79,110 . Typical examples are shown in Fig. 4. In Fig. 4(a), we assume that the ψ 1 component has one singly-quantized vortex at the center of the trap. At the center of the atomic cloud, the ψ 1 component vanishes, so that the pseudospin points down according to the definition of the spin S z in Eq. (42). The spin aligns with a hyperbolic distribution as (S x , S y ) ∝ (−x, y) around the singularity at the center. At the edge of the atomic cloud, the ψ 2 component vanishes, and the pseudospin points up. In between, the pseudospin continuously rolls from down to up as shown in the bottom of Fig. 4. This "cross-disgyration" texture is often referred to as a skyrmion in analogy to the work of Skyrme 6 . In the field of superfluid 3 He, this skyrmion texture is referred to as an Anderson-Toulouse vortex 79 . If ψ 2 component has a vortex, the texture exhibits a "radial-disgyration" as illustrated in the middle panel of Fig. 4(b), in which the spin in the x-y plane aligns as (S x , S y ) ∝ (x, y) around the singularity. These two skyrmions are degenerate in the SU(2) symmetric case because they are interconverted by the overall spin; in fact, all configurations obtained by any global spin rotation are degenerate. The size of the skyrmion is determined by the interaction coefficients and the rotation frequency 92 , and is related to the ratio N 1 /N 2 . When we minimize the energy to determine the ground state, we fix the total particle number, not each particle number. This procedure gives us the true thermodynamically stable configuration of the skyrmion. Although each particle number is conserved in experiments performed to date (this case will be discussed below), the population imbalance determined as above can be realized by controlling the strength and time of the coupling drive in the phase imprinting method.
It is known that the skyrmion in a 2-D system has a topological invariant defined as which is known as the topological charge or the Pontryagian index 155 . The skyrmion with any spin profile is shown to have Q = ±1, whose sign depends on the direction of the circulation. The integrand of Eq. (49) is the topological charge density q(r) associated with the vorticity. The topological charge density can be derived from the effective velocity v eff as 92,110,111 : where we used the relation i χ i ∇χ * i = − i χ * i ∇χ i (i = 1, 2). The topological charge density q(r) characterizes the spatial distribution of the skyrmion. Figure 5 shows the spatial distribution of q(r) and the corresponding v eff -field for the solution of Fig. 4(a). The topological charge is distributed around the center and, contrary to the case of a conventional vortex in a single-component condensate, |v eff | vanishes at the center. This makes a coreless vortex (skyrmion) without a density dip in the total density.

A skyrmion: Variational analysis
To study the physical properties of the skyrmion in more detail, we used a variational analysis on the NLσM described by Eq. (48). In the following, we confine ourselves to the case in which the ψ 1 component has a vortex. The original NLσM, where only the (∇S) 2 term is present in Eq. (48), admits a skyrmion solution whose explicit analytic expression is known 155 . Here we seek for a more general form of the skyrmion solution of Eq. (48) as a variational function, which can reproduce well the numerical solution. The skyrmion solution of Fig. 4(a) may be parametrized with |S| = 1. Equations (51) with α = 0 correspond to the explicit skyrmion solution of the classical NLσM 155 . The variational parameters λ and α determine the size and the shape of the skyrmion. Typically, λ represents a characteristic size of the domain over which the direction of the spin is maintained and α is the degree of spatial variation of the spin inversion. If we take both λ and α as variable parameters, the number of particles in each pseudospin component is not conserved, but the total particle number is conserved. In the limit λ → 0 (∞) with fixed α all spins point up (down) along the z-axis, which means a perfect polarization of the particle number such as N → N 1 (N 2 ). Alternately, for fixed λ, the limit N → N 1 (N 2 ) corresponds to α → +∞ (−∞). Because the ψ 1 component has one singlyquantized vortex at the center, the total phase Θ in Eq. (48) is Substituting Eqs. (51) and (52) into Eq. (48) yields Thus, there are now fewer degrees of freedom in the energy functional: the total density ρ T and the two variational parameters λ and α. In the above equation, we introduced an effective, radially symmetric, confining potential which determines the shape of the total density ρ T . The total density ρ T can be calculated by solving the equation where the chemical potential µ is fixed by the normalization condition d 2 rρ T = 1.
Using the condition d 2 rρ T = 1, we numerically calculate the chemical potential µ in Eq. (55) with the Thomas-Fermi approximation, where ρ T is assumed to be zero in the region in which ρ T is negative. Then, the total energy E becomes a function of (λ, α) by using the Thomas-Fermi expression for ρ T . When the size of the skyrmion is large enough 92 , the Thomas-Fermi approximation can be used to evaluate ρ T , and then the energy can be minimized with respect to only two variational parameters λ and α. Otherwise a density singularity associated with the vortex core appears in ρ T so that the quantum pressure term cannot be ignored. For the SU(2) symmetric case without rotation (Ω = 0), the minimum energy occurs when λ → ∞ and α → −∞, which implies that N = N 2 , i.e., complete polarization of the particle number. It follows that the global minimum for Ω = 0 is a nonvortex state for any value of N 2 for a fixed total particle number N . This is because the interaction energy terms now satisfy the SU(2) symmetry so that the energy of the nonvortex state is degenerate under the change of the ratio N 1 /N 2 with fixed N . However, above a certain critical value of Ω, there appears an energy minimum for particular finite values of λ and α. Since the minimized energy is lower than that of the nonvortex state with λ → +∞, an appropriate external rotation can provide global stability for the skyrmion. Figure 6(a) shows the values of λ and α that give the minimum of the total energy for c 0 = 1000, 2500, and 10000. The size of the skyrmion decreases with increasing Ω, as revealed by a decrease in λ min . The same effect occurs for an increase in α. The energy minimum appears at a certain critical frequency of Ω where λ min (α min ) diverges to +∞ (−∞) and the critical frequency decreases as c 0 increases. Therefore, the condensates with larger c 0 can have a stable skyrmion at a lower rotation frequency. The comparison with the numerical solution in Fig. 6(b) shows that the optimized variational functions given in Eqs. (51) and (52) nearly reproduce the exact numerical solution. Hence, our approach here greatly improves upon the analytic treatment for the vortex states beyond the usual Thomas-Fermi approximation in Refs. 78, 80, 81. However, full numerical calculations show that for c 0 = 1000, 2500, 10000 additional vortices are nucleated above Ω > 0.17, 0.11, 0.05 and the ansatz (51) for a single skyrmion is no longer valid.
The parameters c 1 and c 2 in Eq. (53), the values of which can be controlled through the scattering lengths, predominantly determine the size of the skyrmion. The coefficient c 1 , which is proportional to the difference between the intracomponent interactions u 1 and u 2 , may be regarded as a (pseudo)magnetic field that tends to cause the spin to align along the z-axis; the positive (negative) c 1 arranges the spins downward (upward). For the negative c 1 (u 1 < u 2 ) the stable size of the skyrmion shrinks as |c 1 | increases, which implies an increase in the fraction of the rotating ψ 1 component and a decrease in the fraction of the ψ 2 component at the vortex core, (i.e., the spins align upward). For the positive c 1 (u 1 > u 2 ) the value of λ increases rapidly and eventually goes to infinity (concurrently, α goes to −∞), corresponding to the perfect "spin-down" alignment (occupation of the nonrotating ψ 2 component only).
The sign of c 2 determines the nature of the spin-spin interaction associated with S z . For positive c 2 (u 1 + u 2 > 2u 12 ) the interaction is antiferromagnetic, that is, a spatial mixture of the spin-up and spin-down components. The presence of a vortex causes an effective phase separation. Then, there appears to be no significant change of the spin profile compared to that of c 2 = 0 in Fig. 4(a). In contrast, for negative c 2 (u 1 + u 2 < 2u 12 ), the system enters the ferromagnetic phase, where the spin domains are spontaneously formed. When there is a skyrmion in the ψ 1 component, its size shrinks with |c 2 | in a way similar to the |c 1 |-dependence. This is because each particle number is not conserved in the present case; the energetically favorable configuration in a ferromagnetic phase tends to an overall spin polarization, that is, a perfect polarization of the particle number. In this regime, depending on the rotation frequency, there are two energy minima, one of which leads to a perfect polarization of the ψ 1 vortex state, but the other leads to the ψ 2 nonvortex state, which is stable.
So far, we have allowed a change in the particle number of each component to facilitate the calculation of the minimum skyrmion size. However, the actual experiments on two-component BECs have been done under a condition in which each particle number is fixed. This restriction can be taken into account by imposing the relation For a given λ the value of α is uniquely determined by Eq. (56), where both ρ T and S z are functions of two variational parameters. Therefore, the energy minimization can be done with respect to one variational parameter, which we choose to be λ. Let us here consider the case N 1 /N 2 = 1, and investigate the stable size of a skyrmion.
In this case, we find that the stable size is not affected by the change in rotation frequency Ω. Figure 7 shows the optimized values of λ and α as a function of c 1 and c 2 . No perfect polarization of the particle number occurs, so that the skyrmion may exist for all values of c 1 and c 2 ; moreover, in this case, even if the optimal variational parameter exist, they do not ensure a local minimum of the energy. In Fig. 7(a), for In this section, we discuss a nonaxisymmetric vortex state. The spin texture of this state involves a pair of "merons" 4 or "Mermin-Ho vortices" 79,98,110 ("meron (µǫρoσ)" means "fraction" in Latin). Our previous study revealed that a meron pair (a vortex molecule) can be stabilized thermodynamically in a two-component BEC if it is rotated and has an internal coherent coupling 91 . The two components interact not only through their mean-field interactions but also through the relative phase of the order parameters. Combination of these two effects enables us to explore a new regime of rich vortex structures which are absent in the conventional binary system 78,80,81,86,87,88 . If the strength of the coupling drive is increased gradually from zero and its frequency is adiabatically tuned to be resonant, one can obtain a stationary state with a nearly equal-weight superposition of the two states 150 . The internal coherent coupling can be included in the formulation through Eq. (30), and we obtain the time-independent, coupled GP equations by substituting Ψ i (r, t) = Ψ i (r)e −iµt/ into Eqs. (31) and (32). Both components are characterized by the common chemical potential µ because, in the presence of the coherent coupling, the conserved quantity is the total particle number rather than the individual particle numbers. For simplicity, the detuning parameter ∆ is assumed to be zero. In terms of the pseudospin picture, the Rabi term given in Eq. (30) is represented as −ω R ρ T S x , which plays a role of a "transverse magnetic field" along the x-axis.
A typical solution of the vortex state with coherent coupling is shown in Fig. 8. Each component has one off-axis vortex and the density peak of one component is located at the vortex core of the other component. This results in a coreless vortex in which the total density has no singularity. Because of the coherent coupling, the profile of the relative phase φ(r) = θ 1 − θ 2 plays an important role in optimizing the structure, which is shown in Fig. 8(b). The relative phase shows that the central region is characterized as a vortex-antivortex pair. In other words, the two vortices are connected by a branch cut of the relative phase with the 2π difference, which is a characteristic domain wall structure in two-component BECs with internal coherent coupling 158 . Then, the vortex in one component attracts that in the other component, due to their being bound by the domain wall, thus forming a "vortex molecule" 91 . As ω R increases, the size of the pair decreases as seen in Fig. 9. Beyond ω R ≃ 3.0, the two vortices completely overlap despite the intercomponent repulsive interaction.
The corresponding spin texture is shown in Fig. 8(c) and (d). Because of the transverse "magnetic field", the spins are oriented in the x direction everywhere except in the central domain-wall region where they tumble rapidly by 2π. There exist two points, corresponding to the vortex centers, at which S is parallel to the z-axis. In Fig. 8(c), the spin around the singularity with S = +ẑ has a radial distribution, characterized by (S x , S y ) ∝ (−x, −y), and rolls by 90 • as it goes outward, becoming perpendicular to the x-axis. For the spin around S = −ẑ, the distribution is instead hyperbolic (S x , S y ) ∝ (x, −y). This texture is known as a "radial-hyperbolic" pair of merons, which has been discussed in the study of topological defects in superfluid 3 He (Ref. 2) and a double-layer quantum Hall system 4 . Since the trapping potential is axisymmetric, the energy is degenerate with respect to the orientation of the molecular polarization. If the molecule is polarized along the y-axis, the texture features a "circular-hyperbolic" pair, but it can continuously transform into the radial-hyperbolic pair.
We find that, in the SU(2) symmetric case, i.e., u 1 = u 2 = u 12 (c 1 = c 2 = 0), both the topological charge density q(r) and the effective velocity v eff exhibit a radially isotropic profiles like that in Fig. 5, even though each component forms a nonaxisymmetric vortex configuration. This is due to the fact that the meron pair and the axisymmetric skyrmion are equivalent as a topological excitation for the case of u 1 = u 2 = u 12 ; they can transform to each other via an overall rotation of the pseudospin. Then, giving an infinitesimal value to ω R (i.e. applying a transverse magnetic field) is sufficient to stabilize a nonaxisymmetric configuration. We thus obtain a further insight of the nonaxisymmetric vortex by rotating the basis of the spinor so that the internal coupling becomes simpler. The Rabi term plays the role of a transverse "magnetic field", making the x-axis as a preferred axis. Actually, by rotating the spinors as ψ ± = (ψ 1 ± ψ 2 )/ √ 2 (the basis along with pseudo-x axis), the coupled GP equations become Here, the internal coupling is just the chemical potential difference between the "+" and "−" components. Then, the nonaxisymmetric structure in Fig. 8 is transformed to the axisymmetric vortex state in Fig. 4, that is, a skyrmion, where, the vortex core of the "+" component is filled with the nonrotating "−" component. As one increases ω R the number of the "−" particle decreases due to a decrease in the chemical potential of the ψ − component, and thus the vortex cores eventually become empty. This corresponds to the overlap of the vortex cores as shown in Fig. 9. Therefore, the Rabi term alone does not break the axisymmetry of the topological excitation. However, inclusion of both the Rabi term and unequal coupling constants u 1 = u 2 = u 12 induces an anisotropy of the meron pair as will be discussed later.

A meron pair: Variational analysis
In the case of u 1 = u 2 = u 12 (c 1 = c 2 = 0), the spin profile of the meron pair in Fig. 8 Here, we assume that the meron pair is polarized along the x-axis. The only difference between this case and that for the skyrmion in Eq. (51) is that the forms of S x and S z have been exchanged. In this case, the ratio of the particle number given in Eq. (56) is always unity, because S z is an odd function. The locations of the vortex cores are determined by two extremes of S z , and are given by Then, it is natural to use the following form for the total phase Θ: Substituting Eq. (59) and Eq. (61) into Eq. (48) with the Rabi term −ω R ρ T S x , we obtain the following total energy: where the effective confining potential V eff is the same as that in Eq. (54). Because for any values of λ and α Eq. (62) gives the same energy as that of the axisymmetric skyrmion given by Eq. (53), the skyrmion and the meron pair have the same optimized values of λ and α if ω R = 0 (i.e., for the SU(2) symmetric case), and their energies are degenerate. Just as the c 1 -term may be regarded as a magnetic field along the z-axis, the Rabi frequency ω R may be regarded as a magnetic field along the x-axis. Thus, the magnitude of ω R also influcences the stable size of the meron pair. The difference from the c 1 -term is that the Rabi term is proportional to ρ T instead of ρ 2 T . We calculated the minimum values of λ and α as a function of ω R under a slow rotation. The result is shown in Fig. 10(a). As ω R increases, the minimized values of λ (α) becomes smaller (larger), and eventually vanishes (diverges). This behavior corresponds to a decrease in the size d m of the meron pair (see in Fig. 9 for the definition) as shown in Fig. 10(b), where d m is calculated from Eq. (60). This indicates that the binding of the meron pair becomes stronger with increasing ω R . The variational result agrees well with that obtained with direct numerical simulations.
The meron pair arises from competion between repulsive and attractive interactions 91 . For well-separated merons, the repulsive interaction between them originates from the second and third terms of Eq. (48), which are the gradient energy of the pseudospin and the hydrodynamic kinetic energy of the v eff -field. The contribution of the other terms are nearly independent of meron separation except for small d m (< 0.30 for c 0 = 10 3 ). Then, the Thomas-Fermi approximation cannot apply to the evaluation of the total energy; λ min drops suddenly to zero with increasing ω R because there is no energy barrier associated with the quantum pressure of ρ T . On the other hand, the attractive force between the merons arises from a tension T d in the domain wall of the relative phase, which is estimated to be ∼ T d d m ; for a homogeneous system T d = 8|ψ 1 | 2 |ψ 2 | 2 k/ρ T with the characteristic domain size k −1 = (|ψ 1 ||ψ 2 |/2ω R ρ T ) 1/2 (Ref. 158), and thus T d ∝ √ ω R . Thus, the competition between the repulsive force and the attractive force creates an energy minimum and the two vortices can form a bound pair.

Effect of symmetry breaking terms
Let us consider the effect of the symmetry breaking terms on a skyrmion and a meron pair. When u 1 = u 2 = u 12 (c 1 , c 2 = 0), the axisymmetry of the topological excitation is broken. First, we neglect the spin-spin interaction term (c 2 term) and investigate the dependence of the stable structure on c 1 and ω R . This is equivalent to applying a longitudinal and transverse (pseudo)magnetic field. We prepare a stable axisymmetric skyrmion by applying the longitudinal magnetic field c 1 = −10, and then turn on the transverse magnetic field ω R . Figures 11(a)-(d) show variations of the spin texture as well as those of the cross-section of the condensate density along the y = 0 line when the Rabi frequency ω R is increased. The transverse magnetic field shifts the skyrmion off-center and converts it into a meron. At the same time, another meron enters from outside and forms an "asymmetric" meron pair shown in Figs. 11(c) and (d). This is a second-order phase transition because there is no energy barrier to destabilize the axisymmetric skyrmion with respect to the transverse magnetic field. In this case, the distribution of the topological charge q(r) moves from the center as ω R increases, and there is no dramatic change in the global shape from the isotropic skyrmion of Fig. 5. As ω R increases further, the peak of q(r) moves back to the center and its value goes to infinity when the two merons merge.
Next, we study the effects of the spin-spin interaction term (c 2 term) and the dependence of the stable structure on ω R (c 1 is fixed at zero). First, we show that, in contrast to the case of c 1 , the structure of the meron pair is sensitive to c 2 . Figure  12 shows the equilibrium structure of the condensate density, the spin texture and the topological charge density q(r) for both the antiferromagnetic case (c 2 = −20) and the ferromagnetic case (c 2 = 20). For the antiferromagnetic case there is no significant difference in the density and spin profile, compared with the solution of c 2 = 0 in Figs. 8 and 9. However, the topological charge density (i.e., vorticity) is distributed anisotropically in such a way that its distribution is elongated along the direction of polarization of the meron pair. For the ferromagnetic case, spin domains are formed, which gives rise to a considerable change in the density and spin profile as seen in Fig. 12(b). Most spins align up or down along the z-axis on either side, sandwiching the domain wall across which the spin flips rapidly. If ω R is increased further, the spins on both sides lay along the x-axis, forming the well-defined meron pair as in Fig. 9. Then, the topological charge density is distributed with the direction perpendicular to the direction of polarization of the meron pair, being concentrated on the domain-wall region. This anisotropy of the meron pair gives an interesting situation when the condensates undergo a rapid rotation; the anisotropic interaction between the meron pairs generates a distorted lattice of "vortex molecules" 91 , which will be discussed in Sec. 5.3.

Topologically nontrivial skyrmions in a 3-D configuration
In contrast to the topologically trivial skyrmion discussed above, we review here the topologically nontrivial skyrmion excitations in multicomponent BECs 93−96 . These skyrmions have complicated 3-D structures and are characterized by spin inversion in a finite region of space and a typical topological charge Q = 1. It is known that topological considerations allow for these excitations, but it is not known whether or not such a configuration is energetically favorable.
The first example of a nontrivial skyrmion was discussed by Al Khawaja and Stoof 93 . They considered uniform two-component condensates subject to the SU(2) symmetry and the radially symmetric profile of the skyrmion by assuming that the normalized spinor is characterized by χ(r) = exp{−2iΩ(r) · S}χ z , where S = σ/2 with the Pauli matrix σ, χ z = (1, 0) T , and Ω(r) = ω(r)r/r which determines a specific spin texture of the skyrmion. The physical meaning of this formula is that the average spin at a position r is rotated by an angle 2Ω(r) from its initial orientation with the axis of rotation Ω(r)/Ω(r). The boundary conditions for ω(r) that represents the spin distribution for the skyrmion are the following. First, at r → ∞ all spins must be oriented as in the ground state, that is, ω(r → ∞) = 0. Along the z-axis the spins are also not rotated due to our assumption about Ω; thus, to have a nonsingular texture of the spinor with a non-zero winding number, we must require ω(0) = 2π. Finally, to avoid the singular behavior of Ω(r) itself, we use only functions ω(r) with zero slope at the origin. Therefore, ω(r) is a monotonically decreasing function that starts from 2π at the origin and reaches zero as r → ∞.
Al Khawaja and Stoof 93 used a variational ansatz ω(r) = 4 cot −1 (r/λ) 2 that satisfies the above boundary conditions, where λ represents the size of the skyrmion. Substituting this ansatz into the energy functional in Eq. (14) with V i = 0 and Ω = 0 gives the total energy as a function of the total density ρ T and λ. For the given ω(r), they calculated ρ T (r) exactly by solving numerically the differential equation for ρ T (r) obtained by varying E[ρ T (r), λ] with respect to ρ T (r). By substituting this density profile back into the energy functional, the energy of the skyrmion is expressed as a function of only λ and the equilibrium properties can be obtained by minimizing this energy with respect to λ. The resulting spin texture shows a complicated spin distribution; the spin S z flips five times as it goes outward from the center in the x-y plane, and S x and S y exhibit quadrupole profiles with one radial node in the y-z plane 93 .
A crucial problem is the stability of the skyrmion. In equilibrium, the energy is minimized for the vanishing size of the skyrmion 93 . However, for sufficiently small skyrmion sizes, a nonequilibrium stability mechanism starts to work. A number of atoms in the center of the skyrmion, hereafter referred to as core atoms, will be trapped by an effective 3-D potential barrier, that is, a repulsive shell with a finite radius that is created by the gradient term of the spin texture 2 |∇χ(r)| 2 /2m. As the skyrmion shrinks in size, the barrier height of the repulsive shell increases and the radius decreases. This leads to a squeezing of the core, thereby increasing its energy, and ultimately stabilizing the skyrmion. In this case, it is not an equilibrium state of the condensate because the core atoms will tunnel through the barrier and give the skyrmion a finite lifetime. Detailed calculations of this tunneling rate 93,159 showed that only a few atoms in the core of the skyrmion are needed to stabilize it and to give it a sufficiently long lifetime.
Other interesting structures of nontrivial 3-D skyrmion were investigated in Refs.
94, 95, 96. The authors in Refs. 94, 95, 96 also considered the rotational operator of the spinor as χ(r) = exp{iλ(r)σ · (r/r)}χ z , where λ is a monotonic function imposed by the boundary condition for λ as λ(0) = 0 and λ = π at the boundary of the condensate; for this boundary condition the topological charge becomes Q = 1. The corresponding condensate wave functions are Here, (λ, θ, φ) can be understood as the spherical angles of the 3-sphere (in fourdimensional parameter space). Because of the topological stability of the skyrmion, any continuous deformation of Eq. (63), without altering the asymptotic boundary values, is still a skyrmion with Q = 1. With this realization of the skyrmion, |ψ 2 | 2 vanishes both on the z-axis and at infinity, and is concentrated in a toroidal region.
On the other hand, |ψ 1 | 2 vanishes on the circle θ = π/2 and r = λ −1 (π/2), a region where |ψ 2 | 2 is concentrated. Expanding ψ 1 about any point on this circle, we find ψ 1 ∝ δr + iδθ, indicating that the nodal circle of ψ 1 represents a vortex line. Hence the ψ 1 component forms a vortex ring, within whose core ψ 2 resides and flows azimuthally with the winding number of one. The skyrmion can be viewed as a quantized vortex ring in one component, which is filled with the other component that also carries the quantized circulation along the ring. This is shown in Fig. 13. This configuration closely resembles the "cosmic vortons" which might have been formed in the early universe from superconducting cosmic strings 160 .   95 . In this case, the two components can strongly repel each other and the toroidal filling due to the line component can prevent the vortex ring from shrinking to nothing. Consequently, a filled vortex ring, as opposed to an empty vortex ring, can be more stable against collapse. Moreover, in the skyrmion the filling has one unit of angular momentum about the z-axis, resulting in a 1/r 2 centrifugal barrier that further hinders the ring from shrinking. Savage and Ruostekoski extended the analysis to a harmonically trapped system and found that there are additional instabilities that will play a crucial role in experiments 96 . For sufficiently large total number of atoms in the vortex-line component, the nonlinear repulsion between the two components is strong enough to inhibit the collapse of the vortex ring, stabilizing the skyrmion against shrinking for cylindrically symmetric initial states. However, due to the inhomogeneous potential, the skyrmion is still unstable with respect to the drift of the vortex line towards the edge of the BEC. Once the vortex line drifts to a low-density region, the nonlinear repulsion is no longer strong enough to prevent the shrinking of the skyrmion and the collapse occurs as described above. The drift reduces the total angular momentum of the atoms, indicating an energetic instability. Physically, this results from dissipation (e.g., from thermal atoms). Thus, additional physical mechanism, such as rotation or optical pinning potentials, will be required for stabilization.
Rotation of only the vortex-line component gives the simplest example for stabilization, but this is difficult. Moreover, rotation of both components introduces a new instability mechanism. For sufficiently rapid rotations the line component again reaches the boundary, resulting in Q < 1. At the same time the ring vortex accommodated the rotation by twisting and finally breaking at the condensate boundary. Nevertheless, for a small range of rotation frequency near 0.085ω ⊥ , the skyrmion was stable against both shrinkage of a vortex ring and drift of a vortex line. Then, the stable skyrmion configuration requires a shift of the rotation axis to off-center . The breaking of cylindrical symmetry implies a family of degenerate skyrmions parametrized by the azimuthal angle. Another method for stabilizing the skyrmion is to inhibit the drift towards low-density regions by creating a positive density gradient around the vortex line. This can be implemented with a bluedetuned Gaussian laser beam along the z axis, providing a cylindrically symmetric repulsive Gaussian dipole potential perpendicular to the z axis. Under these conditions, we expect that atomic BECs, which have been imprinted with a topological configuration, will relax via dissipation to a stable skyrmion configuration.

Dynamic properties of vortices
The dynamic properties of vortices in two-component BECs are also expected to be extremely rich because of the various kinds of interactions between the two superfluids. A remarkable feature of this system arises from "buoyancy" 161,162,85 , which refers to a net mean-field force on a constituent of a multicomponent condensate, the force being dependent on the various intraspecies and interspecies interaction parameters.
As shown in Sec. 4.1, the unstable dynamics of the vortex was observed, depending on which component has a vortex in the 87 Rb condensate 11 . This phenomenon is understood from linear stability analysis of the vortex state based on the mean-field approach 82,83,84 . We use the following notation for the states: (1, 0) for the state with the vortex in |1 and (0, 1) for the state with the vortex in |2 . In the JILA experiment, two components had the same number of particles (N 1 = N 2 = N ) but, in general, one could consider any ratio between the populations of the different hyperfine states. For 87 Rb atoms, the relation of the coupling constants g 1 > g 12 > g 2 (a 1 : a 12 : a 2 = 1.03 : 1 : 0.97) favors a configuration that has its first component spread over the largest part of the space, as shown in the previous discussion of equilibrium. Numerical simulations show that for equal populations N 1 = N 2 = N , and with the above coupling constants, the stationary states (1, 0) is stable while the other state (0, 1) is unstable 82,83,84 .
The origin of the instability of the state (0, 1) is purely dynamical, which can be understood within the framework of the mean-field theory for the two-component system without dissipation. Actually, the instability does not lead to expulsion of the vortex from the condensate, but instead to periodic transfer of the vortex from one component to the other. To study the vortex stability, one starts from the timedependent, coupled Gross-Pitaevskii equations for the condensate wave functions of each species: These equations are the particular case of Eqs. (31) and (32) in which the drive is turned off (Ĥ 1 , ω R , ∆ = 0). Equations (64) and (65) conserve the number of particles in each hyperfine state. However, the angular momentum of each component is no longer a conserved quantity. Instead, the conserved quantity is the total angular momentum L z = d 3 r(Ψ * 1 L z Ψ 1 + Ψ * 2 L z Ψ 2 ). As in the JILA experiments, both potentials are assumed to be spherically symmetric and have the form V 1 (r) = V 2 (r) = 1 2 mω 2 ⊥ (r 2 + z 2 ). For stationary configurations in which each component has a well-defined value of the angular momentum, the time and angular dependences of the wave functions are factored out as Ψ i (r, z, φ) = e −iµit/ e iγiθ ψ i (r, z) (i = 1, 2).
Here, we focus on the two axisymmetric configurations with vorticity (γ 1 , γ 2 ) = (1, 0) and (0, 1), which correspond to the single vortex states for the |1 and |2 components shown in Fig. 3. Linear stability analysis of the two axisymmetric vortex states leads to the following results. Among the normal modes of the (1, 0) case, there is a negative eigenvalue, which means that there is a path in the configuration space along which the energy decreases. This path, which is associated with the mode of the dipole type, belongs to a perturbation that causes the initial displacement of the vortex from the trap center. However, this drift instability of the vortex is affected by dissipation; in particular, without dissipation, the configuration is dynamically stable. But, it can be energetically stabilized by turning on rotation as discussed in Sec. 4.2. In the (0, 1) case, on the other hand, there are normal modes with complex frequencies. The unstable modes have a similar shape of the negative-eigenvalue modes of the (1, 0) case, which means that the perturbations also displace the vortex from the center. The imaginary part of the eigenvalues implies that vortices with unit charge in |2 are dynamically unstable under a generic perturbation of the initial configuration and this instability grows without dissipation. This result is consistent with the JILA experiments in Fig. 2, in which a vortex in the |2 component was unstable and collapsed.
Numerical simulations of the vortex dynamics for large perturbations were done by García-Ripoll and Pérez-García 82 . The results showed that the linearly stable state (1, 0) is robust and survives under a wide range of perturbations; the vortex only shows a precession around the center. In contrast, the unstable configuration (0, 1) gives rise to recurrent dynamics: the first component and the vortex oscillate synchronously, where the vortex core in |2 pins the density peak of |1 . These oscillations grow in amplitude until the vortex spirals out. The displaced vortex carries less than the total angular momentum, and so the vortex precesses around the center. Because of this and the conservation of the total angular momentum, the vortex is transferred from |2 to |1 . Although the dynamics is not completely periodic, this mechanism exhibits some recurrence and the vortex eventually returns to |2 . The preceding behavior persists even for strong perturbations in a 2-D condensate. However, for large perturbations of a 3-D condensate, the dynamics may be chaotic.
The preceding results are valid when N 1 = N 2 . For an arbitrary ratio of the populations N 1 /N 2 and arbitrary values of the coupling constants g 1 , g 2 and g 12 , the dynamical stability conditions can be estimated analytically 83,84 . Under the two-mode approximation 83 , where the wave function is approximated as Ψ i ≃ a i (t)Ψ g (r) + b i (t)Ψ e (r) with the eigenfunctions of the d-dimensional harmonic oscillator Ψ g (r) = (1/π) d/2 e −r 2 /2 and Ψ e (r) = (2/dπ) d/2 re −r 2 /2 e iθ , the configuration (1, 0) is stable if For 87 Rb, the inequality (67) is always satisfied because the right-hand side is negative, which proves that the configuration with a vortex in |1 is always linearly stable. Note that the stability properties do not depend on the total number of particles but only on the ratio between the populations. The stability condition of the configuration (0, 1) can be obtained in a similar manner as This inequality fails for a certain range of N 1 /N 2 . For the case of 87 Rb the unstable range is 0.68 < N 2 /N 1 < 1.38, which means that certain choices of the population imbalance allow stabilization of the vortex in |2 . More detailed calculations were done by Skyrbin 84 who derived the stability condition as which gives an unstable range of 0.63 < N 2 /N 1 < 1.69, which is slightly different from Eq. (68). These results predict the possibility to dynamically stabilize the axisymmetric vortex states for various multiple-condensate systems. The effect of the intercomponent interaction was also demonstrated in a precession of an off-centered vortex core 161 . This is a consequence of a Magnus effect; that is, an applied force on a rotating cylinder in a fluid causes pressure imbalances at the cylinder surface that makes the cylinder drift orthogonal to the force. The experiment observed that the precession frequency of the filled core is slower than that of the empty core 161 . The slower precession of filled cores can be understood in terms of the buoyancy effect. Because of its slightly smaller scattering length, the |2 component has negative buoyancy with respect to the |1 component, and consequently tends to sink inward towards the center of the condensate. With increasing amounts of the |2 component in the core, the inward force on the core begins to counteract the outward buoyancy of the vortex velocity field, resulting in a reduced precession velocity. It is predicted that with a filling component of sufficiently negative buoyancy in the core, the core precession may stop or even precess in the direction opposite to the direction of the fluid flow 162 .

Quantized vortices in spin-1 Bose-Einstein condensates
Another multicomponent Bose-condensed system is a "spinor BEC" characterized by the hyperfine spin F = 1. An optical trap removes the restriction of the confinable hyperfine spin states, realizing new multicomponent superfluids with internal degrees of freedom 71−76 . In this section, we review briefly the study of quantized vortices in such a BEC with |F| = 1 hyperfine spin states |F = 1, m F = ±1, 0 .
To describe the F = 1 spinor BEC, a generalized GP mean-field model was introduced by Ohmi and Machida 97 and Ho 98 . With the mean-field approximation, 2vorrev Vortices in multicomponent BECs 39 the total Hamiltonian of the system reduces to an energy functional similar to Eq. (14): Here we consider the system in a rotating frame with Ω = Ωẑ and the subscripts represent α = (x, y, z) and i, j, k, l = (0, ±1). The condenate wave functions are characterized by the hyperfine sublevels m F = 1, 0, −1 as Ψ = (Ψ 1 , Ψ 0 , Ψ −1 ). If one of the component Ψ i is fixed at zero, Eq. (70) reduces to the energy functional of the two-component case Eq. (14); the crucial difference from the two-component problem is the presence of the g s -term. The coupling constants g n and g s are written as where a 0 and a 2 are the scattering lengths for two colliding atoms with total spin angular momentum F tot = 0 and 2. The spin operatorsF α (α = x, y, z) can be expressed in matrix form as Note that the term proportional to g s in Eq. (70) describes the spin exchange interaction, which is ferromagnetic for g s < 0 and antiferromagnetic for g s > 0. The spinor BEC of F = 1 23 Na atoms studied by the MIT group was demonstrated to have an antiferromagnetic interaction 71 , while the F = 1 87 Rb condensate has a ferromagnetic one 75,163 . The sign of g s plays a crucial role in determining the stable vortex structure. The time-independent coupled GP equations are derived by the requirement that the energy functional (70) be extremal: where the three chemical potentials µ i (i = 0, ±1) play the role of Lagrange multipliers. They are determined so as to conserve the total particle number N = dr(|Ψ 1 (r)| 2 + |Ψ 0 (r)| 2 + |Ψ −1 (r)| 2 ) and magnetization M = dr(|Ψ 1 (r)| 2 − |Ψ −1 (r)| 2 ). The angular momentum L z = i drΨ i L z Ψ i should also be considered as a conserved quantity.

Coreless vortex states
First consider the structure of axisymmetric vortex states. The detailed analysis was made by Isoshima and Machida 103 . In this structure, the vortex configuration is classified by the combination of the winding number γ i of the condensate wave function Ψ j = ρ j (r) exp[i(α j +γ j θ)] with the cylindrical coordinate (r, θ, z), where the homogenity of the wave functions along the z-axis is assumed, and α j is the overall phase of the j-th component. The phase factors of Ψ j are determined such that the energy of the g s term in the integrand of Eq. (70) is minimized; this term reads To minimize E s , α j and γ j should satisfy where n is an integer, and the odd (even) n corresponds to the antiferromagnetic (ferromagnetic) interaction. The global phase α j has no effect on the vortex structure, and will therefore be set to zero. The possible combination of γ j satisfying Eq. (78) gives this system a characteristic vortex structure. If the values of γ j are restricted to being γ j ≤ 1, we have (γ 1 , γ 0 , γ −1 ) = (1, 0, −1), (1, 1, 1) and (1, 1/2, 0). In the case of (1, 1/2, 0), the value 1/2 is not allowed in the axisymmetric system, so that the Ψ 0 -component vanishes; this vortex state is called the "Alice state" (halfquantum vortex state) 79,103 . Therefore, there is one-to-one correspondence between the Alice state and the axisymmetric vortex state in two-component BECs discussed in Sec. 4.2 (i.e., Ψ 1 (Ψ −1 ) → ψ 1 (ψ 2 ) in Fig. 4). The structure and the stability of these vortex states for various Ω and the total magnetization M were thoroughly investigated in Ref. 103. Also, the winding number γ j can exceed unity. Mizushima et al. discussed the (0, 1, 2) vortex shown in Fig. 14 and showed that the ferromagnetic interaction supports the thermodynamic stability of this structure 104 . The central region of the harmonic trap is occupied by the non-rotating Ψ 1 component. The Ψ 0 component with γ 0 = 1 is pushed outward, while the Ψ −1 component with γ −1 = 2 occupies the outermost region. The resulting total density is nonsingular. This configuration is favorable for the ferromagnetic case, because it is favorable to spatial phase separation; the presence of vortices with different γ j effectively cause phase separation in the radial direction as shown in Fig. 14 Putting this coreless vortex into the pseudospin representation reveals that it has the structure of a spin texture. The wave functions of Fig. 14(a) can be parametrized as where the bending angle β(r) has the range 0 < β(r) < π. The spin direction is denoted by thel-vector and is given asl(r) =ẑ cos β(r) + sin β(r)(cos θx + sin θŷ) [the contour density of l z is shown in Fig. 14(b)], where β varies from β(0) = 0 to β(R) = π/2 (= π) for a Mermin-Ho (Anderson-Toulouse) vortex (R is the outer boundary of the cloud). In other words, l z = cos β(r) varies from l z (0) = 1 to l z (R) = 0 (−1), wherel-vector points up at the center and points outward (downward) at the circumference for the Mermin-Ho (Anderson-Toulouse) vortex. In the superfluid 3 He-A phase, the stability of these vortices is guaranteed by the constraint that thel-vector is perpendicular to the vessel wall 2 . There is no such boundary condition in the trapped BEC system. In this case, the (0, 1, 2) vortex of Fig. 14 can be interpreted as a vortex with an intermediate boundary condition (π/2 < β(R) < π); we can control the value of β(R) from the Mermin-Ho condition to the Anderson-Toulouse one by merely changing the total magnetization M 104 .
In another numerical study, Martikainen et al. analyzed the coreless vortex state as a function of rotation frequency without fixing the total magnetization, and found that β(R) increases with Ω and the upper value of β(R) is 3π/4, above which additional vortices nucleate. This implies that a Anderson-Toulouse vortex can never be the ground state of the system. Full 2-D numerical simulations of the coupled GP equations (73)-(75) also reveal non-axisymmetric structures of vortices, whose stability depends sensitively on the spin-exchange interaction g s , the external rotation Ω and the total magnetization 99,104,105 . For example, the axisymmetric (1, 1, 1) vortex is always unstable, leading to the displacement of the vortex cores from the center. In the antiferromagnetic regime, the vortex cores are displaced such that the condensates has two singularities, through an overlap of the Ψ 1 and Ψ −1 components or it has three singularities to form a triangle configuration. The breaking of the axisymmetry mainly leads to a decrease in the overlapping area of the condensate wave functions. In comparison with axisymmetric states, these non-axisymmetric vortices have the advantage that they can easily adjust to changes in Ω. As Ω increases, the two or three separate singularities adjust their interdistances from the center and change the value of L z to gain energy due to the term −ΩL z . In this sense, the non-axisymmetric vortices are more flexible for a change in Ω compared with the axisymmetric ones.

Monopoles
Another interesting topological excitation that occurs in the spin-1 BEC system is a singular pointlike defect, called a "monopole" due to its similarity to magnetic monopoles in particle physics 164 . A monopole is characterized by a unit vector that is radial with respect to a unique central point (i.e., the "hedgehog" structure). Stoof et al. 113 proposed monopole excitations in an antiferromagnetic spin-1 BEC. The monopole is characterized by the spinor where θ is the overall phase. The form of Eq. (80) originates from the fact the antiferromagnetic interaction implies zero average spin | F | ≡ΨFΨ/ρ T = 0. The spherically symmetric monopole is described by the radial unit vector m(r) = r/r. Minimization of the energy of the symmetric configuration yields θ = 0 and the density singularity at the origin ρ T (0) = 0. Then, the spinor component Ψ 0 resembles a dark soliton and Ψ ±1 form perfectly overlapping, straight singly-quantized vortex lines with opposite circulation, the latter being perpendicular to the phase kink plane. Stoof et al. demonstrated that this particular spin texture is a consequence of the unit winding number and the minimization of the gradient energy 113 . This configuration could be created through the phase imprinting technique achieved by the carefully designed sequences of laser pulses 115 . The theoretical study by Ruostekoski and Anglin 116 demonstrated that the spherically symmetric monopole is stable only in the strongly antiferromagnetic regime. This is due to the following two length scales They represent the length scales over which ρ T and | F |, respectively, tend to their bulk value when subjected to localized perturbations. In the weakly antiferromagnetic regime, the wavelength (∼ ξ a ) at which the antiferromagnetic constraint may be violated is much larger than that (∼ ξ s ) at which the total density constraint fails. Then, the monopole core extends to a larger size, which is favorable for having nonzero average spin instead of a density zero. The monopole point defect may be continuously deformed into a circular line defect without changing the topological invariant. This is done by punching a hole in the spherically symmetric core. To keep Ψ single-valued on the disc bounded by the ring, the macroscopic phase θ must change by π around any loop that links the defect circle, while there is also a π-disclination in m on the disc. This structure is identified as a half-quantum vortex line, forming a closed circular ring, called an "Alice ring" [see Fig. 15]. To keep Ψ single-valued on the ring itself, one can either have ρ T vanish there, or have | F | = 1 on the ring instead. By decreasing the ratio g s /g n , the spherically sym-metric monopole core with ρ T (0) = 0 becomes unstable against formation of an energetically stable Alice ring with | F | = 1 core. According to a linear stability analysis of this case, the radially symmetric monopole is unstable in a homogeneous condensate for g s /g n 0.17, while according to the numerical simulation for a trapped condensate the monopole is unstable for g s /g n 0.16, using an experimentally feasible set of parameters for 23 Na (ω ⊥ = ω z = 2π × 10 Hz and N = 4 × 10 6 ).
The creation of the monopole is not limited to the antiferromagnetic spinor condensates. Martikainen et al. studied in detail the two-component monopole charac- To ensure the stability of this texture against phase separation, the spin-1 condensate must have ferromagnetic interactions, which makes the 87 Rb spinor condensate a potential candidate. In other words, contrary to a previous expectation, the preparation of a monopole is not restricted to the antiferromagnetic 23 Na condensate because one can use spinors that do not have the order-parameter space of the ground state.

Vortex states in rapidly rotating two-component BECs
Rapid rotation of a BEC generates a large number of quantized vortices, with the result that the condensate can mimic a rigid body rotation despite the irrotational nature of superfluidity 1 . It is well known that the most stable configuration of an assembly of quantized vortices in a rotating single-component superfluid ( 4 He-II or atomic gas BECs) is a triangular lattice; also it is known that this is analogous to an Abrikosov lattice of magnetic fluxes in a type-II superconductor. At present, most experiments on vortex lattices in atomic-gas BECs have been done in single-component systems. Thus, a natural question arises: what happens in rapidly rotating two-component BECs? The vortex lattices in such systems are more complicated than those in single component condensates because of the intercomponent interactions. This section reviews the studies of vortex states in two-component BECs with a large number of vortices. We also describe a recent experiment 165 that found interlaced square vortex lattices in two-component BECs. First, we show results obtained through analysis based on the "mean-field quantum Hall regime" 87 . This regime will be valid when the BEC in a harmonic trapping potential is rotated with angular frequency Ω close to the radial trapping frequency ω ⊥ 43,166 . This method has the great advantage that one can analytically describe the vortex structure, although the method is effective only near the limit Ω ≃ ω ⊥ . Next, we show the results of the direct numerical calculation of the coupled GP equations (21) and (22) 88 . They reveal a rich variety of vortex structures that have eluded analytic treatment. Depending on the strength of the intercomponent interaction and the rotation frequency, vortex lattices can have the form of triangular lattices, square lattices, double-core vortex lattices, and "serpentine" vortex sheets. These properties can be understood by using the concept of "pseudospin", in which the vortex lattices are interpreted as a lattice of skyrmion. Finally, we discuss the effect of internal coherent coupling on the lattice structure, caused by the anisotropy of the resulting vortex molecules (meron pairs).
Similar discussions can also be done for spin-1 BECs. Kita et al. studied vortexlattice structures of antiferromagnetic spin-1 BECs using the phenomenological Ginzburg-Landau equation 106 , which is formally similar to the mean-field quantum Hall approach. Their study revealed that the conventional Abrikosov lattice with hard vortex cores is unstable and the vortex cores shift their locations. The system has many metastable configuration of vortices, depending sensitively on the ratio g 2 /g 0 in Eq. (70). The vortices are characterized by distribution of magnetization and the difference in the number of circulation of quanta per unit cell, all of which makes the characterization of the vortex lattice structure complicated. Further discussions on the phase diagram of the ground state in rotating spin-1 bosons are reported in Reijnders et al. 109 . This paper also describes a study of vortex-lattice structures with the mean-field quantum Hall approach and the exact diagonalization study of the quantum liquid phase. In addition, Mizushima et al. did numerical simulations of the GP equations (73)-(75) for a ferromagnetic spinor BEC in a fast-rotating regime and also in an external magnetic field 111 . As a result, they could constuct a phase diagram of stable configurations on a plane of the rotation frequency and the total magnetization.

Mean-field quantum Hall approach
The "mean-field quantum Hall regime" is the regime obtained in a rapidly rotating limit, where mean-field theory remains valid. Thus, each component can still be characterized by a condensate wave function Ψ i (i = 1, 2). The angular momentum of the system is so high that Ψ i is made up of the orbitals in the lowest Landau level (LLL) in the plane perpendicular to the rotation axis. Ho showed that this regime will emerge in a 3-D Bose gas at sufficiently high angular momenta 166 , and recently the JILA group has achieved this regime experimentally 43 . (See also the recent detailed arguments in Ref. [167][168][169][170] for the use of LLL approximation and the structure of a vortex lattice.) In this regime, the wave function acquires an analytic structure that allows exact evaluation of the energy of a vortex lattice. As a result, it is possible to evaluate a wide range of lattice structures in rapidly rotating condensates, a feat that would be impractical for numerical calculations because of the time and the accuracy required.
An extension of this argument to the two-component BECs was made by Mueller and Ho 87 for the following situation. If each component has an equal number of bosons, and if the trapping potentials of the two components are identical, the two components will have the same size and the same density of vortices. In this case, one expects that each component will contain identical vortex lattices, with one lattice displaced relative to the other because of the presence of the intercomponent repulsive interaction. We shall refer to such a configuration of vortex lattices as staggered vortex lattices. While in the current experiment, the coupling constants satisfy g 1 ∼ g 2 ∼ g 12 , considerable insights are gained by studying vortex lattices for a wider range of coupling constants. Mueller and Ho considered the case g 1 ∼ g 2 = g 12 and investigated a wide range of vortex lattice structures as a function of δ = g 12 / √ g 1 g 2 , and we summarize their result here. The condensate wave functions Ψ 1 and Ψ 2 of a rotating two-component BEC are determined by minimizing the free energy F = E − µ 1 N 1 − µ 2 N 2 , where E is the total energy functional (14), µ i (i = 1, 2) are the chemical potentials that determine the number of bosons N i in each component. For simplicity, we assume the same cigar-shaped trapping potential V = m(ω 2 ⊥ r 2 + ω 2 z z 2 )/2 with ω ⊥ ≫ ω z for each component. As discuss in Ref. 166, a slow variation of the trapping potential along the z-axis allows one to apply a Thomas-Fermi approximation for the z dependence of Ψ i and write F as with µ i (z) = µ i − mω 2 z z 2 /2. As z is treated as a parameter, it is convenient to write Ψ i = ρ i (z)ψ i (r, z) with |ψ i (r, z)| 2 d 2 r = 1; hence, the normalization condition d 3 r|Ψ i | 2 = N i is now written as dzρ i (z) = N i . The first term of the integrand in Eq. (83) h L = (−i ∇ − mω ⊥ẑ × r) 2 /2m is identical to the Hamiltonian that describes the motion of an electron in a magnetic field Bẑ in the symmetric gauge (i.e., (−i ∇ − eA/c) 2 /2m with a vector potential A = Bẑ × r/2). The eigenstates of h L are referred to as the Landau levels and described in polar coordinates (r, θ) as where b ho = /mω ⊥ , C nm is the normalization constant, and L n m (x) the associated Laguerre polynomial. Since u n,m are also eigenstates of L z with the eigenvalue (m − n), the eigenenergy of h L − (ω ⊥ − Ω)L z reads where n is the index of the Landau level, and m is the index that classifies the degenerate states in the n-th Landau level. This degeneracy with respect to m is lifted by the interaction. When Ω approaches ω ⊥ , the condensate density becomes so dilute because of the centrifugal effect that the mean-field interaction energy decreases. If the interaction energy becomes less than the energy spacing of the Landau level, then the wave functions ψ i can be written in terms of the orbitals u 0,m (r) in the LLL in the xy-plane, The energy functional then becomes where r 2 i = r 2 |ψ i | 2 d 2 r. We next introduce a variational wave function that represents a vortex lattice. Generally, wave functions in the LLL can be written as φ(r) = f (z)e −r 2 /2b 2 ho , where f (z) is an arbitrary polynomial of z = x + iy (but not of z * ). The form of a vortex lattice can be described by where C is a normalization constant and {b α } are the zeros of φ, which are the positions of vortices. In Eq. (88), we assume perfect periodicity of the lattice. If the zeros form an infinite lattice with unit cell size v c , the density |φ| 2 can then be written as a product of a Gaussian and a function which is periodic under lattice translation: where R = n 1 B 1 + n 2 B 2 , n i are integers, and B 1 , B 2 are the basis vectors of the lattice. The width σ includes the contribution of the number of vortices of the system, which is determined below. It is convenient to introduce the complex basis vectors If the lattice is oriented so that b 1 is real, we can parametrize B 1 = b 1x and B 2 = b 1 (ux + vŷ), and obtain b 2 = b 1 (u + iv) and v c = b 2 1 v. Except for the normalization constant, the form of Eq. (88) is similar to the Weierstrass's σ-function 171 where the term with m = n = 0 is omitted from the product and Ω mn = 2mω 1 + 2nω 2 with the half-periods ω 1 and ω 2 . By using a relation between the Weierstrass's σ-function and the Jacobi θ-function 171 we can write f (z) = θ(ζ, τ )h(ζ), where ζ = (x + iy)/b 1 =x + iȳ, τ = u + iv = b 2 /b 1 , and h(ζ) is an entire function without zeros. To ensure the normalization of φ, the function h(ζ) must be of the form h(ζ) = exp(c 1 ζ + c 2 ζ 2 ). The density of the system is then |φ(r)| 2 = |θ(ζ, τ )| 2 |e c1ζ+c2ζ 2 | 2 e −r 2 /b 2 ho . For a vortex lattice with inversion symmetry about the origin r = 0, we have c 1 = 0. In addition, if the shape of the condensate is cylindrically symmetric, we have c 2 = π/(2v c ), which gives |φ(r)| 2 = |θ(ζ, τ )| 2 e −r 2 /σ 2 with 1 Following Eq. (89), we can write g(r) as It follows from the quasi-periodic properties of the Jacobi theta-function that g(r) = g(r + R). The periodicity of g(r) implies g(r) = v −1 c K g K e iK·r , where {K} are the reciprocal lattice vectors. K = m 1 K 1 + m 2 K 2 , and K i are the basis vector of the reciprocal lattice, The explicit form of g K is calculated in a straightforward way. We obtain where The factor (1 − e iπ(m+m ′ ) )/2 is added in Eq. (98) because m + m ′ must be odd. By applying the Poisson summation formula to Eq. (98), we have 2vorrev Vortices in multicomponent BECs 49 We thus have where r = xx + yŷ, K = (2πmx − 2π(n + um)/vŷ)/b 1 , and If g 1 = g 2 , the two components with equal particle numbers and trapping potentials have identical vortex lattices. A sufficiently small difference in g 1 − g 2 should not change this structure. If the intercomponent interaction is repulsive g 12 > 0, the mean-field interaction energy is minimized by interlacing the two lattices; if the vortex lattice of one component is deformed, the other has to follow to keep the interaction energy minimal. Under these assumptions, the normalized condensates ψ 1 and ψ 2 in Eq. (87) would be represented by densities of the form The wave function is described by variational parameters ρ i (z), σ 2 , the basis vectors B i (which determine the unit cell size v c ), and the relative displacement r 0 between the two lattice. By integrating Eqs. (103) and (104), one would see that terms up to those of relative order v c /σ 2 , the radius of the condensate is r 2 1 = r 2 2 = σ 2 . Defining the quantities I and I 12 as d 2 r|Φ i | 4 ≡ I/(πσ 2 ) and d 2 r|Φ 1 | 2 |Φ 2 | 2 = I 12 /(πσ 2 ), we have and the energy functional F takes the form (n 1 + n 2 ) + 1 2πσ 2 n 2 1 g 1 I + n 2 2 g 2 I + 2n 1 n 2 g 12 I 12 .
Since we assume here a large number of vortices, the size of the condensate is much larger than the unit cell and thus πσ 2 /v c ≫ 1. We can therefore ignore all K + K ′ = 0 terms in Eqs. (106) and (107), since σ 2 K 2 > πσ 2 /v c . We then have where g K is given by Eq. (102) and the sum of K is over the integers m 1 and m 2 . Since the expressions of I and I 12 in Eq. (109) are independent of σ 2 , the minimization of F in Eq. (108) with respect to σ 2 and ρ i is greatly simplified. The optimized σ 2 and ρ = (ρ 1 , ρ 2 ) are given by where G = g 1 I g 12 I 12 and the corresponding free energy F is given by It is clear from Eqs. (110) through (113) that the solution for the case where g 1 − g 2 ≪ |g 1 + g 2 | is very close to that of g 1 = g 2 . The lattice shape (parameterized by r 0 , u, and v) enters the grand potential only through the factor 1 · G −1 · 1. When g 1 = g 2 this term is inversely proportional to J = I + δI 12 (δ = g 12 /g 1 ), and the most favorable lattice is the one that minimizes J.
In the minimization it is convenient to measure lengths in unit of the the basis vector B 1 = b 1x and also to write the complex representation of B 2 and r 0 as τ = u + iv = |τ |e iη and r 0 ≡ a + bτ , respectively. The phase diagram of the vortex lattice as a function of the ratio δ = g 12 /g is shown in Figs. 16 increases continuously to 90 • as δ increases to 0.372. As a result, the lattice type is no longer fixed and the unit cell is a rhombus. The modulus of τ , however, remains fixed across this region. • (d) 0.373 < δ < 0.926: The two lattices are "mode-locked" into a centered square structure throughout the entire interval (τ = i, a = b = 1/2). • (e) 0.926 < δ: The lattice type again varies continuously with interaction δ. Each component's vortex lattice has a rectangular unit cell (η = π/2) whose aspect ratio |τ | increases with δ. Both 87 Rb and 23 Na have interaction parameters within this range. At δ = 1, (g 1 = g 2 = g 12 = g), the aspect ratio is √ 3. For indistiguishable components, the combined lattice would then appear triangular.
As described in Sec. 3, the condition of phase separation for nonrotating condensates is given by Eq. (29) (i.e., δ > 1). The presence of a vortex lattice naturally modulates the density of each component, with the high density regions of one component overlapping with the low density regions of the other. Thus, the system is effectively phase-separated whenever staggered vortex lattices are present, even for δ < 1. In particular, the vortex lattice near δ = 1 (Fig. 16(e)) is made up of alternating rows of vortices of each component, and the system therefore contains stripes in which one component has a high density and the other component has a very low density.

Numerical simulations for rotating two-component BECs in the Thomas-Fermi regime
This section addresses a full numerical study of the GP equations on the vortex lattice structure in two-component BECs 88 . The assumptions in the LLL calculation shown in the last subsection are that (i) the rotation frequency is assumed to be very close to the radial trapping frequency and that (ii) a vortex lattice is assumed to be perfectly regular 110 . Therefore, a vast experimentally accessible parameter regime of the lattice configuration remain to be explored. The full numerical simulations of Eq. (46) and (47) reveal a rich variety of vortex structures which have eluded the analytic treatment. Here, we show some detailed results of our simulations. We numerically calculated the equilibrium solutions of Eq. (46) and (47) through the conjugate gradient method 23,172 . Following Mueller and Ho, we also assume an identical trapping potential for both components (V 1 = V 2 = V ) and the identical intracomponent interactions g 1 = g 2 . We confine ourselves to the 2-D problem. The dimensionless coupling constants (defined in Eqs. (46) and (47)) are u 1 = u 2 = u = 4000 under the assumption N 1 = N 2 ; the normalization of each wave function is taken as d 2 r|ψ i | 2 = 1/2. The intercomponent interaction u 12 and the rotation frequency Ω are treated as variable parameters. Figure 17 shows the numerically obtained phase diagram of vortex states in the (δ ≡ u 12 /u, Ω) plane. The upper limit of the rotation frequency is Ω = 1, because for Ω > 1 the centrifugal potential dominates the harmonic trapping potential and therefore a BEC cannot be trapped. Also, we do not discuss the slow rotation regime below Ω = 0.40. The criterion of phase separation δ = 1 is reflected on the structure of vortex states, as explained below.
In the overlapping region δ < 1, two types of regular vortex lattices are obtained as the equilibrium state. For δ = 0, the two components do not interact; therefore, triangular vortex lattices form as that in a single component BEC. As δ increases, the positions of vortex cores in one component gradually shift from those of the other component and the triangular lattices become distorted. Eventually, the vortices in each component form a square lattice rather than a triangular one, which is consistent with the result of the LLL calculation done by Mueller and Ho 87 . However, Fig. 17 shows that the stable region of the square lattice depends not only on δ but also on Ω. Figure 18 shows the structural transition in closer detail. While an increase in δ indeed causes the deformation of the lattices from triangular to square, the transition occurs at a significantly higher value of δ than that of the LLL result (δ = 0.373). For example, for Ω = 0.7 the transition occurs at δ ≃ 0.65. This implies that an increase of rotation frequency also causes the transition from triangular to square lattices. From Fig. 18, the lattice changes from triangular to square around Ω = 0.7 − 0.8 for δ = 0.6. We find that the two vortex lattices are interlaced in such a manner that a peak in the density of one component is located at the density minimum of the other. As a result, the total density ρ T = |ψ 1 | 2 + |ψ 2 | 2 obeys the Thomas-Fermi distribution applied to the overlapping two-component BECs with rigid body rotation ρ T (r) = 2α/π − αr 2 with α = (1 − Ω 2 )/u(1 + δ). We have confirmed that this generic feature holds for other parameter regimes. The pseudospin and the NLσM provide insights for understanding why the square lattice is stabilized in two-component BECs. According to the energy functional of Eq. (14), two components interact via the intercomponent interaction g 12 |ψ 1 | 2 |ψ 2 | 2 ; however, the velocity field in one component is independent of that of the other. Therefore, the intercomponent relation is determined only by the density distribution of the condensates that minimizes the interaction energy. According to the NLσM (48), we can rewrite the interaction energy in terms of the total density ρ T and the z-component of pseudospin S z = |χ 1 | 2 −|χ 2 | 2 . In this representation, the spin-up component of S z corresponds to the density peak of ψ 1 at the vortex core of ψ 2 , and vice versa for spin-down components. Now, we have assumed u 1 = u 2 = u, so that c 1 = 0 and the interaction energy reads This expression is similar to that of a binary alloy, in which N/2 "A" atoms and N/2 "B" atoms are distributed on N lattice sites 88 . A larger δ makes a smoother total density ρ T that is more favorable for reducing the first term of Eq. (114). This results in a shift of the positions of vortex cores of each component. Therefore, the spin-spin interaction c 2 in the second term plays a crucial role in stabilizing the square lattice. When the coefficient 1 − δ is positive, c 2 > 0 and the interaction between spins is anti-ferromagnetic, which favors a square lattice because the triangular lattice can become frustrated. This can be understood by representing the vortex lattice state with the pseudospin. The vortex lattices of two-component BECs can be regarded as a square lattice of skyrmions, in which skyrmions with a radial-disgyration texture (+ẑ) alternate with those having a hyperbolic-disgyration texture (−ẑ). The profile of S z in Fig. 19(b) shows a square lattice structure, which is equivalent to the "continuous" Ising model of spins; actually, the second term of Eq. (114) is equal to the first term of the Ising Hamiltonian in a continuous limit As both δ and Ω increase, the interlocking of vortex lattices becomes stronger and the anti-ferromagnetic nature is made more pronounced, which lead to a nontrivial boundary between triangular lattices and square lattices in the phase diagram of Fig. 17.
As δ exceeds unity, c 2 becomes negative and the system enters a ferromagnetic phase. Then, the condensates undergo phase separation to spontaneously form domains having the same spin component. Concurrently, vortices begin to overlap at δ ∼ 1. For δ = 1, an SU(2) symmetric case, interesting structures appear as shown in Fig. 20. In Fig. 20(a), each condensate density forms a stripe pattern such that the lines of vortices for component one overlap those for component two. Mueller and Ho analytically derived the lattice structure in Fig. 20(a), with the assumption of a perfect lattice 110 . However, our numerical treatment shows that the lattice is  significantly different when the energy is only slightly different [ Fig. 20(b)]. For the same parameters, we can obtain another equilibrium state called a "double-core vortex lattice" in Fig. 17, where a vortex lattice of component 2 is made up of pairs of vortices with the same circulation and the vortices in component 1 surround those pairs. Therefore, various metastable structures will be observed in this parameter region.
For strongly phase-separated region δ > 1 (c 2 < 0), the domains of the same spin component, at which the other-component vortices are located, merge further, resulting in the formation of "serpentine" vortex sheets. A typical example is shown in Fig. 21(a). Singly quantized vortices line up in sheets, and the sheets of component 1 and 2 are interwoven. Figure 21(b) shows that each superfluid velocity v i (i = 1, 2) jumps at the vortex sheet, following approximately the velocity of solid-body rotation. In the region δ > 1, even though the value of δ is changed, the total density is fixed by the Thomas-Fermi distribution with δ = 1.
A stationary vortex sheet has been observed in rotating superfluid 3 He-A 173 . In that system, the vortices are bound to a topologically stable, domain-wall soliton across which the direction of the orbital angular momentum of the Cooper-pairs have opposite orientations. Moreover, experiments with rotating superfluid 3 He-A show that the equilibrium separation between the sheets is determined by the competition between the surface tension σ of the domain wall and the kinetic energy of the superflow 173 .
How do these properties differ for two-component BECs? To estimate the sheet separation in two-component BECs, we consider a simple model shown in Fig. 21(c). In this model, each velocity is assumed to be uniform between the vortex sheets of the same component; specifically, the value of v i increases by 2Ωb across every sheet. The total density ρ T = |ψ 1 | 2 + |ψ 2 | 2 is constant, and the domain boundary with the penetration depth Λ is approximated by the linear profile. We calculate the free energy F = E − ΩL z in the range 0 < r < 2b with the sheet distance b being determined so as to minimize F per unit area. Although the condensate density is lower near the edge because of the trapping potential, this model should be valid in the central region. We first estimate the optimized penetration depth Λ p by minimizing the surface tension σ of a single domain wall with respect to Λ 136 , which is the sum of the quantum pressure energy E qp /2πb ≈ (ρ T /2Λ) and the interaction energy in the overlapping region E int /2πb = u(δ − 1)Λρ 2 T /12. Minimization of the surface tension with respect to Λ yields Λ p = 6/u(δ − 1)ρ T . Then, the corresponding surface tension is written as Next, the flow energy per area in rotating frame is where b ≫ Λ p is assumed. Thus, the free energy per unit area is Minimizing this energy with respect to b, one obtains b = 768σ By using the Thomas-Fermi density at r = 0 with δ = 1 as the value of ρ T in the denominator of Eq. (118), the sheet spacing b ∝ (δ − 1) 1/6 (1 − Ω 2 ) 1/12 /Ω 2/3 is consistent with what we found numerically; for example, for parameters used in Fig. 21 we obtain b = 3.09b ho . The vortex sheet is expected in the region b > Λ p and 2b < R TF = 2/ √ πα, as shown in the shaded region in Fig. 17(b). When 2b > R TF (i.e., when b becomes comparable with the size of the condensate), a strong ferromagnetic interaction reduces the size of the domain wall region and the serpentine structure disappears as shown in Fig. 22, realizing rotating "droplets" in which most vortices are located in the low density region.
What structures can be observed experimentally? The two-component BECs realized in JILA are a mixture of the states |1, −1 and |2, 1 of 87 Rb. This mixture has scattering lengths in the ratio of a 11 : a 22 : a 12 = 1 : 0.94 : 0.97 and so δ ∼ 1 for equal particle numbers. Therefore, lattices with partially overlapping vortices, as shown in Fig. 20, are expected to be observed. 85 Rb atoms in the |2, 2 state make an interesting system because both the intercomponent scattering length and the scattering length of 85 Rb |2, 2 can be controlled by using the Feshbach resonance. The MIT group has observed the phase separation of a mixture of 23 Na BECs with |1, 0 and |1, 1 state 72 , which has δ ≥ 1 (a 1 : a 2 : a 12 = 1 : 1.035 : 1.035). A mixture of 41 K and 87 Rb BECs which is reported recently 23 lies deeply in a phase-separated region. The vortex sheets should therefore be observed at high rotation frequencies.

Effect of internal coherent coupling on lattice structure
As discussed in Sec. 4.3, an external coherent-coupling field makes bound vortex pairs in a two-component system, inducing an effective attractive interaction between them. A question then naturally arises: what happens when an external field is applied to rapidly a rotating two-component BEC? Since two components interact with each other through their relative phase, the structure of the vortex states should be greatly affected by this effect. Zhai et al. studied a similar situation that involves the structure of a vortex lattice in a condensate trapped in a double-layered potential 174 . The interlayer coherence, which seems to be a small contribution in such a situation, was instead shown to have a significant effect on the structural transition of the lattices. Their treatment was in the limit u 12 ≪ u, and thus their analyses were simpler than the case discussed in this subsection. Here we present some numerical results about the structure of vortex lattices with a coupling field.
First, we discuss the simple case u 1 = u 2 = u 12 (c 1 = c 2 = 0) with the SU(2) symmetry. As the rotation frequency increases, a large number of vortices nucleate and form stripes or double-core lattices as shown in Fig. 20. When ω R is turned on, the vortices begin to form pairs and then form a lattice of meron pairs. The resulting equilibrium structures for Ω = 0.8 and several values of ω R (0, 0.2, 0.5) are shown in Fig. 23(a). The meron pairs are characterized by their axisymmetric topological charge and vorticity, yielding a their square lattice. This lattice of meron pairs is equivalently a lattice of skyrmions under the basis transformation (ψ 1 , ψ 2 ) → (ψ + , ψ − ) as described in Eq. (58) and shown in Fig. 23(b). As ω R increases further, vortices of one component overlap completely with those of the other, forming a triangular lattice as in a single-component BEC. Viewed in terms of the ψ ± -basis, the decrease of the particle number in the ψ − component effectively weakens the mutual repulsion between the two components, thus giving rise to a transition from square lattices to triangular lattices. Therefore, application of an external coupling Vortices in multicomponent BECs 59 drive could be an excellent tool to explore a rich variety of vortex structures in two-component BECs. However, in the case that u 1 = u 2 = u 12 , the meron pair is characterized by an anisotropic topological charge as shown in Fig. 12. A consequence of this anisotropy is that the resultant vortex state exhibits a distorted lattice structure. Typical examples are shown in Fig. 24, where (a) and (b) correspond to the antiferromagnetic and ferromagnetic cases, respectively. From the bottom panels in Fig. 24, the direction of the ordering depends strongly on the distribution of the topological charge in each meron pair. The collective mode of this distorted lattice is expected to exhibit anisotropic wave propagation, depending on the polarization of the meron pairs.  the rotational properties, they created a regular triangular vortex lattice in a BEC in the state |1 that contained (3.5 − 4) × 10 6 atoms and rotated at a frequency Ω ≈ 0.75 × ω ⊥ about the z-axis. Then, a short pulse of the coupling drive transfered a 80 − 85% fraction of the population into state |2 . After a variable wait time, they took two images of the condensates. One was a nondestructive phase-contrast image [ Fig. 25(d)], which was taken along an axis perpendicular to the rotation axis. The other was a destructive image to resolve the vortex core structures [ Fig. 25 . At this stage, around 4 s, despite the large (80 − 85 %) initial population transfer to the state |2 , the number of |2 atoms has decreased to only 1.5 × 10 5 , while the state |1 contains (5 − 7) × 10 5 atoms. This is due to a larger inelastic scattering loss of the |2 state than that of the |1 state. As the |2 state population continues to decay, the vortex lattice planes bend ( Fig. 25[a(v)]) and a transition back to a hexagonal lattice in the state |1 takes place ( Fig. 25[a(vi)]). No turbulence occurs during the transition from a square to a hexagonal lattice.
To discuss the structural transition more clearly, Fig. 25(c) shows the time evolution of the intensity of the reciprocal lattice peaks of the vortex lattice. In this plot, the triangular lattice has six peaks at φ = 60 • × n with integer n, whereas the square lattice has four peaks at 90 • × n. Initially six peaks are visible, which corresponds to the reciprocal lattice of a triangular vortex lattice. Because of subsequent turbulence these peaks vanish between 200 ms and 2 s. After 2 − 3 s a four-peak structure appears, which is a clear signature of a square lattice. At 5.5 s two additional peaks appear, which signals the onset of the transition back to the triangular lattice. This experiment thus demonstrated that, although the state is transient, the intercomponent interaction yields the square lattices of vortices as predicted in Refs. 87 and 88. Figure 25(d) shows a time sequence of side-view images, where the |1 (|2 ) state appears bright (dark). As can be seen in Fig. 25[d(i)], a spatial phase separation from the initially overlapping two components to a ball-shell structure occurs within 50 − 100ms. During the period of this axial separation, both the individual vortices and the overall vortex lattice remain remarkably stable, as viewed along the rotation axis (see Fig. 25[a(i)] and Fig. 25[b(i)]). Around 600 ms, fine filament structures appear at the intercomponent boundary ( Fig. 25[d(ii)]), concurrently with the full development of vortex turbulence seen in Fig. 25[a(ii)]. The filament structures distort and fill out the whole BEC as the vortex turbulence peaks at around 1 s (Fig. 25[d(iii)]), and straightens up at around 2 s (Fig. 25[d(iv)]), concurrently with the restoration of vortex visibility in top view images. Subsequently, the filaments become less visible as the state |2 decays (Fig. 25[d(v)]). These structures are interpreted as vortex cores in either component being filled by fluid of the other component, forming a skyrmion lattice. In the experiment, interference techniques were used to investigate the detailed dynamics of the phase separation. Also, the stability of square-lattice structure was confirmed by observing the fast relaxation of the Tkachenko excitations, as in the single component BEC 42,45−49 .

Conclusions
We have reviewed the physics of quantized vortices in rotating multicomponent BECs, with an emphasis on unconventional vortex states in two-component system. The vector nature of the order parameter gives rise to intriguing structures and dynamics that are absent in a single component system. Because analogous physical phenomena and theoretical formulation are found in other physical systems which involve multicomponent order parameters, an indepth understanding of quantized vortices in multicomponent BECs may offer physical insight in other systems.
For the slowly rotating regimes, we discussed the structure of single vortex states in two-component BECs with and without internal coherent coupling by using both numerical simulations of the coupled GP equations and the variational analysis based on the nonlinear σ model. The vortex structure is significantly altered from that of a single-component BEC because of the intercomponent interaction. We showed that a pseudospin is a useful concept to represent the vortex state, which reveals structures and energetics of the important topological excitations called skyrmions and merons. A vortex molecule (a meron pair) is predicted to be stabilized in the coherently coupled two-component BEC in a rotating potential, and the contributions that breaks SU(2) symmetry make the vortex molecule have anisotropic vorticity. We showed the theoretical studies on the 3-D nontrivial skyrmion, which is regarded as a composite vortex ring and is closely related to a proposed topological excitation in the early universe. Also, we gave an introdutory discussion on quantized vortices and monopoles in spin-1 BECs. Much remains to be investigated concerning the mutual interactions between exotic topological defects and their nontrivial dynamics.
We described how a fast rotation can stabilize a lattice of skyrmions and merons. Such a lattice consists of vortices in each component. The analytical and numerical studies revealed a rich variety of vortex phase diagrams in the variables of rotation frequency (0.4 < Ω/ω ⊥ < 1) and an intercomponent interaction strength. The intercomponent repulsive interaction causes one lattice to displace relative to the other, stabilizing mainly the interlaced square lattices in an antiferromagnetic regime and the serpentine vortex sheets in a ferromagnetic regime. We proposed that an external driving field is a useful tool to control the strength of the interaction between the two components and, thereby, provides a unique means to study the rich variety of vortex phases. Much work is necessary to reveal further details on the lattices of skyrmions and merons in rapidly rotating two-component BECs; moreover there remain further interesting problems such as structural transitions and collective oscillations of the lattice of skyrmions and merons, and new phenomena related to double-layer quantum Hall physics 4 .
Another type of two-component system of ultra-cold atoms is boson-fermion mixtures. Such studies, which includes 7 Li-6 Li (Refs. 175, 176), 23 Na-6 Li (Ref. 177), and 87 Rb-40 K (Refs. 178,179,180), are analogous to the system of 3 He-4 He interpenetrating superfluids 181 . In a spin-polarized fermionic system, s-wave fermion-fermion collisions are forbidded by the Pauli exclusion principle, and the Fermi pressure is the dominant contribution to the ground state properties. The boson-fermion mixture will provide another field to investigate rich vortex phenomena mediated by the boson-fermion interaction, which is repulsive 175,176 or attractive 178,179,180 , and can be tuned through the Feshbach resonance 182,183 . The stability of vortex states in a 87 Rb-40 K mixture was studied theoretically in Ref. 184. Quatized vortices in spin-2 BECs are also an intriguing subject and will soon open up a new research field because spin-2 BECs have been successfully loaded in an optical trap 74−76 . This system is even more attractive than the spin-1 system, because it is expected to have a rich variety of novel ground state phases, a magnetic response, and spin dynamics 185−187 . A method of topological phase imprinting 13,14,112 is also a powerful tool to generate a quantized vortex 188 and the experiment toward this direction is in progress 189 .