The coupled mode space NEGF method for computing transistor characteristics involves expanding the electronic wavefunction in a mode space basis
$$\Psi(x,y,z) = \sum_n\phi_n(x)\xi_n(y,z;x)$$
where $\xi_n(y,z;x)$ is the nth mode at position x (the x dimension being the direction of transport).
The method involves computing coupling between modes along the transport direction. E.g. coupling coefficients like
$$b_{mn}(x) = \langle\xi_m|\frac{\partial}{\partial x}|\xi_n\rangle$$
The numerical algorithm will divide the transistor into $X$ slices $\Delta x$ apart, and compute the modes on each slice. The issue is my numerical solver will sometimes compute $\xi_n$ as the nth eigenvalue, and sometimes compute $-\xi_n$, so phases might not be consistent across each slice. This changing phase will result in a large $\frac{\partial}{\partial x}|\xi_n\rangle$
My question is: Do the phases need to be consistent across each slice? If they do, how do NEGF algorithms ensure they are consistent?