Pulsing a two-band model to discover topology

In systems with an anomalous quantum Hall effect, the quantized Hall conductivity comes from the integral of a Chern number over some manifold. Usually, this integral is derived via Kubo formula. However, there is geometry involved in how the state evolves, and in fact, we can use the dynamics of the current following a weak pulse in order to find the DC conductivity. The route is easy enough: Say we have a conductivity which when written with respect to time is $\sigma_{xy}(t-t')$, and without loss of generality we apply a pulse $E_x(t) = A_x\delta(t)$, then we can find the current response in the perpendicular direction

This allows us to derive an expression for the DC-conductivity

Geometrically, there is a lot going on with $j_y(t)$ when we have a system with spin-orbit coupling. In particular, take the two band model

where $\mathbf d$ is 3D, $\mathbf p$ is 2D, and $\sigma = (\sigma_x, \sigma_y, \sigma_z)$, the vector of Pauli matrices. The initial states of the system can be represented by where they are on the Bloch sphere $-\hat{\mathbf d}(\mathbf p)$. But once a pulse is supplied, this state will begin to rotate about a different vector $\mathbf d(\mathbf p - e \mathbf A)$. Thus, if we add time dependence to $\hat{\mathbf d}(\mathbf p, t)\equiv -\langle \sigma(t) \rangle$ to represent the state’s location, we can use Heisenberg’s equations of motion to obtain

We can rewrite this equation as

However, this state has an associated current with it, and that can be represented by the operator $j_\mu = -e \partial_\mu \mathbf d(\mathbf p- e \mathbf A) \cdot \sigma$. And the vector $\langle\sigma\rangle = -\hat{\mathbf{d}}(p,t)$. Therefore,

Combining these expressions, we have $\begin{multline} \langle j_\mu \rangle = e^2 \partial_\mu \mathbf d(\mathbf p - e\mathbf A ) \cdot \Bigg[ \hat{\mathbf d}(\mathbf p - e \mathbf A) [\hat{\mathbf d}(\mathbf p - e \mathbf A) \cdot \hat{\mathbf d}(\mathbf p)] \\ - \hbar \frac{\hat{\mathbf d}(\mathbf p-e \mathbf A) \times \frac{\partial \hat{\mathbf d}(\mathbf p, t)}{\partial t}}{2d(\mathbf p - e A)} \Bigg]. \end{multline}$

Or simplified $\begin{multline} \langle j_\mu \rangle = e^2 \partial_\mu d(\mathbf p - e\mathbf A )[\hat{\mathbf d}(\mathbf p - e \mathbf A) \cdot \hat{\mathbf d}(\mathbf p)] \\ - \frac{e^2}{2} \partial_\mu \hat{\mathbf{d}}(\mathbf p - e \mathbf A) \cdot \left[\hat{\mathbf d}(\mathbf p-e \mathbf A) \times \frac{\partial \hat{\mathbf d}(\mathbf p, t)}{\partial t} \right]. \end{multline}$

This is exact. At this point, we make a couple of approximations. First of all, the first term is independent of $t$ so it cannot contribute to the total current if we have a finite DC conductivity. This leaves the second term. We can do the integral over time — there is an order of limits problem but we can get around this by noting that we do not expect the infinite time state to contribute to the energy (or: it averages to something proportional to $\mathbf d(\mathbf p - e \mathbf A)$ anway and so the cross product vanishes), so we discard it and therefore, $\int_0^\infty dt' \, \mathbf d(\mathbf p, t') = - \mathbf d(\mathbf p)$.

Hence, we get the Hall conductivity

At this point, we actually have not expanded in terms of $A_x$ yet. The first term will produce a term that is symmetric in $x$ and $y$—however it drops out due to the cross product $\hat{\mathbf d}(\mathbf p) \times \hat{\mathbf d}(\mathbf p) = 0$ vanishes. Thus, only the second term persists and we see directly that $x$ and $y$ must be different and in fact, we get the well-known formula

This describes the Chern number of some manifold parametrized by $\mathbf p$ (usually the Brillioun zone). It is the number of times the vector $\mathbf d$ wraps the sphere.

To understand why it’s a topological invariant, note that the quantity in the integral looks very much like a Jacobian. In fact, it is; it describes a coordinate transformation from the $\hat {\mathbf d}$ to $(p_x,p_y)$. In this way, the integrand represents an area element on the sphere, and in general $\mathbf p$ is a closed manifold. So $\mathbf d(\mathbf p)$ maps that closed manifold to the sphere, and without any edges or boundaries the area it maps out must be $4 \pi$ times an integer.

This formula is well-known, but this dynamical way of obtaining it is slightly less well-known. We have extended this idea in a paper published last year to handle the out-of-equilibrium case of quenches. In that situation, new phenonmena appear that are quite different from the equilibrium case—terms that we discarded in this calculation become quite relevant.

Integrating org-mode and PushBullet for automated reminders

For a while now, I have been using Emacs’ org-mode as my fully customizable, open-source, plain-text to-do list manager. In order to help me not escape my agenda items though, I wanted some way to receive push notifications, and for that I am using PushBullet. With a little bit of python code, this can easily be accomplished.

Note: This guide uses these services on OS X with Homebrew installed. It should be easily generalizable to other systems.

Set up your machine

Make sure you have the most recent version of emacs. Using Homebrew:

brew install emacs


We are going to interface everything with python, and for that we need to install pushbullet.py with

pip install pushbullet.py


This requires python-magic (which should be automatically installed), but which requires dependencies which we quote here

On Windows, install Cygwin (http://cygwin.com/install.html). To find the libraries, either add /bin to the $PATH or copy cygwin1.dll, cygz.dll, and cygmagic-1.dll to C:\Windows\System32 On OSX: • When using Homebrew: brew install libmagic • When using macports: port install file Now, I assume you have org-mode set up with a list of agenda files hanging around somewhere on your machine (usually identified in your .emacs file or your customizations.el file for Aquamacs). The Python code I cobbled together some python code with help from the org-mode site explaining how to extract agenda information as well as the python implementation I found on this blog in his agenda.py. It’s important to note that the loadfile must define org-agenda-files. Additionally, I had problems if the loadfile tried to do too much upon initialization. If your init file does too much, you may need to make this loadfile something else, something simpler. The push notifications you receive from running this code are specifically tailored to how I wanted them output. One can easily mess around with the section I labeled “PROCESSING Agenda information into notes to send”. Additionally, there are more ways than just pb.push_note to send information (see the documentation here). Set up Automator There are numerous tools to automate tasks on different systems. On OS X, we’ll use Automator. 1. Create a new application. 2. Find the “Run Shell Script” action and drag it to the box on the right. 3. Select shell /bin/bash and insert the following code export PATH=/usr/local/bin:$PATH
python /path/to/org-mode-agenda-push-bullet.py


This is the simplest since the shell /usr/bin/python doesn’t necessarily use the python that we want (located in /usr/local/bin in this case).

From here, we open the Calendar application, create a new event, and insert when we want it to run (frequency, etc.), then under the “alert:” option, we click “Custom…”, then in the first drop down menu “Open file”, in the new drop down menu that appears we click “Other…”, and navigate to where we saved the automator application created above. The rest should be self-explanatory.

If all is well and good, at the specified times, you should get a push bullet reminder at that time (if your computer is on and connected to the internet).

Subtleties in linear response theory

In linear response theory, we consider some small perturbation to a Hamiltonian and look at the response of some observable to that perturbation. In the case considered here, the perturbation is an electric field, and the response is current. The linear response that characterizes these quantities is called the conductance.

There’s a problem though, an electric field accelerates a charge. Consider a classical electron for the time being, then

Or in terms of current $j(t) = -e \dot x(t)$, $\frac{d j}{dt} = \frac{e^2}{m} \mathbf{E}$. From here we can quickly and naively go to frequency space to find $j(\omega) = i \frac{e^2/m}{\omega} \mathbf{E}(\omega)$. Then one might remember that another way to define $\mathbf{E}$ is in terms of a vector potential that is purely time-dependent, so $\mathbf{E}(\omega) = i \omega \mathbf A(\omega)$. Now, if we just plug this into our linear response for the current, we get

All is well and good, right? Well, not quite. In electromagnetism, the constant part of $\mathbf A(t)$ corresponds to the $\omega = 0$ term of $A(\omega)$. This represents what is known as a “pure gauge”. These gauges are physically equivalent to the null field $\mathbf A(\omega=0) = 0$. Thus, whatever linear response is represented above at $\omega = 0$ must be unphysical, right?

Wrong.

Before explainging why this is wrong, let’s give some further context to this linear response theory. The term $- e^2/m$ is actually the single particle term of what is known as the “diamagnetic” response to the conductivity when you add in more electrons (usually distributed in a Fermi distribution). This term persists in quantum mechanics, and no other terms appear to cancel it in the simplest case of $H = \frac{p^2}{2m}$. In fact, while the math becomes more cumbersome, the solution we shall illustate below holds perfectly well for even the non-interacting multi-electron system.

Now, at this point you may have guessed that there’s something strange going on at $\omega = 0$ due to the fact that the electric field accelerates the particle and doesn’t just have a velocity response. At the $\omega = 0$ point, the physical field $E(\omega)$ seems to necessarily be equal to zero in the gauge we have prescribed unless $A(\omega) \sim 1/\omega$ for small $\omega$. This would lead to a divergent $j(\omega)$, restoring our faith that the system is accelerating out of control.

But what about when $\mathbf{A}(\omega=0) = \mathbf{A}_0$? It seems like then we have a true velocity response to an unphysical object. The solution is subtle: At some point in the quick derivation we made an assumption that implied $\mathbf A(t) \rightarrow 0$ at $t \rightarrow -\infty$. This implies that if $\mathbf A(t) = \mathbf A_0$ at any finite time, there had to be some time in between where $d\mathbf A/ dt \neq 0$. Thus, during that “ramp up” time, an electric field was on and it accelerated the charge to a specific velocity resulting in the current $j(\omega) = - \frac{ e^2}{m} \mathbf A(\omega)$.

The assumption is subtle, but the result is rather simple. For now, just assume that $\mathbf A(-\infty) = 0$ and at some $t_0$, $\mathbf A(t_0) = \mathbf A_0$, then we can integrate Eq. \eqref{eq:Newton} to obtain the velocity:

Or, in other words, $\mathbf j(t_0) = -\frac{e^2}{m} \mathbf A_0$, the same as before! This is how a constant $\mathbf A_0$ can be physical: When it represents the change from a different constant vector potential.

Now, to isolate the assumption, let us run through what they were

1. $\frac{d j}{d t} = \frac{e^2}{m} E(t)$.
2. Take the Fourier transform: $j(\omega) = i \frac{e^2/m}{\omega} E(\omega)$.
3. Insert vector potential with $E = -\frac{d}{dt} A$ and assume the Fourier transform exists for $A$: $j(\omega) = - \frac{e^2}m A(\omega)$.
4. Undo the Fourier transform: $j(t) = - \frac{e^2}{m} A(t)$.

Now, let us get the last equation (in #4 above) by a simpler route.

1. $\frac{d j}{dt} = \frac{e^2}{m} E(t)$.
2. Use $E = - \frac{d}{dt} A$ and integrate the above expression from $-\infty$ to $t$: $j(t) - j(-\infty) = -\frac{e^2}{m} ( A(t) - A(-\infty))$.

Two perfectly legitimate calculations resulting in different results. Firstly, this highlights that the first procedure does actually assume $A(-\infty) = 0$. Secondly, the only assumptions that could have given $j(-\infty) = 0$ (an assumption we probably wanted anyway) and $A(-\infty) = 0$ are that they could be given in terms of Fourier transforms. In order for a function to have a Fourier transform it needs to be absolutely integrable—i.e. $% $. Given $A(t)$ as a continuous, piece-wise differentiable function, we need $A(t) \rightarrow 0$ for $t \rightarrow \pm \infty$. This imposes our gauge, and since we are not interested in future times let alone $t \rightarrow +\infty$, we can artificially modify the function as we see fit to accomodate that. But how the function began at $-\infty$ is important, and we must impose that. Hence, we have chosen, at least partially, a gauge.

We are left with a dilemma then about pure plain waves $A(t) = A(\omega) e^{-i \omega t}$. How do those function?

Technically, they are outside of the bounds of the Fourier analysis and we can see that simply by the fact that if we tried to do the above procedure, we couldn’t have a well defined answer as $t \rightarrow -\infty$ (too oscillatory). However, we can approximate the plain wave in terms of an absolutely integrable function $A_\delta(t) = A(\omega) e^{-i (\omega + i \delta) t}$ for any $\delta >0$, and everything works. This shows us explicitly that $t = - \infty$ does have $A_\delta \rightarrow 0$ for all $\delta>0$. And this is the origin of the well known substitution $\omega \rightarrow \omega + i\delta$.

The natural question to ask now is how this works for a real system (with dissipation). Why does such a term not exist at zero frequency?

Unless your system is a superconductor, there is some dissipation in the system. The simplest way to include this is classically: When an electron is going at velocity $\dot x(t)$ it experiences a “drag” that tends to slow it down. Thus, our Newton’s equations become

where $\gamma$ describes how much drag the electron experiences. For more disorder, this would be a larger number. Playing the same Fourier transform game, we can obtain rather quickly that

This is just one step away from the well-known Drude model. We see that if $A(t) = A_0$, then $j(\omega) = 0$. But our gauge choice that we described before is still in place, the only difference is that our “kick” at $t=-\infty$ has an infinite time to dissipate back to rest (the inclusion of $\gamma$ above is critical for $j(\omega=0) =0$). This also suggests a steady state current when $\mathbf E$ is constant: $m\ddot x=0$ implies $j = \frac{e^2}{m \gamma} \mathbf E$. Our current relaxes to zero when there’s nothing around ($\mathbf E = 0$), as we would expect.

When a quantum mechanical description is done—by taking a random disorder potential and averaging over disorder configurations—one obtains similar results. The diamagnetic term for a clean system is real and has a physically well defined explanation.

One may not be surprised that this curious “diamagnetic term” occurs for superconductors, however it is sometimes explained that “gauge symmetry is broken” and that is why such a term exists. This is a misleading statement, but one I will explore in a future post.

Current in Single Particle Quantum Mechanics

For simplicity, I will only use one-dimension in this post, but this can be generalized to higher dimensions rather easily.

Many textbooks on Quantum Mechanics mention current density can be derived from the continuity equation and probability. The usual method for figuring this out is to assume you have some Hamiltonian $H = p^2/2m + V(x)$ where $p$ is the momentum and $x$ is the position. In this way the current density is written in terms of the wave function $\psi(x,t)$ as

This then satisfies the continuity equation

with density $\rho(x,t)=\psi^*(x,t)\psi(x,t)$. It should be noted that if you write the wave function as $\psi(x,t) = \sqrt{\rho(x,t)} e^{i \theta(x,t)}$, then the current is just proportional to the gradient of the phase $j(x,t)= \rho(x,t) \partial_x \theta(x,t)/m$, giving the spatial change in phase a physical significance.

However, there are two lingering questions:

1. Is this current density related to the Heisenberg operator $\dot x(t)$ which tracks the velocity of the system?
2. If so, does it generalize to more arbitrary Hamiltonians?

To answer these questions, we consider the more arbitrary Hamiltonian

where $V(x)$ is some potential and the kinetic energy is some polynomial

We are unworried about bounding the energy, so odd-order Kinetic energy terms are allowed (in the higher dimensional case, the Dirac-like Hamiltonians have linear terms in $p$). At this point, we can take our Heisenberg operator $\dot x(t)$ and find

where $T'$ is the derivative of $T$ with respect to its argument. Now, we would like to obtain a current density from this quantity. We can certainly define the total current at a specific time as

where in the last line we go from the Heisenberg to Schroedinger picture. Now to get density, we need to use a complete set position states, so that

Now, $p$ acts as a derivative on position kets, so that one can verify that

However, there is an ambiguity here since we can write

This ambiguiuty in how to choose the derivatives leaves us with many way to define the current density. Fortunately, only one of these combinations satisfies the continuity equation. To figure out which one that is, let us reverse engineer the continuity equation to obtain a solution. The density is $\rho(x,t) = \psi^*(x,t) \psi(x,t)$, and so using the Schroedinger’s equation, we have

Thus, the continuity equation must become

If we now assume that we have a current density that takes the form

and satisfies the continuity equation, $\partial_t \rho + \partial_x j = 0$, then we can equate operators to obtain

Anticipating the answer, we write the general form of $\vartheta$ as

Then we can take the left hand side Eq. \eqref{eq:diff-ops-cty} and write

On the other hand, we can calculate the right hand side of Eq. \eqref{eq:diff-ops-cty} to be

Equating the left and right sides, we can just read off that $b_{n-1,0} = a_n/n!$, $b_{n-1,n-1}= a_n/n!$ and $b_{n,m+1} = b_{n,m}$, so that $b_{n,m} = a_{n+1}/(n+1)!$.

Thus, we have

Returning all the way to when we were considering $\dot x(t)$ as an integral over position, this suggests that in Eq. \eqref{eq:T-delta}, we want to consider

Given the expression for total current Eq. \eqref{eq:total-current} and integrating the delta function by parts numerous times, we can replace $\partial_x$ with $-\overleftarrow \partial_x$ and $\partial_y$ with $- \overrightarrow\partial_x$, and then the total current is just

which actually integrates the current density! Thus, we have shown that

and that

Indeed, $\dot x(t)$ does track the current of the problem and can even be written as the integral of a current density. Even for the more arbitrary Hamiltonian $H = T(p) + V(x)$.

The δ-function Potential Lattice

I was messing around with some simple problems, and found this simple but illustrative problem. It starts you off in basic quantum mechanics and introduces concepts in a very straightforward way to both get to a more condensed matter perspective while even showing some interesting effects that have experimental consequences (energy bands and band gaps opening).

We look at the relatively simple problem1 of finding the energy spectrum for a particle in the lattice potential

Visualized, it looks like the picture:

The time-independent Schrödinger equation takes the form

where $E$ is the energy.

Since we can solve the problem between the delta-functions quite simply ($U(x) =0$ there), let us restrict our focus to $% $. Here, the wave function takes on the form

where we have $k = \sqrt{2m E}$. Now, we can find an operator that commutes with the Hamiltonian so that we can diagonlize it to help solve the problem — this will be the operator that translates us by $a$.

To make this clear, let us abstract things to operators so that we have a momentum operator $p$ and a position operator $x$, then we have the commutator $[ x, p] = i$. The operator $p$ commutes with functions of $x$ as though it were a derivative $[ p, f( x)] = -i f'( x)$, so considering the translation operator $e^{i a p }$, we can write

and since $U(x)$ is periodic in $a$, we have that $e^{i a p} U( x) e^{ - i a p} = U( x)$. Thus, the operator $T_a = e^{i a p}$ commutes with the Hamiltonian and we can simulataneously diagonalize both it and the Hamiltonian. We say $T_a \lvert \psi\rangle = e^{i a q} \lvert \psi \rangle$ has quasi-momentum $q$. It is important to note that this not the same as real momentum which is not a well-defined quantum number in this problem (that needs translation symmetry).

In other words, we can write our eigenfunctions such that $\psi(x+a) = e^{i q a} \psi(x)$, and this naturally leads us to relate the coefficients for our eigenfunctions above as

Now, we can apply matching conditions at $x = na$ remembering that our wavefunction should be continuous, but that the delta-function will cause the first derivatives to be discontinuous. The equations for the coefficients are

and if we insert the relation between the coefficients at $n-1$ and $n$, we get a matrix equation

This equation has a non-zero solution only if the determinant of the matrix is zero which can be written as

This is an equation which relates the energy to the quasi-momentum. Since $\cos q a$ can only be between -1 and 1, this equation only has a solution when $f(E)$ is between -1 and 1. The oscillatory nature of $f(E)$ means that it should pass in this range multiple (in fact, a countable infinite) number of times.

Armed with this equation, we can use $q$ and an integer to label our energies and we obtain the following energy bands by just solving for $E$ (setting $m = 10$, $a = 1$, and $\alpha = 0.3$)

The dotted lines in this plot represent the spectrum if there were no delta-function potentials (displaced in energy by $\alpha / a$ for clarity). Notice how the introduction of the delta-functions opens up gaps in this energy spectrum, so that there are some energies that are inaccessible. The gaps actually open up when $\lvert f(E)\rvert \geq 1$ since there is no solution to our equation there. This energy gap for small $\alpha$ just goes like $\alpha/a$, vanishing as we’d expect when $\alpha = 0$.

Notice that in this energy spectrum, we see that if $q \rightarrow -q$ we get the same energy.

Additionally, if the energy ever goes negative the solutions turn from plane waves $e^{\pm i k (x-na)}$ into functions localized around the delta functions $e^{\pm\kappa(x-na)}$. In fact, the wave functions look like this for a $q=0$ state:

These only appear when $% $, and are related to the fact that the delta-function potential has a bound state. Additionally, only one band can ever have this state. This is due to the fact that the oscillatory sine and cosine change to their non-oscillatory hyperbolic counterparts. However, these states in the delta-function lattice are spread out throughout the crystal, and can not be said to be “localized” to a specific site — they still all have a definite quasi-momentum.

As with other single particle problems, upon considering the many particle picture, these energy bands get filled up to a set energy level (if we are considering fermions).

1. This is problem 2.53 in Exploring Quantum Mechanics by Galitski, Karnakov, Kogan, and Galitski.