Linear Response

Hide code cell content

 
import mmf_setup; mmf_setup.nbinit()
import os
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 inline
import numpy as np, matplotlib.pyplot as plt
try: from myst_nb import glue
except: glue = None

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.

Linear Response#

General Theory - Sum Rules#

A general formulation of linear response is discussed in [Pitaevskii and Stringari, 2003] and we summarize their results here. The idea is to start consider two linear operators: one \(\op{G}\) that creates the excitation, and another \(\op{F}\) that we measure to detect the excitation. The system evolves with a hermitian Hamiltonian

\[\begin{gather*} \op{H}_{\lambda}(t) = \op{H} + \lambda e^{0^+ t}\op{H}_{\text{pert}}(t), \qquad \op{H}_{\text{pert}} = -(\op{G}e^{\omega t/\I} + \text{h.c.}). \end{gather*}\]

We then consider the linear response in terms of the expectation-value of the operator \(\op{F}\):

\[\begin{gather*} \delta\braket{\op{F}} = \lambda e^{0^+ t}( \chi_{F, G}(\omega)e^{\omega t /\I} + \text{h.c}) = \lambda e^{0^+ t}( \chi_{F, G}(\omega)e^{\omega t /\I} + \underbrace{\chi_{F, G}(-\omega)}_{\chi_{F, G}^*(\omega)}e^{-\omega t /\I}). \end{gather*}\]

If we consider a time-independent Hamiltonian \(\op{H}\) and start with a system in thermal equilibrium at \(t=-\infty\) with temperature \(T\), we obtain

\[\begin{gather*} \chi_{F,G}(\omega) = -\frac{1}{Q\hbar}\sum_{mn} e^{-\beta E_m} \Biggl( \frac{\braket{m|\op{F}|n}\braket{n|\op{G}|m}} {\omega_+ - \omega_{nm}} - \frac{\braket{m|\op{G}|n}\braket{n|\op{F}|m}} {\omega_+ + \omega_{nm}} \Biggr). \tag{7.3} \end{gather*}\]

Do It! Derive this.

Hint: Consider the density matrix.

Recall that a thermal ensemble can be described in terms of the energy eigenstates \(\ket{n}\) by a density matrix

\[\begin{gather*} \op{R} = \sum_{n} \ket{n}e^{-\beta E_{n}}\bra{n}, \qquad \op{H}\ket{n} = \ket{n}E_n, \qquad \beta = \frac{1}{k_B T},\\ \braket{\op{F}} = \Tr \op{R}\op{F} = \sum_{n}\frac{e^{-\beta E_n}}{Q}\braket{n|\op{F}|n},\qquad Q = \sum_{n} e^{-\beta E_n}. \end{gather*}\]

The time-evolution of this state is

\[\begin{gather*} \op{R}(t) = \sum_{n} \ket{n(t)}e^{-\beta E_{n}}\bra{n(t)}, \\ \I\hbar\partial_t\ket{n(t)} = \op{H}_{\lambda}(t)\ket{n(t)}, \qquad \ket{n(-\infty)} = \ket{n}. \end{gather*}\]

Solution

I find it easiest to work with the operator representations of these equations:

\[\begin{gather*} \I \hbar \partial_{t}\op{R}(t) = [\op{H}_{\lambda}(t), \op{R}(t)], \qquad \op{R}(-\infty) = e^{-\beta \op{H}}. \end{gather*}\]

At linear order, there will be no mixing of frequencies, thus, since our perturbation contains only \(e^{\pm \omega t/\I}\), the response will have the same form. We thus include this from the start (but one can check this more generally if desired). Expanding to linear order, we have

\[\begin{gather*} \op{R}(t) = \op{R} + \lambda \underbrace{(\delta\op{R}_{+}e^{\omega t/\I} + \text{h.c.})}_{\delta \op{R}(t)} + O(\lambda^2),\\ \I\hbar \partial_t \delta \op{R}(t) = [\op{H}, \delta\op{R}(t)] + [\op{H}_{\text{pert}}(t), \op{R}] + O(\lambda),\\ \hbar \omega_+ \delta \op{R}_{+} = [\op{H}, \delta\op{R}_{+}] - [\op{G}, \op{R}]\\ -\hbar \omega_- \delta \op{R}_{+}^\dagger = [\op{H}, \delta\op{R}_{+}^\dagger] - [\op{G}^\dagger, \op{R}] \end{gather*}\]

where we have included the convergence factors in

\[\begin{gather*} \omega_{\pm} = \omega \pm \I 0^+. \end{gather*}\]

Now we can compute the response

\[\begin{gather*} Q\delta \braket{\op{F}} = \lambda \Bigl( e^{\omega_+ t/\I}\Tr(\delta\op{R}_+\op{F}) + e^{-\omega_- t/\I}\Tr(\delta\op{R}_+^\dagger\op{F}) \Bigr) \end{gather*}\]

The desired result comes from inserting a complete set of states

\[\begin{multline*} Q\delta \braket{\op{F}} = \lambda \Bigl( e^{\omega_+ t/\I}\sum_{mn}\braket{n|\delta\op{R}_+|m}\braket{m|\op{F}|n} \\ + e^{-\omega_- t/\I} \sum_{mn}\braket{n|\delta\op{R}_+^\dagger|m}\braket{m|\op{F}|n} \Bigr). \end{multline*}\]

Taking the matrix elements of the linear response equation, and using the fact that \(\op{H}\ket{n} = \ket{n}E_n\) and \(\op{R}\ket{n} = \ket{n}e^{-\beta E_n}\) are diagonal, we obtain

\[\begin{gather*} \hbar \omega_+ \braket{n|\delta\op{R}_+|m} = \Bigl( (E_n-E_m)\braket{n|\delta\op{R}_{+}|m} + (e^{-\beta E_n}-e^{-\beta E_m})\braket{n|\op{G}|m} \Bigr),\\ \braket{n|\delta\op{R}_+|m} = \frac{e^{-\beta E_n}-e^{-\beta E_m}}{\hbar \omega_+ - \underbrace{(E_n-E_m)}_{\hbar\omega_{nm}}} \braket{n|\op{G}|m}. \end{gather*}\]

Assembling everything, we have

\[\begin{gather*} \delta \braket{\op{F}} = \frac{\lambda}{\hbar Q} \Biggl( e^{\omega_+ t/\I} \sum_{mn}\frac{(e^{-\beta E_n} - e^{-\beta E_m})\braket{n|\op{G}|m}\braket{m|\op{F}|n}}{\omega_+ - \omega_{nm}} + \text{h.c.}\Biggr), \end{gather*}\]

from which we can read off the stated result, after rearranging some indices.

The sum-rule associated with an operator \(\op{F}\) has the form (see [Wang, 1999])

\[\begin{gather*} \sum_{m}(E_{m} - E_{0})\abs{\braket{m|\op{F}|0}}^2 = \braket{0|[\op{F},\op{H}]\op{F}^\dagger|0}. \end{gather*}\]

These are useful if the final commutators have simple forms.

Do it! Derive this.

Simply insert a complete set of energy eigenstates:

\[\begin{align*} \braket{0|[\op{F},\op{H}]\op{F}^\dagger|0} &= \sum_{m}\braket{0|[\op{F},\op{H}]|m}\braket{m|\op{F}^\dagger|0}\\ &= \sum_{m}(\braket{0|\op{F}\op{H}|m}-\braket{0|\op{H}\op{F}|m})\braket{m|\op{F}^\dagger|0}\\ &= \sum_{m}(\braket{0|\op{F}|m}E_{m} - E_0\braket{0|\op{F}|m})\braket{m|\op{F}^\dagger|0}\\ &= \sum_{m}(E_{m} - E_{0})\braket{0|\op{F}|m}\braket{m|\op{F}^\dagger|0}\\ = \sum_{m}(E_{m} - E_{0})\abs{\braket{m|\op{F}|0}}^2. \end{align*}\]

Here we consider the fluctuations about a stationary state \(\psi_0\) of GPE-like theories to linear order – sometimes called Bogoliubov theory.

Bogoliubov Theory#

We assume that the equations of motion (EoM) arise from a principle of stationary action:

\[\begin{gather*} S[\psi] = \int \left(\int \I\hbar \psi^\dagger\dot{\psi} \d{x} - E[\psi]\right)\d{t},\\ \I\hbar \dot{\psi} = \op{H}[\psi]\psi = \frac{\delta E[\psi]}{\delta \psi^\dagger}. \end{gather*}\]

For standard GPE-like theories, we have

