import mmf_setup; mmf_setup.nbinit()
import os, sys
from pathlib import Path
FIG_DIR = Path(mmf_setup.ROOT) / '../Docs/_build/figures/'
os.makedirs(FIG_DIR, exist_ok=True)
import logging; logging.getLogger("matplotlib").setLevel(logging.CRITICAL)
%matplotlib
import numpy as np, matplotlib.pyplot as plt
try: from myst_nb import glue
except: glue = None

from matplotlib.animation import FuncAnimation
from gpe.contexts import FPS

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

Using matplotlib backend: module://matplotlib_inline.backend_inline

Counterflow Instabilities

Bogoliubov Modes (Phonons)

Recall that \(\xi_{i} = h_{i}/\sqrt{2}\) is a slightly different healing length than we typically use \(\frac{\hbar^2}{2mh_i^2} = \mu_i = g_{ii}n_i\).

As derived in Counterflow Instability, with the simplification considered in [Alperin and Timmermans, 2025] of equal single-particle speeds of sound \(c_{i}\) and healing lengths \(\xi_i\):

\[\begin{gather*} \omega_{B,i}(k) = \pm c_{i}k_i\sqrt{1 + \xi_{i}^2k^2}, \\ \omega_{B,a}(k) = \omega_{B,b}(k) = \omega_{B}(k), \qquad c_{a} = c_{b} = c, \qquad \xi_{a} = \xi_{b} = \xi, \end{gather*}\]

the Bogoliubov modes have the following dispersion

\[\begin{gather*} \omega^2 = \omega_B^2 + v^2k^2 \pm \sqrt{4\omega_B^2v^2k^2 + \gamma^2c^4k^4}, \qquad \gamma^2 = \frac{g_{ab}^2}{g_{a}g_{b}} \end{gather*}\]

where the \(v = v_{a}-v_{b}\) is the relative velocity.

There are three sources of dynamic instability here:

\[\begin{gather*} 1 + \xi^2 k^2 < 0,\\ 4\omega_B^2v^2k^2 + \gamma^2c^4k^4 < 0,\\ \omega_B^2 + v^2k^2 - \sqrt{4\omega_B^2v^2k^2 + \gamma^2c^4k^4} < 0. \end{gather*}\]

The first corresponds to the usual The first only happens when \(\gamma^2 < 0\) – i.e. if one of the modes

Here we explore the generation of counterflow instabilities in miscible two-component BECs by manipulating the velocities. in the case where \(g_{a}n_{a} = g_{b}n_{b} = gn\) with repulsive interactions, the instability will develop when

\[\begin{gather*} \Abs{gn + \frac{\hbar^2 (k^2 - 4q^2)}{4m}} < \abs{g_{ab}}\sqrt{n_a n_b} \end{gather*}\]

where \(k\) is the wavevector of the perturbation, and \(q\) is the relative wave-vector of the two components. The onset of the instability is

\[\begin{gather*} q_{c} = \frac{\sqrt{m(gn - g_{ab}\sqrt{n_a n_b})}}{\hbar}. \end{gather*}\]

As this formula is only strictly valid if \(g_an_a = g_bn_b = gn\), it is useful to assume that \(n_b = n_ag_a/g_b\) or \(n_a = n_b g_b / g_a\). This gives

Recall that the two components have momentum \(\pm \hbar q\), so the relative velocity has a factor of 2: \(v = 2\hbar q/m\).

\[\begin{gather*} \frac{2\hbar q_c}{m} = v_{c} = \underbrace{\sqrt{\frac{gn}{m}}}_{c_{\text{sound}}}\; 2\underbrace{\sqrt{1 - \frac{g_{ab}}{\sqrt{g_ag_b}}}}_{\sqrt{1-\gamma}}. \end{gather*}\]

Thus, the factor on the right \(2\sqrt{1-\gamma}\) expresses the fraction below the Landau critical velocity \(c_{\text{sound}}\) (or speed-of-sound) at which the modulational instability first starts.

Important

Think carefully about this. If \(\gamma = 0\), then the mode is unstable if \(q_c > q_{\text{sound}}\) – i.e. if each fluid exceeds the speed of sound relative to the background frame. However, the relative velocity between the two fluids is \(2c_{\text{sound}}\). This is saying that if \(\gamma < 3/4 = 0.75\), then the fluids can flow past each other at a relative speed faster than the speed of sound, and the flow remains stable. Once \(\gamma > 0.75\), then the instability occurs for subsonic relative flow. We demonstrate this below.

From the imaginary part of the eigenvalues \(\omega\) we obtain the growth rate:

\[\begin{gather*} \delta n(t) \propto e^{t/\tau}, \qquad \abs{\Im \omega} = \frac{1}{\tau}. \end{gather*}\]

