I think that the usual derivation of linear response theory should apply.
Assume that the effect of a field of magnitude $h$ applied
in the $i$ direction ($i$ being $x$ or $y$ or $z$)
can be represented as a perturbation term in the hamiltonian
$$
H = H_0 - h M_i
$$
where
$H_0$ is your original hamiltonian,
and $M_i$ is the total system magnetization in the $i$ direction,
calculated as a sum over all the individual spins.
Then the measured value of a possibly different magnetization component
$M_j$ in the perturbed ensemble is
$$
\langle M_j \rangle_h = \frac{\int M_j\, \exp[-\beta(H_0-h M_i)]}{Q}
\approx \frac{\int M_j\, \bigl(1+\beta h M_i\bigr) \exp[-\beta H_0]}{Q}
\approx \beta h \langle M_i M_j\rangle_0
$$
where $\langle\ldots\rangle_0$
and $\langle\ldots\rangle_h$ represent averages in the
unperturbed and perturbed systems respectively.
Here $Q$ is the unperturbed Hamiltonian and I have assumed for simplicity
that $\langle M_i\rangle_0=0$ for all $i$.
So, in the above equation, we can identify the proportionality coefficient
between $\langle M_j \rangle_h$ and $h$,
and this would lead to a formula very closely analogous to the one you are
using for diagonal elements of $\chi$:
$$
\chi_{ji} = \beta \langle M_i M_j\rangle_0
$$
If the average magnetization in the unperturbed ensemble does not vanish,
then $M_i$ should be replaced by $\Delta M_i=M_i-\langle M_i\rangle_0$ and
similarly for $M_j$, as you have in your expression.
For your classical Heisenberg model, of course, the different magnetization components will be uncorrelated and the off diagonal elements of $\chi$ should be zero.
Only if your Hamiltonian has the appropriate reduced symmetry should you get nonvanishing off diagonal components and, in any case, it is expected that $\chi_{ij}=\chi_{ji}$.
It might be advisable to check directly,
by also doing perturbed simulations
in which small fields $h$ are actually applied in each of the three coordinate directions, separately,
and in each case you measure all three components of the average
magnetization $\langle M_j \rangle_h$,
calculating the relevant three proportionality constants $\chi_{ji}$.
It's prudent to try several different values of $h$
and extrapolate to low $h$,
to confirm that these perturbed simulations are in the linear regime.
Hopefully the results will agree with the fluctuation
expression.
