  Consider an elliptic parameter $k$; we introduce a family of $Z^u$-Dirac
operators $(\mathsf{K}(u))_{u\in\Re(\mathbb{T}(k))}$, relate them to the
$Z$-massive Laplacian of [BdTR17b], and extend to the full $Z$-invariant case
the results of Kenyon [Ken02] on discrete holomorphic and harmonic functions,
which correspond to the case $k=0$. We prove, in a direct statistical mechanics
way, how and why the $Z^u$-Dirac and $Z$-massive Laplacian operators appear in
the $Z$-invariant Ising model, considering the case of infinite and finite
isoradial graphs. More precisely, consider the dimer model on the Fisher graph
${\mathsf{G}}^{\scriptscriptstyle{\mathrm{F}}}$ arising from a $Z$-invariant
Ising model. We express coefficients of the inverse Fisher Kasteleyn operator
as a function of the inverse $Z^u$-Dirac operator and also as a function of the
$Z$-massive Green function; in particular this proves a (massive) random walk
representation of important observables of the Ising model. We prove that the
squared partition function of the Ising model is equal, up to a constant, to
the determinant of the $Z$-massive Laplacian operator with specific boundary
conditions, the latter being the partition function of rooted spanning forests.
To show these results, we relate the inverse Fisher Kasteleyn operator and that
of the dimer model on the bipartite graph
${\mathsf{G}}^{\scriptscriptstyle{\mathrm{Q}}}$ arising from the XOR-Ising
model, and we prove matrix identities between the Kasteleyn matrix of
${\mathsf{G}}^{\scriptscriptstyle{\mathrm{Q}}}$ and the $Z^u$-Dirac operator,
that allow to reach inverse matrices as well as determinants.
