# Vortex dynamics in coherently coupled Bose-Einstein condensates¶

This notebook is the computational appendix of arXiv:1609.03966. We spell out in detail how to obtain the numerical results in the manuscript. We rely on Trotter-Suzuki-MPI, a massively parallel solver for the Gross--Pitaevskii equation, and its Python wrapper. The calculations use only a single computer, but it can take many hours to finish them. For the time-consuming parts, the execution time is mentioned running on an i7-4790 CPU with four physical (eight logical) cores. The total time to run the notebook is about seven hours.

## Table of Contents¶

Preliminaries

Precession of two vortices in a two-component BEC

Transfer of a single vortex between a two-component BEC

Oscillations of the relative density

## Important¶

Throughout the notebook we use a negative value for the Rabi coupling constant. This is needed since the Hamiltonian of the Rabi coupling implemented in trottersuzuki and the Hamiltonian used in the manuscript differ from a sign.

## Preliminaries¶

First we import the necessary modules and ensure that we get identical behaviour in Python 2 and 3.

In [1]:
from __future__ import print_function, division
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
from matplotlib.lines import Line2D
from matplotlib import colors
from matplotlib.collections import LineCollection
import numpy as np
import pandas
from scipy.interpolate import spline, InterpolatedUnivariateSpline
from scipy.optimize import curve_fit
from scipy.integrate import odeint
import time
import trottersuzuki as ts
import seaborn as sns  # prettifies plots
sns.set_context("notebook", font_scale=2., rc={"lines.linewidth": 3})
%matplotlib inline
KERNEL = "cpu"

## Two vortices in a two-component BEC¶

The following function takes two parameters $l_{\Omega}/r_{12}$ and $R$, and returns $\Omega_{rot}/\Omega$. The system is inizialized with two co-rotating vortices located symmetrically across the center of the container. The vortices, and the corresponding domain wall in the relative phase between them, are obtained performing a short imaginary time evolution, which proceeds along the following steps:

1. We start with normalized wavefunctions which take a constant value inside the circular well, $\psi_1 = \psi_2 =\sqrt{1/\pi R^2}$, and vanish outside.
2. We phase-imprint two co-rotating vortices, one per component.
3. We start the imaginary time evolution, in presence of an additional pinning potential (two sharply peaked Gaussians) aimed at keeping the vortex cores stationary (else, they would approach each other during the imaginary time evolution).

Once the gas has stabilized, we remove the pinning potential, and we let the system evolve in real time. The precession frequency $\Omega_{rot}$ is obtained by averaging typically over $\sim5$ full revolutions.

In [2]:
def two_vortices(lOmega_over_r12, radius):
###
#  Definition of variables and functions for the simulation
###

## Fixed variables
length = 2. * radius                           # physical lenght of the lattice's side
initial_vortices_distance = 1.0                # physical distance between vortices cores
n_0 = 4 / (np.pi * length**2)                  # density at the center of the circular box
if lOmega_over_r12 >= 0.5:                     # set r_{12} / \xi
r12_over_Xi12 = 10.
else:
r12_over_Xi12 = 40.

### Lattice parameters
const = 4.                                     # Ensure that one coherence length will be equal to "const" lattice spacings
dim = int(length * const * r12_over_Xi12) + 1  # linear dimension of the lattice

### Hamiltonian parameters
g_a = 10. * r12_over_Xi12**2 / (2. * n_0) # intra-particle coupling (first component)
g_b = 10. * r12_over_Xi12**2 / (2. * n_0) # intra-particle coupling (second component)
g_ab = 0                                  # inter-particle coupling
omega_rabi = -1. / (lOmega_over_r12**2)   # Rabi coupling (it is negative since the library use a different sign convention)

### Solver parameters
delta_t = 2.5e-5  # time for a single step evolution

# Definition of the pinning potential (only for imaginary time evolution)
gaussian_radius = 0.5/r12_over_Xi12
gaussian_height = 100000

def external_potential_imag_a(x, y):
if (x**2 + y**2) >= (length/2.05)**2:
return 1e10
elif np.sqrt((x-1/2)**2+y**2) <= gaussian_radius:
return gaussian_height*np.exp(- ((x-1/2)**2+y**2) / gaussian_radius**2) / (np.pi * gaussian_radius**2)
else:
return 0.

def external_potential_imag_b(x, y):
if (x**2 + y**2) >= (length/2.05)**2:
return 1e10
elif np.sqrt((x+1/2)**2+y**2) <= gaussian_radius:
return gaussian_height*np.exp(- ((x+1/2)**2+y**2) / gaussian_radius**2) / (np.pi * gaussian_radius**2)
else:
return 0.

# Definition of the circular well external potential (only for real time evolution)
def external_potential(x, y):
if (x**2 + y**2) >= (length/2.1)**2:
return 1e10
else:
return 0.

# Normalized initial state's wave function
def const_state(x, y):
return 1./length

# Define the vortex in component 1
def vortex_a(x, y):
z = (x-initial_vortices_distance/2.) + 1j*y
angle = np.angle(z)
return np.exp(1j * angle)

# Define the vortex in component 2
def vortex_b(x, y):
z = (x+initial_vortices_distance/2.) + 1j*y
angle = np.angle(z)
return np.exp(1j * angle)

###
#  Set up of the simulation
###

# Set the geometry of the simulation
grid = ts.Lattice2D(dim, length)

# Set the Hamiltonian
potential_a = ts.Potential(grid)  # Set the pinning potential for the first component
potential_a.init_potential(external_potential_imag_a)
potential_b = ts.Potential(grid)  # Set the pinning potential for the second component
potential_b.init_potential(external_potential_imag_b)
hamiltonian = ts.Hamiltonian2Component(grid, potential_a, potential_b, 1., 1.,
g_a, g_ab, g_b, omega_rabi)

# Set the initial state
state_a = ts.State(grid)        # Initialize the state in the first component
state_a.init_state(const_state)
state_b = ts.State(grid)        # Initialize the state in the second component
state_b.init_state(const_state)

# Imprint the vortices on the corresponding components
state_a.imprint(vortex_a)
state_b.imprint(vortex_b)

# Initialize the solver
solver = ts.Solver(grid, state_a, hamiltonian, delta_t, State2=state_b, kernel_type=KERNEL)

###
#  Set up imaginary time evolution
###

# Calculate initial vortices distance
coord_a = ts.get_vortex_position(grid, state_a, length*0.45)
coord_b = ts.get_vortex_position(grid, state_b, length*0.45)
vortex_distance = np.sqrt((coord_b[0]-coord_a[0])**2 + (coord_b[1]-coord_a[1])**2)

# Decide how many imaginary iterations to perform
# NB: the higher the Rabi coupling, the higher the iterations needed to form the domain wall
iterations = 1000
if lOmega_over_r12 < .2:
num = 1
elif lOmega_over_r12 < .5:
num = 4
elif lOmega_over_r12 < 1.:
num = 15
else:
num = 24

# Start imaginary time evolution
for i in range(num):
solver.evolve(iterations, True)         # Perform imaginary time evolution

# Calculate vortices distance after imaginary time evolution
coord_a = ts.get_vortex_position(grid, state_a, length*0.45)
coord_b = ts.get_vortex_position(grid, state_b, length*0.45)
vortex_distance = np.sqrt((coord_b[0]-coord_a[0])**2 + (coord_b[1]-coord_a[1])**2)

# Check whether the vortices distance has significantly changed and if true abort
if np.abs(vortex_distance - initial_vortices_distance) > grid.delta_x:
return

