# compute change of phase along closed contour

The following image represents the phase of a wavefunction (in radians) on a square lattice, where $m$ and $n$ label the lattice sites. Computationally speaking, it is the density plot of a 41x41 real matrix.

My question is, how can I numerically compute the total change in phase as one goes around a closed loop centered on the central site $m=n=0$?

Let's write the closed loop as a series of lattice sites $(m_1, n_1), (m_2, n_2), \ldots (m_K, n_K)$ where each site in the series should be a neighbor of the preceding site. Denote the phase at the $i$th site as $\phi_i$. If the phase is well approximated as a continuous function, it's trivial to compute the phase accumulated around the contour: $(\phi_2 - \phi_1) + (\phi_3 - \phi_2) + \ldots (\phi_K - \phi_{K-1}) + (\phi_1 - \phi_K) = 0$.

The interesting part comes when the phase is discontinuous. When the phase at a site exceeds $2\pi$ or falls below $0$, your numerical algorithm typically "wraps" the phase by subtracting or adding $2 \pi$. This is clearly happening in your plot quite a bit, for instance when you go from $(m,0)$ to $(m,1)$ with $m \gg 0$.

Concretely, if you see that $\phi_j = 2 \pi - \epsilon_1$ and $\phi_{j+1} = \epsilon_2$ for $\epsilon_1, \epsilon_2 \ll 1$, you probably have a phase wrapping of $+2\pi$. Similarly you could have $-2\pi$ phase wrappings. Between $i = 1$ and $i = K$, sum up all the phase wrappings. This sum is equal to the phase accumulated around the contour.

To find the phase steps in code, you can just test each successive pair of sites in your contour as in the example above. However, the numerical difficulty comes in choosing $\epsilon_1$ and $\epsilon_2$. If your data is too coarse, you will not be able to choose them small enough for reliable unwrapping.

• Ok, so I have a naive question. I see you are talking about phase unwrapping, and I discovered this can be done in 2D as well (eg, see ljmu.ac.uk/external-assets/research/…). However, since this effectively turns my phase function into a continuous one, won't the accumulated phase along the contour be always zero after unwrapping? – Andrei Feb 4 '15 at 17:09
• It's not a naive question at all. Generally, 2D phase unwrapping starts by assuming that $\nabla \times \nabla \phi = 0$ for the "true value" of the phase, i.e., the phase gradient is conservative. If that's true for you, go ahead and try 2D unwrapping. But a lot of current quantum mechanics research involves wavefunctions with vortices, for which the phase gradient is not conservative. I believe 2D phase unwrapping is not applicable in this case. Therefore, I was trying to give you an answer that would work in both cases. – Dave Kielpinski Feb 5 '15 at 0:14
• PS. Apparently there is a Matlab function unwrap for 1D phase and some 2D unwrapping algorithms on the Matlab File Exchange. – Dave Kielpinski Feb 5 '15 at 0:16
• Is this equivalent to computing the winding number (en.wikipedia.org/wiki/Winding_number#Complex_analysis) of my original complex wave function (from which the phase was extracted) for a closed path around the origin? – Andrei Feb 5 '15 at 13:00
• @Andrei Yes, it's equivalent. I do not have a definite opinion on which method will be numerically more stable, though. – Dave Kielpinski Feb 6 '15 at 1:15