\[\begin{gather*} E[\psi] = \int \d{x}\; \left( \psi^{\dagger}K(\op{p})\psi + \mathcal{E}(n) \right), \qquad n = \abs{\psi}^2,\\ \I\hbar \dot{\psi} = \Bigl(K(\op{p}) + \mathcal{E}'(n)\Bigr)\psi, \end{gather*}\]

where \(\mathcal{E}(n)\) is the equation of state and \(K(\op{p})\) is the single-particle dispersion relation (kinetic energy). More generally, we might have an external potential, chemical potential, etc. We roll all of these into \(\mathcal{E}(n)\) for this discussion and keep this general. Thus, for the standard GPE we would have

\[\begin{gather*} K_{GPE}(\op{p}) = \frac{\op{p}^2}{2m} = \frac{-\hbar^2\nabla^2}{2m},\\ \mathcal{E}_{GPE}(n) = g\frac{n^2}{2} + V(x, t)n,\\ \I\hbar \dot{\psi} = \Biggl( \underbrace{\frac{-\hbar^2\nabla^2}{2m}}_{K(\op{p})} + \underbrace{gn + V(x, t)}_{\mathcal{E}'_{GPE}(n)} \Biggr)\psi. \end{gather*}\]

We now replace the spatial dependence with the usual Hilbert-space structure, and express the equations of motion in terms of kets:

\[\begin{gather*} \I\hbar \dot{\ket{\psi}} = \op{H}[\psi]\ket{\psi}, \qquad \psi(x) = \braket{x|\psi}. \end{gather*}\]

The idea of Bogoliubov theory is to consider fluctuations about a stationary state \(\ket{\psi_0}\):

\[\begin{gather*} \overbrace{\op{H}[\psi_0]}^{\op{H}_0}\ket{\psi_0} = \ket{\psi_0}\mu,\qquad \ket{\psi} = e^{\mu t/\I\hbar}\Bigl( \ket{\psi_0} + \overbrace{\underbrace{\ket{u_+}}_{\ket{u}}e^{\omega t/\I} + \underbrace{\ket{u_-}}_{-\ket{v^*}}e^{-\omega^* t/\I}}^{\ket{u}} \Bigr), \end{gather*}\]

finding \(\ket{u_{\pm}}\) that satisfy the equations of motion to linear order. The presence of the non-linear term \(n(x) = \braket{x|\psi}\braket{\psi|x}\) mixes the positive and negative frequency components. Specifically

\[\begin{alignat*}{2} n(x) &= n_0(x) &&+ \overbrace{u(x)\psi_0^*(x) + \text{h.c.}}^{\delta n(x)} \\ &= n_0(x) && + u_{+}(x)e^{\omega t/\I}\psi_0^*(x) + u_{-}(x)e^{-\omega^* t/\I}\psi_0^*(x)\\ &&& + u_{+}^*(x)e^{-\omega^* t/\I}\psi_0(x) + u_{-}^*(x)e^{\omega t/\I}\psi_0(x),\\ &= n_0(x) &&+ e^{\omega t/\I}\Bigl(u_+(x)\psi_0^*(x) + u_-^*(x)\psi_0(x)\Bigr)\\ &&&+ e^{-\omega^* t/\I}\Bigl(u_-(x)\psi_0^*(x) + u_+^*(x)\psi_0(x)\Bigr) + O(u^2). \end{alignat*}\]

Expanding the Hamiltonian to linear order, we thus have

\[\begin{gather*} \op{H}[\psi] = \op{H}_0 + \mathcal{E}''(n_0)\overbrace{(u\psi_0^* + u^*\psi_0)(\op{x})}^{ u(\op{x})\psi_0^*(\op{x}) + u^*(\op{x})\psi_0(\op{x})} + O(u)^2. \end{gather*}\]

Finally, expanding the equation of motion and multiplying through by \(e^{-\mu t /\I\hbar}\), we have

\[\begin{gather*} \mu \overbrace{\ket{\psi}e^{-\mu t /\I\hbar}} ^{\ket{\psi_0} + \ket{u}} + \I\hbar\dot{\ket{u}} = \overbrace{\op{H}_0\ket{\psi_0}}^{\ket{\psi_0}\mu} + \op{H}_0\ket{u} + \mathcal{E}''(n_0)(u\psi_0^* + u^*\psi_0)(\op{x})\ket{\psi_0} + O(u^2),\\ \ket{u}\mu + \I\hbar\dot{\ket{u}} = \op{H}_0\ket{u} + \mathcal{E}''(n_0)(u\psi_0^* + u^*\psi_0)(\op{x})\ket{\psi_0} + O(u^2). \end{gather*}\]

We now expand \(\ket{u} = \ket{u_+}e^{\omega t/\I} + \ket{u_-}e^{-\omega^* t/\I}\), and collect the positive and negative frequency terms, demanding that both vanish independently.

\[\begin{gather*} u_+\mu + \hbar\omega u_+ = \braket{x|\op{H}_0|u_+} + \mathcal{E}''(n_0)(u_+\psi_0^* + u_{-}^*\psi_0)\psi_0,\\ u_-\mu - \hbar\omega^* u_- = \braket{x|\op{H}_0|u_-} + \mathcal{E}''(n_0)(u_-\psi_0^* + u_{+}^*\psi_0)\psi_0. \end{gather*}\]

These can be packaged as a generalized eigenvalue equation after conjugating the second equation

\[\begin{gather*} \begin{pmatrix} \op{H}_0 - \mu + n_0\mathcal{E}''(n_0) & \psi_0^2\mathcal{E}''(n_0)\\ (\psi^*_0)^2\mathcal{E}''(n_0) & \op{H}_0^* - \mu + n_0\mathcal{E}''(n_0) \end{pmatrix} \begin{pmatrix} u_+\\ u_-^* \end{pmatrix} = \hbar\omega \begin{pmatrix} 1\\ & -1 \end{pmatrix} \begin{pmatrix} u_+\\ u_-^* \end{pmatrix}. \end{gather*}\]

On Conjugate Pairs of Eigenvalues: \(\omega_n\), \(-\omega_n^*\).

This is a generalized eigenvalue problem of the form \(\mat{A}\ket{n} = \mat{B}\ket{n}\omega_{n}\) with Hermitian matrices \(\mat{A} = \mat{A}^\dagger\) and \(\mat{B} = \mat{B}^\dagger\). Numerically this generalized eigenvalue problem is easy to work with if \(\mat{B}\) is positive definite – i.e. has positive eigenvalues (see e.g. scipy.linalg.eigh()). In this case we can use its diagonalization to write

\[\begin{gather*} \mat{B} = \mat{U}^\dagger\mat{D}\mat{U} =\mat{U}^\dagger\sqrt{\mat{D}}\mat{U} \underbrace{\mat{U}^\dagger\sqrt{\mat{D}}\mat{U}}_{\mat{b}} = \mat{b}^\dagger\mat{b}, \end{gather*}\]

where \(\mat{b}\) is invertable since \(\mat{D}\) is diagonal with positive entries.

Then we have a hermitian problem \(\tilde{\mat{A}}\tilde{\ket{n}} = \tilde{\ket{n}}\omega_n\) with

\[\begin{gather*} \mat{A}\mat{b}^{-1}\mat{b}\ket{n} = \mat{b}^\dagger\mat{b}\ket{n}\omega_{n}\\ \underbrace{(\mat{b}^{-1})^\dagger\mat{A}\mat{b}^{-1}}_{\tilde{\mat{A}}} \underbrace{\mat{b}\ket{n}}_{\tilde{\ket{n}}} = \underbrace{\mat{b}\ket{n}}_{\tilde{\ket{n}}}\omega_{n}. \end{gather*}\]

I.e., we can transform the problem into a standard Hermitian eigenvalue problem.

In our case, unfortunately, \(\mat{B}\) is not positive definite, making the problem somewhat challenging. There are some nice properties of our system, however (see e.g. [Wu and Niu, 2003]). Specifically, we can write the as \(\mat{B} = \mat{Z}\) where \(\mat{Z} = \mat{\sigma}_{z}\otimes \mat{1}\). These following from the property that

\[\begin{gather*} \mat{X}\mat{A}\mat{X} = \mat{A}^*, \qquad \mat{X}\mat{B}\mat{X} = -\mat{B}^* \end{gather*}\]

where \(\mat{X} = \mat{\sigma}_{x}\otimes \mat{1}\). Hence, if \(\ket{n}\) is an eigenvector with eigenvalue \(\omega_{n}\), then \(\mat{X}\ket{n}^*\) is an eigenvector with eigenvalue \(-\omega_{n}^*\):

\[\begin{gather*} \mat{A}\ket{n} = \mat{B}\ket{n}\omega_n\\ \mat{X}\mat{A}\mat{X}\mat{X}\ket{n} = \mat{X}\mat{B}\mat{X}\mat{X}\ket{n}\omega_n\\ \mat{A}^*\mat{X}\ket{n} = -\mat{B}^*\mat{X}\ket{n}\omega_n\\ \mat{A}\mat{X}\ket{n}^* = \mat{B}\mat{X}\ket{n}^*(-\omega_n^*). \end{gather*}\]

Thus, all eigenvalues with non-zero real parts appear in conjugate pairs \(\omega_n\) and \(-\omega_n^*\). Purely imaginary eigenvalues may in principle be isolated, but will have eigenvectors which satisfy \(\mat{X}\ket{n}^* = \ket{n}\).

In many cases – for example, if the dispersion relationship is quadratic \(K(p) = p^2/2m\), the matrix \(\mat{A}=\mat{A}^*\) is self-conjugate and \(\mat{B}^* = -\mat{B}\). (This may not be the case with more complicated dispersions that break parity for example.) In these cases, we can further simplify the equations. Introduce the following projectors:

\[\begin{gather*} \mat{P}_{\pm} = \frac{\mat{1} \pm \mat{Z}}{2}, \qquad \mat{X}\mat{P}_{\pm}\mat{X} = \mat{P}_{\pm}.\\ \ket{n} = \frac{\mat{X} + \mat{Z}}{\sqrt{2}}\ket{\tilde{n}},\\ \mat{A}\ket{n} = \frac{\mat{A}\mat{X} + \mat{A}\mat{Z}}{\sqrt{2}}\ket{\tilde{n}},\\ \end{gather*}\]
\[\begin{gather*} \mat{P}_{+}\mat{A}\mat{P}_{+}\mat{P}_{-}\mat{A}\mat{P}_{-}\ket{\tilde{n}} =\ket{\tilde{n}}\omega_n \end{gather*}\]

HAH = P()P + M()M

HAH = (X+Z)A(X+Z)/2 = (XAX + ZAZ + XAZ + ZAX)/2 = (4A + ZAZ + AXZ + ZXA)/2

On the \(\psi_0^2\) off-diagonals: \(\tilde{K}(\op{p})\).

The factor of \(\psi_0^2\) and conjugate in the off-diagonals is a nuisance, since the corresponding factor of \(n_0\) appears on the diagonals. In many cases, one can take the stationary state \(\psi_0 = \psi_0^*\) to be real, in which case \(\psi_0^2 = n_0\), but there are examples where this cannot be done (states with twisted boundary conditions for example, as discussed below).

In these cases, the terms differ by a phase \(\psi_0^2 = e^{2\I\phi}n_0\) which can be absorbed in a redefinition of \(u_{-} \rightarrow e^{2\I\phi} u_{-}^*\). If the phase is spatially dependent, however, then this operation will generally not commute with the kinetic energy but can be convenient. The argument goes like this (and shows a more streamlined calculation):

First include the chemical potential in \(\op{H}_0 = \op{H}[\psi_0] - \mu\op{1}\) so that \(\op{H}_0\psi_0 = 0\), and consider the state

\[\begin{gather*} \op{H}_0\psi_0 = 0, \qquad \psi = \psi_0(1 + \epsilon) = \psi_0 + \psi_0\epsilon, \qquad \epsilon = \epsilon_+ e^{\omega t/\I} + \epsilon_- e^{-\omega t/\I}. \end{gather*}\]

The expansion of the density now contains only \(n_0\): no factors of \(\psi_0\) or \(\psi_0^*\) remain:

\[\begin{gather*} n = n_0(1 + \epsilon + \epsilon^*) + O(\epsilon^2). \end{gather*}\]

Inserting this into the EoM and keeping terms of order \(\epsilon\), we have

\[\begin{gather*} \I \hbar \psi_0\dot{\epsilon} = \op{H}_0(\psi_0\epsilon) + \mathcal{E}''(n_0)n_0(\epsilon + \epsilon^*)\psi_0. \end{gather*}\]

The last term – from the expansion of \(n\) – contains both \(\epsilon\) and \(\epsilon^*\). This is a consequence of the non-linear nature of the GPE-like equation and gives rise to a coupled set of equations for the linear response. In this case, there are no terms like \(\psi_0^2\) as in the original derivation. The factored perturbation simply gives \(n_0\) factors with the complication that \(\op{H}_0\) now acts on the product \(\psi_0\epsilon\). To proceed, we must define a new operator \(\op{\tilde{H}}_0\) that satisfies

\[\begin{gather*} \op{H}_0(\psi_0 \epsilon) = \psi_0\op{\tilde{H}_0}(\epsilon). \end{gather*}\]

Only the kinetic energy is a problem, so we really just need to an operator \(\tilde{K}(\op{p})\) such that

\[\begin{gather*} K(\op{p})(\psi_0\epsilon) = \psi_0\tilde{K}(\op{p})\epsilon. \end{gather*}\]

Then we can express these equations as

\[\begin{gather*} \begin{pmatrix} \op{\tilde{H}} + n_0\mathcal{E}'' - \hbar \omega & n_0\mathcal{E}''\\ n_0\mathcal{E}'' & \op{\tilde{H}}^* + n_0\mathcal{E}'' + \hbar \omega^* \end{pmatrix} \begin{pmatrix} \epsilon_+\\ \epsilon_-^* \end{pmatrix} = \vect{0},\\ \op{\tilde{H}} = \tilde{K}(\op{p}) + \mathcal{E}'(n_0) + V(x, t) - \mu. \end{gather*}\]

This modified approach will be useful when we can explicitly form \(\tilde{K}(\op{p})\).

Consider how this works for various powers of \(p\) as a consequence of the General Leibniz rule

\[\begin{align*} \op{p}(\psi_0\epsilon) &= (\op{p}\psi_0) \epsilon + \psi_0(\op{p}\epsilon)\\ \op{p}^2(\psi_0\epsilon) &= (\op{p}^2\psi_0) \epsilon + 2(\op{p}\psi_0)(\op{p}\epsilon) + \psi_0(\op{p}^2\epsilon),\\ \op{p}^n(\psi_0\epsilon) &= \sum_{m=0}^{n}\binom{n}{m}(\op{p}^{m}\psi_0)(\op{p}^{n-m} \epsilon). \end{align*}\]

For the usual quadratic dispersion \(K(\op{p}) = \op{p^2}/2m\) and expanding about a state \(\psi_0\) with momentum \(p_0\), we have

\[\begin{align*} \tilde{K}(\op{p}) = \frac{(p_0+\vect{p})^2}{2m} = \frac{p_0^2}{2m} + \frac{p_0\op{p}}{m} + \frac{\op{p}^2}{2m} = E_0 + v_0\op{p} + \frac{\op{p}^2}{2m},\\ \tilde{K}^*(\op{p}) = \frac{(p_0-\vect{p})^2}{2m} = \frac{p_0^2}{2m} - \frac{p_0\op{p}}{m} + \frac{\op{p}^2}{2m} = E_0 - v_0\op{p} + \frac{\op{p}^2}{2m}. \end{align*}\]

Note that \(\op{p}^* = -\op{p}\). This is exactly the transformation required to restore Galilean invariance. I.e., the normal modes about the state \(\psi_0\) with momentum \(p_0\) (and velocity \(v_0=p_0/m\)) should be the same as the normal modes in a co-moving frame where \(\psi_0\) has zero momentum. More on this below.

Solving the Bogoliubov Equations.

For this discussion, we will consider the following form of the Bogoliubov equations:

\[\begin{gather*} \begin{pmatrix} \mat{A} & \mat{B}\\ \mat{B} & \mat{A}^* \end{pmatrix} \begin{pmatrix} u_{+}\\ u_{-}^* \end{pmatrix} = \mat{z} \begin{pmatrix} u_{+}\\ u_{-}^* \end{pmatrix} \omega, \qquad \mat{z} = \begin{pmatrix} \mat{1}\\ & -\mat{1} \end{pmatrix}. \end{gather*}\]

This can be written as a non-hermitian eigenvalue problem by letting \(u=u_{+}\) and \(v = -u_{-}^*\), recovering the form of the Bogoliubov-de Gennes (BdG) equations seen in the Hartree-Fock-Bogoliubov (HFB) theory of superconductivity:

\[\begin{gather*} \begin{pmatrix} \mat{A} & -\mat{B}\\ \mat{B} & -\mat{A}^* \end{pmatrix} \begin{pmatrix} u\\ v \end{pmatrix} = \begin{pmatrix} u\\ v \end{pmatrix} \omega. \end{gather*}\]

To make further progress, consider

\[\begin{gather*} \begin{pmatrix} \alpha\\ \beta \end{pmatrix} = \frac{1}{\sqrt{2}} \begin{pmatrix} \mat{1} & \mat{1}\\ \mat{1} & -\mat{1} \end{pmatrix} \begin{pmatrix} u\\ v \end{pmatrix},\\ \frac{1}{2} \begin{pmatrix} \mat{1} & \mat{1}\\ \mat{1} & -\mat{1} \end{pmatrix} \begin{pmatrix} \mat{A} & -\mat{B}\\ \mat{B} & -\mat{A}^* \end{pmatrix} \begin{pmatrix} \mat{1} & \mat{1}\\ \mat{1} & -\mat{1} \end{pmatrix} \begin{pmatrix} \alpha\\ \beta \end{pmatrix} = \begin{pmatrix} \alpha\\ \beta \end{pmatrix} \omega,\\ \frac{1}{2} \begin{pmatrix} \mat{A}-\mat{A}^* & \mat{A} + \mat{A}^* + 2\mat{B}\\ \mat{A} + \mat{A}^* - 2\mat{B} & \mat{A} - \mat{A}^* \end{pmatrix} \begin{pmatrix} \alpha\\ \beta \end{pmatrix} = \begin{pmatrix} \alpha\\ \beta \end{pmatrix} \omega. \end{gather*}\]
\[\begin{gather*} \begin{pmatrix} \mat{A} & \I\mat{B}\\ \I\mat{B} & -\mat{A}^* \end{pmatrix} \begin{pmatrix} u_{+}\\ \I u_{-}^* \end{pmatrix} = \begin{pmatrix} u_{+}\\ \I u_{-}^* \end{pmatrix} \omega, \end{gather*}\]
I = np.eye(2)
X = np.array([[0, 1], [1, 0]])
Z = np.array([[1, 0], [0, -1]])
H = (X+Z)/np.sqrt(2)
P, M = (I+Z)/2, (I-Z)/2
A = np.array([[11+1j, 0.12], [0.12, 11-1j]])
#A = np.array([[11, 0.12], [0.12, 11]])
assert np.allclose(X@A@X, A.conj())
H@A@H

array([[ 1.11200000e+01-2.23711432e-17j, -2.42933074e-16+1.00000000e+00j],
       [-4.71276680e-16+1.00000000e+00j,  1.08800000e+01-2.23711432e-17j]])

Homogeneous Matter#

To see how this works, consider homogeneous matter where \(\psi_0\) has momentum \(\hbar k_0\). Translation invariance implies that we can use plane-waves \(u_{\pm} = \sqrt{n_0}\epsilon_{\pm}e^{\I (k_0 \pm k) x}\).

\[\begin{gather*} \psi = \sqrt{n_0}e^{\I (k_0 x - \mu t/\hbar)}\Bigl( 1 + \epsilon_+ e^{\I(kx - \omega t)} + \epsilon_- e^{-\I(kx - \omega t)} \Bigr),\\ = e^{\mu t/\hbar \I}\Bigl( \sqrt{n_0}e^{\I k_0 x} + \sqrt{n_0}\epsilon_+ e^{\I\bigl((k+k_0)x - \omega t\bigr)} + \sqrt{n_0}\epsilon_- e^{-\I\bigl((k-k_0)x - \omega t\bigr)} \Bigr). \end{gather*}\]

In this case \(\op{H}\) and \(\op{H}^*\) can be treated as a numbers

\[\begin{gather*} H = \frac{\hbar^2(k_0 + k)^2}{2m} + \mathcal{E}'(n_0) - \mu, \qquad H^* = \frac{\hbar^2(k_0 - k)^2}{2m} + \mathcal{E}'(n_0) - \mu. \end{gather*}\]

We start by expanding about the ground state with \(k_0=0\) so that \(H = H^*\). Taking the determinant, we have

\[\begin{gather*} H^2 + 2Hn_0\mathcal{E}''(n_0) - H \hbar (\omega - \omega^*) - \hbar^2\abs{\omega}^2 = 0. \end{gather*}\]

If \(\omega = \omega^*\), then we find the usual phonon dispersion relation:

\[\begin{gather*} \hbar\omega(k) = \pm\sqrt{H(H + 2n_0\mathcal{E}'')} = \pm\sqrt{H(H + 2mc^2)},\qquad c = \sqrt{\frac{n_0\mathcal{E}''(n_0)}{m}}, \end{gather*}\]

where (as we shall see shortly) \(c\) is the long wavelength speed of sound, expressed in terms of the compressibility \(\mathcal{E}''(n_0)\).

The chemical potential \(\mu\) cancels the constant pieces so that \(\op{H}\psi_0 = 0\) and:

\[\begin{gather*} H = \frac{\hbar^2 k^2}{2m}. \end{gather*}\]

Phonon Dispersion#

This gives the familiar dispersion of phonons about the ground state \(p_0 = 0\):

\[\begin{align*} \omega(k) &= \pm ck \sqrt{1 + \frac{\hbar^2k^2}{4m^2c^2}}, & c &= \sqrt{\frac{n_0\mathcal{E}''(n_0)}{m}},\\ E(p) = \hbar\omega &= \pm cp \sqrt{1 + \frac{p^2}{2m} \frac{1}{2mc^2}}. \end{align*}\]

This is linear in the long wavelength limit \(k \rightarrow 0\) with the advertised speed of sound \(c\). Note that when expressed in terms of the energy and momentum, this is a classical result without any explicit factors of \(\hbar\).

More generally, expanding about state \(\psi_0\) with momentum \(p_0\),

\[\begin{gather*} H = \overbrace{\frac{\hbar^2k^2}{2m}}^{H_0} + v_0k = \frac{\hbar^2(k+k_0)^2}{2m} - E_0,\\ H^* = \frac{\hbar^2k^2}{2m} - v_0k = \frac{\hbar^2(k-k_0)^2}{2m} - E_0,\\ \begin{vmatrix} H_0 - \hbar(\omega - v_0k) + n_0\mathcal{E}''(n_0) & n_0\mathcal{E}''(n_0)\\ n_0\mathcal{E}''(n_0) & H_0 + \hbar(\omega - v_0k)^* + n_0\mathcal{E}''(n_0) \end{vmatrix} = 0 \end{gather*}\]

where \(v_0 = p_0/m\) and \(E_0 = mv_0^2/2\). Taking the determinant, we have the slightly more cumbersome

\[\begin{gather*} H_0^2 + 2H_0n_0\mathcal{E}''(n_0) - H_0 \hbar (\omega - \omega^*) - \hbar^2\abs{(\omega-v_0k)}^2 = 0. \end{gather*}\]

And if \(\omega = \omega^*\), we find the boosted phonon dispersion:

\[\begin{gather*} \hbar\omega(k) - v_0k = \pm\sqrt{H_0(H_0 + 2n_0\mathcal{E}'')}. \end{gather*}\]
hbar = m = 1
c = 1
k = np.linspace(-2, 2, 101)
k0 = 0.5

def get_w(k,  k0=k0, c=c):
    v0 = hbar*k0/m
    p = hbar*k
    H0 = p**2/2/m
    w0 = np.sqrt(H0*(H0 + 2*m*c**2))/hbar
    return w0 + v0*k, -w0 + v0*k
    
fig, ax = plt.subplots()
for w in get_w(k):
    ax.plot(k, w)
../_images/f7fa4fb9863fc4a856caba1530759a2063e0f74a48a8f19fd61ee2973a2da16b.png

Counterflow Instability#

As an application, consider a miscible (\(g_{ab}^2 < g_{aa}g_{bb}\)) two-component superfluid. The homogeneous system becomes unstable if the two components flow past each other at a sufficiently large relative velocity \(\abs{v_a - v_b} > v_c\) – a counterflow instability. Here we use the Bogoliubov theory to see this instability.

The Hamiltonian is derived by varying the following self-energy:

\[\begin{gather*} \mathcal{E}(n_a, n_b) = \cdots + \frac{1}{2}\vect{n}^\dagger\mat{g}\vect{n}, \qquad \mat{g} = \begin{pmatrix} g_{aa} & g_{ab}\\ g_{ab} & g_{bb} \end{pmatrix}, \qquad \vect{n} = \begin{pmatrix} n_a\\ n_b \end{pmatrix}. \end{gather*}\]

We consider perturbing the following state

\[\begin{gather*} \ket{\psi_0} = \begin{pmatrix} \sqrt{n_a}e^{\I q x + \mu_a t/\I\hbar}\\ \sqrt{n_b}e^{-\I q x + \mu_b t/\I\hbar} \end{pmatrix}, \qquad \mu_a = \frac{\hbar^2q^2}{2m} + g_{aa}n_a + g_{ab}n_b, \qquad \mu_b = \frac{\hbar^2q^2}{2m} + g_{ab}n_a + g_{bb}n_b. \end{gather*}\]

To this we add a perturbation

\[\begin{gather*} \ket{u} = \begin{pmatrix} u_a\sqrt{n_a}e^{\I q x + \mu_a t/\I\hbar}\\ u_b\sqrt{n_b}e^{-\I q x + \mu_b t/\I\hbar} \end{pmatrix}e^{\I (kx - \omega t)} - \begin{pmatrix} v_a^*\sqrt{n_a}e^{\I q x + \mu_a t/\I\hbar}\\ v_b^*\sqrt{n_b}e^{-\I q x + \mu_b t/\I\hbar} \end{pmatrix}e^{-\I (kx - \omega t)},\\ \delta \vect{n} = \begin{pmatrix} n_a(u_a - v_a)\\ n_b(u_b - v_b) \end{pmatrix} e^{\I(kx-\omega t)} + \begin{pmatrix} n_a(u_a^* - v_a^*)\\ n_b(u_b^* - v_b^*) \end{pmatrix} e^{-\I(kx-\omega t)} + \cdots \end{gather*}\]

These give rise to the following Bogoliubov equations

\[\begin{gather*} \begin{pmatrix} \overbrace{K(q\mat{z} + k\mat{1}) - K(q\mat{z})} ^{\frac{(\hbar k)^2}{2m}\mat{1} + \hbar k v \mat{z}} + \mat{gn} - \hbar\omega \mat{1} & \mat{gn}\\ \mat{gn} & \underbrace{K(q\mat{z} - k\mat{1}) - K(q\mat{z})} _{\frac{(\hbar k)^2}{2m}\mat{1} - \hbar k v \mat{z}} + \mat{gn} + \hbar\omega \mat{1} \\ \end{pmatrix} \begin{pmatrix} u_{a}\\ u_{b}\\ -v_{a}\\ -v_{b} \end{pmatrix} = 0, \end{gather*}\]

where \(K(k) = \hbar^2k^2/2m\), and

\[\begin{gather*} \mat{gn} = \begin{pmatrix} g_{aa}n_{a} & g_{ab}n_b \\ g_{ab}n_{a} & g_{bb}n_b \end{pmatrix},\\ K(q\mat{z} \pm k\mat{1}) - K(q\mat{z}) = \frac{\hbar^2}{2m}k\Bigl(k\mat{1} \pm 2 q\mat{z}\Bigr) = \frac{(\hbar k)^2}{2m}\mat{1} \pm \hbar k v \mat{z},\\ v = \frac{\hbar q}{m} = v_a = -v_b. \end{gather*}\]

The structure of the equations is then

\[\begin{gather*} \mat{M} = \begin{pmatrix} \frac{\hbar^2 k^2}{2m}\mat{1} + \hbar k v \mat{z} + \mat{gn} & -\mat{gn}\\ \mat{gn} & -\frac{\hbar^2 k^2}{2m}\mat{1} + \hbar k v \mat{z} - \mat{gn}\\ \end{pmatrix},\\ (\mat{M} - \hbar \omega \mat{1}) \begin{pmatrix} u_{a}\\ u_{b}\\ v_{a}\\ v_{b} \end{pmatrix} = 0. \end{gather*}\]

Nonlinear Terms: \(\mat{gn}\).

Consider the equation of \(\psi_a\). It contains the non-linear term

\[\begin{gather*} \I\hbar \dot{\psi}_a = \cdots + (g_{aa}n_a + g_{ab}n_b)\psi_a. \end{gather*}\]

Expanding this will give three contributions – two from expanding the densities, and one from expanding \(\psi_a\) itself. The piece obtained by expanding \(\psi_a\) will cancel will the chemical potential \(\mu_a\), so we will only have the other two terms. Considering the equation with \(e^{-\I\omega t\) for \(u_a\), we will have which have the form

\[\begin{gather*} \hbar \omega \sqrt{n_a}u_a = \cdots + g_{aa}n_a (u_a - v_a)\sqrt{n_a} + g_{ab}n_b(u_b - v_b)\sqrt{n_a}. \end{gather*}\]
\[\begin{gather*} \begin{pmatrix} g_{aa}n_a & g_{ab}n_b & g_{aa}n_a & g_{ab}n_b\\ g_{aa}n_a & g_{ab}n_b & g_{aa}n_a & g_{ab}n_b\\ \vdots \end{pmatrix} \end{gather*}\]

Determinant of \(\mat{M}\) (incomplete).

To compute the determinant of \(\mat{M}\), we first write the tensor structure in terms of the 2×2 matrices

\[\begin{gather*} \mat{1} = \begin{pmatrix} 1 & 0\\ 0 & 1 \end{pmatrix}, \qquad \mat{z} = \begin{pmatrix} 1 & 0\\ 0 & -1 \end{pmatrix}, \qquad \mat{o} = \begin{pmatrix} 1 & -1\\ 1 & -1 \end{pmatrix}: \end{gather*}\]
\[\begin{gather*} \mat{M} = \frac{\hbar^2 k^2}{2m}\mat{z}\otimes\mat{1} + \hbar k v \mat{1}\otimes\mat{z} + \mat{o}\otimes\mat{gn}. \end{gather*}\]

Swapping the middle two rows and columns corresponds to reordering these tensors:

\[\begin{gather*} \mat{N} = \frac{\hbar^2 k^2}{2m}\mat{1}\otimes\mat{z} + \hbar k v \mat{z}\otimes\mat{1} + \mat{gn}\otimes\mat{o}, \qquad (\mat{N} - \hbar \omega \mat{1}) \begin{pmatrix} u_{a}\\ v_{a}\\ u_{b}\\ v_{b} \end{pmatrix} = 0. \end{gather*}\]

Recalling that \(v = v_{a} = -v_{b}\), the block structure is:

\[\begin{gather*} \mat{N} = \begin{pmatrix} \frac{\hbar^2 k^2}{2m}\mat{z} + \hbar k v_{a}\mat{1} + g_{aa}n_a\mat{o} & g_{ab}n_{b}\mat{o}\\ g_{ab}n_{a}\mat{o} & \frac{\hbar^2 k^2}{2m}\mat{z} + \hbar k v_{b}\mat{1} + g_{bb}n_{b}\mat{o}\\ \end{pmatrix}. \end{gather*}\]
\[\begin{gather*} \mat{z}\mat{o} = \begin{pmatrix} 1 & -1\\ -1 & 1 \end{pmatrix}, \qquad \mat{o}\mat{z} = \begin{pmatrix} 1 & 1\\ 1 & 1 \end{pmatrix}. \end{gather*}\]

Can one use properties of determinants like those in [Silvester, 2000] or [Abadir and Magnus, 2005]to simplify this?

\[\begin{gather*} (\hbar\omega_{i})^2 = \tilde{k}^2(2g_i n_i + \tilde{k}^2)\\ \det{\mat{M}} = -4g_ag_bn_an_b \gamma^2 \tilde{k}^4 + \Bigl((\hbar\omega - 2\tilde{k}\tilde{q})^2 - (\hbar \omega_a)^2\Bigr) \Bigl((\hbar\omega + 2\tilde{k}\tilde{q})^2 - (\hbar \omega_b)^2\Bigr)\\ = (\hbar\omega)^{4} - 2 \tilde{k}^{2}(g_an_a + g_bn_b + \tilde{k}^{2} + 4 \tilde{q}^{2})(\hbar\omega)^{2} + 8\tilde{k}^{3}\tilde{q}(g_bn_b - g_an_a)(\hbar\omega) +\\ + 4 g_a g_b n_a n_b \tilde{k}^{4} + 2 (g_a n_a + g_b n_b) \tilde{k}^{4} (\tilde{k}^{2} - 4\tilde{q}^{2}) + \tilde{k}^{8} - 8 \tilde{k}^{4} \tilde{q}^{2} (\tilde{k}^{2} - 2\tilde{q}^{2}) - 4 g_{ab}^{2}n_a n_b \tilde{k}^{4},\\ \end{gather*}\]

Solving this gives

\[\begin{gather*} \det{\mat{M}} = -4g_ag_bn_an_b \gamma^2 \frac{\hbar^4 k^4}{4m^2} + \Bigl((\hbar\omega + v_a k)^2 - (\hbar \omega_{B,a})^2\Bigr) \Bigl((\hbar\omega + v_b k)^2 - (\hbar \omega_{B,b})^2\Bigr),\\ (\hbar\omega_{B,i})^2 = \frac{\hbar^2}{2m}k^2\left(2g_i n_i + \frac{\hbar^2}{2m}k^2\right)\\ \end{gather*}\]

[Alperin and Timmermans, 2025] express this in terms of the individual speeds of sound \(c_{i}\) and modified healing lengths \(\xi_{i} = h_{i}/\sqrt{2}\):

\[\begin{gather*} \gamma^2 = \frac{g_{ab}^2}{g_{a}g_{b}}, \qquad c_{i}^2 = \frac{g_{i}n_i}{m_i},\qquad \xi_{i}^2 = \frac{\hbar^2}{4m_{i}^2 c_{i}^2}, \qquad \omega_{B,i}^2 = k^2 c_{i}^2(1 + \xi_{i}^2k^2),\\ c_a^2 c_b^2 \gamma^2 k^4 = \prod_{i\in\{a, b\}} \Bigl((\omega + \vect{v}_{i}\cdot\vect{k})^2 - \omega_{B,i}^2\Bigr). \end{gather*}\]

If we now assume that \(c_{a} = c_{b}\) and \(\xi_{a} = \xi_{b}\) so that \(\omega_{B,a} = \omega_{B,b} = \omega_{B}\), and take \(v = v_b = -v_a\), then the linear terms cancel and we have

\[\begin{gather*} c^4 \gamma^2 k^4 = \Bigl((\omega - vk)^2 - \omega_{B}^2\Bigr) \Bigl((\omega + vk)^2 - \omega_{B}^2\Bigr)=\\ = \Bigl(\omega^2 + v^2k^2 - \omega_{B}^2\Bigr)^2 - 4v^2k^2\omega^2. \end{gather*}\]

This has the solution

\[\begin{gather*} \omega^4 - 2(\omega_{B}^2 + v^2k^2)\omega^2 + \Bigl(v^2k^2 - \omega_{B}^2\Bigr)^2 -c^4\gamma^2k^4 = 0,\\ \omega^2 = \omega_{B}^2 + v^2k^2 \pm \sqrt{4\omega_{B}^2 v^2k^2 + \gamma^2 c^4k^4}. \end{gather*}\]

For repulsive interactions, \(\gamma^2 > 0\), hence the only way for an instability to manifest is if the r.h.s. is negative:

\[\begin{gather*} \Bigl(\omega_{B}^2(k) + v^2k^2\Bigr)^2 < 4\omega_{B}^2 v^2k^2 + \gamma^2 c^4k^4\\ \Bigl(\omega_{B}^2(k) - v^2k^2\Bigr)^2 < \gamma^2 c^4k^4\\ \Bigl(1 - \frac{v^2}{c^2} + \xi^2k^2\Bigr)^2 < \gamma^2. \end{gather*}\]

  1. For repulsive interactions, we can always induce an instability by setting

    \[\begin{gather*} v = c\sqrt{1 + \xi^2k^2} \end{gather*}\]

    at which point the l.h.s vanishes.

  2. This instability happens at \(k=0\) for \(v=c\): This is the usual Landau criterion.

  3. In the presence of interactions \(\gamma^2>0\), the actual \(k=0\) critical velocity is lower:

    \[\begin{gather*} v_{c} = c\sqrt{1 - \gamma} \end{gather*}\]

Hide code cell source

 
fig, axs = plt.subplots(3, 1, figsize=(3,6))
kmax = 6.5
k = np.linspace(0, kmax, 500)
gamma2 = 1.0
m = 1/2
c = xi = 1.0
hbar = 2*m*c*xi
gn = m*c**2
w_unit = 2*m*c**2/hbar
for v, ax in zip([1, 2, 5], axs):
    q = np.sqrt(2*m) * v/2
    D = 2*k**2*np.sqrt((2*gn + k**2)*4*q**2 + gamma2*gn**2)
    P = k**2*(2*gn + k**2 + 4*q**2)
    ax.plot(k*xi, np.sqrt(P+D) / w_unit, 'g-')
    with np.errstate(invalid='ignore'):
        ax.plot(k*xi, np.sqrt(P-D) / w_unit, 'b-')
        ax.plot(k*xi, np.sqrt(-(P-D)) / w_unit, 'b:')
    ax.set(ylim=(0, 10), xticks=np.arange(7), ylabel="$\omega(k)$ [$2mc^2/\hbar$]")
    ax.text(4.5, 8, f"${v=}c$")
ax.set(xlabel=r"$k$ [$\xi$]");
../_images/193bb59fc36c7f53ae579d663681b346f9a06a657d9f529e879572335a4ac23d.png

  1. In the first version of [Alperin and Timmermans, 2025] they claimed that allowing \(\gamma^2 < 0\) – i.e. one of the components is attractive – one can realize a roton minimum in the quasi-particle dispersion relationship, however, the approximation \(c_{a}^2 = c_{b}^2\) and \(\xi_a^2 = \xi_b\)^2 cannot be realized if \(\gamma^2 < 0\). In particular, if we take \(g_{a} <0\) to be attractive, then

    \[\begin{gather*} c_{a}^2 = -c^2, \qquad \xi_a^2 = -\xi^2,\\ \omega_{B,a}^2 = -c^2k^2(1-\xi^2k^2), \qquad \omega_{B,b}^2 = c^2k^2(1+\xi^2k^2). \end{gather*}\]

    The dispersion relation now satisfies

    \[\begin{align*} \abs{\gamma^2}c^4 k^4 &= \bigl((\omega - vk)^2 - \omega_{B,a}^2\bigr) \bigl((\omega + vk)^2 - \omega_{B,b}^2\bigr),\\ &= \bigl((\omega - vk)^2 + c^2k^2(1-\xi^2k^2)\bigr) \bigl((\omega + vk)^2 - c^2k^2(1+\xi^2k^2)\bigr),\\ &= (\omega^2 + v^2k^2 - c^2\xi^2k^4)^2 - (c^2k^2 - 2\omega v k)^2. \end{align*}\]

    This has the form

    \[\begin{gather*} \omega^4 + p\omega^2 + q \omega + r = 0, \qquad p = -2(c^2\xi^2k^4 + v^2k^2), \qquad q = 4vc^2k^3, \\ r = (v^2k^2 - c^2\xi^2k^4)^2 - c^4k^4 = c^4\xi^4k^8 + (v^4-c^4 - 2v^2c^4\xi^2)k^4 \end{gather*}\]

    In the limit of large \(k\) we have

    \[\begin{gather*} p \rightarrow -2 c^2\xi^2k^4, \qquad q = 4vc^2k^3, \qquad r \rightarrow c^4\xi^4 k^8,\\ \Delta \rightarrow -c^{12}k^{12}( 512 k^{12}\xi^{12} + 4096k^{6}v^{2}\xi^{6} + 6912v^{4}) < 0\\ \end{gather*}\]

    Hence, there are always imaginary frequencies, and the state is unstable.

Old discussion

If we assume \(g_an_a = g_bn_b = gn\), then the linear term vanishes and we have a quadratic equation in terms of \((\hbar\omega)^2\) that can be easily solved to obtain the following condition for an instability:

\[\begin{gather*} (2gn + \tilde{k}^2 - 4\tilde{q}^2)^2 < 4 g_{ab}^{2}n_a n_b. \end{gather*}\]

We note a couple of points about this result:

  1. The r.h.s. is manifestly non-negative, so there is always an instability since we can set the l.h.s. to zero by inducing a counterflow with momentum

    \[\begin{gather*} \tilde{q} = \sqrt{\frac{gn}{2} + \frac{\tilde{k}^2}{4}}. \end{gather*}\]

    The physical interpretation is simple: The minimum such \(\tilde{q}\) for \(\tilde{k} = 0\) is the speed of sound

    \[\begin{gather*} \frac{\hbar q}{m} = \sqrt{\frac{gn}{m}} = c. \end{gather*}\]

    Thus, this simply reduces to the Landau criterion for the stability of the superfluid.

  2. The actual critical velocity is lower due to the interactions \(g_{ab}\):

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

  3. The quasiparticle dispersion has the form:

    \[\begin{gather*} \gamma^2 = \frac{g_{ab}^2}{g_ag_b}, \qquad g_{ab}^2 n_a n_b = (\gamma gn)^2,\\ (\hbar\omega)^{2} = \tilde{k}^{2}(2gn + \tilde{k}^{2} + 4 \tilde{q}^{2}) \pm 2\tilde{k}^2\sqrt{(2 gn + \tilde{k}^2)4\tilde{q}^2 + \gamma^2(gn)^2}. \end{gather*} \]

    In terms of \(\gamma\), the condition of instability is

    \[\begin{gather*} \Abs{1 + \frac{\tilde{k}^2 - 4\tilde{q}^2}{2gn}} < \abs{\gamma} \end{gather*}\]

Details

The standard GPE dispersion relationship is

\[\begin{gather*} \omega^2 = k^2 c^2\Bigl(1 + \frac{\hbar^2 k^2}{2m}\frac{1}{2mc^2}\Bigr) \end{gather*}\]

which they express below Eq. (7) as \(\omega^2_{B,i} = k^2c_{i}^2(1+\xi_i^2 k^2)\), which requires

\[\begin{gather*} c = \sqrt{\frac{gn}{m}}, \qquad \xi = \frac{\hbar}{2mc}. \end{gather*}\]

This differs from the usual healing length \(h = \sqrt{2} \xi\) where

\[\begin{gather*} \frac{\hbar^2}{2m h^2} = gn. \end{gather*}\]
\[\begin{gather*} (\hbar\omega)^{4} - 2\tilde{k}^{2}(2gn + \tilde{k}^{2} + 4 \tilde{q}^{2})(\hbar\omega)^{2}+\\ + 4 (gn)^2 \tilde{k}^{4} + 4 gn \tilde{k}^{4} (\tilde{k}^{2} - 4\tilde{q}^{2}) + \tilde{k}^{8} - 8 \tilde{k}^{4} \tilde{q}^{2} (\tilde{k}^{2} - 2\tilde{q}^{2}) - 4 g_{ab}^{2}n_a n_b \tilde{k}^{4} = 0,\\ (\hbar\omega)^{2} = \tilde{k}^{2}(2gn + \tilde{k}^{2} + 4 \tilde{q}^{2}) \pm 2\tilde{k}^2\sqrt{8 gn \tilde{q}^2 + 4 \tilde{k}^2 \tilde{q}^2 + g_{ab}^{2}n_a n_b}. \end{gather*}\]

If the interactions are repulsive, then the discriminant is positive, so we only have an instability if the \(\omega^2 < 0 \) is negative:

\[\begin{gather*} (2gn + \tilde{k}^{2} + 4 \tilde{q}^{2}) < 2\sqrt{8 gn \tilde{q}^2 + 4 \tilde{k}^2 \tilde{q}^2 + g_{ab}^{2}n_a n_b},\\ (2gn + \tilde{k}^2 - 4\tilde{q}^2)^2 < 4 g_{ab}^{2}n_a n_b,\\ \end{gather*}\]

To compare with :

\[\begin{align*} \vect{V}_{j} = \frac{\hbar k_j}{m} &\equiv \pm \frac{\hbar q}{m}\\ \vect{V}_R &\equiv \frac{2\hbar q}{m},\\ \vect{P} &\equiv \hbar k,\\ \gamma &\equiv \frac{g_{ab}}{g} = \frac{g_{ab}}{\sqrt{g_ag_b}}. \end{align*}\]
\[\begin{gather*} 4\left(\frac{\hbar^2k^2}{4m^2} + \frac{gn}{m}(1-\abs{\gamma}\right) < \frac{4\hbar^4 k^2q^2}{m^2} < 4\left(\frac{\hbar^2k^2}{4m^2} + \frac{gn}{m}(1+\abs{\gamma}\right)\\ \left(\frac{\hbar^4 k^2q^2}{m^2} - \frac{\hbar^2k^2}{4m^2} - \frac{gn}{m}\right)^2 < \frac{g_{ab}^2n_an_b}{m^2}. \end{gather*}\]

For fixed \(\tilde{q}\) we can find the wavevector \(\tilde{k}\) that maximizes the instability by minimizing the (negative) \((\hbar\omega)^2\):

\[\begin{gather*} (2gn + 2\tilde{k}^2 + 4\tilde{q}^2) - 2\sqrt{\cdots} - \frac{2\tilde{k}^2}{\sqrt{\cdots}}(4\tilde{q}^2) = 0\\ (2gn + 2\tilde{k}^2 + 4\tilde{q}^2)\sqrt{\cdots} = 16 \tilde{q}^2(gn + \tilde{k}^2) + 2g_{ab}^{2}n_a n_b\\ \end{gather*}\]

To simplify these, we define

\[\begin{gather*} \gamma = \frac{g_{ab}}{\sqrt{g_ag_b}}, \qquad g_{ab}^2 n_a n_b = (\gamma gn)^2,\\ (\hbar\omega)^{2} = \tilde{k}^{2}(2gn + \tilde{k}^{2} + 4 \tilde{q}^{2}) \pm 2\tilde{k}^2\sqrt{(2 gn + \tilde{k}^2)4\tilde{q}^2 + \gamma^2(gn)^2}. \end{gather*}\]
\[\begin{gather*} w^2 = ak + k^2 - 2k\sqrt{qk + b}, \qquad a = 2gn + 4\tilde{q}^2, \qquad q = 4\tilde{q}^2, \qquad b = 8gn\tilde{q}^2 + g_{ab}n_{a}n_{b},\\ (a + 2k) - 2\sqrt{qk + b} - \frac{k}{\sqrt{qk + b}}q = 0\\ (a + 2k)\sqrt{qk + b} = 3qk + 2b\\ (a + 2k)^2(qk + b) = (3qk + 2b)^2\\ 4ak^3 + (4b + 2aq - 9q^2)k^2 + (2ab + a^2q - 6qb)k + ba^2 - 4b^2 = 0\\ \end{gather*}\]

This is a rather nasty cubic…

Hide code cell source

 
import numpy as np
k, q = 1.2, 3.4
ga, gb, gab = 1.1, 2.2, 3.3
gn = 1.3
na, nb = gn/ga, gn/gb
m, hbar = 1.2, 1.4
Z = np.array([[1, 0], [0, -1]])
O = np.array([[1, 0], [0, 1]])
Gn = np.array([
    [ga*na, gab*nb],
    [gab*na, gb*nb]])

def K(k):
    return (hbar*k)**2/2/m
    
A = K(q*Z + k*O) - K(q*Z) + Gn
B = -Gn
C = Gn
D = -(K(q*Z - k*O) - K(q*Z) + Gn)
M = np.bmat([[A, B], [C, D]])
hw = np.linalg.eigvals(M)
gp, gm = (ga+gb)/2, (ga-gb)/2
n_p, n_m = na+nb, na-nb
a, c = (hbar*k)**2/2/m, hbar**2 * k * q / m
nesc = 2*gn*a + a**2 + c**2
desc = 2*np.sqrt((a*c)**2 + a**2*gab**2*na*nb + 2*c**2*gn*a)
hw_ = np.sqrt([nesc + desc,  nesc - desc])
hw_ = np.sort(np.concatenate([hw_, -hw_]))
assert np.allclose(np.sort(hw), hw_)

Hide code cell content

 
import sympy
from sympy import var, Matrix, Identity, latex, det
k, q, w, ga, gb, gab, na, nb = var('k, q, w, g_a, g_b, g_ab, n_a, n_b')
M = Matrix([
  [k*(k+2*q)+ga*na, gab*nb, -ga*na, -gab*nb],
  [gab*na, k*(k-2*q)+gb*nb, -gab*na, -gb*nb],
  [ga*na, gab*nb, -(k*(k-2*q)+ga*na), -gab*nb],
  [gab*na, gb*nb, -gab*na, -(k*(k+2*q)+gb*nb)]])
print(latex(det(M-w*sympy.eye(4)).simplify().collect(w)))

4 g_{a} g_{b} k^{4} n_{a} n_{b} + 2 g_{a} k^{6} n_{a} - 8 g_{a} k^{4} n_{a} q^{2} - 4 g_{ab}^{2} k^{4} n_{a} n_{b} + 2 g_{b} k^{6} n_{b} - 8 g_{b} k^{4} n_{b} q^{2} + k^{8} - 8 k^{6} q^{2} + 16 k^{4} q^{4} + w^{4} + w^{2} \left(- 2 g_{a} k^{2} n_{a} - 2 g_{b} k^{2} n_{b} - 2 k^{4} - 8 k^{2} q^{2}\right) + w \left(- 8 g_{a} k^{3} n_{a} q + 8 g_{b} k^{3} n_{b} q\right)

Dynamics Response#

We now consider driving the system where \(V(x, t) = V_0 \cos(\omega t)\cos(kx) = \tfrac{1}{2}V_0(e^{\omega t/\I} + e^{-\omega t/\I})\cos(kx)\) is small. This adds an inhomogeneous term which we treat as the same order a \(\epsilon_{\pm}\):

\[\begin{gather*} \begin{pmatrix} \op{H} + n_0\mathcal{E}'' - \hbar \omega & n_0\mathcal{E}''\\ n_0\mathcal{E}'' & \op{H}^* + n_0\mathcal{E}'' + \hbar \omega^* \end{pmatrix} \begin{pmatrix} \epsilon_+(x) \\ \epsilon_-^*(x) \end{pmatrix} = -\frac{V_0}{2}\cos(kx) \begin{pmatrix} \psi_0\\ \psi_0^* \end{pmatrix}, \end{gather*}\]

the formal solution of which can be found by inverting the matrix.

Homogeneous Matter#

The full linear response follows from inverting the matrix, but we need to be a bit careful. We start by setting \(k_0 = 0\) so that parity is conserved. Then we expand:

\[\begin{gather*} \psi = \sqrt{n_0}\Bigl( 1 + (\epsilon_{+}e^{-\I\omega t} + \epsilon_{-}e^{\I\omega t})\cos(kx)\Bigr),\\ n = n_0 + \sqrt{n_0}\cos(kx)2\Re\Bigl( \epsilon_{+}e^{-\I\omega t} + \epsilon_{-}e^{\I\omega t}\Bigr) \end{gather*}\]

Collecting all terms, we have two independent sets of equations:

\[\begin{gather*} \begin{pmatrix} A - \hbar\omega e^{\I\eta} & B\\ B & A + \hbar\omega e^{-\I\eta} \end{pmatrix} \begin{pmatrix} \epsilon_{+}\\ \epsilon_{-}^* \end{pmatrix} = -\frac{V_0\sqrt{n_0}}{4} \begin{pmatrix} 1\\ 1 \end{pmatrix},\\ A = H_0 + B, \qquad B = n_0\mathcal{E}''(n_0). \end{gather*}\]
\[\begin{gather*} \begin{pmatrix} \epsilon_{+}\\ \epsilon_{-}^* \end{pmatrix} = -\frac{V_0\sqrt{n_0}}{4\Bigl( (A - \hbar\omega e^{\I\eta})(A + \hbar\omega e^{-\I\eta}) - B^2 \Bigr)} \begin{pmatrix} A + \hbar\omega e^{-\I\eta} & -B\\ -B & A - \hbar\omega e^{\I\eta} \end{pmatrix} \begin{pmatrix} 1\\ 1 \end{pmatrix},\\ = -\frac{V_0\sqrt{n_0}}{4\Bigl( A^2 - \hbar^2\omega^2 - B^2 - 2\I A\hbar\omega \sin(\eta) \Bigr)} \begin{pmatrix} A + \hbar\omega e^{-\I\eta} -B\\ A - \hbar\omega e^{\I\eta} - B \end{pmatrix},\\ \epsilon_{+} = -\frac{V_0\sqrt{n_0}(A + \hbar\omega e^{-\I\eta} -B)}{4\Bigl( A^2 - \hbar^2\omega^2 - B^2 - 2\I A\hbar\omega \sin(\eta)\Bigr)}\\ \epsilon_{-} = -\frac{V_0\sqrt{n_0}(A - \hbar\omega e^{\I\eta} - B)}{4\Bigl( A^2 - \hbar^2\omega^2 - B^2 - 2\I A\hbar\omega \sin(\eta)\Bigr)},\\ \epsilon_{+} = -\frac{V_0\sqrt{n_0}(A + \hbar\omega e^{-\I\eta} -B)}{4\Bigl( A^2 - \hbar^2\omega^2 - B^2 - 2\I A\hbar\omega \sin(\eta)\Bigr)}\\ \epsilon_{-} = -\frac{V_0\sqrt{n_0}(A - \hbar\omega e^{\I\eta} - B)}{4\Bigl( A^2 - \hbar^2\omega^2 - B^2 - 2\I A\hbar\omega \sin(\eta)\Bigr)} \end{gather*}\]
\[\begin{split} \begin{pmatrix} \op{H} + g n_0 & g\psi_0^2\\ g \bar{\psi}_0^2 & \bar{\op{H}} + g n_0 \end{pmatrix}\cdot \begin{pmatrix} u(x)\\ v(x) \end{pmatrix} = \omega \begin{pmatrix} \mat{1} \\ & -\mat{1} \end{pmatrix}\cdot \begin{pmatrix} u(x)\\ v(x) \end{pmatrix}, \end{split}\]

where \(n_0 = \abs{\psi_0}^2\) and \(\op{H} = -\hbar^2\nabla^2/2m + gn_0 + \op{V}_{\text{ext}}\) is the single-particle Hamiltonian for the ground state. To solve this numerically, we write this as \(\mat{A}\cdot\vect{q} = \omega \mat{B}\cdot\vect{q}\).

These matrices have the following properties: \(\mat{A} = \mat{A}^\dagger\) and \(\mat{B} = \mat{B}^\dagger\) are Hermitian, and the matrix \(\mat{C} = \mat{C}^{-1} = \bigl(\begin{smallmatrix}&\mat{1}\\\mat{1}\end{smallmatrix}\bigr)\) conjugates \(\mat{A}\):

\[\begin{split} \mat{C}\cdot\mat{A}\cdot\mat{C} = \bar{\mat{A}}\\ \mat{C}\cdot\mat{B}\cdot\mat{C} = -\bar{\mat{B}}. $$. Thus, if we have one eigenvalue $\omega_{+}$ and eigenvector $\vect{q}_{+}$, then, we must have another pair $\omega_{-} = \bar{\omega}_{+}$ and $\vect{q}_{-} = \mat{C}\cdot\bar{\vect{q}}_{+}$:\end{split}\]

\mat{A}\vect{q}{+} = \omega+ \mat{B}\vect{q}{+}\ \bar{\mat{A}}\mat{C}\vect{q}{+} = \omega_+ \bar{\mat{B}}\mat{C}\vect{q}_{+}. $$

Furthermore, if \(\psi_0\) is a stationary state, then \(\op{H}\psi_0 = \mu\psi_0\).

One has two choices: solve the non-symmetric eigenvalue problem \((\mat{B}^{-1}\cdot\mat{A})\cdot\vect{q} = \omega\vect{q}\), or try to massage this into a form where \(\mat{B}\) is positive definite.

\[ (\mat{A} + 2\mat{1})\cdot\vect{q} = (\omega\mat{B} + 2\mat{I})\cdot\vect{q} \]

Fixed Particle Number#

Let \(\psi = \psi_0 + \d{\psi}\). The change in density is:

\[ \dot{n} = \dot{\psi}^\dagger \psi_0 + \psi_0^\dagger \dot{\psi} + \dot{\psi}^\dagger \dot{\psi} \]
\[ n - n_0 = e^{-\I\omega t}(u \psi_0^* + v\psi_0) + e^{\I\omega t}(u^* \psi_0 + v^*\psi_0^*) + e^{-2\I\omega t} u^*v^* + e^{2\I\omega t} uv + (u^*u + v^*v) \]

particle number is:

\[\begin{split} \d{N} = \int\left( \d{\psi}^\dagger \psi_0 + \psi_0^\dagger \d{\psi} + \d{\psi}^\dagger \d{\psi} \right)\d{x}\\ = e^{-\I\omega t}\int(u \psi_0^* + v\psi_0)\d{x} + e^{\I\omega t}\int(u^* \psi_0 + v^*\psi_0^*)\d{x} + e^{-2\I\omega t} \int u^*v^*\d{x} + e^{2\I\omega t}\int uv\d{x} + \int(u^*u + v^*v)\d{x} \end{split}\]
from pytimeode.evolvers import EvolverABM, EvolverSplit
from mmfutils.contexts import NoInterrupt

[I 03:37:54 numexpr.utils] NumExpr defaulting to 2 threads.

import bec;reload(bec)
from bec import State, u
s = State(N=1e1, Nxyz=(2**7,), Lxyz=(20*u.micron,))
s.pre_evolve_hook()
s._N = 1e1
s.normalize()
s.plot()
s.cooling_phase = 1j

e = EvolverSplit(s, dt=1.0, normalize=True)

E_tol = 1e-8
E1 = e.y.get_energy()
E0 = E1 + 2*E_tol

log = True

with NoInterrupt(ignore=True) as interrupted:
    while not interrupted and abs(E1 - E0) > E_tol:
        e.evolve(100)
        E0, E1 = E1, e.y.get_energy()
        plt.clf()
        e.y.plot(log=log)
        display(plt.gcf())
        clear_output(wait=True)

e.dt = 0.01
E0 = E1 + 2*E_tol        
with NoInterrupt(ignore=True) as interrupted:
    while not interrupted and abs(E1 - E0) > E_tol/100:
        print abs(E1 - E0)
        e.evolve(100)
        E0, E1 = E1, e.y.get_energy()
        plt.clf()
        e.y.plot(log=log)
        display(plt.gcf())
        clear_output(wait=True)

fact = 100.0
e = EvolverABM(e.get_y(), dt=1.0/fact)
E0 = E1 + 2*E_tol        
with NoInterrupt(ignore=True) as interrupted:
    while not interrupted and abs(E1 - E0) > E_tol/fact/10:
        print abs(E1 - E0)
        e.evolve(100)
        E0, E1 = E1, e.y.get_energy()
        plt.clf()
        e.y.plot(log=log)
        display(plt.gcf())
        clear_output(wait=True) 
s = e.get_y()

  Cell In[8], line 31
    print abs(E1 - E0)
    ^
SyntaxError: Missing parentheses in call to 'print'. Did you mean print(...)?

Normal Modes#

Here we consider the fluctuations about a stationary state \(\psi_0\) of the GPE:

\[ \psi = \psi_0 + u(x)e^{\I\omega t} + v^*(x) e^{-\I\omega t}. \]

This gives the following generalized eigenvalue problem for the modes:

\[\begin{split} \begin{pmatrix} \op{H} + g n_0 & g\psi_0^2\\ g \bar{\psi}_0^2 & \bar{\op{H}} + g n_0 \end{pmatrix}\cdot \begin{pmatrix} u(x)\\ v(x) \end{pmatrix} = \omega \begin{pmatrix} \mat{1} \\ & -\mat{1} \end{pmatrix}\cdot \begin{pmatrix} u(x)\\ v(x) \end{pmatrix}, \end{split}\]

where \(n_0 = \abs{\psi_0}^2\) and \(\op{H} = -\hbar^2\nabla^2/2m + gn_0 + \op{V}_{\text{ext}}\) is the single-particle Hamiltonian for the ground state. To solve this numerically, we write this as \(\mat{A}\cdot\vect{q} = \omega \mat{B}\cdot\vect{q}\).

These matrices have the following properties: \(\mat{A} = \mat{A}^\dagger\) and \(\mat{B} = \mat{B}^\dagger\) are Hermitian, and the matrix \(\mat{C} = \mat{C}^{-1} = \bigl(\begin{smallmatrix}&\mat{1}\\\mat{1}\end{smallmatrix}\bigr)\) conjugates \(\mat{A}\):

\[\begin{split} \mat{C}\cdot\mat{A}\cdot\mat{C} = \bar{\mat{A}}\\ \mat{C}\cdot\mat{B}\cdot\mat{C} = -\bar{\mat{B}}. $$. Thus, if we have one eigenvalue $\omega_{+}$ and eigenvector $\vect{q}_{+}$, then, we must have another pair $\omega_{-} = \bar{\omega}_{+}$ and $\vect{q}_{-} = \mat{C}\cdot\bar{\vect{q}}_{+}$:\end{split}\]

\mat{A}\vect{q}{+} = \omega+ \mat{B}\vect{q}{+}\ \bar{\mat{A}}\mat{C}\vect{q}{+} = \omega_+ \bar{\mat{B}}\mat{C}\vect{q}_{+}. $$

Furthermore, if \(\psi_0\) is a stationary state, then \(\op{H}\psi_0 = \mu\psi_0\).

One has two choices: solve the non-symmetric eigenvalue problem \((\mat{B}^{-1}\cdot\mat{A})\cdot\vect{q} = \omega\vect{q}\), or try to massage this into a form where \(\mat{B}\) is positive definite.

\[ (\mat{A} + 2\mat{1})\cdot\vect{q} = (\omega\mat{B} + 2\mat{I})\cdot\vect{q} \]

Fixed Particle Number#

Let \(\psi = \psi_0 + \d{\psi}\). The change in density is:

\[ \dot{n} = \dot{\psi}^\dagger \psi_0 + \psi_0^\dagger \dot{\psi} + \dot{\psi}^\dagger \dot{\psi} \]
\[ n - n_0 = e^{-\I\omega t}(u \psi_0^* + v\psi_0) + e^{\I\omega t}(u^* \psi_0 + v^*\psi_0^*) + e^{-2\I\omega t} u^*v^* + e^{2\I\omega t} uv + (u^*u + v^*v) \]

particle number is:

\[\begin{split} \d{N} = \int\left( \d{\psi}^\dagger \psi_0 + \psi_0^\dagger \d{\psi} + \d{\psi}^\dagger \d{\psi} \right)\d{x}\\ = e^{-\I\omega t}\int(u \psi_0^* + v\psi_0)\d{x} + e^{\I\omega t}\int(u^* \psi_0 + v^*\psi_0^*)\d{x} + e^{-2\I\omega t} \int u^*v^*\d{x} + e^{2\I\omega t}\int uv\d{x} + \int(u^*u + v^*v)\d{x} \end{split}\]

To linear order, conservation of particle number thus implies that

\[ \int(u \psi_0^\dagger + v\psi_0)\d{x} = 0. \]
import scipy.linalg
sp = scipy
H = s.get_H()
n_0 = s.get_density().ravel()
psi_0 = 1*s[...].ravel()
A = np.bmat(
    [[H + np.diag(s.g*n_0), np.diag(s.g*psi_0**2)],
     [np.diag(s.g*psi_0.conj()**2), H.conj() + np.diag(s.g*n_0)]])
A = np.asarray(A)
i = np.eye(len(H))
z = np.zeros_like(H)
B = np.asarray(np.bmat([[i, z], [z, -i]]))
C = np.asarray(np.bmat([[z, i], [i, z]]))
I = np.eye(2*len(H))

assert np.allclose(C.dot(C), I)
assert np.allclose(C.dot(A).dot(C), A.conj())
ws, uvs = sp.linalg.eig(np.linalg.inv(B).dot(A))
inds = np.argsort(abs(ws.real))
ws = ws[inds]
uvs = uvs[:, inds]
us, vs = uvs.reshape((2, len(H), 2*len(H)))

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[9], line 3
      1 import scipy.linalg
      2 sp = scipy
----> 3 H = s.get_H()
      4 n_0 = s.get_density().ravel()
      5 psi_0 = 1*s[...].ravel()
      6 A = np.bmat(

NameError: name 's' is not defined

ws[:10], s.ws[0]

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[10], line 1
----> 1 ws[:10], s.ws[0]

NameError: name 'ws' is not defined

s.cooling_phase = 1.0
s.pre_evolve_hook()
s.cooling_phase = 1.0
s.t = 0.0

mode = 4
d = 0.001
u_, v_ = us[:, mode], vs[:, mode]
s[...] = psi_0 + d*(u_ + v_.conj())
w = ws[mode].real
T = 2*np.pi / w 
dt = T/100/100
e = EvolverSplit(s, dt=dt, normalize=False)
log = False
x = s.xyz[0]
with NoInterrupt(ignore=True) as interrupted:
    while not interrupted:
        e.evolve(100)
        plt.clf()
        #e.y.plot(log=log)
        t = e.y.t
        psi = psi_0 + d*(u_*np.exp(1j*w*t) + v_.conj()*np.exp(-1j*w*t))
        plt.plot(x, abs(e.y[...])**2 - abs(psi_0)**2)
        plt.plot(x, abs(psi)**2 - abs(psi_0)**2)

        #plt.ylim(0, 2)
        plt.xlim(-6, 6)
        plt.title("N = {}".format(e.y.get_N()))
        display(plt.gcf())
        clear_output(wait=True)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[11], line 1
----> 1 s.cooling_phase = 1.0
      2 s.pre_evolve_hook()
      3 s.cooling_phase = 1.0
      4 s.t = 0.0

NameError: name 's' is not defined

np.allclose(np.linalg.inv(B).dot(A).dot(uvs), ws[None,:]*uvs)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[12], line 1
----> 1 np.allclose(np.linalg.inv(B).dot(A).dot(uvs), ws[None,:]*uvs)

NameError: name 'uvs' is not defined

Hydrodynamics#

Consider the following GPE-like theory

\[\begin{gather*} \I\hbar e^{\I\eta}\dot{\psi} = \left( \frac{-\hbar^2\nabla^2}{2m} + \mathcal{E}'(n) + \beta \dot{n} \right)\psi. \end{gather*}\]

Now apply the Madelung transform \(\psi = \sqrt{n}e^{\I\phi}\). We have the following derivatives (similarly with spatial derivatives):

\[\begin{gather*} \dot{\psi} = \left(\frac{\dot{n}}{2n} + \I\dot{\phi}\right)\psi, \\ \psi'' = \left( \frac{n''}{2n} - \frac{(n')^2}{2n^2} + \I\phi'' \right)\psi + \left(\frac{n'}{2n} + \I\phi'\right)^2\psi = \Biggl( \overbrace{\frac{\sqrt{n}''}{\sqrt{n}}} ^{\frac{2nn'' - (n')^2}{4n^2}} - (\phi')^2 + \I\frac{(n\phi')'}{n} \Biggr)\psi,\\ \I\hbar e^{\I\eta}\left(\frac{\dot{n}}{2n} + \I\dot{\phi}\right) = \frac{-\hbar^2}{2m} \left( \frac{\sqrt{n}''}{\sqrt{n}} - (\phi')^2 + \I\frac{(n\phi')'}{n} \right) + \mathcal{E}'(n) + \beta\dot{n}. \end{gather*}\]

Consider first \(\eta=0\). Collecting the imaginary and real parts, we have

\[\begin{gather*} \dot{n} = -\vect{\nabla}\cdot\underbrace{ \Bigl(n\overbrace{\frac{\hbar \vect{\nabla}\phi}{m}}^{\vect{v}}\Bigr)}_{\vect{j}},\\ - \hbar\dot{\phi} = \frac{-\hbar^2}{2m} \left( \frac{\nabla^2\sqrt{n}}{\sqrt{n}} - (\vect{\nabla}\phi)^2 \right) + \mathcal{E}'(n) + \beta\dot{n}. \end{gather*}\]

Introducing the velocity \(\vect{v} = \hbar \vect{\nabla}\phi / m\) and current \(\vect{j} = n\vect{v}\), we have the traditional conservation of particle number and Euler equations:

\[\begin{gather*} \dot{n} + \vect{\nabla}\cdot\vect{j} = 0, \qquad \hbar\dot{\phi} + \frac{1}{2} m v^2 = \frac{\hbar^2}{2m} \frac{\nabla^2\sqrt{n}}{\sqrt{n}} - \mathcal{E}'(n) - \beta\dot{n}. \end{gather*}\]

One typically takes the gradient of the last equation (divided by \(m\)) to obtain:

\[\begin{gather*} \frac{\hbar}{m}\vect{\nabla}\dot{\phi} + \vect{\nabla}\frac{v^2}{2} = \dot{\vect{v}} + \vect{\nabla}\frac{v^2}{2} = \underbrace{ \frac{-\vect{\nabla}}{m} \overbrace{ \left( \frac{-\hbar^2}{2m} \frac{\nabla^2\sqrt{n}}{\sqrt{n}} + \mathcal{E}'(n) + \beta\dot{n} \right) }^{V} }_{\vect{F}/m = -\vect{\nabla}V/m} \end{gather*}\]

Convective Derivative

We have one complication here – the typical convective derivative is:

\[\begin{gather*} D_t\vect{A} = \dot{\vect{A}} + \vect{v}\cdot\vect{\nabla}\vect{A}, \end{gather*}\]

but the second term is not generally equivalent to \(\vect{\nabla}v^2/2\). Specifically, in index form, we have

\[\begin{gather*} \dot{v}_i + \nabla_{i}\frac{v^2}{2} = \dot{v}_i + v_j\nabla_{i}v_j = D_t v_i + v_j\nabla_{i}v_j - v_j\nabla_{j}v_i = D_t v_i + \underbrace{v_j(\nabla_{i}v_j - \nabla_{j}v_i)}_{\propto \vect{v}\cdot(\vect{\nabla}\times\vect{v})}. \end{gather*}\]

This correction, however, vanishes for potential flow \(v \propto \vect{\nabla}\phi\) since the divergence of a curl vanishes. With indices, this is obvious:

\[\begin{gather*} \nabla_{i}v_j - \nabla_{j}v_i \propto \nabla_{i}\nabla_{j}\phi - \nabla_{j}\nabla_{i}\phi = 0. \end{gather*}\]

Summarizing, we have

\[\begin{gather*} \dot{n} + \vect{\nabla}\cdot\vect{j} = 0,\qquad mD_t\vect{v} = -\vect{\nabla}V,\\ \end{gather*}\]

where

\[\begin{gather*} D_t A = \dot{A} + (\vect{v}\cdot\vect{\nabla})A, \qquad V = \frac{-\hbar^2}{2m}\frac{\nabla^2\sqrt{n}}{\sqrt{n}} + \mathcal{E}'(n) + \beta\dot{n}. \end{gather*}\]

Hi#