Hamiltonians which commute both as operators and as connections This is something which I suspect is written up in introductory books on mathematical physics if I knew where to look. Suppose I have some parameters $t_1$, ..., $t_k$ ranging over a neighborhood in $\mathbb{R}^k$. I also have $k$ matrix-valued functions of the $t$'s: $H_1(t_1, \ldots, t_k)$, ... $H_k(t_1, \ldots, t_k)$. These obey both
$$[H_i, H_j]=0 \quad (\ast)$$
and
$$[\partial_i+H_i, \partial_j+H_j] =0 \quad (\dagger).$$
For those who don't like the language of connections, we can expand $(\dagger)$ as $\partial H_i/\partial t_j - \partial H_j/\partial t_i + [H_i, H_j]=0$ or, in the presence of $(\ast)$, as 
$$\frac{\partial H_i}{\partial t_j} = \frac{\partial H_j}{\partial t_i}.$$
Equation $(\ast)$ tells us that, assuming the $H_i$ are individually diagonalizable, we can find $u(t)$ a simultaneous eigenvector for all the $H_i$:
$$H_i(t) u(t) = \lambda_i(t) u(t). \quad (\ast\ast)$$
Equation $(\dagger)$ tells us that the vector-valued PDE 
$$\frac{\partial v}{\partial t_i} + H_i v=0 \quad (\dagger \dagger)$$
will have a unique solution $v(t_1, t_2, \ldots, t_n)$ for any initial value.
I'm pretty sure there is supposed to be a relation between the solutions to  $(\ast \ast)$ and $(\dagger \dagger)$. What is the right statement, and what is the keyword to read about this situation?
Motivation: I'm trying to work through the papers of Varchenko, Scherbak and others on the KZ equation. I think it would really clear my head to just see this scenario described abstractly without all the details of which operators they are thinking about.
$\def\mg{\mathfrak{gl}_n}$ Edit to spell out the relation. Let $V_1$, $V_2$, ..., $V_n$ be representations of $\mg$. So $U(\mg)^{\otimes n}$ acts on $V_1 \otimes V_2 \otimes \cdots \otimes V_n$. Let $\Omega \in U(\mg) \otimes U(\mg)$ be the Casimir. (Note: The element I learned to call the Casimir was a central element $c$ in $U(g)$. In terms of that element, $\Omega = \Delta(c) - c \otimes 1 - 1 \otimes c$.) Let $\Omega_{ij}$ be $\Omega$ acting in positions $i$ and $j$.
For generic parameters $z_1$, ..., $z_n$, define $H_i = \sum_{j \neq i} \Omega_{ij}/(z_i-z_j)$. Then, as I understand it, the KZ equation is $(\partial_i + H_i) v(z_1, \ldots, z_n)=0$, where $v$ is a function valued in $V_1 \otimes V_2 \otimes \cdots \otimes V_n$. The $H_i$'s obey both $(\ast)$ and $(\dagger)$ (a nice exercise). And people seem to be very interested in solving both $(\ast \ast)$ "diagonalizing the action of the Gaudin subalgebra" and $(\dagger \dagger)$ "solving the KZ equation". So I was hoping to understand how they relate, and why.
 A: Hi David,
