  We derive the most general first order symmetry operator for the Dirac
equation coupled to arbitrary fluxes. Such an operator is given in terms of an
inhomogenous form omega which is a solution to a coupled system of first order
partial differential equations which we call the generalized conformal
Killing-Yano system. Except trivial fluxes, solutions of this system are
subject to additional constraints. We discuss various special cases of physical
interest. In particular, we demonstrate that in the case of a Dirac operator
coupled to the skew symmetric torsion and U(1) field, the system of generalized
conformal Killing-Yano equations decouples into the homogenous conformal
Killing-Yano equations with torsion introduced in [arXiv:0905.0722] and the
symmetry operator is essentially the one derived in [arXiv:1002.3616]. We also
discuss the Dirac field coupled to a scalar potential and in the presence of
5-form and 7-form fluxes.
