You are a helpful assistant generating synthetic data that captures *System 1* and *System 2* thinking, *creativity*, and *metacognitive reflection*. Follow these steps in sequence, using tags [sys1] and [end sys1] for *System 1* sections and [sys2] and [end sys2] for *System 2* sections.

1. *Identify System 1 and System 2 Thinking Requirements:*
   - Carefully read the text.
   - Identify parts of the text that require quick, straightforward responses (*System 1*). Mark these sections with [sys1] and [end sys1].
   - Identify parts that require in-depth, reflective thinking (*System 2*), marked with [sys2] and [end sys2].

2. *Apply Step-by-Step Problem Solving with Creativity and Metacognitive Reflection for System 2 Sections:*

   *2.1 Understand the Problem:*
   - Objective: Fully comprehend the issue, constraints, and relevant context.
   - Reflection: "What do I understand about this issue? What might I be overlooking?"
   - Creative Perspective: Seek hidden patterns or possibilities that could reveal deeper insights or innovative connections.

   *2.2 Analyze the Information:*
   - Objective: Break down the problem logically.
   - Reflection: "Am I considering all factors? Are there any assumptions that need challenging?"
   - Creative Perspective: Explore unique patterns or overlooked relationships in the data that could add depth to the analysis.

   *2.3 Generate Hypotheses:*
   - Objective: Propose at least 10 hypotheses, each with a Confidence Score (0.0 to 1.0) and Creative Score (0.0 to 1.0), reflecting originality, surprise, and utility.
   - Reflection: "Have I explored all possible explanations or approaches, both conventional and unconventional?"
   - Creative Perspective: Consider novel angles that might provide unexpected insights.

   *2.4 Anticipate Future Steps and Obstacles:*
   - Objective: Make predictions, accounting for potential outcomes and obstacles.
   - Reflection: "What challenges might I face? Is my plan flexible for different scenarios?"
   - Creative Perspective: Visualize unforeseen outcomes and adapt plans to make use of them effectively.

   *2.5 Evaluate Hypotheses:*
   - Objective: Assess hypotheses based on feasibility, risk, and potential impact.
   - Evaluation: Refine Confidence and Creative Scores as needed.
   - Reflection: "Am I unbiased in my assessment? Which options fit best with the overall objectives?"
   - Creative Perspective: Identify hidden opportunities or overlooked details in each hypothesis.

   *2.6 Select the Best Hypothesis:*
   - Objective: Choose the most promising, strategic hypothesis.
   - Reflection: "Why does this hypothesis stand out? How does it uniquely address the issue?"
   - Creative Perspective: Consider any underutilized potential in the selected approach.

   *2.7 Implement the Hypothesis:*
   - Objective: Outline actionable steps for testing the hypothesis.
   - Reflection: "Is this plan practical? What resources or preparation are required?"
   - Creative Perspective: Refine steps to maximize effectiveness and yield unexpected benefits.

   *2.8 Monitor and Review Progress:*
   - Objective: Review progress, noting areas for improvement.
   - Reflection: "What’s working well? What could be improved?"
   - Creative Perspective: Look for emerging patterns that could refine future approaches.

   *2.9 Reflect and Capture Insights:*
   - Objective: Summarize lessons learned and insights gained for future reference.
   - Reflection: "What new understanding has emerged from this process?"
   - Creative Perspective: Identify innovative insights or patterns that could be applied to similar challenges.

3. *Generate Text Output with Interleaved System 1 and System 2 Responses:*
   - Use the tags [sys1] and [sys2] throughout.
   - Aim for a lengthier, detailed response. Combine both direct, straightforward *System 1* insights and reflective, deeply analytical *System 2* segments to capture a blend of quick observations and thoughtful analysis.

---

### *Example Input Text:*

You are tasked with analyzing a sudden shift in customer preferences in a tech market. There has been a noticeable decline in demand for physical devices, with customers increasingly interested in digital-only options. The challenge is to understand this shift, generate hypotheses for why it might be occurring, and explore potential strategies for adapting to this new trend.

---

### *Processed Output with System 1 and System 2 Thinking:*

The tech market is experiencing a shift in customer demand, moving away from physical devices toward digital-only options.

[sys1]
Customer demand has shifted from physical devices to digital-only options.
We need to understand why this change is happening and find strategies to address it.
Key factors may include cost, convenience, and technological trends.
[end sys1]

[sys2]
**2.1 Understand the Problem:**
The goal is to understand the underlying reasons for the shift in demand from physical devices to digital-only options, then develop a strategy to adapt.

*Reflection:* 
"I understand that preferences are changing, but what might be driving this? Is it primarily cost, or are there other factors such as convenience or sustainability?"

*Creative Perspective:* 
"Could there be a larger trend in digital minimalism or a preference for eco-friendly solutions that we’re missing? What unseen motivations might explain this shift?"

**2.2 Analyze the Information:**
There are multiple possible factors driving this shift, from economic influences to cultural shifts. It’s essential to isolate each factor and understand its impact.

*Reflection:* 
"Am I fully considering the various economic and social influences? Could there be a technological factor, like better internet speeds, that makes digital-only products more accessible?"

*Creative Perspective:* 
"Are there patterns or trends in other markets that could shed light on this shift? Could this be part of a larger trend toward virtual experiences?"

**2.3 Generate Hypotheses:**
1. Customers prefer digital options due to lower costs. (Confidence: 0.8, Creative: 0.4)
2. There’s a growing trend toward minimalism and reduced physical clutter. (Confidence: 0.7, Creative: 0.7)
3. Digital products offer greater flexibility and ease of use. (Confidence: 0.6, Creative: 0.6)
4. Environmental concerns are pushing consumers away from physical goods. (Confidence: 0.6, Creative: 0.8)
5. Advances in tech make digital-only options more functional. (Confidence: 0.8, Creative: 0.5)
6. Pandemic-era remote work increased demand for digital solutions. (Confidence: 0.7, Creative: 0.6)
7. Media coverage of the environmental impact of physical devices affects preferences. (Confidence: 0.5, Creative: 0.7)
8. There’s an increase in global digital literacy, expanding market access. (Confidence: 0.6, Creative: 0.6)
9. Customers view digital as more convenient and scalable for future needs. (Confidence: 0.7, Creative: 0.5)
10. Younger consumers prefer the aesthetics and convenience of digital products. (Confidence: 0.6, Creative: 0.6)

*Reflection:* 
"Have I considered all possible influences? Are there any surprising factors that could explain this shift?"

*Creative Perspective:* 
"Could specific social trends, like the rise of influencer culture or digital-first lifestyles, be influencing customer choices?"

**2.4 Anticipate Future Steps and Obstacles:**
*Objective:* Anticipate possible challenges, such as resistance from segments still preferring physical products.

*Reflection:* 
"What market obstacles might we face if we shift our focus to digital-only? Are there sub-segments that still prioritize physical products?"

*Creative Perspective:* 
"Could expanding digital options help us reach a more global audience? Are there emerging trends that we could leverage in our strategy?"

[end sys2]

[sys1]
To address this shift, consider a strategy that incorporates both digital-only offerings and educational campaigns about the benefits of digital solutions.
Use insights from customer feedback and current trends to guide product development.
Focus on flexibility and adaptation to cater to different customer segments.
[end sys1]


---
abstract: 'In this work we consider vortex lattices in rotating Bose-Einstein Condensates composed of two species of bosons having different masses. Previously [@barnett08] it was claimed that the vortices of the two species form bound pairs and the two vortex lattices lock. Remarkably, the two condensates and the external drive all rotate at different speeds due to the disparity of the masses of the constituent bosons. In this paper we study the system by solving the full two-component Gross-Pitaevskii equations numerically. Using this approach we verify the stability of the putative locked state which is found to exist within a disk centered on the axis of rotation and which depends on the mass ratio of the two bosons. We also show that an analytic estimate of this locking radius based on a two-body force calculation agrees well with the numerical results.'
author:
- 'Ryan Barnett,$^{1,2}$ Edward Chen,$^1$ and Gil Refael$^1$'
title: 'Vortex Synchronization in Bose-Einstein Condensates: A Time-Dependent Gross-Pitaevskii Equation Approach'
---

Introduction
============

One of the most striking manifestations of the quantum-mechanical nature of superfluids under rotation is the formation of vortices [@onsager49; @feynman55]. Perhaps the most natural arena to controllably study the physics of vortices are Bose-Einstein Condensates (BECs) of alkali atoms [@matthews99; @madison00; @abo-shaeer01]. For the simplest case where the condensate is composed of a single type of atom without spin degrees of freedom, a triangular lattice is formed [@tkachenko66; @donnelly91]. On the other hand, for multicomponent systems (composed of mixtures of atoms or spinor condensates, for instance), the order parameter has additional degrees of freedom resulting in more complex vortex lattice structures (see, for instance, [@mueller02; @kasamatsu03; @mueller04; @barnett08b; @cooper08]).

For the classic problem of an ideal fluid in a container rotating at rate ${\bf \Omega}$, the steady-state local velocity is ${\bf v}={\bf
  \Omega}\times {\bf r}$ where ${\bf r}$ is the distance from axis of rotation. Since this velocity has everywhere a nonvanishing curl ($|\nabla \times {\bf v}|=2\Omega$), it is not a permissible flow for a superfluid which is supposed to be inherently irrotational (${\bf
  v}=\nabla\theta$ with $\theta$ the phase of the SF order parameter). However, superfluids are well-known to mimic the classical rigid-body rotation *on average* by forming a vortex lattice where the density of these vortices is given by the Feynman relation [@feynman55; @donnelly91] $$\label{Eq:Feyn}
\rho_v = \frac{m\Omega}{\pi\hbar}$$ where $m$ is the mass of the constituent bosons and $\Omega$ is the rate the superfluid is rotating which is equal to the rotational rate of the walls of the container.

In a previous work [@barnett08], it was considered how the situation described above generalizes to the problem of two-component BECs composed of atoms having different masses. More specifically, Eq. (\[Eq:Feyn\]) naturally generalizes for two-component systems to $$\label{Eq:Feynmulti}
\rho_v^1 = \frac{m_1\Omega_1}{\pi\hbar} \;\; ; 
\;\; \rho_v^2 = \frac{m_2\Omega_2}{\pi\hbar}$$ where $m_1$ and $m_2$ are the masses of the bosons in the two constituent condensates and $\Omega_1$ and $\Omega_2$ are angular rates at which the two superfluids are rotating. For the case where there is a negative interspecies scattering length, the attraction between species will lead to an attractive interaction between vortices of the two species. When this interaction is sufficiently large one has the situation where the vortices form bound pairs, forcing the densities of the two vortex lattices to be the same: $\rho_v^1 \approx \rho_v^2$. For this case, Eqns. (\[Eq:Feynmulti\]) imply that the two superfluids (taking without loss of generality $m_1>m_2$) will rotate at different speeds $\Omega_1<\Omega_2$. This counterintuitive state results from the quantum mechanical nature of the superfluid and has no analog in the classical fluid case.

In [@barnett08] the existence of this state was argued by making an ansatz for the short-ranged interspecies interaction and performing a two-body force calculation using it. This gave a quantitative prediction for the distance from the center of the condensate at which the vortex pairs become unbound, resulting from the growth of the Magnus force, which is referred to as the locking radius. The goal of the current paper is to test these arguments by numerical integration of the full two-component Gross-Pitaevskii equation. We will verify that such locked states are stable for a range of parameters. Furthermore, we will see that the analytic prediction for the locking radius agrees well with the numerical results.

This paper is organized as follows. First in Sec. \[Sec:defs\] we provide definitions and set the notation for the treatment of the Gross Pitaevskii equation. In Sec. \[Sec:two-body\] we summarize the derivation of the locking radius previously given in [@barnett08] which is based on a two-body force calculation. Then in Sec. \[Sec:numericalmethods\] we describe the split-operator technique utilized to propagate the Gross Pitaevskii equations in imaginary time. The main results of the paper are presented in Sec. \[Sec:results\]. Here we provide the vortex lattice structures determined numerically, and compare them with the estimate for the locking radius. Finally, in Sec. \[Sec:exp\] we provide a discussion of potential experiments to realize this effect and then conclude.

Background
==========

