Adiabatic evolution of superposition of states I am supposed to find a specific superposition of eigenstates of a time-dependent hamiltonian $H(t)$. The hamiltonian is of the form : $H(t) =\sum_i \left(J\vec{\sigma}_i\cdot \vec{\sigma}_{i+1} + h(t) {\sigma}_i^z \right)$. Periodic Boundary Condition.
I have taken the following state at $t=0$,
$$\vert \Psi(0) \rangle = \vert 0(0) \rangle + e^{i\phi} \vert 1(0)\rangle$$
Now, the hamiltonian is changing with time, I know that at $t = T$, this state will evolve to:
$$\vert \Psi(T) \rangle = e^{i\Phi_0}\vert 0(T) \rangle + e^{i\phi}e^{i\Phi_1} \vert 1(T)\rangle$$
where $\Phi_i$ are the sum of dynamical and Berry phases.
The trouble I'm having is that the states $\vert 1(T) \rangle$ and $\vert 0(T) \rangle$ that I'm obtating through some numerics (Lanczos algorithm) are having an arbitrary phase that is not same for all diagonalisations (I have to do multiple diagonalizations for different values of $t$).
Because of this, I'm not exactly getting the evolution of the same state for all the times $t$.
I want the superposition state for all times $t$.
 A: I'll try a partial answer to get he ball rolling.
Under the adiabatic approximation, the phases are calculated as follows:

Under a slowly changing Hamiltonian $H ( t )$
with instantaneous eigenstates $ | n (t)\rangle $
and corresponding energies $E_n ( t )$, a quantum system
evolves from the initial state
$ |\psi(0)\rangle =\sum _{n}c_{n}(0)|n(0)\rangle $
to the final state
$ |\psi(t)\rangle =\sum _{n}c_{n}(t)|n(t)\rangle$
where the coefficients undergo the change of phase
$ c_{n}(t)=c_{n}(0)e^{i\theta _{n}(t)}e^{i\gamma _{n}(t)}$
with the dynamical phase
$ \theta_{m}(t)={\frac {-1}{\hbar }}\int _{0}^{t}E_{m}(t')dt'$
and geometric phase
$ \gamma _{m}(t)=i\int _{0}^{t}\langle m(t')| \partial_{t'} | {m}(t')\rangle dt'.$

The issues OP has identified with the random phase of numerical diagonalisation ruin the last expression. Specifically, when a random phase is introduced,
$\partial_t |m(t)\rangle$ is no longer defined. This is not just a numerical problem, it's fundamental - it is always possible (in the adiabatic limit) to introduce a U(1) gauge transformation $|n(t)\rangle \mapsto e^{-i\beta(t)} |n(t)\rangle$ that refers to the same physical state. $\beta(t)$ is any slowly-varying smooth function. In other words, the phase is not unique at a given "long time" point.
To proceed with Berry's argument, we need to specialise to the case that the Hamiltonian depends on a set of parameters $R$. In our case, we will slightly generalise the Hamiltonian:
$$ \mathcal{H}(\vec{h}) = \sum_i J \sigma_i \cdot \sigma_{i+1} + \vec{h}\cdot \sigma_i $$
The spectrum and eigenstates of $H$ are dependent only on the 3-vector $\vec{h}$, so are relabeled $E(h)$, $|n(h)\rangle$ respectively.
The Berry phase is $\gamma_n = \int_0^T i\langle n | \partial_{t} | n(t)\rangle dt = \int_\mathcal{C} i\langle n | \partial_{h^i} | n(t)\rangle dl^i$, where the right-hand expression is a line integral over in $h$-space between $h(t=0)$ and $h(t=T)$.
Only in the special case that $h$ returns to its original value at some later time does the Berry phase become gauge invariant, which is to say, experimentally relevant/measurable.
I now summarise the argument given here. Suppose $\mathcal{C}$ is a closed curve bounding an area $S$.
Define the Berry potential $A_i = i\langle n(h) | \frac{\partial}{\partial h^i} | n(h) \rangle $, such that $\gamma_n = \oint_\mathcal{C} \vec{A}\cdot d\vec{h}$. Apply Stokes' theorem:
$$\oint_\mathcal{C} \vec{A}\cdot d\vec{h} = \iint_S i d\left[\langle n(h) | \partial_i n(h) \rangle  dh^i\right]$$
$$ = \iint_S i d\left[\langle n(h) |\right]  \partial_i n(h) \rangle  dh^i + i \langle n(h) |d\left[ \partial_i n(h) \rangle \right]  dh^i $$
$$ = \sum_m \iint_S i \langle \partial_i n(h) |m(h)\rangle\langle m(h) | \partial_j n(h) \rangle dh^i \wedge dh^j = \iint_S \Omega_{ij}  dh^i \wedge dh^j $$
Then, using
$$ \langle n | \partial_j H | m \rangle + \langle n | H |\partial_j m \rangle = \delta_{mn} \partial_j E_m + E_m \langle n| \partial_j m \rangle $$
to eliminate dependence on $\partial_j m$, one can rewrite the Berry curvature as
$$\Omega_{ij} = i \sum_{m\neq n} \langle \partial_i n(h) |m(h)\rangle\langle m(h) | \partial_j n(h) \rangle$$
$$ = i \sum_{m\neq n} \frac{\langle n | \partial_i H | m \rangle \langle m | \partial_j H | n \rangle - (j\leftrightarrow i) }{(E_n-E_m)^2}$$
which is gauge invariant! $^\text{up to an overall, irrelevant factor of $2\pi$}$
Returning to your original situation, which corresponds to the case that the 3D $h(t)$ curve is constrained to the $z$ axis, we find the following -

*

*The area of $S$ is zero.

*The Berry curvature is given by $ \Omega^{\mu \nu} = i \sum_{m\neq n} \sum_i \frac{\langle n | \sigma_i^\mu | m \rangle \langle m | \sigma_i^\nu | n \rangle - ( \mu \leftrightarrow \nu) }{(E_n-E_m)^2}$.

*Unless the Berry curvature exhibits a singularity somewhere, the Berry phase is zero. However, it's likely that there will be a singularity somewhere in $h$ space.

Please don't take these calculations as gospel, there may be missing factors of 2 or $i$.
A: Your states have an arbitrary global phase due to $U(1)$ gauge symmetry. To remove it, you need to fix a common gauge across all your states. A convenient way to do this is to iteratively impose parallel transport, i.e. modifying the phase of your $|n_i\rangle$ state such that the condition $\text{Im}\langle n_{i-1}|n_i\rangle = 0$ is fulfilled. Note that in this case the (discretized) Berry phase will not be given as a sum of phase differences between your subsequent states along your closed loop in Hilbert space, but simply as the phase difference between your initial and final states, $\phi=\text{Im}\ln\langle n_{0}|n_N\rangle$.