Of course, the growth will saturate at some point, thus, to measure this, we must consider only the linear regime.

For our experiment, we use the following parameters:

• healing_length: \(\hbar/\sqrt{2m g n}\) where \(gn = (g_{aa}n_a + g_{bb}n_b)/2\).

• na: \(n_a\).

• nb: \(n_b\).

• dgn: \(gna - gnb\). This could probably be more physical.

• gamma: \(g_{ab}/\sqrt{g_{a}g_{b}}\): Miscibility.

From these, we can define the initial background densities \(n_a\) and \(n_b\), and then the coupling constants.

_TINY = np.finfo(float).tiny
from gpe import bec, bec2
from gpe.utils import ExperimentBase, AsNumpyMixin, StateWithExperimentMixin, _GPU
from gpe.minimize import MinimizeState

u = bec2.u

class State(StateWithExperimentMixin, bec2._StateBase):
    def get_ns_TF(self, Vs_TF):
        return self.experiment.get_ns_TF(state=self, Vs_TF=Vs_TF)
    
    @property
    def t0(self):
        return self.experiment.t_unit


class Experiment(bec2.V0Mixin, bec2.GPEMixin, AsNumpyMixin, ExperimentBase):
    hbar = m = 1.0
    Nx = 64
    healing_length_micron = 1.0  # Defines gn = (g_an_a + g_bn_b)/2
    dx_healing_length = 0.5
    na = nb = 1.0
    dgn = 0.0
    gamma = 0.9
    State = State
    
    # We allow for a small bump for species a to test the Landau criterion.
    V0 = 0.0
    sigma = 1.0
    
    ######################################################################
    # Methods required by IExperiment
    def init(self):
        self.healing_length = self.healing_length_micron * u.micron
        self.dx = self.dx_healing_length * self.healing_length
        self.Lx = self.dx * self.Nx
        na, nb = ns = np.array([self.na, self.nb])
        gn = (self.hbar / self.healing_length)**2 / 2 / self.m
        ga, gb = gn / ns
        gab = self.gamma * np.sqrt(ga*gb)
        self.gs = (ga, gb, gab)
        self.mus = (gn, gn)
        
        # Estimate of the critical relative wavevector for the instability:
        self._qc = np.sqrt(self.m*(gn - gab*np.sqrt(na*nb)))/self.hbar
        
        # This is not obvious: can we somehow automate this?  If not, it should
        # be described in an interface.
        self.shape = (2, self.Nx)
        
        self.state_args = dict(
            Nxyz=(self.Nx,),
            Lxyz=(self.Lx,),
            ms=(self.m, self.m),
            mus=self.mus,
        )
        super().init()

    def get_Vext_GPU(self, state=None):
        """Return (Va, Vb, Vab), the external potentials."""
        if state is None:
            state = self
        zeros = state.xp.zeros(state.shape)
        mu_a = mu_b = 0.0
        if state.initializing and getattr(self, "mus", None):
            mu_a, mu_b = self.mus
        Va, Vb, Vab = zeros - mu_a, zeros - mu_b, zeros
        Va += self.V0 * np.exp(-(state.xyz[0]/self.sigma)**2/2)
        return Va, Vb, Vab
    
    def get_state(self, initialize=True):
        """Quickly return a valid `State` object."""
        return self.State(experiment=self, **self.state_args)

    def get_initial_state(self, _E_tol=1e-12, _psi_tol=1e-12):
        state = self.get_state()
        minimizer = MinimizeState(state, fix_N=True)
        psi0 = minimizer.minimize(E_tol=_E_tol, psi_tol=_psi_tol).get_psi()
        state = self.get_state()
        state.set_psi(psi0)
        return state

    def get_initialized_state(self, state):
        """Return a valid state initialized from `state`.

        This is used in chained simulations where a specified state of one
        simulation is used to initialize a state for further use.  For example,
        for expansion."""
        state_ = self.get_state()
        state_.set_psi(state.get_psi())
        return state_
    # End of methods required by IExperiment
    ######################################################################
   
    ######################################################################
    # Helper methods
    def get_t_growth(self, qn, kn=None):
        """Return the timescale for exponential growth of an unstable mode.
        
        Arguments
        ---------
        qn : int
            Which relative momentum `q = state.basis.kx[qn]`.
        kn : int, None
            Which mode to consider `k = state.basis.kx[kn]`.  If `None`, then return the
            shortest timescale (most-unstable mode).
        """
        kx = self.get_state().basis.kx
        ga, gb, gab = self.gs
        na, nb = self.na, self.nb
        f = self.hbar / np.sqrt(2*self.m)
    
        @np.vectorize
        def get_T(qn, kn):
            q = f * kx[qn]
            k = f * kx[kn]
            M = [[k**2 + 2*k*q + ga*na, gab*nb, -ga*na, -gab*nb],
                 [gab*na, k**2 - 2*k*q + gb*nb, -gab*na, -gb*nb],
                 [ga*na, gab*nb, -(k**2 - 2*k*q + ga*na), -gab*nb],
                 [gab*na, gb*nb, -gab*na, -(k**2 + 2*k*q + gb*nb)]]
            ws, uvs = np.linalg.eig(M)
            with np.errstate(divide='ignore'):
                T = 1 / (abs(ws.imag).max() / self.hbar)
            return T
        
        if kn is None:
            kn = np.arange(self.Nx//2)
            return min(get_T(qn, kn))
        return get_T(qn, kn)
        
    def plot_bogoliubov(self, ax=None):
        """Show the most unstable Bogoliubov modes."""
        if ax is None:
            fig, ax = plt.subplots()
        else:
            fig = ax.figure
        
        qn = np.arange(self.Nx//2)[1:20]
        kn = np.arange(self.Nx//2)[1:40]
        
        ga, gb, gab = self.gs
        na, nb = self.na, self.nb
        Ts = self.get_t_growth(qn[:, None], kn[None, :])
        mesh = ax.pcolormesh(qn, kn, Ts.T/self.t_unit)
        ax.set(xlabel="$n_q$", ylabel="$n_k$")
        fig.colorbar(mesh, label=rf"Growth time $t_{{\text{{growth}}}}$ [{self.t_name}]")
        return ax
            
e = Experiment()
s = e.get_state()
e.plot_bogoliubov();

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

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[2], line 17
     13     def t0(self):
     14         return self.experiment.t_unit
     15 
     16 
---> 17 class Experiment(bec2.V0Mixin, bec2.GPEMixin, AsNumpyMixin, ExperimentBase):
     18     hbar = m = 1.0
     19     Nx = 64
     20     healing_length_micron = 1.0  # Defines gn = (g_an_a + g_bn_b)/2

AttributeError: module 'gpe.bec2' has no attribute 'V0Mixin'

To see the development of the counterflow instability, we give the components a relative velocity above the critical velocity, and add a small amount of noise to seed the instability. From the eigenvalues, we can determine the timescale for the growth of various modes \(t_{\text{growth}}\).

from gpe.contexts import FPS
from pytimeode.evolvers import EvolverABM
rng = np.random.default_rng(seed=2)
qn = 5

e = Experiment(Nx=64)
s = e.get_state()
q = s.basis.kx[qn]
s.set_psi(
    s.get_psi() 
    * np.exp(1j * np.array([q, -q])[:, None] * s.x[None,:])
    * np.exp(1j*(rng.random(size=s.shape)-0.5)*0.001))

# Evolve to t_growth
t_growth = e.get_t_growth(qn=qn)
T = 10*t_growth
dt = 0.4*s.t_scale
Nt = 20
dT = T/Nt
steps = int(np.ceil(dT/dt))
dt = dT/steps

ev = EvolverABM(s, dt=dt)
states = [s]
for n in range(Nt):
    ev.evolve(steps)
    states.append(ev.get_y())

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[3], line 6
      2 from pytimeode.evolvers import EvolverABM
      3 rng = np.random.default_rng(seed=2)
      4 qn = 5
      5 
----> 6 e = Experiment(Nx=64)
      7 s = e.get_state()
      8 q = s.basis.kx[qn]
      9 s.set_psi(

NameError: name 'Experiment' is not defined

fig, ax = plt.subplots()
for s in FPS(states, display=True, fig=fig, embed=True):
    ax.cla()
    s.plot()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[4], line 2
      1 fig, ax = plt.subplots()
----> 2 for s in FPS(states, display=True, fig=fig, embed=True):
      3     ax.cla()
      4     s.plot()

NameError: name 'states' is not defined

fig, ax = plt.subplots()
n0 = np.array([e.na, e.nb])[:, None]
ns = np.array([s.get_density() - n0 for s in states])
n_max = np.max(ns)
for s in FPS(states, display=True, fig=fig, embed=True):
    ax.cla()
    na, nb = s.get_density() - n0
    ax.plot(s.x, na, label='$n_a-n_0$')
    ax.plot(s.x, nb, label='$n_b-n_0$')
    ax.set(ylim=(-n_max, n_max))
    ax.legend()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[5], line 2
      1 fig, ax = plt.subplots()
----> 2 n0 = np.array([e.na, e.nb])[:, None]
      3 ns = np.array([s.get_density() - n0 for s in states])
      4 n_max = np.max(ns)
      5 for s in FPS(states, display=True, fig=fig, embed=True):

NameError: name 'e' is not defined

ts = [s.t for s in states]
fig, ax = plt.subplots()
n_maxs = ns.max(axis=2)
ax.semilogy(ts, n_maxs)
ax.plot(ts, 9e-5*np.exp(ts/t_growth))
ax.set(ylim=(1e-4,1))
inds = slice(1, -1, None)
Ps = np.polyfit(ts[inds], np.log(n_maxs)[inds], deg=1).T
plt.plot(ts, np.exp(np.polyval(Ps[0], ts)), ':C0')
plt.plot(ts, np.exp(np.polyval(Ps[1], ts)), ':C1')
1/Ps[:, 0], t_growth

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[6], line 1
----> 1 ts = [s.t for s in states]
      2 fig, ax = plt.subplots()
      3 n_maxs = ns.max(axis=2)
      4 ax.semilogy(ts, n_maxs)

NameError: name 'states' is not defined

fig, ax = plt.subplots()
nts = np.fft.fftshift(np.fft.fft(ns, axis=-1), axes=[-1])
kn = np.fft.fftshift(s.basis.kx) / (2*np.pi / e.Lx)
assert np.allclose(kn, np.round(kn, 0))
kn = np.round(kn, 0).astype(int)
mesh = ax.pcolormesh(ts/t_growth, kn, np.log10(abs(nts[:, 0, :])).T, vmin=-3)
fig.colorbar(mesh, label=r"$\log_{10}\delta n$")
ax.set(ylabel="Mode $n$ of perturbation $k_n$.",  xlabel="$t/t_{\mathrm{growth}}$");

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[7], line 2
      1 fig, ax = plt.subplots()
----> 2 nts = np.fft.fftshift(np.fft.fft(ns, axis=-1), axes=[-1])
      3 kn = np.fft.fftshift(s.basis.kx) / (2*np.pi / e.Lx)
      4 assert np.allclose(kn, np.round(kn, 0))
      5 kn = np.round(kn, 0).astype(int)

NameError: name 'ns' is not defined

This plot shows the growth of the various modes as a function of time. We see that, up until time \(t \approx 7t_{\text{growth}}\), only the unstable modes grow. After this time, the non-linearities become significant, and amplitude is transferred to other modes (primarily the harmonics).

fig, axs = plt.subplots(10, 1, figsize=(4, 8), sharex=True, gridspec_kw=dict(hspace=0))
i0 = 2
n_max = nts[i0 + len(axs), 0, :].max()
for ax, i in zip(axs, range(i0, len(ts), 2)):
    ax.plot(kn, abs(nts[i, 0, :]))
    ax.plot(kn, 1.1e-3*np.exp(ts[i]/e.get_t_growth(qn, kn=kn)), '--')
    ax.set(xlim=(0, kn.max()))
    ax.text(0.7, 0.7, rf"$t={ts[i]/t_growth:.1f} t_{{\rm growth}}$",
           transform=ax.transAxes)
ax.set(xlabel="Mode index $n$ of $k_n$.");

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[8], line 3
      1 fig, axs = plt.subplots(10, 1, figsize=(4, 8), sharex=True, gridspec_kw=dict(hspace=0))
      2 i0 = 2
----> 3 n_max = nts[i0 + len(axs), 0, :].max()
      4 for ax, i in zip(axs, range(i0, len(ts), 2)):
      5     ax.plot(kn, abs(nts[i, 0, :]))
      6     ax.plot(kn, 1.1e-3*np.exp(ts[i]/e.get_t_growth(qn, kn=kn)), '--')

NameError: name 'nts' is not defined

Growth of modes (solid blue) compared with

\[\begin{gather*} 1.1\times 10^{-3}e^{t/t_{\text{growth}}(k_n)} \end{gather*}\]

where \(t_{\text{growth}} = 1/\Im(\omega_n)\) is the exponential growth rate of the corresponding modes. We see that by \(t/t_{\text{growth}} \approx 8\), the linear growth has saturated, and non-linear interactions start transferring amplitude to the other modes.

Landau Criterion

The main point of [Huynh et al., 2023] is that there are particular parameter regimes where stationary supersonic flow past a barrier exists. See also [Leboeuf and Pavloff, 2001] and [Paris-Mandoki et al., 2017].

Here we explore the Landau criterion using the same simulation. [Huynh et al., 2023] has a thorough analysis of stationary flow past a barrier in the 1D GPE, which clearly demonstrates several points:

1. If the barrier is attractive, then there exists stable flow for \(v < c\) (and only special cases for stationary supersonic flow \(v > c\)).

2. If the barrier is repulsive, then the threshold is strictly less than \(c\). For a wide barrier, this corresponds with the suppression in density at the height of the barrier. This suppression is modified for narrow barriers.

We briefly demonstrate the first point below by turning off the inter-species coupling (\(\gamma = 0\)), and inducing a flow slightly above and below the speed of sound with a small attractive potential.

Show code cell source

Hide code cell source

%%time
rng = np.random.default_rng(seed=2)
e = Experiment(Nx=64, gamma=0.0, V0=-0.01)  # Slight attractive barrier - no g_ab
qs = np.sqrt(e.m*e.gs[0]*e.na)   # q for sound.
s0 = e.get_state()
s1 = e.get_state()
qn = np.where(s.basis.kx > qs)[0][0]
q0, q1 = s.basis.kx[qn-1:qn+1]    # Modes just below and above threshold
for s, q in zip((s0, s1), (q0, q1)):
    s.set_psi(s.get_psi() * np.exp(1j * np.array([q, 0])[:, None] * s.x[None,:]))

t_growth = 10
T = 10*t_growth
dt = 0.4*s.t_scale
Nt = 20
dT = T/Nt
steps = int(np.ceil(dT/dt))
dt = dT/steps

ev0 = EvolverABM(s0, dt=dt)
ev1 = EvolverABM(s1, dt=dt)
evs = [ev0, ev1]

fig, ax = plt.subplots()
for n in FPS(Nt, fig=fig, embed=True):
    [ev.evolve(steps) for ev in evs]
    ax.cla()
    n0, n1 = [ev.y.get_density()[0] for ev in evs]
    for (n, q) in zip((n0, n1), (q0, q1)):
        ax.plot(s.x, n, label=f"$q={q/qs:.2f}q_s$")
    ax.set(ylim=(0.6,1.2))
    ax.legend(loc="lower right")

CPU times: user 85 μs, sys: 0 ns, total: 85 μs
Wall time: 88.7 μs

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[9], line 1
----> 1 get_ipython().run_cell_magic('time', '', 'rng = np.random.default_rng(seed=2)\ne = Experiment(Nx=64, gamma=0.0, V0=-0.01)  # Slight attractive barrier - no g_ab\nqs = np.sqrt(e.m*e.gs[0]*e.na)   # q for sound.\ns0 = e.get_state()\ns1 = e.get_state()\nqn = np.where(s.basis.kx > qs)[0][0]\nq0, q1 = s.basis.kx[qn-1:qn+1]    # Modes just below and above threshold\nfor s, q in zip((s0, s1), (q0, q1)):\n    s.set_psi(s.get_psi() * np.exp(1j * np.array([q, 0])[:, None] * s.x[None,:]))\n\nt_growth = 10\nT = 10*t_growth\ndt = 0.4*s.t_scale\nNt = 20\ndT = T/Nt\nsteps = int(np.ceil(dT/dt))\ndt = dT/steps\n\nev0 = EvolverABM(s0, dt=dt)\nev1 = EvolverABM(s1, dt=dt)\nevs = [ev0, ev1]\n\nfig, ax = plt.subplots()\nfor n in FPS(Nt, fig=fig, embed=True):\n    [ev.evolve(steps) for ev in evs]\n    ax.cla()\n    n0, n1 = [ev.y.get_density()[0] for ev in evs]\n    for (n, q) in zip((n0, n1), (q0, q1)):\n        ax.plot(s.x, n, label=f"$q={q/qs:.2f}q_s$")\n    ax.set(ylim=(0.6,1.2))\n    ax.legend(loc="lower right")\n')

File <timed exec>:2
      1 'Could not get source, probably due dynamically evaluated source code.'

NameError: name 'Experiment' is not defined

We now run the same simulation, without the barrier, but with the two fluids flowing past each other with the same relative velocity. If \(\gamma < 3/4 = 0.75\), then we don’t see an instability, even at a relative flow velocity of \(c\), but if we increase \(\gamma\) slightly, we see the instability:

%%time
from gpe.contexts import FPS
from pytimeode.evolvers import EvolverABM
rng = np.random.default_rng(seed=2)
e0 = Experiment(Nx=64, gamma=0.6, V0=-0.0)
e1 = Experiment(Nx=64, gamma=0.76, V0=-0.0)
qs = np.sqrt(e0.m*e.gs[0]*e.na)
s = e0.get_state()
qn = np.where(s.basis.kx > qs/2)[0][0]  # Just above the speed of sound.
q = s.basis.kx[qn]
ss = [e.get_state() for e in [e0, e1]]
for s in ss:
    s.set_psi(
        s.get_psi() 
        * np.exp(1j * np.array([q, -q])[:, None] * s.x[None,:])
        * np.exp(1j*(rng.random(size=s.shape)-0.5)*0.01))

# Evolve to t_growth
t_growth = e1.get_t_growth(qn=qn)
T = 10*t_growth
dt = 0.4*s.t_scale
Nt = 20
dT = T/Nt
steps = int(np.ceil(dT/dt))
dt = dT/steps

evs = [EvolverABM(s, dt=dt) for s in ss]

fig, ax = plt.subplots()
for n in FPS(Nt, fig=fig, embed=True):
    [ev.evolve(steps) for ev in evs]
    ax.cla()
    for ev in evs:
        na, nb = ev.y.get_density()
        e = ev.y.experiment
        l, = ax.plot(s.x, na, label=f"$\gamma={e.gamma:.2f}$")
        ax.plot(s.x, nb, ls='--', c=l.get_c())
        
    ax.set(ylim=(0,2.0), title=f"$q={q/qs:.2f}q_s$")
    ax.legend(loc="lower right")

CPU times: user 62 μs, sys: 0 ns, total: 62 μs
Wall time: 64.8 μs

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[10], line 1
----> 1 get_ipython().run_cell_magic('time', '', 'from gpe.contexts import FPS\nfrom pytimeode.evolvers import EvolverABM\nrng = np.random.default_rng(seed=2)\ne0 = Experiment(Nx=64, gamma=0.6, V0=-0.0)\ne1 = Experiment(Nx=64, gamma=0.76, V0=-0.0)\nqs = np.sqrt(e0.m*e.gs[0]*e.na)\ns = e0.get_state()\nqn = np.where(s.basis.kx > qs/2)[0][0]  # Just above the speed of sound.\nq = s.basis.kx[qn]\nss = [e.get_state() for e in [e0, e1]]\nfor s in ss:\n    s.set_psi(\n        s.get_psi() \n        * np.exp(1j * np.array([q, -q])[:, None] * s.x[None,:])\n        * np.exp(1j*(rng.random(size=s.shape)-0.5)*0.01))\n\n# Evolve to t_growth\nt_growth = e1.get_t_growth(qn=qn)\nT = 10*t_growth\ndt = 0.4*s.t_scale\nNt = 20\ndT = T/Nt\nsteps = int(np.ceil(dT/dt))\ndt = dT/steps\n\nevs = [EvolverABM(s, dt=dt) for s in ss]\n\nfig, ax = plt.subplots()\nfor n in FPS(Nt, fig=fig, embed=True):\n    [ev.evolve(steps) for ev in evs]\n    ax.cla()\n    for ev in evs:\n        na, nb = ev.y.get_density()\n        e = ev.y.experiment\n        l, = ax.plot(s.x, na, label=f"$\\gamma={e.gamma:.2f}$")\n        ax.plot(s.x, nb, ls=\'--\', c=l.get_c())\n\n    ax.set(ylim=(0,2.0), title=f"$q={q/qs:.2f}q_s$")\n    ax.legend(loc="lower right")\n')

File <timed exec>:4
      1 'Could not get source, probably due dynamically evaluated source code.'

NameError: name 'Experiment' is not defined

Here we choose \(q\) slightly greater than \(q_s/2\) so that the relative velocity is above the Landau critical velocity. The counterflow instability develops only if \(\gamma > 0.75\).

Boosting

Another way to realize the counterflow instability is to consider independent coordinate transformations for the two species. Consider a single component gas in a time-dependent force (linear potential):

\[\begin{gather*} \I\hbar \dot{\psi}(x, t) = H(\op{p}, x)\psi(x, t) - a(t)x\psi(x, t). \end{gather*}\]

The linear potential breaks the periodicity of the solution, frustrating attempts to use a Fourier basis. However, we can restore periodicity by adding a phase:

\[\begin{align*} \psi(x, t) =& e^{\I\phi(x, t)}\Psi(x, t),\\ \dot{\psi} =& \I \dot{\phi}e^{\I\phi}\Psi + e^{\I\phi}\dot{\Psi}\\ \psi_{,x} =& e^{\I\phi}(\nabla + \I\phi_{,x})\Psi,\\ \psi_{,xx} =& e^{\I\phi}(\nabla + \I\phi_{,x})^2\Psi. \end{align*}\]

The last pattern allows us to write

\[\begin{gather*} H(\op{p}, x)\psi(x, t) = e^{\I\phi}H\Bigl(\op{p} + \hbar\phi_{,x}, x\Bigr)\Psi. \end{gather*}\]

Hence, we have the following equation of motion for \(\Psi(X, t)\):

\[\begin{gather*} \I\hbar\dot{\Psi} = \Biggl(H\Bigl(\op{p} + \hbar\phi_{,x}, X + x_0(t)\Bigr) + \hbar \dot{\phi}\Biggr)\Psi. \end{gather*}\]

Thus, we have three effects

1. The dispersion relationship gets modified:

   \[\begin{gather*} E(\op{p}) \rightarrow E(\op{p} + \hbar \phi_{,x}). \end{gather*}\]

   This is done in the kinetic energy – we do not change the meaning of \(k\) in our basis, which should simply be thought of as specifying which basis functions are being used for the periodic function \(\Psi(x,t)\).

2. We have a new potential:

   \[\begin{gather*} V(x, t) \rightarrow V(x, t) + \hbar \dot{\phi}. \end{gather*}\]

3. The physical wavefunction has an additional phase. This affects the computation of quantities like the current density. Specifically

   \[\begin{gather*} j = \Re (\psi^\dagger \op{p} \psi) = \Re \bigl(\Psi^\dagger (\op{p} + \hbar \phi_{,x})\Psi\bigr). \end{gather*}\]

Going back to our original problem, we can cancel the linear term by defining

\[\begin{gather*} \hbar \dot{\phi} = a(t)x, \qquad \hbar\phi(x,t) = A(t) x, \qquad \dot{A}(t) = a(t). \end{gather*}\]

The periodicity of the problem is now restored with a modified dispersion. Note that the wavefunction \(\psi(x, t) = e^{\I\phi(x, t)}\Psi(x, t)\) is not periodic anymore unless \(\phi(0, t) = \phi(L, t) + 2\pi n\), but the density \(n(x, t)\) remains periodic. This is a form of dynamical twisted boundary conditions which allows for phase to continuously build up in \(\psi(x, t)\) without discontinuous phase slips. For practical reasons, we must provide the integrated \(A(t)\) to the code, not the actual force \(F = m a(t)\).

For a two-component system, we can do this independently to the two components, dynamically inducing the counterflow instability.

Details.

Originally I included a shift \(X = x - x_0(x)\), but realized after that this was not needed. Here we keep the shift for generality. This implies that

\[\begin{gather*} \pdiff{X}{x} = 1. \end{gather*}\]

Hence, \(\partial_{X} \equiv \partial_{x} \equiv \nabla\). Now expand the derivatives and apply both product and chain rules:

\[\begin{align*} \psi_{,x} =& \I \phi_{,x}e^{\I\phi}\Psi + e^{\I\phi}\Psi_{,X} = e^{\I\phi}(\nabla + \I\phi_{,x})\Psi,\\ \psi_{,xx} =& (\I \phi_{,xx} - \phi_{,x}^2) e^{\I\phi}\Psi + 2\I\phi_{,x}e^{\I\phi}\Psi_{,X} + e^{\I\phi}\Psi_{,XX},\\ (\partial_{X} + \I\phi_{,x})^2\Psi =& (\partial_{X} + \I\phi_{,x})(\I\phi_{,x}\Psi + \Psi_{,X})\\ =& \I\phi_{,xx}\Psi + 2\I\phi_{,x}\Psi_{,X} + \Psi_{,XX} -\phi_{,x}^2\Psi + \I\phi_{,x}\Psi_{,X}. \end{align*}\]

The transformation effected on the momentum is:

\[\begin{gather*} \op{p} = -\I\hbar \nabla \rightarrow -\I\hbar (\nabla + \I \phi_{,x}) = -\I\hbar \nabla + \hbar\phi_{,x} \end{gather*}\]

Note: this only strictly works for Hamiltonians of the form

\[\begin{gather*} H(\op{p}, \op{x}) = E(\op{p}) + V(\op{x}). \end{gather*}\]

Terms like spin-orbit coupling, or angular momenta with products of \(\op{x}\) and \(\op{p}\) will have more complicated transformation properties.

Putting everything together

\[\begin{gather*} \I\hbar\Bigl(\I \dot{\phi}e^{\I\phi}\Psi + e^{\I\phi}\dot{\Psi} - e^{\I\phi}\dot{x}_0(t)\Psi_{,X}\Bigr) = e^{\I\phi}H\Bigl(\op{p} + \hbar\phi_{,x}, X + x_0(t)\Bigr)\Psi\\ \I\hbar\dot{\Psi} = H\Bigl(\op{p} + \hbar\phi_{,x}, X + x_0(t)\Bigr)\Psi +\hbar \dot{\phi}\Psi + \I\hbar\dot{x}_0(t)\Psi_{,X}\\ = H\Bigl(\op{p} + \hbar\phi_{,x}, X + x_0(t)\Bigr)\Psi +\hbar \dot{\phi}\Psi - \dot{x}_0(t)\op{p}\Psi_{,X}. \end{gather*}\]

As noted above, we don’t need the shift \(x_0(t)\). However, if this were included, we would then we need to interpolate from one abscissa to the other to compute the non-linear terms. This can be done by adding a phase to the FFT:

\[\begin{gather*} \psi(x) = \sum_{k} \psi_{k} e^{\I k x}\\ \psi(x-x_0) = \sum_{k} \psi_{k} e^{\I k (x - x_0)} = \sum_{k} \Bigl(e^{-\I k x_0}\psi_{k}\Bigr) e^{\I k}. \end{gather*}\]

Example

Here we consider an example of applying a constant force to the two components (in opposite directions):

\[\begin{gather*} a(t) = a, \qquad A(t) = at, \qquad \phi_{\pm}(x, t) = \pm a \hbar t x. \end{gather*}\]

from gpe.utils import _GPU

@_GPU.add_non_GPU_methods
class StateCounterflow(State):
    """State implementing the counterflow protocol.
    
    In terms of the notes here, `self.data` stores the periodic $\Psi(x, t)$.
    
    Follows the structure of gpe.bec2.StateScaleBase with the idea
    of eventual unification.
    """
    def init(self):
        super().init()
        for key in ["kx2", "k2"]:
            # These need to be recomputed as we scale, so we must delete any cached values.
            if hasattr(self, key):
                delattr(self, key)

    def _get_kx2_GPU(self, kx=None, Lx=None):
        """Return the effective `kx**2` and `twist_phase_x`  for the kinetic
        energy matrix.  This will be passed as the `kx2` and `twist_phase_x`
        arguments to the `basis.lagrangian` function.
        """
        x = self.basis.xyz[0]
        phi_x = self.experiment.m * self.experiment.a * self.t / self.hbar
        if kx is None:
            kx = self.xp.array([self.basis.kx + phi_x, 
                                self.basis.kx - phi_x])

        if Lx is None:
            Lx = self.basis.Lx

        kx2 = kx**2
        twist_phase_x = self.xp.asarray([self.xp.exp(1j*phi_x * x), 
                                         self.xp.exp(-1j*phi_x * x)])
        return kx2, twist_phase_x

    def get_psi_GPU(self):
        # Includes scaling factor and phase
        kx2, twist_phase_x = self._get_kx2_GPU()
        return self.data * twist_phase_x

    def set_psi(self, psi):
        # Includes scaling factor and phase
        kx2, twist_phase_x = self._get_kx2()
        self.set_data(psi / twist_phase_x)

    ######################################################################
    # Required by interface IStateDFT
    #
    def apply_laplacian(self, factor, exp=False, **_kw):
        """Apply the laplacian multiplied by `factor` to the state.

        Arguments
        ---------
        factor : array-like
           The result will be multiplied by this factor.
        exp : bool
           If `True` then `exp(factor*laplacian)(y)` will be computed instead.
        """
        # Here we need to add the scale factors.  We do this with the
        # k2 argument.
        if "k2" not in _kw:
            kx2, twist_phase_x = self._get_kx2_GPU()
            k2 = (self.xp.asarray(sum(k**2 for k in self.basis._pxyz[1:]))[None, None, ...] 
                  + kx2)
            _kw["k2"] = k2

        super().apply_laplacian(factor=factor, exp=exp, **_kw)


class ExperimentCounterflowBase(Experiment):
    a = 1.0  # Relative acceleration
    State = StateCounterflow
    
class ExperimentCounterflow(ExperimentCounterflowBase):
    tc = 1.0   # Time when acceleration should give critical velocity
    def init(self):
        super().init()
        vc = self.hbar * self._qc / self.m
        self.a = vc / self.tc

def get_angle(psi):
    N = len(psi)
    angles = [np.angle(psi[0])]
    for n in range(1, len(psi)):
        th = angles[-1]
        angles.append(th + np.angle(psi[n]/psi[n-1]))
    return np.array(angles) - angles[N//2]

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[11], line 72
     68 
     69         super().apply_laplacian(factor=factor, exp=exp, **_kw)
     70 
     71 
---> 72 class ExperimentCounterflowBase(Experiment):
     73     a = 1.0  # Relative acceleration
     74     State = StateCounterflow
     75 

NameError: name 'Experiment' is not defined

e = ExperimentCounterflow(tc=30, Nx=256)
s = e.get_state()
s.set_psi(s.get_psi() * np.exp(1j*(rng.random(size=s.shape)-0.5)*0.001) 
          * (1 + 0.01*np.exp(-s.x**2)))

T = 6*e.tc
dt = 0.4*s.t_scale
Nt = 100
dT = T/Nt
steps = int(np.ceil(dT/dt))
dt = dT/steps

ev = EvolverABM(s, dt=dt)
states = [s]

fig, axs = plt.subplots(2, 1)
for n in FPS(Nt, display=True, fig=fig, embed=True):
    ev.evolve(steps)
    states.append(ev.get_y())
    s = states[-1]
    [ax.cla() for ax in axs]
    plt.sca(axs[0])
    s.plot()
    axs[1].plot(s.x, get_angle(s.get_psi()[0]))
    axs[1].plot(s.x, get_angle(s.get_psi()[1])) 

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[12], line 1
----> 1 e = ExperimentCounterflow(tc=30, Nx=256)
      2 s = e.get_state()
      3 s.set_psi(s.get_psi() * np.exp(1j*(rng.random(size=s.shape)-0.5)*0.001) 
      4           * (1 + 0.01*np.exp(-s.x**2)))

NameError: name 'ExperimentCounterflow' is not defined