The Problem

import mmf_setup;mmf_setup.nbinit()

This cell adds /home/docs/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src to your path, and contains some definitions for equations and some CSS for styling the notebook. If things look a bit strange, please try the following:

  • Choose "Trust Notebook" from the "File" menu.
  • Re-execute this cell.
  • Reload the notebook.

Here we derive some numerical techniques for imprinting a vortex dipole in a periodic box after Billam:2014.

The Problem#

In 2D infinite matter, a vortex at position \(z_0 = x_0 + \I y_0\) can be described, up to an overall phase, by

\[ \psi(z) = f(\abs{z-z_0}) (z-z_0) \propto e^{\I\phi(z)} \]

where \(f(r) \sim r^2\) for small \(r \ll \xi\) describes the core structure. For our purposes, we will neglect this dependence and focus on the form of the phase \(\phi(z)\):

\[ \phi(z) = \tan^{-1}\frac{y-y_0}{x-x_0}. \]

This structure is generally incompatible with periodic boundary conditions due to the global rotation (but see…).

If we consider istead a vortex dipole with an additional anti-vortex at \(z_1 = x_1 + \I y_1\), then we can produce a periodic lattice by considering all images:

\[ \psi_{D}(z) \propto (z-z_0)(z-z_1)^{*}, \qquad \psi(z) \propto \prod_{n_x, n_y} \psi_{D}\Bigl(z-L(n_x,n_y)\Bigr), \qquad L(n_x, n_y) = n_xL_x + \I n_yL_y. \]
%matplotlib inline
import numpy as np, matplotlib.pyplot as plt
from tqdm import tqdm
from itertools import product
Lx = 2*np.pi
Ly = 2*np.pi
d = 2.0
z0 = Lx/4 + 0.25j*Ly
z1 = Lx/4+d + (0.25*Ly+d)*1j

def psi_D(z):
    psi = (z-z0)*(z-z1).conjugate()
    return psi/abs(psi)

def get_psi(z, Nx=12, Ny=12):
    psi = np.prod([
        psi_D(z-(nx*Lx+1j*ny*Ly)) 
        for nx, ny in (product(
            range(-Nx//2, Nx//2+1),
            range(-Ny//2, Ny//2+1)))], axis=0)
    return psi/abs(psi)

N = 100
abs(get_psi(np.pi + 2j*np.pi, Nx=N) - get_psi(np.pi + 0j*np.pi, Nx=N))

np.float64(1.5893967749154707)

x = np.linspace(0, Lx, 100)
y = np.linspace(0, Ly, 102)
z = x[:, None] + 1j*y[None, :]
psi = get_psi(z, Nx=N, Ny=N)
plt.pcolormesh(x, y, np.angle(psi).T)

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[4], line 4
      1 x = np.linspace(0, Lx, 100)
      2 y = np.linspace(0, Ly, 102)
      3 z = x[:, None] + 1j*y[None, :]
----> 4 psi = get_psi(z, Nx=N, Ny=N)
      5 plt.pcolormesh(x, y, np.angle(psi).T)

Cell In[2], line 16, in get_psi(z, Nx, Ny)
     15 def get_psi(z, Nx=12, Ny=12):
---> 16     psi = np.prod([
     17         psi_D(z-(nx*Lx+1j*ny*Ly))
     18         for nx, ny in (product(
     19             range(-Nx//2, Nx//2+1),

Cell In[2], line 12, in psi_D(z)
     11 def psi_D(z):
---> 12     psi = (z-z0)*(z-z1).conjugate()
     13     return psi/abs(psi)

KeyboardInterrupt: 

Ns = 2**np.arange(2, 8)
psis = np.array([get_psi(z, Nx=N, Ny=N) for N in Ns])
plt.loglog(1/Ns, abs(psis - psis[-1]).max(axis=(-1)).max(axis=-1), '+')
\[ \psi = \psi_0 + a\frac{1}{N^2} \]
shape = (psis.shape[0], )
N = 5
psi0 = np.polyfit(1/Ns[:N]**2, psis[:N].reshape((N, np.prod(psis.shape[1:]))), 4)[-1].reshape(psis.shape[1:])
plt.plot(np.log(Ns), np.log(abs(psis-psi0).max(axis=(-1)).max(axis=-1)), '+')
plt.plot(np.log(Ns), -2.3-2*np.log(Ns))

Ns = 2**np.arange(2, 9)
psis = np.array([get_psi(z, Nx=N, Ny=N) for N in Ns])
plt.loglog(1/Ns, abs(psis - psi0).max(axis=(-1)).max(axis=-1), '+')

def A(x):
    return np.arctan(x)

_x = 0.1
_a = 0.2
sum(A(_x/((2*k-1)*np.pi/2 - _a)) - A(_x/((2*k-1)*np.pi/2 + _a)) for k in range(1, 100)), A(np.tanh(_x)*np.tan(_a))

x, y = z.real, z.imag
x0, y0, x1, y1 = z0.real, z0.imag, z1.real, z1.imag
X0, Y0, X1, Y1 = x-x0, y-y0, x-x1, y-y1


Eq5 = sum(A(Y1/((2*m - 1)*np.pi + (X1 - np.pi))) 
        - A(Y1/((2*m - 1)*np.pi - (X1 - np.pi)))
          for m in range(1, 1000))
Eq6 = -A(np.tanh(Y1/2) * np.tan((X1 - np.pi)/2))
plt.imshow(Eq6-Eq5)
plt.colorbar()

# Formula from [Billam:2014]
def get_theta(z, z0, z1, N=5):
    # X^+ = X0, X^- = X1
    x, y = z.real, z.imag
    x0, y0, x1, y1 = z0.real, z0.imag, z1.real, z1.imag    
    X0, Y0, X1, Y1 = x-x0, y-y0, x-x1, y-y1
    res = (
        sum(np.arctan(np.tanh((Y1 + 2*np.pi*n)/2) * np.tan((X1 - np.pi)/2))
            - np.arctan(np.tanh((Y0 + 2*np.pi*n)/2) * np.tan((X0 - np.pi)/2))
            + np.pi * (1.0*(X0>0) - (X1>0))
            for n in range(-N, N+1))
        - (x0 - x1)*y/2/np.pi)
    return res

th = get_theta(z, z0=z0, z1=z1)
inds = (50,50)
fig, axs = plt.subplots(1,2, figsize=(10,5))
axs[0].imshow((th - th[inds] + np.pi) % (2*np.pi) - np.pi)
N = 2
psi0 = get_psi(z, Nx=N, Ny=N)
axs[1].imshow(np.angle(psi0/psi0[inds]))

z[inds]

N = 20
inds = (50,50)
psi0 = get_psi(z, Nx=N, Ny=N)
plt.plot((f[:, 50]/f[inds]).real)
plt.plot((psi0[:, 50]/psi0[inds]).real)

inds = (25,25)
th = get_theta(z, z0=z0, z1=z1)
f = np.exp(1j*th)
plt.imshow(abs(f/f[inds] -  psi0/psi0[inds]))
plt.colorbar()

1.*(abs(z) < 1) - (abs(z) < 2)