###
#  Set up real time evolution
###

# Change external potential to the circular wall without the pinning potential
potential_a.init_potential(external_potential)
potential_b.init_potential(external_potential)
solver.update_parameters()

# Iterations to perform to get a tot_angle
max_it_real_time = 20
if lOmega_over_r12 < .5:
iterations = 30
elif lOmega_over_r12 < 1.:
iterations = 8000
else:
iterations = 16000

angles_12 = []   # List of precession angles
times = []       # List of times
const_angle = 0. # Angle offset

# Start real time evolution
for cont in range(max_it_real_time):
solver.evolve(iterations)

# Calculate the precession angle of the vector r_12
coord_a = ts.get_vortex_position(grid, state_a, length*0.45)  # Get position of the vortex in the first component
coord_b = ts.get_vortex_position(grid, state_b, length*0.45)  # Get position of the vortex in the second component

angles = np.angle(coord_a[0] - coord_b[0] + 1j * (coord_a[1] - coord_b[1])) # Calculate precession angle

# Add an offset
if cont > 2:
last_angle = angles_12[len(angles_12)-2]
else:
last_angle = angles
if abs(angles + const_angle - last_angle) > np.pi:
const_angle += 2. * np.pi * np.sign(last_angle)

angles_12.append(angles + const_angle)
times.append(delta_t * iterations * (cont+1))

###
#  Calculate precession frequency
###

# fit angles curve
fit_par = np.polyfit(times, angles_12, 1)  # fit_par[0] is the precession frequency

return fit_par[0] / np.abs(omega_rabi)  # return \Omega_{rot} / \Omega

It follows a function for the simulation of two counter-rotating vortices. The precession frequency and the translational velocity are calculated.

In [3]:
def two_counter_rotating_vortices(lOmega_over_r12, radius):
###
#  Definition of variables and functions for the simulation
###

## Fixed variables
length = 2. * radius                           # physical lenght of the lattice's side
initial_vortices_distance = 1.0                # physical distance between vortices cores
n_0 = 4 / (np.pi * length**2)                  # density at the center of the circular box
if lOmega_over_r12 >= 0.3:                     # set r_{12} / \xi
r12_over_Xi12 = 20.
else:
r12_over_Xi12 = 40.

### Lattice parameters
const = 2.
dim = int(length * const * r12_over_Xi12) + 1      # linear dimension of the lattice

### Hamiltonian parameters
g_a = r12_over_Xi12**2 / (2. * n_0) # intra-particle coupling (firt component)
g_b = r12_over_Xi12**2 / (2. * n_0) # intra-particle coupling (second component)
g_ab = 0                                  # inter-particle coupling
omega_rabi = -1. / (lOmega_over_r12**2)   # Rabi coupling (it is negative since the library use a different sign convention)

### Solver parameters
delta_t = 2.5e-5  # time for a single step evolution

# Definition of the pinning potential (only for imaginary time evolution)
gaussian_radius = 0.5/r12_over_Xi12
gaussian_height = 100000

def external_potential_imag_a(x, y):
if (x**2 + y**2) >= (length/2.05)**2:
return 1e10
elif np.sqrt((x-1/2)**2+y**2) <= gaussian_radius:
return gaussian_height*np.exp(- ((x-1/2)**2+y**2) / gaussian_radius**2) / (np.pi * gaussian_radius**2)
else:
return 0.

def external_potential_imag_b(x, y):
if (x**2 + y**2) >= (length/2.05)**2:
return 1e10
elif np.sqrt((x+1/2)**2+y**2) <= gaussian_radius:
return gaussian_height*np.exp(- ((x+1/2)**2+y**2) / gaussian_radius**2) / (np.pi * gaussian_radius**2)
else:
return 0.

# Definition of the circular well external potential (only for real time evolution)
def external_potential(x, y):
if (x**2 + y**2) >= (length/2.05)**2:
return 1e10
else:
return 0.

# Normalized initial state's wave function
def const_state(x, y):
return 1./length

# Define the vortex in component 1
def vortex_a(x, y):
z = (x-initial_vortices_distance/2.) + 1j*y
angle = np.angle(z)
return np.exp(1j * angle)

# Define the vortex in component 2
def vortex_b(x, y):
z = (x+initial_vortices_distance/2.) + 1j*y
angle = np.angle(z)
return np.exp(-1j * angle)

###
#  Set up of the simulation
###

# Set the geometry of the simulation
grid = ts.Lattice2D(dim, length)

# Set the Hamiltonian
potential_a = ts.Potential(grid)  # Set the pinning potential for the first component
potential_a.init_potential(external_potential_imag_a)
potential_b = ts.Potential(grid)  # Set the pinning potential for the second component
potential_b.init_potential(external_potential_imag_b)
hamiltonian = ts.Hamiltonian2Component(grid, potential_a, potential_b, 1., 1.,
g_a, g_ab, g_b, omega_rabi)

# Set the initial state
state_a = ts.State(grid)        # Initialize the state in the first component
state_a.init_state(const_state)
state_b = ts.State(grid)        # Initialize the state in the second component
state_b.init_state(const_state)

# Imprint the vortices on the corresponding components
state_a.imprint(vortex_a)
state_b.imprint(vortex_b)

# Initialize the solver
solver = ts.Solver(grid, state_a, hamiltonian, delta_t, State2=state_b)

###
#  Set up imaginary time evolution
###

# Calculate initial vortices distance
coord_a = ts.get_vortex_position(grid, state_a, length*0.45)
coord_b = ts.get_vortex_position(grid, state_b, length*0.45)
vortex_distance = np.sqrt((coord_b[0]-coord_a[0])**2 + (coord_b[1]-coord_a[1])**2)

# Decide how many imaginary iterations to perform
# NB: the higher the Rabi coupling, the higher the iterations needed to form the domain wall
iterations = 50
num = 20

# Start imaginary time evolution
for i in range(num):
solver.evolve(iterations, True)         # Perform imaginary time evolution
time = iterations * (i + 1) * delta_t

# Calculate vortices distance after imaginary time evolution
coord_a = ts.get_vortex_position(grid, state_a, length*0.45)
coord_b = ts.get_vortex_position(grid, state_b, length*0.45)
vortex_distance = np.sqrt((coord_b[0]-coord_a[0])**2 + (coord_b[1]-coord_a[1])**2)

###
#  Set up real time evolution
###

# Change external potential to the circular wall without the pinning potential
potential_a.init_potential(external_potential)
potential_b.init_potential(external_potential)
solver.update_parameters()

# Iterations
max_it_real_time = 14
if lOmega_over_r12 <= .1:
iterations = 30
elif lOmega_over_r12 <= .2:
iterations = 200
elif lOmega_over_r12 <= .4:
iterations = 500
else:
iterations = 1000

angles_12 = []   # List of precession angles
times = []       # List of times
const_angle = 0. # Angle offset
time = 0.
position_x_a = []
position_y_a = []
position_x_b = []
position_y_b = []

# Start real time evolution
for cont in range(max_it_real_time):
solver.evolve(iterations)

# Calculate the precession angle of the vector r_12
coord_a = ts.get_vortex_position(grid, state_a, length*0.45)  # Get position of the vortex in the first component
coord_b = ts.get_vortex_position(grid, state_b, length*0.45)  # Get position of the vortex in the second component
position_x_a.append(coord_a[0])
position_y_a.append(coord_a[1])
position_x_b.append(coord_b[0])
position_y_b.append(coord_b[1])

