Here we demonstrate the use of the Bogoliubov-de Gennes (BdG) equations to look at the stability of fluctuations in quantum systems. The examples here are motivated by Gross-Pitaevskii (GP) and related equations, but aimed at explaining qualitatively the meaning and interpretation of the solutions.

In [1]:

import mmf_setup;mmf_setup.nbinit()
import numpy as np
from matplotlib import pyplot as plt
import holoviews as hv
hv.notebook_extension(logo=False)

from holoviews import Element2D
from holoviews.plotting.mpl import ElementPlot

class Figure(Element2D):
    """This Figure class allows you to wrap a matplotlib figure 
    but use it as a Holoviews Element.  Eventually everything should
    be done in Holoviews, but until some bugs are worked out, this
    allows you to use all of the matplotlib functionality.
    """
    def __init__(self, fig, *v, **kw):
        plt.close('all')
        Element2D.__init__(self, fig, *v, **kw)

class FigurePlot(ElementPlot):

    def initialize_plot(self, ranges=None):
        element = self.hmap.last
        self.handles['fig'] = element.data
        return self.handles['fig']

    def update_handles(self, key, axis, element, ranges, style):
        self.handles['fig'] = element.data

    def update_frame(self, key, ranges=None, element=None):
        element = self._get_frame(key)
        self.handles['fig'] = element.data
        return self.handles['fig']

hv.Store.register({Figure: FigurePlot}, 'matplotlib')

