Eigenvalues in Floquet theory

It sounds to me like you are getting non-orthogonal eigenvectors. Are you sure that your Floquet Hamiltonian is exactly Hermitian? Are you calling an eigensolver that is specifially for Hermitian matrices that is guaranteed to produce eigenvectors that are very close to orthongonal (roughly within machine precision).

I think in the eigen package, you want to use SelfAdjointEigenSolver.

The other obvious issues would be you are extracting rows instead of columns, or indexing starting at 1 and not at 0.

By the way, be careful how you computing the Floquet Hamiltonian. If you use a generic matrix logarithm the result might be very far from Hermitian.

I should give a disclaimer that I'm referencing my own work, but it pretty much addresses this exact point. Half of the problems you want to solve is with the code implementation and half with the theory itself.

I will assume that you are using a Fourier space representation, and a large enough cutoff that the Floquet Hamiltonian solutions converge, i.e.: $$H^F=begin{bmatrix}ddots&vdots ldots&H^{(0)}+omega&H^{(-1)}&H^{(-2)} &H^{(1)}&H^{(0)}&H^{(-1)} &H^{(2)}&H^{(1)}&H^{(0)}-omega&ldots &&&vdots&ddotsend{bmatrix}$$ , where: $$H(t)=sum_kH^{(k)}e^{-ikomega t}$$

Then I recommend you use a Sparse matrix solver and find the eigenpairs of a Brillouin-zone in the middle of the quasi-energy region (usually $[-frac{omega}{2},frac{omega}{2}]$ will do just fine). Otherwise you will diverge from the real Floquet eigenstate solutions. I have more experience using MATLAB and it provides most of the tools you need to do so, but for C++ you can check Intel MKL or OneMKL for a sparse eigen solver. Then each eigenvector and eigevalue are ordered accordingly, just make sure you get the whole eigenvector. If you don't have any resonant solutions or don't care about precise eigenstate definition and ordering, than you can just proceed with using those, e.g. calculating the transition probability for each Fourier component.

Otherwise you should use the average energy to re-order and redefine the eigenstates (for the theory you can check the work I linked at the beginning). But in short you need to calculate the average energy matrix of the resonant eigenstates: $$H_{mn}=int_0^T<Psi_m(t)|H(t)|Psi_n(t)>dt$$ It's basically the same matrix as the Floquet Hamiltonian without the $komega$. Then calculate the eigenpairs again and chain them to get the proper eigenstates. Outside the resonant region, just calculate the expectation value. Combine these solutions and you have well-order eigenstates.