angles = np.angle(coord_a[0] - coord_b[0] + 1j * (coord_a[1] - coord_b[1])) # Calculate precession angle

# Add an offset
if cont > 2:
last_angle = angles_12[len(angles_12)-2]
else:
last_angle = angles
if abs(angles + const_angle - last_angle) > np.pi:
const_angle += 2. * np.pi * np.sign(last_angle)

angles_12.append(angles + const_angle)
times.append(delta_t * iterations * (cont+1))
time = iterations * (cont + 1) * delta_t

###
#  Calculate precession frequency
###

# fit angles curve
times = np.array(times)
angles_12 = np.array(angles_12)
times = times[~np.isnan(angles_12)]
angles_12 = angles_12[~np.isnan(angles_12)]
fit_par = np.polyfit(times[0:20], angles_12[0:20], 1)  # fit_par[0] is the precession frequency

###
#  Calculate translational velocity
###
T = 2. * np.pi / np.fabs(omega_rabi)
time = np.arange(1,max_it_real_time+1) * iterations * delta_t / T
fit_par_translationalVelocity_a = np.polyfit(time, position_y_a, 1)
fit_par_translationalVelocity_b = np.polyfit(time, position_y_b, 1)
mean_translationalVelocity = 0.5*(fit_par_translationalVelocity_a[0]+fit_par_translationalVelocity_a[0])

return fit_par[0] / np.fabs(omega_rabi), mean_translationalVelocity

Simulation for co-rotating vortices. We explore the parameter space $l_{\Omega} / r_{12} \in [0.1, 10]$ and $R = 2$. The manuscript also has results for $R=3.5, 5$, but the calculations take a long time, so here only one value is evaluated.

In [4]:
interval = np.linspace(-1, 1, 21)
parameters = [2]
# parameters = [2, 3.5, 5]
df = pandas.DataFrame(index=range(len(interval)*len(parameters)),
columns=["length_x", "l_omega_r12", "omega_rot/omega"], dtype=np.float64)
i = 0
time0 = time.time()
for lOmega_over_r12 in np.power(10, interval):
for radius in parameters:
time_i = time.time()
df.iloc[i] = [2*radius, lOmega_over_r12, two_vortices(lOmega_over_r12, radius)]
print("Execution time:", time.time()-time_i)
i += 1
print("Total time:", time.time()-time0)
Execution time: 124.72676301002502
Execution time: 115.6088182926178
Execution time: 113.28329038619995
Execution time: 107.32662296295166
Execution time: 170.9677221775055
Execution time: 168.91597509384155
Execution time: 171.0564215183258
Execution time: 578.9320180416107
Execution time: 549.2432420253754
Execution time: 516.2964489459991
Execution time: 998.1947374343872
Execution time: 1012.504718542099
Execution time: 1160.1199383735657
Execution time: 1554.923291683197
Execution time: 1018.2498126029968
Execution time: 976.794002532959
Execution time: 978.2081477642059
Execution time: 978.4028224945068
Execution time: 990.4037802219391
Execution time: 1100.2374074459076
Execution time: 1572.6919741630554
Total time: 14957.08951497078

It follows the simulation for counter-rotating vortices. We explore the parameter space $l_{\Omega} / r_{12} \in [0.1, 0.5]$ and $R = 5$.

In [5]:
radius = 5.
LOmega_over_r12 = np.array([0.5, 0.4, 0.3, 0.2, 0.1])
Omega_rot = []
TranslationalVelocity = []
time0 = time.time()
for lOmega_over_r12 in LOmega_over_r12:
time_i = time.time()
omega_rot, translationalVelocity = two_counter_rotating_vortices(lOmega_over_r12, radius)
Omega_rot.append(omega_rot)
TranslationalVelocity.append(translationalVelocity)
print("Execution time:", time.time()-time_i)

print("Total time:", time.time()-time0)
Execution time: 241.59611558914185
Execution time: 139.57053637504578
Execution time: 136.5041525363922
Execution time: 218.66592621803284
Execution time: 135.8393702507019
Total time: 872.1766188144684

Fig. 1a plots the precession frequency of two co-rotating vortices, one per component, along with the ones for the counter-rotating vortices. Colors indicate results obtained for different radii $R$ of the circular container. The dash-dotted line is the strong-coupling limit.The following script produces the figure:

In [6]:
plt.figure(figsize=(9, 6))
mycolors = {4:"#6EE1F7", 7:"#43CA2E", 10:"#861010"}

for radius in parameters:
length_x = int(2*radius)
ordered_left = df[(df['length_x'] == length_x) &\
(df['l_omega_r12'] <= 0.50118723) &\
(df['l_omega_r12'] != 0.79432823) &\
(df['l_omega_r12'] != 0.63095734)\
].sort('l_omega_r12')
ordered_right = df[(df['length_x'] == length_x) &\
(df['l_omega_r12'] >= 0.50118723) &\
(df['l_omega_r12'] != 0.79432823) &\
(df['l_omega_r12'] != 0.63095734)\
].sort('l_omega_r12')
x_l = ordered_left['l_omega_r12']
y_l = ordered_left['omega_rot/omega']
x_r = ordered_right['l_omega_r12']
y_r = ordered_right['omega_rot/omega']
x_rm = x_r.as_matrix()
y_rm = y_r.as_matrix()
x_r2 = np.linspace(x_rm.min(), x_rm.max(), 300)
y_r2 = spline(x_rm, y_rm, x_r2)
plt.plot(x_l, y_l, "--", dashes = [12, 6], label="", color=mycolors[length_x])
plt.plot(x_r2, y_r2, label='$R/r_{12}=%.1f$' % radius, color=mycolors[length_x])

# add the strong coupling limit
x = np.linspace(0, 2, 100)
y = -4*np.sqrt(2)*x/np.pi
plt.plot(x, y, '-.', label="$\mathrm{strong} \,\,\, \mathrm{coupling}$");

# add counter-rotating vortices
plt.plot(LOmega_over_r12,Omega_rot,'D', color=mycolors[4], lw=2, label="$\mathrm{counter-rotating} \,\,\, \mathrm{vortices}$");

# set labels and limits
plt.xscale('log')
plt.xlabel("$l_\Omega/r_{12}$")
plt.ylabel("$\Omega_\mathrm{rot}/\Omega$")
plt.xlim(0.1, 10)
plt.ylim(-2, 4)
plt.yticks([-2, -1, 0, 1, 2, 3, 4], ['$-2$', '$-1$', '$0$', '$1$', '$2$', '$3$', '$4$'])
plt.xticks([0.1, 1, 10], ['$10^{-1}$', '$10^{0}$', '$10^{1}$', ])
plt.axhline(0, color='black', linewidth=1)
lgd = plt.legend(bbox_to_anchor=(0, 1), loc=2, borderaxespad=0.);
plt.show()
/export/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:6: FutureWarning: sort(columns=....) is deprecated, use sort_values(by=.....)
/export/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:7: FutureWarning: sort(columns=....) is deprecated, use sort_values(by=.....)

The same data are plotted with different axes, to highlight the behavior at weak-coupling. Here $\omega_{12} \equiv \hbar/M r_{12}^2$, and the dots are the result expected for a single vortex in a single component BEC inside a cylinder. This scripts generates Fig. 1b in the manuscript:

In [7]:
plt.figure(figsize=(9, 6))

