  This work pioneers the quantization of primordial fermion perturbations in
hybrid Loop Quantum Cosmology (LQC). We consider a Dirac field coupled to a
spatially flat, homogeneous, and isotropic cosmology, sourced by a scalar
inflaton, and treat the Dirac field as a perturbation. We describe the
inhomogeneities of this field in terms of creation and annihilation variables,
chosen to admit a unitary evolution if the Dirac fermion were treated as a test
field. Considering instead the full system, we truncate its action at quadratic
perturbative order and construct a canonical formulation. In particular this
implies that, in the global Hamiltonian constraint of the model, the
contribution of the homogeneous sector is corrected with a quadratic
perturbative term. We then adopt the hybrid LQC approach to quantize the full
model, combining the loop representation of the homogeneous geometry with the
Fock quantization of the inhomogeneities. We assume a Born-Oppenheimer ansatz
for physical states and show how to obtain a Schr\"odinger equation for the
quantum evolution of the perturbations, where the role of time is played by the
homogeneous inflaton. We prove that the resulting quantum evolution of the
Dirac field is indeed unitary, despite the fact that the underlying homogeneous
geometry has been quantized as well. Remarkably, in such evolution, the fermion
field couples to an infinite sequence of quantum moments of the homogeneous
geometry. Moreover, the evolved Fock vacuum of our fermion perturbations is
shown to be an exact solution of the Schr\"odinger equation. Finally, we
discuss in detail the quantum backreaction that the fermion field introduces in
the global Hamiltonian constraint. For completeness, our quantum study includes
since the beginning (gauge-invariant) scalar and tensor perturbations, that
were studied in previous works.
