--- execution: timeout: 60 jupytext: cell_metadata_filter: -all formats: md:myst,ipynb main_language: python text_representation: extension: .md format_name: myst format_version: 0.13 jupytext_version: 1.17.1 kernelspec: name: python3 display_name: Python 3 (ipykernel) language: python --- ```{code-cell} 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 ``` (sec:Counterflow)= # Counterflow Instabilities ## Bogoliubov Modes (Phonons) :::{margin} 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 {ref}`sec:CounterflowLinearResponse`, with the simplification considered in {cite}`Alperin: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 :::{margin} 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. ```{code-cell} _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(); ``` 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}}$. ```{code-cell} 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()) ``` ```{code-cell} fig, ax = plt.subplots() for s in FPS(states, display=True, fig=fig, embed=True): ax.cla() s.plot() ``` ```{code-cell} 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() ``` ```{code-cell} 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 ``` ```{code-cell} 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}}$"); ``` 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). ```{code-cell} 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$."); ``` :::{margin} 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 :::{margin} The main point of {cite}`Huynh:2023` is that there are particular parameter regimes where stationary **supersonic** flow past a barrier exists. See also {cite}`Leboeuf:2001` and {cite}`Paris-Mandoki:2017`. ::: Here we explore the Landau criterion using the same simulation. {cite}`Huynh: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. ```{code-cell} :tags: [hide-input] %%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") ``` 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: ```{code-cell} %%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") ``` 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. :::{admonition} Details. :class: dropdown 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*} ```{code-cell} 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] ``` ```{code-cell} 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])) ```