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