I want to get the eigenvectors of the Floquet operator and then derive them with respect to 'k' to be able to do the inner product that is in the Zak phase formula and then integrate over 'k' from -pi to pi. As I will get three eigenvectors, I expect to get three Zak phases.

I am writing code, and I expect help.

import numpy as np
from scipy.integrate import quad
from scipy.linalg import expm

# Define parameters
J1 = 1
J2 = 2
a = 0.25
T = 2

def floquet_operator(k):
    A1 = np.array([[0, J1 + a, 0],
                          [J1 + a, 0, J1],
                          [0, J1, 0]])
    A2 = np.array([
        [0, J2 * np.exp(-1j * k), 0],
        [J2 * np.exp(1j * k), 0, J2 * np.exp(-1j * k)],
        [0, J2 * np.exp(1j * k), 0]
    ], dtype=complex)
    
    exp_A1 = expm(-1j * A1 * (T / 2))
    exp_A2 = expm(-1j * A2 * (T / 2))
    return exp_A1 @ exp_A2