I think there is indeed a relation, which I learned precisely from papers of Varchenko among others. All of this is rather classical and can be found e.g. in Etingof-Frenkel-Kirilov book "Lectures on Representation Theory and Knizhnik-Zamolodchikov Equations". 
The fact that the $H_i$ satisfies this stronger condition is equivalent to say that for any parameter $\kappa$ the operators $\kappa \partial_i+H_i$ alos satisfies ($\dagger$).
Hence you can take asymptotic expansion of solutions at $\kappa \rightarrow 0$ on some neighbourhood $D$ of some $z_0$, of the form
$$e^{S(z)/\kappa} (f_0(z)+O(\kappa))$$
where $S$ is a scalar valued function. Then you can show that, assuming the $H_i(z)$ are simultaneously digonalizable then $f_0$ is a common eigenvector of them, with eigenvalues $\partial_i S$. Conversly given a common eigenvector at some $z_0$ you can construct an asymptotic solution. So the usual trick, widely used in the study of the KZ equation, is to also take some asymptotic limit w.r.t. the variable $z_i$ in such a way that eigenvectors are "easy" to find. The standard example in the KZ case is the asymptotic zone 
$$|z_i-z_1| \ll |z_j -z_1|\quad if\quad i < j $$
for which, up to some change of variable, the equation can be written
$$\kappa \partial_i f= \left ( \Omega_i/u_i +reg\right)f\quad i=1\dots n-1$$
where $\Omega_i=\sum_{k < i} \Omega_{k,i+1}$ and $reg$ is regular at $u=0$. 
Then given some common eigenvector $v$ of the $\Omega_i$ with eigenvalues $\mu_i$ there exists a unique solution of the form
$$(\prod u_i^{\mu_i/\kappa})(v+r(u))$$
where $r(u)$ is regular at $u=0$ and $r(0)=0$.
I'm not very familiar with D-modules (and by the way I would be happy is someone extends on this), but you can rephrase it as follows: viewing $\kappa$ as a formal variable leads to a filtration on the algebra of differential operators on $V$ (the vector space acted on by the $H_i$) which in turn is nothing but the usual filtration by the degree of differential operators. Taking the associated graded turns the equation
$$(\partial_i+H_i)f=0$$
into the equation
$$(y_i+H_i)f=0$$
Whose solutions are clearly commons eigenvectors of $H_i$. So I'm rather confident that you can say that the spectrum of the $H_i$ for all common eigenvectors is the characteristic variety of the D-module of solutions of the differential equation you started with.
A: I can't comment on the case of several operators $H_i$, but for a single operator, the eigenvector equation $(**)$
$$ H(τ) \psi_τ = λ(τ) \psi_τ $$
and the time-dependent Schrödinger equation $(\dagger\dagger)$
$$ (i\frac{\partial}{\partial t} - H(t)) ψ(t) = 0 $$
are related by the adiabatic theorem.
Not sure if that's what you are looking for, but I would be very surprised if your setting didn't have a similar intuition.

Essentially, the idea of the adiabatic theorem is the following: the eigenvector equation describes, for each parameter $\tau$, an instantaneous eigenvector $\psi_{\tau}$. This gives a solution $\psi_{\tau}(t) = e^{-itλ(τ)} ψ_τ$ to the "instantaneous" Schrödinger equation
$$ (i\frac{\partial}{\partial t} - H(τ)) ψ_τ(t) = 0 $$
where the Hamiltonian $H$ is considered at a fixed time $τ$.
Now, if the Hamiltonian $H(t)$ varies very "slowly" in time, then it is reasonable to expect that the full Schrödinger equation will essentially follow the solutions to the "instantaneous" Schrödinger equation(s). First it evolves like a solution of the instantaneous equation with $H(0)$, then for $H(\Delta t)$ a small time step after, and so on.
This can be made precise by rescaling time to $τ=t/T$ and obtaining an asymptotic expansion
$$ ψ(t) = e^{-i∫λ(τ)dt} ψ_τ + \mathcal O(1/T) $$
in the limit $T\to ∞$ and in the $L^2$ sense. More details can be found wherever you can find details about the adiabatic theorem.
A: To expand on Greg's answer regarding the adiabatic theorem. You are looking for situations where the adiabatic evolution is exact. This is the case for a Hamiltonian of the form 
$H = i\left[\frac{\partial P}{\partial t},P\right]$
where $P$ is a projector onto your chosen instantaneous eigenstate. This comes from T, Kato, J. Phys. Soc. J. Jpn. 5, 435 (1950).
A: Without further assumptions, it does not seem that much can be said. Consider the case k=1. You are asking for a connection between the eigenvalue problem for H(t) and the equation dv/dt+Hv=0. But time-dependent linear ODE systems cannot in general be related to the eigenvalue problem.
On the other hand, the condition $\partial H_i/\partial t_j=\partial H_j/\partial t_i$ implies that
$H_i=\partial K/\partial t_i$ for some $K$. Let us now strengthen your assumptions and assume that the $H_i$ commute not only with each other, but also with $K$. Then solutions of ($\dagger\dagger$) can be written as $\exp(-K(t))w$ for fixed $w$, and solutions of ($**$) can be written as $u=\partial v/\partial t_i$, where $v$ is an eigenfunction of $K$.
