---
execution:
  timeout: 30
jupytext:
  formats: md:myst,ipynb,py
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.16.6
kernelspec:
  display_name: Python 3 (ipykernel)
  language: python
  name: python3
---

```{code-cell} ipython3
:init_cell: true

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

# SOC Supersolids

+++

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

$$
  \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,\\
$$

*(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

```{code-cell} ipython3
%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))
```

## 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:

$$
  \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.
$$

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$.


[Martone:2015]: https://doi.org/10.1140/epjst/e2015-02386-x 'G. I. Martone, "Visibility and stability of superstripes in a spin-orbit-coupled Bose-Einstein condensate", Eur. Phys. J. Spec. Top. 224(3), 553--563 (2015)'

+++

### Demler

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

$$
  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}}.
$$

We thus have:

$$
  \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
$$


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$.

```{code-cell} ipython3
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:

$$
  \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},
$$

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

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

The effective Hamiltonian for $\Psi_R$ is:

$$
  \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}.
$$

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

```{code-cell} ipython3
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
```

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


$$
  \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}\\
$$


[Hamner:2013]: http://dx.doi.org/10.1103/PhysRevLett.111.264101 'C. Hamner, Yongping Zhang, J. J. Chang, Chuanwei Zhang, and P. Engels, "Phase Winding a Two-Component BEC in an Elongated Trap: Experimental Observation of Moving Magnetic Orders and Dark-Bright Solitons", prl 111(26), 264101 (2013) [1306.6102](http://arXiv.org/abs/1306.6102)'

+++

## References
* [Approach for making visible and stable stripes in a spin-orbit-coupled Bose-Einstein superfluid](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.90.041604): 

* [Martone:2015]: Characterization of the phase diagram.

* [Kasper:2020]: [Simulating a quantum commensurate-incommensurate phase transition using two Raman coupled one dimensional condensates](https://arxiv.org/abs/2002.12311): 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](https://journals.aps.org/prl/pdf/10.1103/PhysRevLett.124.045702): Different context, but similar supersolid physics.
* [Hamner:2013]: Peter's phase-winding experiment.

[Hamner:2013]: http://dx.doi.org/10.1103/PhysRevLett.111.264101 'C. Hamner, Yongping Zhang, J. J. Chang, Chuanwei Zhang, and P. Engels, "Phase Winding a Two-Component BEC in an Elongated Trap: Experimental Observation of Moving Magnetic Orders and Dark-Bright Solitons", prl 111(26), 264101 (2013) [1306.6102](http://arXiv.org/abs/1306.6102)'

[Martone:2015]: https://doi.org/10.1140/epjst/e2015-02386-x 'G. I. Martone, "Visibility and stability of superstripes in a spin-orbit-coupled Bose-Einstein condensate", Eur. Phys. J. Spec. Top. 224(3), 553--563 (2015)'

+++

# 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.

```{code-cell} ipython3
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)
```

```{code-cell} ipython3
E = e.get_dispersion()
ks = np.linspace(-3,3,100)
plt.plot(ks, E(ks)[0])
E.get_k0()
```

```{code-cell} ipython3
s = e.get_initial_state()
clear_output()
s.plot(show_momenta=True);
history = [s]
dhistory = [s - s]
```

```{code-cell} ipython3
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')
```

```{code-cell} ipython3
import mmfutils.plot.colors
cm = mmfutils.plot.colors.Colormaps.diverging
```

```{code-cell} ipython3
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

```{code-cell} ipython3
str(b"dijf".decode()
```

```{code-cell} ipython3

```