SOC Supersolids

$\newcommand{\vect}[1]{\mathbf{#1}} \newcommand{\uvect}[1]{\hat{#1}} \newcommand{\abs}[1]{\lvert#1\rvert} \newcommand{\norm}[1]{\lVert#1\rVert} \newcommand{\I}{\mathrm{i}} \newcommand{\ket}[1]{\left|#1\right\rangle} \newcommand{\bra}[1]{\left\langle#1\right|} \newcommand{\braket}[1]{\langle#1\rangle} \newcommand{\Braket}[1]{\left\langle#1\right\rangle} \newcommand{\op}[1]{\mathbf{#1}} \newcommand{\mat}[1]{\mathbf{#1}} \newcommand{\d}{\mathrm{d}} \newcommand{\D}[1]{\mathcal{D}[#1]\;} \newcommand{\pdiff}[3][]{\frac{\partial^{#1} #2}{\partial {#3}^{#1}}} \newcommand{\diff}[3][]{\frac{\d^{#1} #2}{\d {#3}^{#1}}} \newcommand{\ddiff}[3][]{\frac{\delta^{#1} #2}{\delta {#3}^{#1}}} \newcommand{\floor}[1]{\left\lfloor#1\right\rfloor} \newcommand{\ceil}[1]{\left\lceil#1\right\rceil} % The following allows KaTeX to render these for significantly improved % performance on CoCalc. \newcommand{\DeclareMathOperator}[2]{\newcommand{#1}{#2}} \DeclareMathOperator{\Tr}{Tr} \DeclareMathOperator{\erf}{erf} \DeclareMathOperator{\erfi}{erfi} \DeclareMathOperator{\sech}{sech} \DeclareMathOperator{\sn}{sn} \DeclareMathOperator{\cn}{cn} \DeclareMathOperator{\dn}{dn} \DeclareMathOperator{\sgn}{sgn} \DeclareMathOperator{\order}{O} \DeclareMathOperator{\diag}{diag} % \newcommand{\mylabel}[1]{\label{#1}\tag{#1}} \newcommand{\degree}{\circ}$

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.

[I 03:35:37 numexpr.utils] NumExpr defaulting to 2 threads.
[I 03:35:38 root] Patching zope.interface.document.asReStructuredText to format code

SOC Supersolids#

Here we consider a two-component BEC with SOC of equal Rashba and Dresselhaus couplings.

\[\begin{split} \op{H} = \begin{pmatrix} K(p) - \mu + \frac{\delta}{2} & \frac{\Omega}{2} e^{2\I k_Rx}\\ \frac{\Omega}{2} e^{-2\I k_Rx} & K(p) - \mu - \frac{\delta}{2}\\ \end{pmatrix}, \qquad \mat{R} = \begin{pmatrix} e^{\I k_Rx}\\ & e^{-\I k_Rx} \end{pmatrix},\\ \mat{R}^\dagger\cdot\op{H}\cdot\mat{R} = \mat{H}_R = \begin{pmatrix} K(p+k_R) - \mu + \frac{\delta}{2} & \frac{\Omega}{2}\\ \frac{\Omega}{2}& K(p-k_R) - \mu - \frac{\delta}{2}\\ \end{pmatrix}\\ \I\hbar \dot{\Psi} = \op{H}\Psi, \qquad \Psi = \mat{R}\Psi_R,\\ \I\hbar \dot{\Psi}_R = (\op{H}_R - \I\hbar\mat{R}^\dagger\dot{\mat{R}})\Psi_R,\\ \end{split}\]

(In this case, $\dot{\mat{R}} = 0$, but we will consider a more general case below.)

Natural scales are expressed in terms $E_R$ and $k_R$ with the dimensionless parameters $k$, $d$, and $w$:

\[ 4E_R = 4\frac{\hbar^2k_R^2}{2m}, \qquad k = \frac{p}{\hbar k_R}, \qquad d = \frac{\delta}{4E_R}, \qquad w = \frac{\Omega}{4E_R}. \]

One can realize this by choosing units so that $\hbar = m/2 = k_R = 1$. In terms of these, the single-particle dispersion bands are:

\[ \frac{E_{\pm}(p)}{2E_R} = \frac{k^2 + 1}{2} \pm \sqrt{(k - d)^2 + w^2}. \]

In addition to these single-particle properties, one has non-linear couplings with an interaction energy of the form:

\[ E_{\mathrm{int}} = g_{aa} \frac{n_a^2}{2} + g_{bb} \frac{n_a^2}{2} + g_{ab}n_a n_b. \]

The relevant dimensionless quantities are:

\[ e_1 = \frac{\bar{n}}{4E_R}\frac{g_{aa} + g_{bb} + 2g_{ab}}{8}, \qquad e_2 = \frac{\bar{n}}{4E_R}\frac{g_{aa} + g_{bb} - 2g_{ab}}{8}, \qquad e_3 = \frac{\bar{n}}{4E_R}\frac{g_{bb} - g_{aa}}{4}. \]

Homogeneous States#

Excluding the super

%pylab inline --no-import-all
k = np.linspace(-1.1,1.1,100)
def E_m(k, d=0, w=1.0):
    return (k**2+1)/2 - np.sqrt((k-d)**2+w**2)

plt.plot(k, E_m(k, d=0.01, w=0.9))

%pylab is deprecated, use %matplotlib inline and import the required libraries.
Populating the interactive namespace from numpy and matplotlib

[<matplotlib.lines.Line2D at 0x7a61b5f50c20>]

../_images/5b00d1bd69e28c875d84f88ff36cc8b309f4a3015ea4d4a5d89d358032f301fc.png

Supersolid Stripe Phase#

As discussed in Martone:2015, a supersolid stripe phase can be realized if $\delta = 0$ and $w$ is sufficiently small if the superfluid is miscible $g=\sqrt{g_{aa}g_{bb}} > g_{ab}$. (The miscible $^{87}$Rb hyperfine states $a=\ket{1, -1}$ and $b=\ket{1,0}$ states have $a_{aa} = 100.40a_B$, $a_{bb}=100.86a_B$, and $a_{ab}= 100.41a_B$. The relevant quantity is the geometric mean $a = 100.63a_B$, so the miscibility criterion is satisfied.)

They use the following Anzatz:

\[\begin{split} \Psi = \sqrt{\bar{n}}\left[ C_+ \begin{pmatrix} \cos\theta_+\\ -\sin\theta_+ \end{pmatrix} e^{\I k_+ x} + C_- \begin{pmatrix} \sin\theta_-\\ -\cos\theta_- \end{pmatrix} e^{-\I k_- x} \right], \qquad \abs{C_+}^2+\abs{C_-}^2 = 1. \end{split}\]

The supersolid phase occurs when $\abs{C_\pm} = 1/\sqrt{2}$ with $k_{\pm} = \pm k_R$ with density

\[ \frac{n(x)}{\bar{n}} = 1 + \frac{\Omega}{2(2E_R + G_1)}\cos(2k_R x + \phi) = 1 + \frac{w}{1 + 2e_1}\cos(2k_R x + \phi) \]

(Note: in their notations $e_i = G_i/4E_R$.) .

In the limits of small $w$ and small $d$, the conditions for a stripe phase is:

\[ \lim_{w\rightarrow 0}: \quad \abs{d + 2e_3} \leq 4e_2, \qquad \lim_{d \rightarrow 0}: \quad w < \sqrt{\frac{1}{1+\frac{e_1}{e_2}}} \approx 0.033. \]

The actual region is approximately triangular with these vertices. The bound on $w$ is for the $a=\ket{1, -1}$ and $b=\ket{1,0}$ states of $^{87}$Rb, which also sets an upper bound on the maximum relative density contrast is $\delta n / \bar{n} < w/(1 + 2e_1) < w$.

Demler#

To compare with [Kasper:2020], we note the following correspondence:

\[\begin{split} g \equiv g_{aa} = g_{bb}, \qquad 0 \equiv g_{ab},\qquad 0 \equiv \delta, \qquad \frac{Q}{2} \equiv k_R, \qquad M \equiv m, \\ 2\hbar J \equiv \Omega, \qquad x_0 = \sqrt{\frac{\hbar}{2MJ}} \equiv \frac{\hbar}{\sqrt{m\Omega}},\qquad E_0 = \frac{\hbar^2}{2Mx_0^2} \equiv \frac{\Omega}{2},\\ \bar{Q} = x_0Q \equiv \sqrt{\frac{8 E_R}{\Omega}} = \sqrt{\frac{2}{w}}, \qquad \bar{\mu} = \frac{\mu}{E_0} \equiv \frac{2 \mu}{\Omega}, \qquad \bar{g} = \frac{g}{x_0E_0} \equiv \frac{2g\sqrt{m}}{\hbar\sqrt{\Omega}}. \end{split}\]

We thus have:

\[\begin{split} \mu = g \bar{n} = \frac{\Omega}{2}\bar{\mu}\\ e_1 = e_2 = \frac{w}{8}\bar{\mu}, \qquad e_3 = 0, \qquad w = \frac{2}{\bar{Q}^2}, \qquad \end{split}\]

In Figure 2. they have $\bar{\mu}=30$ and $\bar{g}=0.6$. Since $d=e_3=0$, and $e_1=e_2$, there should be a SS phase for $w < \sqrt{1/2} = 0.707$ corresponding to $\bar{Q} > \sqrt[4]{8} = 1.68$. This is consistent with their homogeneous (commensurate) phase for $\bar{Q}=0.7$ and a inhomogeneous (incommensurate) SS-phase for $\bar{Q}=1.9$.

np.random.seed(2)
g, m, k_R, Omega, hbar, n = np.random.random(6)
E_R = hbar**2*k_R**2/2/m
w = Omega/4/E_R
mu = g*n

Q = 2*k_R
J = Omega/2/hbar
x0 = np.sqrt(hbar/2/m/J)
E0 = hbar**2/x0**2/2/m
Qbar = x0*Q
mubar = mu/E0
gbar = g/x0/E0

assert np.allclose(E0, Omega/2)
assert np.allclose(Qbar, np.sqrt(2/w))
assert np.allclose(mubar, 2*mu/Omega)
assert np.allclose(gbar, 2*g*np.sqrt(m/Omega)/hbar)

Hamner Phase Winding#

In Hamner:2013, Peter looks at phase winding from a coupled system with the following Hamiltonian:

\[\begin{split} \op{H} = \begin{pmatrix} K(\op{p}) + \Delta_0 + \delta x& \Omega_0\\ \Omega_0 & K(\op{p}), \end{pmatrix}, \qquad K(\op{p}) = \frac{\op{p}^2}{2m}, \end{split}\]

with the usual non-linear couplings and $\Omega_0/2\pi \hbar = 7.4$kHz. Here they couple the $\ket{1,-1}$ and $\ket{2,0}$ states with $a_{aa} = 100.40a_B$, $a_{bb}=94.57a_B$, and $a_{ab}= 98.13a_B$. These states are weakly immiscible, so supersolidity is not relevant, but perhaps there is an analogy?

Consider now the transformation

\[\begin{split} \mat{R}_{a} = \begin{pmatrix} e^{-\I \delta t x /\hbar}\\ & 1 \end{pmatrix}, \qquad \Psi = \mat{R}_a\Psi_r. \end{split}\]

The effective Hamiltonian for $\Psi_R$ is:

\[\begin{split} \I\hbar\dot{\Psi}_R = \mat{H}_R\dot{\Psi}_R, \qquad \mat{H}_R = \mat{R}_a^\dagger\cdot\op{H}\cdot\mat{R}_a - \I\hbar \mat{R}_a^\dagger\dot{\mat{R}}_a = \begin{pmatrix} K(\op{p}_x - \delta t) + \Delta_0 & \Omega_0 e^{\I \delta t x/\hbar}\\ \Omega_0 e^{-\I \delta t x/\hbar} & K(\op{p}_x) \end{pmatrix}. \end{split}\]

This looks somewhat like a SOC Hamiltonian with a time-dependent $k_R = \delta t x/2\hbar$.

ga, gb, gab = (100.4, 94.57, 98.13)
#ga, gb, gab = (100.4, 100.86, 100.41)
print(np.sqrt(1/(1+(ga+gb + 2*gab)/(ga+gb - 2*gab))))
from gpe.bec import u
delta = u.hbar * 2*np.pi *2.3*u.kHz
wc = 0.033
Omega0 = u.hbar * 2*np.pi * u.kHz
t = 2*np.sqrt(u.m * Omega0/wc)/delta
t/u.ms

nan

/tmp/ipykernel_4142/2068196133.py:3: RuntimeWarning: invalid value encountered in sqrt
  print(np.sqrt(1/(1+(ga+gb + 2*gab)/(ga+gb - 2*gab))))

np.float64(2.2339557556563774)

If we identify this with an SOC system, then we have a time-varying $k_R = \delta t/2\hbar$, hence we might obtain a transition to the supersolid phase when:

\[ w = \frac{m\Omega_0}{\hbar^2 k_R^2} < 0.033,\qquad t > \frac{2}{\delta}\sqrt{\frac{m\Omega_0}{0.033}}. \]

We consider this system now

\[\begin{split} \op{H} = \begin{pmatrix} K(p-\kappa t) & \frac{\Omega}{2} e^{2\I \kappa t x}\\ \frac{\Omega}{2} e^{-2\I \kappa t x} & K(p+\kappa t)\\ \end{pmatrix}, \qquad \mat{R} = \begin{pmatrix} e^{\I \kappa t x}\\ & e^{-\I \kappa t x} \end{pmatrix},\qquad -\I\hbar\mat{R}^\dagger\dot{\mat{R}} = \begin{pmatrix} \hbar\kappa x \\ & -\hbar\kappa x \end{pmatrix}\\ (\op{H}_R - \I\hbar\mat{R}^\dagger\dot{\mat{R}}) = \begin{pmatrix} K(p) + \hbar \kappa x & \frac{\Omega}{2}\\ \frac{\Omega}{2}& K(p) - \hbar \kappa x\\ \end{pmatrix}\\ \end{split}\]

References#

Approach for making visible and stable stripes in a spin-orbit-coupled Bose-Einstein superfluid:
Martone:2015: Characterization of the phase diagram.
[Kasper:2020]: Simulating a quantum commensurate-incommensurate phase transition using two Raman coupled one dimensional condensates: Paper by Demler’s group that produces a similar Hamiltonian with a time-varying field coupling two tubes that can tunnel. A SOC-like Hamiltonian emerges (they call it a Pokrovsky-Talapov (PT) model). They do not mention the SOC case.
Rotating a Supersolid Dipolar Gas: Different context, but similar supersolid physics.
Hamner:2013: Peter’s phase-winding experiment.

Trap#

Here we use our code to construct the supersolid strip phase in a small external trap. We use this trap to break the degeneracy of the ground state to help our minimizer find the solution.

import mmf_setup;mmf_setup.nbinit()
from gpe.imports import *

import soc_supersolid;reload(soc_supersolid)
from soc_supersolid import u
e = soc_supersolid.Supersolid(
    #lattice_height_V_TF=0.1*0.5, T__x=np.inf,
    lattice_height_V_TF=0.0, T__x=np.inf, 
    #xi_micron=0.125/2,
    xi_micron=0.1,
    #dx=0.01 * u.micron,
    #w=2.7/4.0, d=0.0,
    #w=2.7/4.0, d=0./16,lattice_k_k_r=0.78,
    w=0.1/4.0, d=-0.001223419561528528,
    lattice_k_k_r=1.0,
    #lattice_height_V_TF=0.0, T__x=2.0, 
    barrier_x=0*u.micron,
    cells_x=20)
s = e.get_state()
s = e.get_initial_state()
clear_output(wait=True)

res = e.plot(s)
e1, e2, e3 = res['es']
w, d = res['w'], res['d']
w_c = np.sqrt(1/(1+e1/e2))
print(f"SS window: w={w:.4f} < {w_c:.4f}, d={d:.4f} in {-2*e3:.4f}+-{4*e2:.4f}")
if not(w < w_c):
    print(f"w={w:.4f} > {w_c:.4f} too large for SS")
if not(abs(d + 2*e3) < 4*e2):
    print(f"d = {d:.4f} outside window ({-2*e3-4*e2:.4f}, {-2*e3+4*e2:.4f}) for SS")
n = s.get_density().sum(axis=0)
(n.max() - n.min())/(n.max() + n.min()), w/(1+2*e1)

SS window: w=0.0250 < 0.0331, d=-0.0012 in -0.0012+-0.0009

(np.float64(0.01260954569538818), np.float64(0.01786075625567671))

../_images/2e5ff47d763fce646f0259ef44fcab38fa4a8e214d8d73a2a492f6815bd9e413.png

E = e.get_dispersion()
ks = np.linspace(-3,3,100)
plt.plot(ks, E(ks)[0])
E.get_k0()

np.float64(0.9996882148689227)

../_images/cbe5f42b3dcba173eed6369e1b69dd86ff6d917d9583c1fb30437caab2bffaca.png

s = e.get_initial_state()
clear_output()
s.plot(show_momenta=True);
history = [s]
dhistory = [s - s]

../_images/beed35dff35d2b637d17a442db2bd7b8d588df2432ba8f9f133e6985d5350ba6.png

ev = EvolverABM(history[-1], dt=0.3*s.t_scale)
pe = None
with NoInterrupt() as interrupted:
    while not interrupted:
        for n in range(10):
            ev.evolve(100)
            history.append(ev.get_y())
            dhistory.append(history[-1] - history[0])
        pe = ev.y.plot(plot_elements=None, history=dhistory)
        display(pe.fig)
        clear_output(wait=True)
        plt.close('all')

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[8], line 3
      1 ev = EvolverABM(history[-1], dt=0.3*s.t_scale)
      2 pe = None
----> 3 with NoInterrupt() as interrupted:
      4     while not interrupted:

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pytimeode/evolvers.py:411, in EvolverABM.evolve(self, steps)
    410 for kt in range(steps):
--> 411     self.do_step()
    412 assert np.allclose(self.y.t, t0 + steps * self.dt)

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pytimeode/evolvers.py:427, in EvolverABM.do_step(self)
    425 else:
    426     # self.do_step_ABM()
--> 427     self.do_step_ABM_apply()

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pytimeode/evolvers.py:594, in EvolverABM.do_step_ABM_apply(self)
    593 dy = dys.pop()
--> 594 dy = self.get_dy(y=y, dy=dy)
    596 ys.insert(0, y)

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pytimeode/evolvers.py:106, in EvolverBase.get_dy(self, y, t, dy)
    105 with y.lock:
--> 106     dy = y.compute_dy_dt(dy=dy)
    107     if dy is None:

File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/bec2.py:614, in _StateBase.compute_dy_dt(self, dy, subtract_mu)
    613 Ky = y.copy()
--> 614 Ky.apply_laplacian(factor=self.K_factor)
    615 Va, Vb, Vab = self.get_V_GPU()

File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/bec2.py:710, in _StateBase.apply_laplacian(self, factor, exp, **_kw)
    709         _kw[_k] = _v
--> 710 self.data[...] = self.basis.laplacian(self.data, factor=factor, exp=exp, **_kw)

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/mmfutils/math/bases/bases.py:431, in PeriodicBasis.laplacian(self, y, factor, factors, exp, kx2, k2, kwz2, **_kw)
    430 # Apply K
--> 431 yt = self.fftn(y)
    432 laplacian_y = self.ifftn(
    433     self._apply_K(yt, kx2=kx2, k2=k2, exp=exp, factor=factor, factors=factors)
    434 )

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/mmfutils/math/bases/bases.py:500, in PeriodicBasis.fftn(self, x)
    499 axes = self.axes % len(x.shape)
--> 500 return self._fftn(x, axes=axes)

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/mmfutils/performance/fft.py:458, in fftn(a, s, axes)
    457     _FFT_CACHE[key] = get_fftn(a=a.copy(), s=s, axes=axes)
--> 458 res = _FFT_CACHE[key](a)
    459 if _COPY_OUTPUT and not res.flags["OWNDATA"]:

KeyboardInterrupt: 

During handling of the above exception, another exception occurred:

AttributeError                            Traceback (most recent call last)
Cell In[8], line 3
      1 ev = EvolverABM(history[-1], dt=0.3*s.t_scale)
      2 pe = None
----> 3 with NoInterrupt() as interrupted:
      4     while not interrupted:
      5         for n in range(10):
      6             ev.evolve(100)

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/mmf_contexts/contexts.py:537, in NoInterrupt.__exit__(self, exc_type, exc_value, traceback)
    535 kernel = getattr(IPython.get_ipython(), "kernel", None)
    536 if kernel:
--> 537     del kernel.pre_handler_hook
    538     del kernel.post_handler_hook

AttributeError: 'IPythonKernel' object has no attribute 'pre_handler_hook'

import mmfutils.plot.colors
cm = mmfutils.plot.colors.Colormaps.diverging

s = history[0]
n0 = s.get_density()
xs = s.xyz[0].ravel()
ts = [_s.t/_s.t_unit for _s in history]
plt.subplot(211)
ns = np.array([(_s.get_density()-n0)[0] for _s in history])
imcontourf(ts, xs, ns, cmap=cm, vmax=0.4, vmin=-0.4);plt.xlabel('t [ms]')
plt.colorbar()
plt.subplot(212)
ns = np.array([(_s.get_density()-n0)[1] for _s in history])
imcontourf(ts, xs, ns, cmap=cm, vmax=0.4, vmin=-0.4);plt.xlabel('t [ms]')
plt.colorbar()

Lattice#

str(b"dijf".decode()