Finding the correct order of eigenvectors of a parameter-dependent Hermitian matrix

Question

so, I have a symmetric, analytic matrix $\mathbf{H}(x)$ ($x$ is real). Because $\mathbf{H}(x)$ is analytic and $x$ is real, it is possible to find analytic functions for the eigenvectors and the eigenvalues of $\mathbf{H}(x)$. at $x=0$, each eigenvalue $\vec v_i(0)$ is associated with an eigenvalue $\lambda_i(0)$, where the eigenvalues are ordered $\lambda_{i}\leq\lambda_{i+1}$.

As $x$ changes, some eigenvector $\vec v_i(x)$ "switch place" because the eigenvalues cross. This is illustrated in the following example. The green and grey lines switch place at $x\approx1.2$, because I used a diagonalization algorithm that orders the eigenvalues by their values:

Algorithms such as numpy's og scipy's $\texttt{eigh}$, of course, do not have information about the function $\mathbf{H}(x)$ and simply order the eigenvalues by their value for each $x$.

So, my questions are:

1. Does anyone know of a good algorithm to re-order the eigenvector-eigenvalue pair in such a way that they are analytic functions (that is, to "swap green and grey" at $x\approx1.2$ etc.?)
2. In my application, I will only have the eigenvalues and eigenvectors for few $x$ (I can generate them for more values, but that is expensive), but I still need them in such an order that they are associated with the right function. Hence, the function $\mathbf{H}(x)$ is not explicitely available.
3. Some eigenvalues are degenerate for all $x$ (in addition to them switching places with other eigenvalues). So, numpy/scipy will return a "random" orthogonal pair of vectors. Of course, the subspace spanned by these eigenvectors will be analytic, but some type of algorithm needs to be applied to make the individual eigenvectors analytic.

What I have been thinking so far:

1. Because the eigenvectors are analytic, when no crossings occur, we have that $\vec v_i(x)^T \vec v_j(x+\delta x)\approx\delta_{ij}$. Based on this observation, we can "regain" the correct order. This seems to work pretty well, unlike (Continuity of eigenvectors of parametric matrix), where it seems to have failed.
2. The same answer as 1., where some threshold is used for when more samples need to be obtained.
3. One can use the orthogonal Procrustes algorithm to "similarize" the pair of vectors as much as possible to some arbitrary reference,for example $x=0$, though I am not sure whether this will yield analytic eigenvectors.

Before I implement any more of this, it would be nice to hear if someone has a more elegant solution or possibly knows of an implementation in some library (ideally Python).

Answer 1

As the commenters have noted, this is a common issue in condensed matter physics. There are a few different solutions. The one I currently use (pretty good) is explained in this blog post.

The rough idea is that when these eigenvalue crossings occur, there is some sense in which you want to track 'the same eigenstate' you had before. This is not necessarily accurate, as when the two eigenvalues cross, the eigenstates can mix their character. However, this is a good heuristic. What you do is try to match the eigenvectors between neighboring x, by reordering eigenvalues such that the associated eigenfunctions have maximal overlap between x and x'. The blog post I linked has a lot more detail, and what's better, code implementations in Python.