Basis state of non-interacting fermions

Question

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$?