$\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{\op}[1]{\mathbf{#1}} \newcommand{\mat}[1]{\mathbf{#1}} \newcommand{\d}{\mathrm{d}} \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}}} \DeclareMathOperator{\erf}{erf} \DeclareMathOperator{\Tr}{Tr} \DeclareMathOperator{\order}{O} \DeclareMathOperator{\diag}{diag} \DeclareMathOperator{\sgn}{sgn}$

This cell 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.

HoloViewsJS successfully loaded in this cell.

1. Model¶

We consider a 1D model that followsing from the following Energy density functional:

$$ E[\psi(x)] = \braket{\psi|\op{K}|\psi} + \int \mathcal{E}(n(x))\d{x}, \qquad n(x) = \abs{\psi(x)}^2 = \braket{x|\psi}\braket{\psi|x}, \qquad \psi(x) = \braket{x|\psi}. $$\begin{align} \frac{\delta E}{\delta \psi^\dagger} &= \op{K}\psi + \mathcal{E}'(n)\psi, & \frac{\delta E}{\delta \psi} &= \psi^\dagger\op{K} + \mathcal{E}'(n)\psi^\dagger\\ \frac{\delta^2 E}{(\delta \psi)^2} &= \mathcal{E}''(n)\psi^\dagger{}^2, & \frac{\delta^2 E}{(\delta \psi^\dagger)^2} &= \mathcal{E}''(n)\psi^2, & \frac{\delta^2 E}{\delta \psi^\dagger \delta \psi} &= \op{K} + \mathcal{E}'(n) + n\mathcal{E}''(n)\\ \end{align}

Dynamics in this system follow from finding the stationary solutions of the action for the Schrödinger field

$$ S[\psi(x,t)] = \int\d{t}\;\left[ \int \d{x}\; \psi^\dagger(x, t)\I \hbar\partial_t\psi(x,t) - E[\psi(x, t)]. \right] $$

Stationary points satisfy the GP equation:

$$ \I\hbar\partial_t \ket{\psi} = \frac{\delta E[\psi]}{\delta\psi^\dagger} = \Bigl(\op{K} + \ket{x}\mathcal{E}'\bigl(n(x)\bigr)\bra{x}\Bigr)\ket{\psi} = \op{H}[n]\psi $$

with a definite chemical potential $\mu$ which is the eigenvalue of the Hamiltonian:

$$ \mu = \hbar\omega, \qquad \op{H}[n]\psi = \mu\psi. $$

1.1 BdG¶

Now consider adding a perturbation $\psi_1$ to the stationary state $\psi_0$:

$$ \psi = \psi_0 + \psi_1, \qquad \op{H}[n_0] \psi_0 = \hbar\omega_0\psi_0 = \mu\psi_0. $$

The energy of this perturbation is:

$$ E[\psi] = E[\psi_0] + \mu \int \d{x} (\psi_0^\dagger \psi_1 + \mathrm{h.c.}) + \frac{1}{2}\int \d{x} \begin{pmatrix} \psi_1^\dagger & \psi_1^T \end{pmatrix}\cdot \mat{M} \cdot \begin{pmatrix} \psi_1 \\ \psi_1^* \end{pmatrix} + \order(\psi_1^3) \\ \mat{M} = \begin{pmatrix} \frac{\delta^2 E}{\delta \psi \delta \psi^\dagger} & \frac{\delta^2 E}{(\delta \psi)^2} \\ \frac{\delta^2 E}{(\delta \psi^\dagger)^2} & \frac{\delta^2 E}{\delta \psi^\dagger \delta \psi} \end{pmatrix} = \begin{pmatrix} \op{K} + \mathcal{E}'(n_0) + n_0\mathcal{E}''(n_0) & \mathcal{E}''(n_0)\psi_0^2\\ \mathcal{E}''(n_0)(\psi_0^\dagger)^2& \op{K} + \mathcal{E}'(n_0) + n_0\mathcal{E}''(n_0) \end{pmatrix} $$

The object of the BdG analysis or linear response is to find perturbations that keep their shape as they evolve in time. In order to close the solution, we need to add both positive and negative frequency states:

$$ \psi(x, t) = e^{-\I\omega_0 t}\left(\psi_0(x) + u(x)e^{-\I\omega t} + v^*(x)e^{\I\omega t}\right). $$

Applying the GP equation and keeping terms linear in the perturbation yields the equations:

$$ (\mat{M} - \hbar\omega_0)\begin{pmatrix} u(x)\\ v(x) \end{pmatrix} = \hbar\omega \begin{pmatrix} u(x)\\ -v(x) \end{pmatrix}\\ = \begin{pmatrix} \op{K} + \mathcal{E}'(n_0) + n_0\mathcal{E}''(n_0) - \hbar\omega_0 & \mathcal{E}''(n_0)\psi_0^2(x)\\ \mathcal{E}''(n_0)\psi_0^\dagger{}^2(x) & \op{K} + \mathcal{E}'(n_0) + n_0\mathcal{E}''(n_0) - \hbar\omega_0 \end{pmatrix} \cdot \begin{pmatrix} u(x)\\ v(x) \end{pmatrix}. $$

Inserting these into the energy, we find that the energy of the perturbation is:

$$ E[\psi_0 + \psi_1] = E[\psi_0] + \mu \delta N + \Re\omega\int(\abs{u}^2 - \abs{v}^2)\d{x} $$

where $\delta N = \int (\psi^\dagger\psi - \psi_0^\dagger\psi_0)\d{x}$.

If we consider homogeneous states, we can write:

$$ \psi(x, t) = e^{-\I(\omega_0 t - kx)}\left( \sqrt{n} + ue^{-\I(\omega t-qx)} + v^*e^{\I(\omega t-qx)}\right), \qquad \omega_0 = \epsilon(k) + n\mathcal{E}'(n),\\ \begin{pmatrix} \epsilon(k+q) - \epsilon(k) + n\mathcal{E}''(n) & n\mathcal{E}''(n)\\ n\mathcal{E}''(n) & \epsilon(k-q) - \epsilon(k) + n\mathcal{E}''(n) \end{pmatrix} \cdot \begin{pmatrix} u\\ v \end{pmatrix} = \omega \begin{pmatrix} u\\ -v \end{pmatrix},\\ \omega = \epsilon_- \pm \sqrt{\epsilon_+\bigl(\epsilon_+ + 2n\mathcal{E}''(n)\bigr)}, \qquad \epsilon_- = \frac{\epsilon(k+q) - \epsilon(k-q)}{2}, \qquad \epsilon_+ = \frac{\epsilon(k+q) + \epsilon(k-q) - 2\epsilon(k)}{2},\\ \psi_+ = \begin{pmatrix} 1, \epsilon_- \end{pmatrix} $$

Since this is not a Hermitian eigenproblem, $\omega$ can be complex. When $\omega$ is complex $\abs{u} = \abs{v}$ and the perturbation has zero energy since $\abs{u}^2 - \abs{v}^2 = 0$. When $\omega$ is real, the energy of the state is.

In [40]:

A = 1.0
B = 1.0
C = 2.0
ep = (A+B)/2
em = (A-B)/2
D2 = ep*(ep+2*C)
d2 = abs(D2)
wp = em 
M = np.array([[A+C, C], [C, B+C]])
Z = np.array([[1, 0], [0, -1]])
w, uv = np.linalg.eig(Z.dot(M))
w = w[0]; u, v = uv[:,0]
(v/u)**2, (A+C-w)/(B+C+w)

Out[40]:

(0.14589803375031546, 0.14589803375031543)

1.2 Phonons¶

Consider the density contrast of these pure phonon excitations to order $\order(u^2)$:

$$ \frac{n(x, t) - n}{\sqrt{n}} \approx \Re ue^{-\I(\omega t - qx)} + \Re v^*e^{\I(\omega t - qx)} = a \cos(\omega t - qx + \phi). $$

Thus, the pure phonons have a wavelength $1/q$ and phase velocity $v=\omega/q$. Note, however, that if you just give an excitation like $\delta\psi = \delta e^{\I qx}$ as a perturbation, then this will in general be a linear combination of the four modes with $\pm q$ and $\pm \omega$ so one will often see interference patterns.

2. Numerical Example¶

In [5]:

hbar = m = 1

class SpinOrbitCoupledBEC(object):
    """Represents a 1D spin orbit coupled BEC in the lowest band."""
    def __init__(self, d=-0.1, w=0.2, g=0.3,
                L=100.0, N=256, **kw):
        self.d = d
        self.w = w
        self.g = g
        self.L = L
        self.N = N
        self.__dict__.update(kw)
        self.init()
        
    def init(self):
        """Initialize the object."""
        N, L = self.N, self.L
        dx = self.dx = L/N
        x = self.x = np.arange(N)*dx - L/2.0
        k = self.k = 2*np.pi * np.fft.fftfreq(N, dx)
        
        # Typically we use the FFT, but to be explicit,
        # we form the actual matrices here
        FT = np.exp(-1j*k[:, None]*(x[None,:] + L/2.0))
        IFT = np.exp(1j*(x[:, None] + L/2.0)*k[None, :])/N
        assert np.allclose(IFT.dot(FT), np.eye(N))
        np.random.seed(1)
        psi = np.random.random(N) + np.random.random(N)*1j - 0.5 - 0.5j
        assert np.allclose(np.fft.fft(psi), FT.dot(psi))

        # Here is the kinetic energy matrix
        self.K = IFT.dot(np.diag(self.get_K(k)).dot(FT))
        
    def get_K(self, k):
        """Return the dispersion relationship"""
        return (hbar*k)**2/2.0/m

    def get_K(self, k, s=-1):
        """Return the dispersion relationship.
        
        This version is for a Spin-Orbit coupled BEC
        """
        return (k**2 + 1.0)/2.0 + s*np.sqrt((k+self.d)**2 + self.w**2)

    def get_E(self, n, d=0):
        """Return the d'th derivative of the equation of state."""
        if d == 0:
            return self.g*n**2/2.0
        elif d == 1.0:
            return self.g*n
        elif d == 2:
            return self.g
        else:
            return 0.0

    def get_BdG(self, k, q, n):
        em = (self.get_K(k+q) - self.get_K(k-q))/2.0
        ep = (self.get_K(k+q) + self.get_K(k-q) - 2*self.get_K(k))/2.0
        desc = ep*(ep+2*n*self.get_E(n, d=2)) + 0j
        return em - np.sqrt(desc), em + np.sqrt(desc)
    
    def apply_H(self, psi, n=None):
        """Return `H*psi` with density n."""
        if n is None:
            n = abs(psi)**2
        return self.K.dot(psi) + self.get_E(n, d=1)*psi

In [6]:

import ipywidgets

@ipywidgets.interact(
    k=(-3.0,3.0, 0.01),
    gn=(0.0, 2.0),
    d=(-1.0, 1.0, 0.01),
    w=(0.0, 1.0, 0.01),
)
def draw(k=-0.9, gn=0.3, d=-0.1, w=0.2):
    s = SpinOrbitCoupledBEC(d=d, w=w, g=1.0)
    ks = np.linspace(-6, 6, 1200)
    wm, wp = s.get_BdG(k=k, q=ks, n=gn/s.g)
    fig = plt.figure()
    plt.plot(ks, s.get_K(ks), '-b')
    plt.plot(s.k, s.get_K(s.k), '+b')
    plt.plot(ks, s.get_K(ks, s=1), '-b')
    plt.plot(ks+k, wm.real, '-g')
    plt.plot(ks+k, wp.real, '-g')    
    plt.plot(ks+k, wm.imag, '--g')
    plt.xlabel('k')
    plt.axvline(k, c='y')
    plt.axis([-3,3, -1, 1])
    return Figure(fig)

By using the FFT here, we enforce the boundary conditions of a periodic box. The homogeneous solutions are plane-wave states:

$$ \psi_k = \sqrt{n}e^{-(\I\hbar\omega t - k x)}, \qquad k = \frac{2\pi j}{L}, \qquad \hbar\omega = \mu = \frac{\hbar^2 k^2}{2m} + \mathcal{E}'(n). $$

(The lattice restricts $j$ to be integer in order to assure an exact solution.)

In [7]:

s = SpinOrbitCoupledBEC()

def get_psi_k(k, n=1.0):
    return np.sqrt(n)*np.exp(1j*k*s.x)

for k in s.k:
    for n in [1.0, 2.0]:
        mu = s.get_K(k) + s.get_E(n, d=1)
        psi = get_psi_k(k, n=n)
        assert np.allclose(s.apply_H(psi), mu*psi)

Now we consider adding a small perturbation to these solutions. We add

In [8]:

import holoviews as hv
hv.notebook_extension()

from pytimeode import interfaces, evolvers

class State(interfaces.ArrayStateMixin, SpinOrbitCoupledBEC):
    interfaces.implements(interfaces.IStateForABMEvolvers)

    def __init__(self, psi=None, cooling_phase=1.0, **kw):
        SpinOrbitCoupledBEC.__init__(self, **kw)
        if psi is None:
            psi = 0 * self.x
        self.data = psi + 0j
        self.cooling_phase = cooling_phase

    def compute_dy(self, dy, t=0.0):
        psi = self.data
        n = abs(self.data)**2
        Hpsi = self.apply_H(psi, n=n)
        mu = psi.conj().dot(Hpsi)/psi.conj().dot(psi)
        dy[...] = (Hpsi - mu * psi)/1j/self.cooling_phase
        return dy
    
    def plot(self, k=None):
        """Nice plot of the current state."""
        psi = self[...]
        xs = self.x
        ks = self.k
            
        n = abs(psi)**2
        
        # Compute population ratio K from velocity
        j_quasi = hbar*(psi.conj()*np.fft.ifft(1j*self.k*np.fft.fft(psi))).imag
        k_quasi = np.ma.divide(j_quasi, n).filled(0)
        desc = np.sqrt((k_quasi+self.d)**2 + self.w**2)
        K = (k_quasi+self.d)/desc
        na = (1+K)*n/2.0
        nb = (1-K)*n/2.0
        
        N = (psi.conj().dot(psi)*self.dx).real
        E = ((psi.conj().dot(self.K.dot(psi)) 
              + self.get_E(n).sum())*self.dx).real
        msg = "N={:0.2f}, E={:0.2f}".format(N, E)
        
        psi_k = np.fft.fft(psi)/np.sqrt(self.N)

        fig = plt.figure(figsize=(20, 5))
        plt.subplot(121)
        plt.plot(xs, n, 'b-')
        plt.plot(xs, na, 'b-.')
        plt.plot(xs, nb, 'b:')        
        plt.ylim(0, n.max())
        plt.xlabel('x')
        plt.ylabel('n(x)')
        
        plt.subplot(122)
        plt.plot(np.fft.fftshift(ks), np.fft.fftshift(np.log10(abs(psi_k)**2)))
        plt.xlabel('k')
        plt.ylabel('log(n(k))')

        plt.twinx()
        ks = np.linspace(-6, 6, 400)
        plt.plot(ks, self.get_K(ks), '-b')
        plt.plot(ks, self.get_K(ks, s=1), '-b')
        plt.plot(self.k, self.get_K(self.k), '+b')

        if k is not None:
            wm, wp = self.get_BdG(k=k, q=ks, n=n.mean())
            plt.plot(ks+k, wm.real, '-g')
            plt.plot(ks+k, wp.real, '-g')    
            plt.plot(ks+k, wm.imag, '--g')
            plt.axvline(k, c='y')

        plt.axis([-3,3, -1, 1])
        plt.suptitle(msg)
        return Figure(fig)
    
    def demonstrate_stability(self, dpsi=None,
                              perturbation='x',
                              pos=0, width=0, amplitude=0.1, 
                              k=0, n=1.0,
                              steps=200, dt_Emax=0.1):
        """This sets up a homogeneous state with specified k and perturbs
        it with dpsi.  It the evolves this state so that the evolution
        of the instabilities can be seen.
        
        Arguments
        ---------
        perturbation : 'x', 'k'
            If dpsi is None, then the perturbation is a packet in position
            or momentum space at `pos` with width `width'.
        dpsi : array-like
            Perturbation.  If None, then use a delta function.
        k : float
            Momentum of homogeneous state.  Chosen to be the closest k
            supported by the lattice.
        n : float
            Background density of initial state.
        cooling_phase : complex
            Choose something like 1 + 0.01j to add some cooling.
        """
        k = self.k[np.argmin(np.abs(k - self.k))]
        if dpsi is None:
            if perturbation == 'x':
                X = self.x
            elif perturbation == 'k':
                X = self.k
            else:
                raise NotImplementedError
            ind = np.argmin(abs(X - pos))
            pos = X[ind]
            if width == 0:
                dpsi = 0*X
                dpsi[ind] = 1.0
            else:
                dpsi = np.exp(-(X-pos)**2/width**2/2.0)

            if perturbation == 'k':
                dpsi = np.fft.ifft(dpsi)*np.sqrt(self.N)

            dpsi *= amplitude
            
        psi0 = np.sqrt(n)*np.exp(1j*k*self.x)
        self[...] = psi0 + dpsi
        return hv.DynamicMap(self.get_data(steps=steps, dt_Emax=dt_Emax, k=k))

    def get_data(self, steps=200, dt_Emax=0.1, k=None):
        Emax = abs(self.get_K(self.k)).max()
        evolver = evolvers.EvolverABM(self, dt=dt_Emax/Emax)
        while True:
            yield evolver.y.plot(k=k)
            evolver.evolve(steps)
        
# Silly check to make sure things work.
hv.DynamicMap(State(np.exp(-s.x**4)).get_data())

HoloViewsJS successfully loaded in this cell.

Out[8]:

In [9]:

%%output fps=20
s = State(d=0, cooling_phase=1+0.1j)
s.demonstrate_stability(k=-0.47)

Out[9]:

In [35]:

%%output fps=20
s = State(d=-0.2, cooling_phase=1+0.5j)
s.demonstrate_stability(k=-2)

Out[35]:

In [10]:

%%output fps=20
s = State(d=-0.2, cooling_phase=1+0.01j)
s.demonstrate_stability(k=-1, perturbation='k', pos=0.3, width=0.1, amplitude=0.2)

Out[10]:

In [25]:

%%output fps=20
s = State(d=-0.2, cooling_phase=1+0.01j)
s.demonstrate_stability(k=-1, perturbation='k', pos=-1.1, width=0, amplitude=0.2)

Out[25]:

In [14]:

%%output fps=20
s = State(d=-0.2, cooling_phase=1+0.01j)
s.demonstrate_stability(k=-1, perturbation='x', pos=0, width=2, amplitude=0.2)

Out[14]:

In [18]:

s = State(d=0, cooling_phase=1+1j, steps=1000)
r = 5
na = np.exp(-(s.x/r)**2)
nb = np.exp(-(s.x/r)**2)
n = na + nb
K = (nb - na)/n
k = s.w/np.sqrt(1./K**2 - 0.9999999999)

phi = np.cumsum(k)*s.dx

phase_slip = (phi[-1]-phi[0])/2./np.pi
phase_slip -= int(phase_slip)
phi -= phase_slip * np.arange(len(phi))/len(phi) * 2*np.pi
phase_slip = (phi[-1]-phi[0])/2./np.pi
s[...] = np.sqrt(n)*np.exp(1j*phi)
hv.DynamicMap(s.get_data())

Out[18]:

In [78]:

s[...] = np.exp(-(s.x/5)**2/2)*np.exp(1j*s.k[16]*s.x)
hv.DynamicMap(s.get_data())

Out[78]:

In [ ]: