Susceptibilities and response functions

response function = susceptibility = (pure or mixed) second derivative of a (Helmholtz, Gibbs, etc.) free energy.

Magnetization (a first, not second derivative of a free energy) is not a response function as the free energy is not observable, so one cannot observe its response to a change of some variable.

It's been a while since this question was asked, but I think there's a big picture missing in these answers. The connection to probability theory will provide a robust framework to understand why @Arnold's first statement makes sense. Further, linear response theory (I discuss here Differentiating Propagator, Green's function, Correlation function, etc) is a formal way to work backward up the derivative chain. (i.e. predict magnetization $m$ given your response function/susceptibility $\chi$).

Part 1 - Average Magnetization The average magnetization is defined as, $$\left< m \right> = \sum_i m_i \times p_i $$. But what is $p_i$?? Let's use our thermodynamics class and consider the partition function. Which is the sum of all the Boltzmann factors and defines the probability measure for a system in thermal equilibrium, $$\mathcal{Z} = \sum_i e^{-\beta E(s_i)}$$ where, $\beta = \frac{1}{T}$, and $E(s_i)=E_i$ is the energy of the $i^{th}$ state. This is just an isolated system. The probability of having energy $E_i$ is defined by this probability measure as, $$p_i=\frac{e^{-\beta E_i}}{\mathcal{Z}}$$ Now let's turn on some external magnetic field $H_j$ which interacts with each of the $j$ particles in state $i$ via their magnetization $m_j(s_i)$. Here the magnetic field interacts with each particle magnetization so we use $m$. $$\mathcal{Z} = \sum_i e^{-\beta \left(E_i - \sum_j H_j m_j(s_i) \right)}$$ Note, that this is the right sign because $E - mH$ looks like the first law of thermo $U-pV$. We'll suppress the state dependence in $m_j(s_i)=m_j$ cause it's not relevant to our upcoming math. Anywho, The probability of being in state $s_i$ is defined given as, $$p_i=\frac{e^{-\beta E_i + \beta \sum_j H_j m_j}}{\mathcal{Z}}$$ NOW we can compute the average magnetization. $$\left< m \right> = \sum_i m_i \times p_i = \sum_i m_i \times \left(\frac{e^{-\beta E_i + \beta \sum_j H_j m_j}}{\mathcal{Z}} \right)$$ Rearrange terms, $$ \left< m \right> = \frac{1}{\mathcal{Z}} \sum_i m_k e^{-\beta E_i + \beta \sum_j H_j m_j} = \frac{T}{\mathcal{Z}} \frac{\partial}{\partial H_k} \sum_i e^{-\beta E_i + \beta \sum_j H_j m_j} = \frac{T}{\mathcal{Z}} \frac{\partial}{\partial H_k} \mathcal{Z}$$ So in this case we recognize that the partition function is the Moment Generating Function (MGF) (https://en.wikipedia.org/wiki/Moment-generating_function#Definition) of probability distribution function. The average value of the magnetization $\left< m \right>$ is the first moment (i.e the first derivate of $\mathcal{Z}$).

Part 2 - The Free Energy We can got one step further by using $\frac{1}{f[x]}\partial_x f(x) = \partial_x \log[f(x)]$. $$ \left< m \right> = T \frac{\partial}{\partial H_k} \log[\mathcal{Z}]$$ Recognize that $F=-T \log [Z]$ and all spins have the same coupling $H_j = H, \, \, \forall j$ $$ \left< m \right> = - \frac{\partial F}{\partial H} $$ If we plus in this value of $m$ into your final equation we recover that $$ -\frac{\partial}{\partial H} m = \frac{\partial}{\partial H} \frac{\partial F}{\partial H} = \frac{\partial^2 F}{\left( \partial H \right)^2} $$

This is not surprising!... "Why Dude!? Cause it seems pretty cool!!"... well I'll tell you! The partition function is the Moment Generating Function. The log of the MGF is the Cumulant Generating Function (https://en.wikipedia.org/wiki/Cumulant#Definition). This tells us that the second derivative of $F$ is the variance of the distribution / correlation in the system. $$\frac{\partial^2 F}{\left( \partial H \right)^2} = \left< m^2 \right> - \left< m_i \right>^2 $$ The susceptibility is by definition the correlation in the system (https://phys.libretexts.org/Bookshelves/Thermodynamics_and_Statistical_Mechanics/Book%3A_Statistical_Mechanics_%28Styer%29/09%3A_Strongly_Interacting_Systems_and_Phase_Transitions/9.04%3A_Correlation_Functions_in_the_Ising_Model)

$$\frac{\partial^2 F}{\left( \partial H \right)^2} = \left< m^2 \right> - \left< m_i \right>^2 = \chi$$

In total, I have just shown that (1) the partition function $\mathcal{Z}$ is the moment generating function and (2) the free energy, $F$ is the cumulant generating function. And this explains why the derivative look the way they do! Furthermore, an open question here is "how to interpret @Arnold's second statement", because a physicist can certainly measure the free energy in their simulations and experiments. A quick example from my field is that we use the free energy to study a phase transition in the nuclear strong force http://faculty.washington.edu/srsharpe/int07/Petreczky2.pdf While that heuristic may work for him, I'm not sure that's a broadly applicable.