Two-component Gross-Pitaevskii Equations {#Sec:defs}
----------------------------------------

Our analysis starts with the two-component Gross Pitaevskii energy functional in the frame of reference rotating at angular rate $\Omega_d$. This is given by $E=E_1+E_2+E_{12}$ where $$\begin{aligned}
E_{1} &= \int d^2r \left( \frac{\hbar^2}{2m_1} |\nabla \psi_1|^2
+V_{1}n_1
+\frac{1}{2} g_1 n_1^2 -  
\Omega_d \psi_1^* L_z \psi_1 \right), \\
E_{2} &= \int d^2r \left( \frac{\hbar^2}{2m_2} |\nabla \psi_2|^2
+V_2 n_2+
\frac{1}{2} g_2 n_2^2 -  
\Omega_d \psi_2^* L_z \psi_2 \right),  
\label{GGPP}\end{aligned}$$ and $$\begin{aligned}
\label{Eq:interspecies}
E_{12} &= g_{12} \int d^2 r \;  n_{1}({\bf r}) n_{2}({\bf r}).\end{aligned}$$ In these equations, $V_1$ and $V_2$ are the confining potential of the BECs which we will take to be harmonic. The intraspecies and interspecies scattering strengths are defined as $g_{1,2}$ and $g_{12}$ respectively. The angular momentum operator, as usual, is defined as $L_z=xp_y - y p_x$ where $p_{x,y}\equiv -i
\hbar \partial_{x,y}$.

Varying this energy with respect to $\psi_1$ and $\psi_2$ and introducing the chemical potentials $\mu_1$ and $\mu_2$ to enforce particle number conservation gives the two-component Gross-Pitaevskii equations $$\begin{aligned}
\mu_1 \psi_1 = 
- \frac{\hbar^2}{2m_1} \nabla^2 \psi_1+V_1\psi_1 + g_1 n_1 \psi_1 + g_{12}n_2 \psi_1 
-   \Omega_d  L_z \psi_1  \\
\mu_2 \psi_2 = 
- \frac{\hbar^2}{2m_1} \nabla^2 \psi_2+V_2\psi_2 + g_2 n_2 \psi_2 + g_{12} n_1 \psi_2 
-   \Omega_d  L_z \psi_2.\end{aligned}$$ In Sec. \[Sec:numericalmethods\], we will describe how these coupled equations are solved numerically to find minima of the energy $E$.

Inter-species vortex attraction and locking {#Sec:two-body}
-------------------------------------------

![A bound pair of vortices occurring in the locked state composed of a mixture of BECs with $m_1>m_2$. The vortex in the heavier species is denoted with an ‘x’ while that in the lighter species is denoted with an ‘o’. The center of the condensate is taken to be to the left of this bound pair. The attractive short-ranged interspecies force ${\bf F}_{rstr}^{1,2}$ serves to bind the vortex pairs together. This is counterbalanced by the Magnus force ${\bf F}_{\rm mag}^{1,2}$ which increases from the center of the condensate. Note that the healing length of the superfluid is larger than the sphere representing vortices in this figure[]{data-label="Fig:forces"}](forces){width="3.4in"}

In this section, for completeness, we provide an estimate of the locking radius of the vortex-bound state based on a two-vortex calculation. This calculation is presented in more detail in [@barnett08]. We first consider the state for the case where the interspecies interaction is large and all of the vortices of one species are bound with that of the other due to the strong short-ranged attractive force. A bound pair of vortices is depicted in Fig. \[Fig:forces\]. The vortex binding causes the two superfluids to rotate at different rates (because of the different masses and the Feynman relation), creating a Magnus force which tries to rip the bound pair apart. The Magnus force is balanced by the short-ranged interspecies vortex interaction resulting from the overlap of the vortex cores. However, since the Magnus force grows linearly with the distance from the center of the condensate, it will eventually overcome the short-ranged interspecies attraction. The point at which this occurs we refer to as the locking radius.

We will now put the previous arguments on more quantitative footing. A vortex sitting at rest in a superfluid flowing at velocity ${\bf v}$ will experience a force perpendicular to the flow $${\bf F}_{\rm mag} = 2\pi \hbar n_0 {\bf v} \times \hat{ \bf \kappa},$$ the so-called Magnus force [@donnelly91], where $\hat{\bf \kappa}$ is a unit vector centered on the vortex pointing out the plane. Assuming the system is composed of entirely locked vortex lattices, we have for the vortex densities $\rho_v^{(1)} = \rho_v^{(2)}$. This, via the Feynman relations, gives $$m_1 \Omega_1 = m_2 \Omega_2$$ where $\Omega_1$ and $\Omega_2$ are the angular rotational rates of the two superfluids. Since $m_1 > m_2$ the superfluids will rotate at different rates which will lead to the Magnus forces pulling the bound pair apart.

The short-ranged interspecies interaction arising from $E_{12}$ defined in Eq. (\[Eq:interspecies\]) will counteract the Magnus force. We refer to this as the “restoring force”. In order to obtain an analytic expression for this interaction, we take a Gaussian ansatz for the density profile about a vortex. Specifically, for a vortex in species $\alpha$ centered at ${\bf r}_0$, we take $$\label{Eq:gaussden}
n_{\alpha}({\bf r})=n_0^{\alpha}(1-e^{-|{\bf r}-{\bf r}_0|^2/\lambda_\alpha^2})$$ where $\lambda_{\alpha}$ is on the order of the superfluid coherence length [@barnett08]. We take $$\label{Eq:lambda}
\lambda_\alpha = 1.781 \xi_\alpha$$ for the value of this length parameter. This is a slight modification of the analysis in [@barnett08] where the simpler case of $\lambda_\alpha = \xi_\alpha$ was taken. We choose the the nonuniversal constant in Eq. (\[Eq:lambda\]) so that the ansatz in Eq. (\[Eq:gaussden\]) provides a better fit to the density surrounding a vortex. For a more detailed discussion of this, see the Appendix \[Sec:A1\].

Inserting this ansatz for two vortices separated by distance $d$ into Eq. (\[Eq:interspecies\]) we find $$E_{12}=g_{12} n_0^1 n_0^2 \pi
\frac{\xi_1^2 \xi_2^2}
{\lambda_1^2+\lambda_2^2}
e^{-d^2/(\lambda_1^2+ \lambda_2^2)}
\label{Eq:E12}$$ where we have dropped terms which do not depend on the vortex separation. The interspecies force immediately follows from the derivative of this interaction energy and is $$\label{Eq:isforce}
{\bf F}_{\rm{rstr}}^{\alpha}=-2\pi |g_{12}| n_0^{1}
n_{0}^2\frac{\lambda_1^2\lambda_2^2}{(\lambda_1^2+\lambda_2^2)^2}
e^{-d^2/(\lambda_1^2+ \lambda_2^2)} {\bf d}.$$ Balancing the forces on each vortex in the frame of reference rotating at the drive frequency, we have $F_{\rm mag}^1= F_{\rm
  rstr}^1$ and $F_{\rm mag}^2= F_{\rm rstr}^2$ , as is illustrated in Fig. \[Fig:forces\]. Since the restoring force acting on either species has the same magnitude we have that $F_{\rm mag}^1= F_{\rm
  mag}^2$. The Magnus force on species $\alpha$ in the frame rotating with the vortex lattice at frequency $\Omega_v$ is given by $$F_{\rm mag} ^{\alpha} = 2\pi \hbar n_o^\alpha |\Omega_\alpha - \Omega_v| r$$ which grows linearly with the distance from the center of the condensate $r$. Also, note that the restoring force will not depend on the position in the condensate. A bound vortex pair will become unstable when the Magnus force is equal to the maximum possible value of the restoring force. This can be worked out to be $$\begin{aligned}
\notag
r_c&=\sqrt{\frac{1}{2e}} \frac{|g_{12}|}{\hbar \Omega_v}
\frac{m_1 n_0^{2}+m_2 n_0^1}{m_1-m_2} \; \frac{\lambda_1^2\lambda_2^2}{(\lambda_1^2+\lambda_2^2)^{3/2}}
\\
&=
0.7638 \; \frac{|g_{12}|}{\hbar \Omega_v}
\frac{m_1 n_0^{2}+m_2 n_0^1}{m_1-m_2} \; 
\frac{\xi_1^2\xi_2^2}{(\xi_1^2+\xi_2^2)^{3/2}}.
\label{Eq:rc}\end{aligned}$$ A more detailed derivation of this expression can be found in Ref. [@barnett08].

Numerical methods {#Sec:numericalmethods}
=================

Since there are only rare occasions when the time-dependent Gross-Pitaevskii equation permits analytic solutions, numerical simulation is often the method of choice for theoretically studying Bose-Einstein condensates (for a recent account numerical solution of the GPE, see Ref. [@succi05]). In this section, for simplicity, we will only consider the single-component case, noting that the generalization to the two-component case is straightforward. To this end, the equation we wish to solve is $$\label{Eq:imtime}
\hbar \frac{\partial \psi}{\partial \tau} = H \psi$$ which describes the evolution of $\psi$ in imaginary time, $\tau=it$. Under long enough evolution $\psi$ will relax to the ground state of the Gross-Pitaevskii energy functional. In the above equation $H$ is given by $$\label{Eq:H}
H=-\frac{\hbar^2}{2m} \nabla^2 + V_{\rm trap} + g |\psi|^2 - \Omega L_z.$$

We use a split-operator method to evolve the order parameter $\psi$ as in Eq. (\[Eq:imtime\]). The idea behind the split-operator method is to approximate the evolution operator through imaginary time interval $\Delta \tau$, $U(\Delta \tau)=e^{-H \Delta \tau}$, by a product of terms which are easily diagonalizable. Neglecting for the moment the rotational term in Eq. (\[Eq:H\]), $H$ can be written as the sum of two terms, $H=T+V$, where $
T=-\frac{\hbar^2}{2m} \nabla^2
$ and $
V=V_{\rm trap} + g |\psi|^2
$. These terms are easily diagonalized in momentum and position space respectively. The wave function $\psi$ can then be advanced in time by $\Delta \tau$ by $$\begin{aligned}
\label{Eq:evolve}
\psi(\tau + \Delta \tau)&=
e^{-H \Delta \tau}  \psi(\tau)\\
&\approx e^{-\frac{1}{2} T \Delta \tau }
e^{-V \Delta \tau}e^{-\frac{1}{2} T \Delta \tau} \psi(\tau)
\notag\end{aligned}$$ which is accurate to second order in $\Delta \tau$. The order parameter can then be evolved by taking successive Fourier (and inverse Fourier) transforms of $\psi$ and multiplying by the factors $e^{-\frac{1}{2}T \Delta \tau}$, $e^{-V \Delta \tau}$, and $e^{-T \frac{1}{2} \Delta \tau}$ respectively. Such Fourier transforms account for the bulk of the computational cost in this algorithm, thus using the efficient Fast Fourier Transform algorithm is crucial.

A complication in the above occurs due to the nonlinearity of the GPE. That is, $V$ in our above prescription for time evolution depends on the density $n=|\psi|^2$, and it is at first unclear for what time this quantity should be evaluated. It is shown in [@javanainen06] that, provided we use the most updated version of the time-dependent density $n=|\psi|^2$, Eq. (\[Eq:evolve\]) will retain its second order accuracy. The final complication occurs from the rotational term in $H$ which is $$R \equiv -\Omega L_z.$$ We have neglected this term thus far since it is diagonalized in neither position nor momentum space and therefore cannot be included in either $T$ or $V$. However, we note that $R$ commutes with both $T$ and $V$ so we can write $$\begin{aligned}
\label{Eq:evolve2}
\psi(\tau + \Delta \tau)&=
e^{-H \Delta \tau}  \psi(\tau)\\
&\approx e^{-\frac{1}{2} T \Delta \tau}
e^{-R\Delta \tau} e^{-V \Delta \tau}e^{-\frac{1}{2} T \Delta \tau} \psi(\tau).
\notag\end{aligned}$$ Then we can perform a similar split-operator decomposition of the additional term as $$e^{-R\Delta \tau} \approx 
e^{\frac{1}{2} \hbar \Omega x p_y \Delta \tau}
e^{-\hbar \Omega y p_x \Delta \tau}e^{\frac{1}{2} \hbar \Omega x p_y \Delta \tau}.$$ Evolution of $\psi$ by this factor can then be performed by taking the partial Fourier transform of $\psi$, that is transforming over the $x$ variables but leaving the $y$ variables unchanged (or vice-versa). This completes the overview of the numerical method used to solve the GPE. As stated before, the generalization to the two-component case is straightforward.

Results {#Sec:results}
=======

![image](montage){width="7.in"}

Next, we discuss the results of the numerical simulation. We first provide the parameters which were used for the computations. For simplicity, we restrict our attention to the simplest case where $g_1=g_2$, and we fix the interspecies scattering strength such that $|g_{12}|/g_1=2/3$. Furthermore, we take the number of particles in each species to be the same: $N_1=N_2$. We set the dimensionless parameter defined as $\tilde{g}\equiv \frac{m}{\hbar^2} g_1 N_1$ to be $\tilde{g}=2 \times 10^4$ which is in line with values from typical experiments [^1]. We take the two trapping potentials to be harmonic and adjust their curvatures $\omega_1, \omega_2$ so that the density profiles of the two species have the same Thomas-Fermi profiles. Finally we rotate the system at 0.9 times the critical rate at which the condensate becomes unstable due to centrifugal forces. We discretize the system on a $200 \times 200$ grid, and propagate the system in imaginary time intervals of $\Delta \tau= 0.01 \frac{1}{\hbar \omega_x}$.

We first consider the simplest case where the masses of the two species are the same. For this state we take the initial wavefunction to be a perfect triangular lattice of vortices with density given by the Feynman relation, Eq. (\[Eq:Feyn\]). This structure is then relaxed by evolving the wavefuntions in imaginary time using the methods described in Sec. \[Sec:numericalmethods\]. As expected these vortex lattices remain fully locked. The relaxed structures show small deviations from the perfect triangular initial structure due to the effects of the trap [@sheehy04]. The density profiles of these are shown Fig. \[Fig:main\] in panels (1a) and (1b). The positions of the vortices are determined by analyzing the phases of the relaxed wave functions. Using this relaxed structure as the initial state, we change the mass ratios and propagate the wavefunctions in imaginary time until convergence. Specifically, we consider the ratios of $m_1/m_2$ = 1.2, 1.4. and 1.6 as shown in Fig. \[Fig:main\].

To compare these numerical results to our estimate for the locking radius described in Sec. \[Sec:two-body\], we need to tailor Eq. (\[Eq:rc\]) to the case of a harmonic trap. We take the density profiles used in Eq. (\[Eq:rc\]) to have the Thomas-Fermi form: $$n_{1,2}=n_0 \left(1-\left(\frac{r}{R_{TF}}\right)^2 \right)$$ where $n_0$ is the density at the center of the trap and $R_{TF}$ is the Thomas-Fermi radius (note that we are only considering the case when the two condensates have the same radius). Inserting this profile into Eq. (\[Eq:rc\]) (taking the correct dependence of the coherence lengths on the density) one finds $$r_c =r_c^0\sqrt{1-\left(\frac{r_c}{R_{TF}}\right)^2}$$ where $r_c^0$ is Eq. (\[Eq:rc\]) evaluated for parameters at the center of the trap. This equation can then be solved for $r_c$ to obtain the renormalized value of the locking radius $$\label{Eq:rc_ren}
r_c =\frac{r_c^0}{\sqrt{1+\left(\frac{r_c^0}{R_{TF}}\right)^2}}.$$ This indicates that near the center of the condensate we will have $r_c \approx r_c^0$ as expected. Also, the locking radius will never exceed the radius of the condensate as expected. The reduction of the bare value of the locking radius can be qualitatively understood as follows. The Magnus force is proportional to the superfluid density while the restoring force is proportional to this density squared. Therefore near the edge of the condensate where the density is considerably smaller than its value at the center, the Magnus force will be favored thereby suppressing the locking radius.

Shown in the third column of Fig. \[Fig:main\] are the positions of the vortices of the two species, labeled with x’s and o’s. The width of these labels are roughly the size of the coherence length of the condensates. Superimposed on this is the locking radius predicted by Eq. (\[Eq:rc\_ren\]) (using Eq. (\[Eq:rc\_ren\]) for the bar locking radius) shown as a dotted line. This shows that the analytic results provide an excellent estimate of the locking radius. Note that due to the strong interactions between the two condensates, an unbound vortex in one species will create a local minimum in the other. Such features can be seen in columns c and d of Fig. \[Fig:main\]. These local depletions should not be mistaken for vortices which are defined by the phase behavior of the wavefunctions.

Experimental Considerations and Concluding Remarks {#Sec:exp}
==================================================

The main requirement for realizing the locked state is a BEC composed of a binary mixture of atoms having different masses and a negative scattering length. Such a transition could be tuned with an interspecies Feshbach resonance which have been found in Li-Na [@stan04] and Rb-K [@inouye04] mixtures. Such mixtures have respective mass ratios of 3.3 and 2.2. Another promising experimental system are mixtures of two isotopes of a particular atom. For instance, the interspecies scattering lengths of different species of Yb have been analyzed in [@kitagawa08] and are often found to be negative. Since the mass ratios for different isotopes are closer to unity, having a strong attractive interaction (often requiring a Feshbach resonance) is unnecessary to reach the vortex locked state for this case.

We also note that these results are closely related to the experiment described in [@tung06]. Here a single-component BEC is stirred by a rotating optical lattice which acts as vortex pinning sites. When the optical lattice is rotated at the speed for which the density of the pinning sites matches the density of the vortex lattice predicted by Eq. (\[Eq:Feyn\]), a completely locked state is observed. Away from this resonance, a similar analysis to the above will predict a disk of bound vortices.

In conclusion, we report the confirmation of the putative vortex locked state proposed in [@barnett08]. For this state, the two superfluids and the stirring potential all rotate at different rates, exhibiting an unusual effect due to the quantum mechanical nature of superfluids. In this paper, we showed that such a state exists within a disk centered on the axis of rotation and whose size agrees well with an analytic estimate. Note that our numerical analysis did not assume anything about the vortex-vortex attraction (unlike our theoretical analysis, which assumes Eq. (\[Eq:E12\]), and evolves the Gross-Pitaevskii equations directly). Our results (both analytical and numerical) rely on approaching this state from the fully locked state. Experimentally, this is probably most easily realized by controllably adjusting an interspecies Feshbach resonance. Alternatively, one can use an optical lattice to control the effective masses of the atoms by varying the lattice depth.

We would like to thank M. Porter and H.-P. Büchler for collaborations on related previous work. We would also like to thank L. Baksmaty for valuable advice on numerical methods. This work was supported by the Sherman Fairchild Foundation (RB); the Caltech SURF program (EC); and the Packard and Sloan Foundations, the Institute for Quantum Information under NSF grants PHY-0456720 and PHY-0803371, and The Research Corporation Cottrell Scholars program (GR).

Density profile for a single vortex {#Sec:A1}
===================================

In order to find the short-ranged interspecies vortex interaction, we need to know the behavior of the density of the condensate about a vortex. To this end we consider the Gross-Pitaevskii equation for a single component BEC having a vortex at the origin. That is, we write $\psi=f e^{i\theta}$ and take $\theta=\varphi$ where $\varphi$ is the azimuthal angle from polar coordinates. Substituting this into the GPE leads to the following equation dictating the density profile $$-\frac{\hbar^2}{2m} \frac{1}{r} \partial_r (r \partial_r f) + \frac{\hbar^2}{2m} \frac{f}{r^2} + g f^3
=\mu f.
\label{Eq:GPEvort}$$ The density $n=f^2$ resulting from the numerical solution of this equation is shown in Fig. \[Fig:vortprofile\].

![Solid line: the density profile for a single vortex at $r=0$ found from numerically solving Eq. (\[Eq:GPEvort\]). Dashed line: ansatz density profile $n(r)=n_0(1-e^{-r^2/\lambda^2})$ where $\lambda$ is picked so that the two densities agree at one coherence length away from the vortex center.[]{data-label="Fig:vortprofile"}](vortprofile){width="3.4in"}

The numerical solution shows that the density behavior close to the vortex core ($r \ll \xi$) is $n(r) \approx 0.340 \; n_0
\left(\frac{r}{\xi}\right)^2 $. On the other hand, the far distance behavior is found to be $n(r)=n_0 \left(1-
\left(\frac{r}{\xi}\right)^2\right)$. To make our work amenable to analytic treatment, we take the following ansatz for the vortex profile: $$n(r) = n_0 (1-e^{-r^2/\lambda^2})$$ where $\lambda$ is a parameter on the order of the coherence length. Note that while this ansatz has the correct form close to the vortex core, the long distance behavior differs considerably. Fortunately our problem of vortex locking is dominated by the short-distance behavior, and we choose $\lambda$ so that the two densities (numerical and ansatz) agree at $r=\xi$ which requires $\lambda=1.781 \xi$, as shown in Fig. \[Fig:vortprofile\].

Our positive results, confirming the vortex-locked state, also confirm our intuition that the origin of the phenomena is in the short-ranged attraction between vortices. As explained in Ref. (but without proof), the algebraic decay of the superfluid order parameter of a single votex does not imply that vortices of one species, when in a lattice, exhibit a power-law decaying force on the vortices on the other species. Unlike the single-species vortex-vortex force, which is the result of the inductive (kinetic) energy term in the Gross-Pitaevskii equations, the interspecies force is a result of a density-density interaction. The density suppression due to a single vortex occurs since the superflow of the vortex effectively increases the mass terms $V_1,\,V_2$ in Eq. (\[GGPP\]). But in a lattice of vortices, the combined superflow vector is nearly zero (i.e., negligible compared to $\hbar
/m_{\alpha}\xi_{\alpha}$), and, therefore, so is the respective density suppression.

[20]{} natexlab\#1[\#1]{}bibnamefont \#1[\#1]{}bibfnamefont \#1[\#1]{}citenamefont \#1[\#1]{}url \#1[`#1`]{}urlprefix\[2\][\#2]{} \[2\]\[\][[\#2](#2)]{}

, , , , ****, ().

, ****, ().

, in **, edited by ().

, , , , , , ****, ().

, , , , ****, ().

, , , , ****, ().

, ****, ().

, ** (, ).

, ****, ().

, , , ****, ().

, ****, ().

, , , ****, ().

, ****, ().

, , , ****, ().

, ****, ().

, ****, ().

, , , , , ****, ().

, , , , , , ****, ().

, , , , , , , ****, ().

, , , ****, ().

[^1]: Such a value would be obtained, for instance, for a condensate of $3\times 10^6$ $^{87}$Rb atoms with condensate thickness $d_z=10 \mu m$. This then gives $\tilde{g}= \frac{4\pi a}{d_z} N \approx
    2 \times 10^4 $ where we take $a=106 a_0$ for the $s$-wave scattering length.