for radius in parameters:
length_x = int(2*radius)
ordered = df[df['length_x'] == length_x].sort('l_omega_r12')
x = 1. / (ordered['l_omega_r12'])**2
y = x * ordered['omega_rot/omega']
if length_x == 4:
label='$R/r_{12}=2$'
if length_x == 7:
label='$R/r_{12}=3.5$'
if length_x == 10:
label='$R/r_{12}=5$'

plt.plot(x, y, label='$R/r_{12}=%.1f$' % radius, color=mycolors[length_x])
# positions to inter/extrapolate
xi = np.linspace(0, .25, 10)
# spline order: 1 linear, 2 quadratic, 3 cubic ...
order = 1
# do inter/extrapolation   (the data must be monotonously increasing)
data_x = (np.array(x)[::-1])[:7]
data_y = (np.array(y)[::-1])[:7]
s = InterpolatedUnivariateSpline(data_x, data_y, k=order)
yi = s(xi)
# add the result for a single off-centered vortex in a cylindrical container; see Pethick&Smith, Chap.(9.4), Eq.(9.46)
b = 0.5;
plt.plot([0], [1 / (((length_x/2)**2 - b**2))], 'o', label="", markersize=10, color=mycolors[length_x])

# set labels and limits
plt.xlabel("$\Omega/\omega_{12}$")
plt.ylabel("$\Omega_\mathrm{rot}/\omega_{12}$")
plt.xlim(-0.02, .5)
plt.xticks([0., .25, 0.5], ['$0$', '$0.25$', '0.5'])
plt.ylim(-.8, .4)
plt.yticks([-.6, -.3, 0, .3], ['$-0.6$', '$-0.3$', '$0$', '$0.3$']);
plt.axhline(0, color='black', linewidth=1);
plt.axvline(0, color='black', linewidth=1);
lgd = plt.legend(bbox_to_anchor=(0.65, 1), loc=2, borderaxespad=0.);
plt.show()
/export/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:5: FutureWarning: sort(columns=....) is deprecated, use sort_values(by=.....)

The translational velocity of two counter-rotating vortices is shown and compared with the theoretical result for strong coupling. See Fig. 2 in the manuscript.

In [8]:
plt.figure(figsize=(9,6));
plt.plot(LOmega_over_r12, np.fabs(TranslationalVelocity),'D',label='$R/r_{12}=5$', color="#861010")

# theoretical curve
y = 4.*np.sqrt(2)*np.array([0,1])
plt.plot(np.array([0,1]),y,'-', dashes=[15, 4, 2, 4], lw=2, label="$\mathrm{strong} \, \mathrm{coupling}$")

# set labels and limits
plt.xlabel("$l_{\Omega}/r_{12}$")
plt.ylabel("$v_{\mathrm{pair}} T / r_{12}$")
plt.xlim(0.05,.52)
plt.ylim(0,3)
plt.yticks([0,.5,1,1.5,2,2.5,3], ['$0$','$0.5$','$1$','$1.5$','$2$','$2.5$','$3$'])
plt.xticks([0.1, 0.2, 0.3, 0.4, 0.5],['$0.1$', '$0.2$', '$0.3$','$0.4$','$0.5$'])
plt.axhline(0, color='black',linewidth=1);
plt.axvline(0, color='black',linewidth=1);
lgd = plt.legend(bbox_to_anchor=(0.02, 0.95), loc=2, borderaxespad=0.);
plt.show()
sns.set_style("white")  # set white style for the following plots

## Transfer of a single vortex between a two-component BEC¶

We consider a two-component BEC with equal populations ($N_1 = N_2 = 1$), in a harmonic trap of frequency $\omega$ to illustrate the coherent oscillations of vorticity induced by the Rabi coupling $\Omega$. First we define some helper functions to compute the theoretical trajectory of the vortex:

In [12]:
def dynamic_equations(y, t, Omega, kappa, R, xi):
u, phi = y
c = 4./3. - 128./(45.*np.pi)
uDot = -.25*Omega*(4./3.-c*u**2) * np.sin(phi-kappa*t) / (1.-u**2)
phiDot = (1./(R**2 * (1.-u**2))*(np.log(R/xi) + .5*np.log(1-u**2) - .5)\
- .25*Omega /(u*(1-u**2)) * (4./3. - 3.*c*u**2)* np.cos(phi-kappa*t))
return [uDot, phiDot]

def dynamic_equations_ret(y, t, Omega, kappa, R, xi):
u, phi = y
c = 4./3. - 128./(45.*np.pi)
uDot = .25*Omega*(4./3.-c*u**2) * np.sin(phi+kappa*t) / (1.-u**2)
phiDot = -(1./(R**2 * (1.-u**2))*(np.log(R/xi) + .5*np.log(1-u**2) - .5)\
- .25*Omega /(u*(1-u**2)) * (4./3. - 3.*c*u**2)* np.cos(phi+kappa*t))
return [uDot, phiDot]

def get_theoretical_trajectory(cycle1, time1, cycle2, time2, Rabi_period, Omega_Rabi, kappa, TF_radius, xi_A):
u = np.array([])
phi = np.array([])
# first component
if cycle1 < 1.:
y0 = [0.0001, 0.]                           # initial conditions
time_series = np.arange(0., Rabi_period/2., 0.01)           # time of evolution
# integrate dynamic equations
r = odeint(dynamic_equations, y0, time_series, args=(Omega_Rabi, kappa, TF_radius, xi_A))
time1_series = time_series
if cycle1 >= 1:
y0 = [0.0001, 0.]                                     # initial conditions
time_series = np.arange(0., .4*Rabi_period, 0.01) - time1        # time of evolution
# integrate dynamic equations
r = odeint(dynamic_equations_ret, y0, time_series, args=(Omega_Rabi, kappa, TF_radius, xi_A))
u = r[:,0]
phi = r[:,1]
y0 = [0.0001, 0.]                                    # initial conditions
time_series = np.arange(0., .4*Rabi_period, 0.01) + time1        # time of evolution
# integrate dynamic equations
r = odeint(dynamic_equations, y0, time_series, args=(Omega_Rabi, kappa, TF_radius, xi_A))
time1_series = np.concatenate((time_series - 0.4*Rabi_period, time_series), axis=0)
u = np.concatenate((u[::-1], r[:,0]))
phi = np.concatenate((phi[::-1], r[:,1]))
vTraj_x1 = []
vTraj_y1 = []
for idx in range(len(u)):
vTraj_x1.append(u[idx] * np.cos(phi[idx]))
vTraj_y1.append(u[idx] * np.sin(phi[idx]))

# second component
y0 = [0.0001, 0]                                     # initial conditions
time_series = np.arange(0., .4*Rabi_period, 0.01) - time2        # time of evolution
# integrate dynamic equations
r = odeint(dynamic_equations_ret, y0, time_series, args=(Omega_Rabi, kappa, TF_radius, xi_A))
u = r[:,0]
phi = r[:,1]
y0 = [0.0001, 0]                                    # initial conditions
time_series = np.arange(0., .4*Rabi_period, 0.01) + time2        # time of evolution
# integrate dynamic equations
r = odeint(dynamic_equations, y0, time_series, args=(Omega_Rabi, kappa, TF_radius, xi_A))
time2_series = np.concatenate((time_series - 0.4*Rabi_period, time_series), axis=0)

u = np.concatenate((u[::-1], r[:,0]))
phi = np.concatenate((phi[::-1], r[:,1]))
vTraj_x2 = []
vTraj_y2 = []
for idx in range(len(u)):
vTraj_x2.append(u[idx] * np.cos(phi[idx]))
vTraj_y2.append(u[idx] * np.sin(phi[idx]))

