I am trying to calculate the periodic dynamics of many-body systems (spin-$1/2$ $XY$) Hamiltonian, where,

\begin{equation}

H_1 = \sum_{i=1}^{N-1}(\sigma^{x}_{i}\sigma^{x}_{i+1}+\sigma^{y}_{i}\sigma^{y}_{i+1})+\sum_{i}h_{i}\sigma^{z}

\end{equation}

and

\begin{equation}

H_2 = \sum_{i}h_{i}\sigma^{z}

\end{equation}

where, $h_i$ is the disorder field and $\sigma$'s are the Pauli spin-$1/2$ matrices. The dynamics take as

\[

H(t) =

\begin{cases}

H_1 & \text{if $0<T/2$} \\

H_2 & \text{if $T/2<t<T$}.

\end{cases}

\]

Here $T$ is the time-period of the dynamics. The Floquet operator of evolution is given by

$U_{F}=e^{-iH_2 T/2}e^{-i H_1 T/2}$.

The Above Hamiltonian can be written as the fermionic operators $c^{\dagger}(c)$ where $c^{\dagger}(c)$ are fermionic creation (annihilation) operator as

\begin{equation}

H(t)=\sum_{mn}C^{\dagger}_{m}M_{mn}(t)C_{n}

\end{equation}

where $C^{\dagger}=[c^{\dagger} c]$.

The basis state is given by

\begin{equation}

c^{\dagger}_{i_1}c^{\dagger}_{i_2}\ldots c^{\dagger}_{i_N}|0\rangle.

\end{equation}

Where $i_1,i_2,\ldots,i_N$ are set of integers such that $1\leq i_1, i_2,\ldots, i_N\leq N$.

The time evolved state is given by

\begin{equation}

|\psi(t)\rangle=e^{-iC^{\dagger}M C T}c^{\dagger}_{i_1}c^{\dagger}_{i_2}\ldots c^{\dagger}_{i_N}|0\rangle

\end{equation}

My aim is to calculate the average enegrgy $e(T)=\frac{1}{2}\langle \Psi(t)|(H_1+H_2)|\Psi(t)\rangle$.From the Hamiltonian, $M$ can be obtained. How to understand the basis state and how to do numerical with the state $c^{\dagger}_{i_1}c^{\dagger}_{i_2}\ldots c^{\dagger}_{i_N}|0\rangle$?