return np.array(vTraj_x1), np.array(vTraj_y1), np.array(vTraj_x2), np.array(vTraj_y2), time1_series, time2_series

We prepare the system by phase-imprinting a single vortex with positive circulation in the center of the first component, and we find the corresponding ground state by performing a short evolution in imaginary time in absence of Rabi coupling, which allows the formation of a vortex core with the suitable profile at the center of the first component. After equilibration, we switch on the Rabi coupling, and let the system evolve in real time for a variable time $t>0$. The simulation is performed by the following function:

In [9]:
def single_vortex(dim, length, Rabi_period, coupling, Rabi_periods_to_be_simulated, gAB_over_g,
densities1, densities2, phases1, phases2, coord1x, coord1y, coord2x, coord2y):

TF_radius = (4*coupling/np.pi)**0.25
grid = ts.Lattice2D(dim, length)

## Define the Hamiltonian
potential = ts.HarmonicPotential(grid, 1., 1.)

g1, g2 = coupling, coupling
g12 = gAB_over_g * coupling
omega = 0
hamiltonian = ts.Hamiltonian2Component(grid, potential, potential, 1., 1., g1, g12, g2, omega)

## Define the initial state
width = max(TF_radius, 1)
w = 0.5 / width**2
state1, state2 = ts.GaussianState(grid, w), ts.GaussianState(grid, w)

## Imprint vortex phase
def vortex_phase(x, y):
return np.exp(1j * np.angle(x + 1j*y))

state1.imprint(vortex_phase)

## Initialize the solver
frames_per_Rabi_period = 30
real_iterations = 700  # number of iterations performed by each instance of ts
deltaT = Rabi_period/(real_iterations * frames_per_Rabi_period)

solver = ts.Solver(grid, state1, hamiltonian, deltaT, State2=state2, kernel_type=KERNEL)

#################################################################
## IMAGINARY EVOLUTION
iterations = 2000
NumEvo = 4
for i in range(0, NumEvo):
solver.evolve(iterations, True)

#################################################################
## REAL-TIME EVOLUTION
omega = -2*np.pi / Rabi_period
hamiltonian.omega_r = omega  # Switch on the Rabi coupling
solver.update_parameters()

iterations = real_iterations
NumEvo = int(frames_per_Rabi_period * Rabi_periods_to_be_simulated)

coord1 = ts.get_vortex_position(grid, state1, width).reshape((2, 1))
coord2 = ts.get_vortex_position(grid, state2, width).reshape((2, 1))
time_coord = np.array([0.])

### coherence length of BEC_A
xi_B = 1./np.sqrt(2 * g2 * state2.get_particle_density().max()/(grid.delta_x**2))
xi_A = xi_B
kappa = .485

## Calculate theoretical trajectories
cycle1 = 0.
time1 = 0.
cycle2 = 0.
time2 = .5 * Rabi_period
vTraj_x1, vTraj_y1, vTraj_x2, vTraj_y2,\
time1_s, time2_s = get_theoretical_trajectory(cycle1, time1, cycle2, time2, Rabi_period,
np.abs(omega), kappa, TF_radius, xi_A)

## Collects all theoretical trajectories
vTrajTot_x1 = np.array([])
vTrajTot_y1 = np.array([])
vTrajTot_x2 = np.array([])
vTrajTot_y2 = np.array([])
time1Tot_s = np.array([])
time2Tot_s = np.array([])

vTrajTot_x1 = np.concatenate((vTrajTot_x1, vTraj_x1), axis=0)
vTrajTot_y1 = np.concatenate((vTrajTot_y1, vTraj_y1), axis=0)
vTrajTot_x2 = np.concatenate((vTrajTot_x2, vTraj_x2), axis=0)
vTrajTot_y2 = np.concatenate((vTrajTot_y2, vTraj_y2), axis=0)
time1Tot_s = np.concatenate((time1Tot_s, time1_s), axis=0)
time2Tot_s = np.concatenate((time2Tot_s, time2_s), axis=0)

## Collects particle densities and phases
densities1.append(state1.get_particle_density())
densities2.append(state2.get_particle_density())
phases1.append(state1.get_phase())
phases2.append(state2.get_phase())

## Collects simulated trajectories
coord1x.append(coord1[0])
coord1y.append(coord1[1])
coord2x.append(coord2[0])
coord2y.append(coord2[1])

for i in range(0, NumEvo):
solver.evolve(iterations)
if (i+1 == frames_per_Rabi_period//4 or i+1 == frames_per_Rabi_period//2 or\
i+1 == 3*frames_per_Rabi_period//4 or i+1 == frames_per_Rabi_period):
densities1.append(state1.get_particle_density())
densities2.append(state2.get_particle_density())
phases1.append(state1.get_phase())
phases2.append(state2.get_phase())
coord1x.append(coord1[0])
coord1y.append(coord1[1])
coord2x.append(coord2[0])
coord2y.append(coord2[1])

_coord1 = ts.get_vortex_position(grid, state1, width).reshape((2, 1))
_coord2 = ts.get_vortex_position(grid, state2, width).reshape((2, 1))
coord1 = np.concatenate((coord1, _coord1.reshape((2, 1))), axis=1)
coord2 = np.concatenate((coord2, _coord2.reshape((2, 1))), axis=1)
time = deltaT*(i+1)*iterations
time_coord = np.concatenate((time_coord, np.array([time])), axis=0)

if ((not np.isnan(coord1[0][len(coord1[0])-1])) and np.isnan(coord1[0][len(coord1[0])-2])) or\
((not np.isnan(coord2[0][len(coord2[0])-1])) and np.isnan(coord2[0][len(coord2[0])-2])):
cycle1 = (time + .5*Rabi_period) // Rabi_period
time1 = cycle1 * Rabi_period
cycle2 = time // Rabi_period
time2 = (cycle2 + .5) * Rabi_period
vTraj_x1, vTraj_y1, vTraj_x2, vTraj_y2,\
time1_s, time2_s = get_theoretical_trajectory(cycle1, time1, cycle2, time2, Rabi_period,
np.abs(omega), kappa, TF_radius, xi_A)
vTrajTot_x1 = np.concatenate((vTrajTot_x1, vTraj_x1), axis=0)
vTrajTot_y1 = np.concatenate((vTrajTot_y1, vTraj_y1), axis=0)
vTrajTot_x2 = np.concatenate((vTrajTot_x2, vTraj_x2), axis=0)
vTrajTot_y2 = np.concatenate((vTrajTot_y2, vTraj_y2), axis=0)
time1Tot_s = np.concatenate((time1Tot_s, time1_s), axis=0)
time2Tot_s = np.concatenate((time2Tot_s, time2_s), axis=0)

return coord1, coord2, time_coord, vTrajTot_x1, vTrajTot_y1, vTrajTot_x2, vTrajTot_y2, time1Tot_s, time2Tot_s

Fig. 3 plots the contours of equal dimensionless energy $\tilde E$, for $N_1 = N_2$, $\bar\phi=0$, and $R/\xi = 5$. From left to right, $M\Omega R^2/\hbar = 0,1,10$. The figure is generated by the following code:

In [10]:
fig = plt.figure(figsize=(14.5, 5.9))
fontsize = 32

TF_radius = 2.24638933492  # T-F radius for g=20
xi_A = 0.423384837924      # and coherence length

for indx in [0, 1, 2]:
plt.subplot(1, 3, indx+1)
if indx == 0:
Omega_Rabi = 0.
elif indx == 1:
Omega_Rabi = 1./TF_radius**2
else:
Omega_Rabi = 10./TF_radius**2

# vectors containing the physical coordinates
points = 1000
y_vec = np.linspace(-TF_radius, TF_radius, points, endpoint=False)
x_vec =np.linspace(-TF_radius, TF_radius, points, endpoint=False)

## Compute vortex energy for all the position in the mesh
X, Y = np.meshgrid(x_vec, y_vec)
U = np.sqrt(X**2 + Y**2)/TF_radius
N_A, N_B = 1, 1
c = 4./3 - 128./(45*np.pi)
Z_A = (1-U**2) * (2 * np.log(TF_radius/xi_A) + np.log(1-U**2)) + 2*U**2 - 1\
+  Omega_Rabi * TF_radius**2 * np.sqrt(N_B/N_A) * np.absolute(4*U/3 - c*U**3) * np.cos(np.angle(X+1.0j*Y))

#################################################################
## Make plot

# set ticks and labels
ticks_1, label_ticks_1 = [0., points/4., points/2., points*3/4., points],\
['$'+str(-1)+'$', '$'+str(-0.5)+'$', '$0$', '$'+str(0.5)+'$', '$'+str(1)+'$']
ticks_2, label_ticks_2 = [0., points/4., points/2., points*3/4., points],\
['', '', '', '', '']
if indx == 0:
plt.yticks(ticks_1, label_ticks_1, fontsize=fontsize)
plt.ylabel('$y/R$', fontsize=fontsize)
else:
plt.yticks(ticks_2, label_ticks_2, fontsize=fontsize)
plt.xticks(ticks_1, label_ticks_1, fontsize=fontsize)
plt.xlabel('$x/R$', fontsize=fontsize)

# set levels and position of the colorbars
if indx == 0:
levels = np.linspace(0.8, 2.6, 9, endpoint=False)
cbarpos = [0.1, 0.4, 0.5, 0.015]
labelcbar = '0'
elif indx == 1:
levels = np.linspace(0., 3.2, 8, endpoint=False)
cbarpos = [0.5, 0.5, 0.69, 0.015]
labelcbar = '1'
else:
levels = np.linspace(-9, 15, 8, endpoint=False)
cbarpos = [0.1, 1.1, 0.3, 0.015]
labelcbar = '2'
barticks = []
for jj in range(len(levels)):
if indx != 2:
barticks.append('$'+str(levels[jj])+'$')
else:
barticks.append('$%4.0f$' % (levels[jj]))

cont_ = plt.contourf(Z_A, levels, cmap='Blues')    # plot the energy levels

# set ticks of the colorbars
if indx == 0:
cbarpos = [0.13, .89, 0.245, 0.015]
labelcbar = '0'
cbaxes0 = fig.add_axes(cbarpos, label=labelcbar)
cbar = plt.colorbar(cont_, cax = cbaxes0, orientation='horizontal')
cbar.ax.tick_params(labelsize=fontsize-8)
cbar.ax.xaxis.set_ticks_position('top')
cbar.ax.set_xticklabels(barticks, rotation=45)
elif indx == 1:
cbarpos = [0.42, 0.89, 0.245, 0.015]
labelcbar = '1'
cbaxes1 = fig.add_axes(cbarpos, label=labelcbar)
cbar = plt.colorbar(cont_, cax = cbaxes1, orientation='horizontal')
cbar.ax.xaxis.set_ticks_position('top')
cbar.ax.tick_params(labelsize=fontsize-8)
cbar.ax.set_xticklabels(barticks, rotation=45)
else:
cbarpos = [0.73, .89, 0.245, 0.015]
labelcbar = '2'
cbaxes2 = fig.add_axes(cbarpos, label=labelcbar)
cbar = plt.colorbar(cont_, cax = cbaxes2, orientation='horizontal')
cbar.ax.xaxis.set_ticks_position('top')
cbar.ax.tick_params(labelsize=fontsize-8)
cbar.ax.set_xticklabels(barticks, rotation=45)

plt.tight_layout(rect=[0, 0, 1, 0.9])
plt.show()
/export/anaconda3/lib/python3.5/site-packages/matplotlib/figure.py:1742: UserWarning: This figure includes Axes that are not compatible with tight_layout, so its results might be incorrect.
warnings.warn("This figure includes Axes that are not "

Coherent quantum phase slip is shown in the following simulation. This is a time-consuming step running for about three minutes.

In [13]:
densities1, densities2 = [], []
phases1, phases2 = [], []
coord1x, coord1y = [], []
coord2x, coord2y = [], []

# linear size of the lattice
dim = 150
length = 9.

# number of Rabi periods to be simulated
Rabi_periods_to_be_simulated = 1

# start evolution
gAB_over_g =1.
Rabi_period = np.pi
coupling = 20.
time0 = time.time()
coord1, coord2, _, _, _, _, _, _, _ = single_vortex(dim, length, Rabi_period, coupling, Rabi_periods_to_be_simulated,
gAB_over_g, densities1, densities2, phases1, phases2, coord1x,
coord1y, coord2x, coord2y)
print("Total time:", time.time()-time0)
Total time: 88.06991386413574

We plot the densities and phases of the two components and superimpose the trajectory of the vortex's core (Fig. 4 in the manuscript): transfer of vorticity observed in a harmonically trapped two-component BEC, in presence of a Rabi frequency $\Omega = 2\omega$. In this simulation we used $gN=g_{12}N=40\hbar^2/M$, so that $R^2/l_\Omega^2 \approx10$ and $Na/a_\perp \approx 3$. From left to right, the columns show, respectively, the particle densities $n_1$ and $n_2$, the phase $S_1$ of $\psi_1$, and the phase difference $S_1-S_2$ at given times: from top to bottom, $t/T=\{0,\,0.25,\,0.5, \, 0.75,\,1\}$, with $T=2\pi/\Omega$ the Rabi period. Markers indicate the trajectory of the vortex core, up to the time at which the screenshot is taken (circles and squares for vortex core in first and second component respectively), while the continuous lines shown the theoretical trajectory. The color of the markers indicate at what (past) time the core was at that specific position.

In [14]:
############ Helper Functions ############
TF_radius_ = 2.24638933492
TF_radius = 2.24638933492
Rabi_period = np.pi
Omega_Rabi = 2.*np.pi/Rabi_period
g_A = 20.
g_B = g_A
g_AB = g_A
xi_A =  0.434474936715

kappa = .485

def dynamic_equations_plot(y, t, Omega, kappa, R, xi):
u, phi = y
c = 4./3. - 128./(45.*np.pi)
uDot = -.25*Omega*(4./3.-c*u**2) * np.sin(phi-kappa*t) / (1.-u**2)
phiDot = (1./(R**2 * (1.-u**2))*(np.log(R/xi) + .5*np.log(1-u**2) - .5)\
- .25*Omega /(u*(1-u**2)) * (4./3. - 3.*c*u**2)* np.cos(phi-kappa*t))
return [uDot, phiDot]

def dynamic_equations_ret_plot(y, t, Omega, kappa, R, xi):
u, phi = y
c = 4./3. - 128./(45.*np.pi)
uDot = .25*Omega*(4./3.-c*u**2) * np.sin(phi+kappa*t) / (1.-u**2)
phiDot = -(1./(R**2 * (1.-u**2))*(np.log(R/xi) + .5*np.log(1-u**2) - .5)\
- .25*Omega /(u*(1-u**2)) * (4./3. - 3.*c*u**2)* np.cos(phi+kappa*t))
return [uDot, phiDot]

def get_theoretical_trajectory_plot(cycle1, time1, time2):

u = np.array([])
phi = np.array([])
# first component
y0 = [0.0001, 0.]                           # initial conditions
time_series = np.arange(0., Rabi_period/2., 0.01)           # time of evolution
# integrate dynamic equations
r0 = odeint(dynamic_equations, y0, time_series, args=(Omega_Rabi, kappa, TF_radius, xi_A))
time1_series = time_series

if cycle1 >= 1:
y0 = [0.0001, 0.]                                     # initial conditions
time_series = np.arange(0., .4*Rabi_period, 0.01) - time1        # time of evolution
# integrate dynamic equations
r = odeint(dynamic_equations_ret, y0, time_series, args=(Omega_Rabi, kappa, TF_radius, xi_A))
u = r[:,0]
phi = r[:,1]
y0 = [0.0001, 0.]                                    # initial conditions

u = np.concatenate((u[::-1], r0[:,0]))
phi = np.concatenate((phi[::-1], r0[:,1]))
vTraj_x1 = []
vTraj_y1 = []
for idx in range(len(u)):
vTraj_x1.append(u[idx] * np.cos(phi[idx]))
vTraj_y1.append(u[idx] * np.sin(phi[idx]))

# second component
y0 = [0.0001, 0]                                     # initial conditions
time_series = np.arange(0., .4*Rabi_period, 0.01) - time2        # time of evolution
# integrate dynamic equations
r = odeint(dynamic_equations_ret, y0, time_series, args=(Omega_Rabi, kappa, TF_radius, xi_A))
u = r[:,0]
phi = r[:,1]
y0 = [0.0001, 0]                                    # initial conditions
time_series = np.arange(0., .4*Rabi_period, 0.01) + time2        # time of evolution
# integrate dynamic equations
r = odeint(dynamic_equations, y0, time_series, args=(Omega_Rabi, kappa, TF_radius, xi_A))
time2_series = np.concatenate((time_series - 0.4*Rabi_period, time_series), axis=0)

u = np.concatenate((u[::-1], r[:,0]))
phi = np.concatenate((phi[::-1], r[:,1]))
vTraj_x2 = []
vTraj_y2 = []
for idx in range(len(u)):
vTraj_x2.append(u[idx] * np.cos(phi[idx]))
vTraj_y2.append(u[idx] * np.sin(phi[idx]))

return np.array(vTraj_x1), np.array(vTraj_y1), np.array(vTraj_x2), np.array(vTraj_y2)

############ Plot ############
fig = plt.figure(figsize=(18, 20.5));

for row in range(0, 5):
# Calculate theoretical trajectories
if row<3:
cycle1 = 0
time1 = 0
time2 = Rabi_period / 2
else:
cycle1 = 1
time1 = Rabi_period
time2 = Rabi_period / 2

vTraj_x1, vTraj_y1, vTraj_x2, vTraj_y2 = get_theoretical_trajectory_plot(cycle1, time1, time2)

cont = row*10
fontsi = 23
rescale = 0.65
markersize = 90
imag_time =  False

length_x_ = 9.
delta_x = length_x_ / dim

# rescaling of the pcolor plots
dim_x = int(dim * rescale / 2.) *2
rescale_ = dim_x / dim
length_x = length_x_ * rescale_
TF_radius = TF_radius_ * rescale_

density_A_ = densities1[row]
density_B_ = densities2[row]
phase_A_ = phases1[row]
phase_B_ = phases2[row]

edge = (dim - dim_x) // 2
density_A = density_A_[edge:dim-edge, edge:dim-edge]
density_B = density_B_[edge:dim-edge, edge:dim-edge]
phase_A = phase_A_[edge:dim-edge, edge:dim-edge]
phase_B = phase_B_[edge:dim-edge, edge:dim-edge]

_coord1x = coord1x[row]
_coord1y = coord1y[row]
_coord2x = coord2x[row]
_coord2y = coord2y[row]

# vectors containing the physical coordinates (in h.o. units, and centered in the lattice)
x_vec = (np.arange(dim_x) - dim_x*0.5) * delta_x
y_vec = x_vec

# draw a mask around the phase pattern, with a radius slightly larget than the TF radius
# set the first pixel to -pi, and the second to +pi
# (so the range of plotting is fixed, and the mask doesn't change color)

approx_cloud_radius = max(TF_radius_, 1)
disk_mask = np.zeros((dim_x, dim_x))
for ix in range(0, dim_x):
for iy in range(0, dim_x):
if ((ix-dim_x*0.5)**2 + (iy-dim_x*0.5)**2) > (1.25 * (2*approx_cloud_radius/length_x) * dim_x / 2)**2:
disk_mask[ix, iy] = 0.99
disk_mask = np.ma.masked_where(disk_mask < 0.5, disk_mask)

phase_A[0, 0] = -np.pi
phase_B[0, 0] = -np.pi
phase_A[0, 1] = np.pi
phase_B[0, 1] = np.pi

################ Draw theoretical trajectories ################
vTraj_x1 /= rescale_
vTraj_y1 /= rescale_
vTraj_x2 /= rescale_
vTraj_y2 /= rescale_

for idx in range(len(vTraj_x1)):
vTraj_x1[idx] = (vTraj_x1[idx] * 2*TF_radius/length_x + 1.)*dim_x*0.5
vTraj_y1[idx] = (vTraj_y1[idx] * 2*TF_radius/length_x + 1.)*dim_x*0.5

for idx in range(len(vTraj_x2)):
vTraj_x2[idx] = (vTraj_x2[idx] * 2*TF_radius/length_x + 1.)*dim_x*0.5
vTraj_y2[idx] = (vTraj_y2[idx] * 2*TF_radius/length_x + 1.)*dim_x*0.5
# drop nan values
vTraj_x1 = vTraj_x1[~np.isnan(vTraj_x1)]
vTraj_y1 = vTraj_y1[~np.isnan(vTraj_y1)]
vTraj_x2 = vTraj_x2[~np.isnan(vTraj_x2)]
vTraj_y2 = vTraj_y2[~np.isnan(vTraj_y2)]

if cycle1 == 0:
c_ind1 = 1-np.linspace(0., .4, len(vTraj_y1))
else:
c_ind1 = np.linspace(0.6, 1, len(vTraj_y1)//2)
c_ind1 = 1-np.concatenate([c_ind1, np.linspace(0., .4, len(vTraj_y1)-len(vTraj_y1)//2)], axis=0)
points1 = np.array([vTraj_x1,vTraj_y1]).T.reshape(-1,1,2)
segments1 = np.concatenate([points1[:-1], points1[1:]], axis=1)

lc1 = LineCollection(segments1, cmap=plt.get_cmap('Greys'), norm=plt.Normalize(0,1))
lc1.set_array(c_ind1)
lc1.set_linewidth(5)
lc11 = LineCollection(segments1, cmap=plt.get_cmap('Greys'), norm=plt.Normalize(0,1))
lc11.set_array(c_ind1)
lc11.set_linewidth(5)
lc111 = LineCollection(segments1, cmap=plt.get_cmap('Greys'), norm=plt.Normalize(0,1))
lc111.set_array(c_ind1)
lc111.set_linewidth(5)

c_ind2 = 1-np.linspace(0., 1, len(vTraj_y2))
points2 = np.array([vTraj_x2,vTraj_y2]).T.reshape(-1,1,2)
segments2 = np.concatenate([points2[:-1], points2[1:]], axis=1)

lc2 = LineCollection(segments2, cmap=plt.get_cmap('Greys'), norm=plt.Normalize(0.,1.3))
lc2.set_array(c_ind2)
lc2.set_linewidth(7)
lc21 = LineCollection(segments2, cmap=plt.get_cmap('Greys'), norm=plt.Normalize(0,1.3))
lc21.set_array(c_ind2)
lc21.set_linewidth(7)

#################################################################
# Make plot
###################

data=[density_A, density_B, phase_A, np.remainder(phase_A-phase_B+np.pi, 2*np.pi)]

x_centers = [_coord1x, _coord2x, _coord1x, _coord2x]
y_centers = [_coord1y, _coord2y, _coord1y, _coord2y]
cmaps = ['afmhot', 'afmhot', 'hsv', 'hsv']

norm = matplotlib.colors.Normalize(vmin = 0., vmax = 1.)

for i in np.arange(len(data)):
ax = plt.subplot(5, 4, i+1+row*4)
if i == 0:
marker = 'o'
subplot_title = "$n_1$"
elif i == 1:
marker = 'D'
subplot_title = "$n_2$"
elif i == 2:
marker = 'o'
subplot_title = "$S_1$"
else:
marker = 'D'
marker2 = 'o'
subplot_title = "$S_1 - S_2 \, \mathrm{mod} \, 2 \pi$"
if row == 0:
ax.set_title(subplot_title, fontsize=fontsi, va='bottom')

x_ticks_1 = [dim_x//2 - int(TF_radius_/delta_x),
dim_x//2 - int(TF_radius_/(2*delta_x)),
dim_x//2,
dim_x//2 + int(TF_radius_/(2*delta_x)),
dim_x//2 + int(TF_radius_/(delta_x))]
x_ticks_2 = ['', '', '', '', '']
y_ticks_1 = [dim_x//2 - int(TF_radius_/delta_x),
dim_x//2 - int(TF_radius_/(2*delta_x)),
dim_x//2,
dim_x//2 + int(TF_radius_/(2*delta_x)),
dim_x//2 + int(TF_radius_/(delta_x))]
y_ticks_2 = ['$'+str(-1)+'$', '$'+str(-0.5)+'$', '$0$', '$'+str(0.5)+'$', '$'+str(1)+'$']
if row == 4:
plt.xlabel('$x/R$', fontsize=fontsi)
plt.xticks(y_ticks_1, y_ticks_2, fontsize=fontsi)
if i==0:
plt.ylabel('$y/R$', fontsize=fontsi)
plt.yticks(y_ticks_1, y_ticks_2, fontsize=fontsi)
else:
plt.yticks(x_ticks_1, x_ticks_2, fontsize=fontsi)
else:
plt.xticks(x_ticks_1, x_ticks_2, fontsize=fontsi)
if i==0:
plt.ylabel('$y/R$', fontsize=fontsi)
plt.yticks(y_ticks_1, y_ticks_2, fontsize=fontsi)
else:
plt.yticks(x_ticks_1, x_ticks_2, fontsize=fontsi)
##############

plt.pcolormesh(data[i], cmap=cmaps[i], linewidth=0, rasterized=True)

# Plot theoretical trajectories
if i==0:
ax.add_collection(lc1)
if i==2:
ax.add_collection(lc11)

if i==1:
plt.scatter(2,2, s=1)
if row != 0:
ax.add_collection(lc2)
if i==3:
ax.add_collection(lc111)
if row != 0:
ax.add_collection(lc21)

#plot mask
if i >= 2:
plt.pcolor(disk_mask, cmap='Greys_r', linewidth=0, rasterized=True)
#plot trajectories
color_indices = np.linspace(0, 1, _coord1x.size) * _coord1x.size / 40.
scatter_plot = plt.scatter(x_centers[i] / delta_x + 0.5*dim_x, y_centers[i] / delta_x + 0.5*dim_x,\
c=color_indices, cmap='Greys_r', norm=norm, s=markersize, edgecolors='Green', marker=marker)
if i == 3:
plt.scatter(x_centers[2] / delta_x + 0.5*dim_x, y_centers[2] / delta_x + 0.5*dim_x,\
c=color_indices, cmap='Greys_r', norm=norm, s=markersize, edgecolors='Green', marker=marker2)

cbaxes = fig.add_axes([0.93, 0.135, 0.015, 0.74])
cbar = plt.colorbar(scatter_plot, cax = cbaxes, ticks = [0, 0.5, 1], format='$%1.1f$')
cbar.ax.tick_params(labelsize=fontsi)
cbar.set_label(r'$t/T$', labelpad=5, y=0.5, fontsize=fontsi)
cbar.ax.invert_yaxis()

plt.tight_layout(rect=(0, 0, 0.92, 0.995))
plt.show()
/export/anaconda3/lib/python3.5/site-packages/matplotlib/figure.py:1742: UserWarning: This figure includes Axes that are not compatible with tight_layout, so its results might be incorrect.
warnings.warn("This figure includes Axes that are not "

Top row: Initial condition, with a vortex at the center of the first component.

Second row: The system imaged after one quarter of a Rabi period ($t=T/4$). The vortex core follows the trajectory marked by the small circles and squares: from black ($t=0$) to gray ($t=T/4$), it travels until the edge of the first component; just before $t=T/4$, while the first vortex is still inside the first cloud, a second vortex starts entering the second cloud, and a domain wall in the relative phase is clearly visible in the fourth column.

Third row: After half a Rabi period ($t=T/2$), the first vortex has completely disappeared from the first component, while the new one gradually migrates to the center of the second component.

Fourth row: At $t=3T/4$, two vortices are again visible, one in each component.

Bottom row: After a full Rabi period ($t=T$), the vortex leaves the second component, reappears inside the first, and returns back to its center. Continuous lines show the theoretical trajectories for $\kappa = 0.485 \omega_\perp$.

In the following, we compare the time dependence of the vortex energy given by the simulated and predicted trajectories. The simulation takes about five minutes to finish.

In [15]:
densities1, densities2 = [], []
phases1, phases2 = [], []
coord1x, coord1y = [], []
coord2x, coord2y = [], []

# linear size of the lattice
dim = 150
length = 15.

# number of Rabi periods to be simulated
Rabi_periods_to_be_simulated = 5

# start evolution
gAB_over_g =1.
Rabi_period = 5*np.pi
coupling = 20.
time0 = time.time()
coord1, coord2, time_coord, vTrajTot_x1, vTrajTot_y1,\
vTrajTot_x2, vTrajTot_y2, time1_th, time2_th = single_vortex(dim, length, Rabi_period, coupling,\
Rabi_periods_to_be_simulated, gAB_over_g,\
densities1, densities2, phases1, phases2,\
coord1x, coord1y, coord2x, coord2y)
print("Total time:", time.time()-time0)
Total time: 251.6512804031372

This is Fig. 7 in the manuscript:

In [16]:
#sns.reset_origin()

R = 2.24638933492        # T-F radius for g=20
xi = 0.423384837924      # and coherence length
k = 0.485
Omega  = 2.*np.pi / Rabi_period

c = 4./3. - 128./(45.*np.pi)
T = 2*np.pi/Omega

u = np.sqrt(vTrajTot_x1**2 + vTrajTot_y1**2)
phi = np.angle(vTrajTot_x1 + 1j * vTrajTot_y1)
t = time1_th
ThEnergy_1 = (1-u**2)*(2.*np.log(R/xi) + np.log(1-u**2)) + 2.*u