I want to write down a spin hamiltonian for $N$ sites (let's say for starters each site is spin = 1/2) that yields particular spin correlation functions, $$C_{ij} = \langle \hat{\vec{S}}_{i} \cdot \hat{\vec{S}}_{j} \rangle - \langle \hat{\vec{S}}_{i} \rangle \cdot \langle \hat{\vec{S}}_{j} \rangle$$ and a particular spectral gap; that is, if $E_0$ and $E_1$ are the first non-equal eigenvalues, I want to reproduce also the difference $$\varepsilon = E_1 - E_0.$$ This is sort of an inverse problem, I guess. The question then, which is maybe more mathematical than physical, is: Can I find coefficients $J_{ij}$ such that the hamiltonian $$\hat{H_1} = \sum_{i,j} J_{ij}\hat{\vec{S}}_{i} \cdot \hat{\vec{S}}_{j}$$ has a ground state that yields the prescribed values $C_{ij}$ and has the spectral gap $\varepsilon$? If such simple hamiltonian is not enough, can I achieve it with paramaters $J_{ij}$ and $\delta$ in a hamiltonian that looks like $$\hat{H_2} = \sum_{i,j} J_{ij}\hat{\vec{S}}_{i} \cdot \hat{\vec{S}}_{j} + \delta \hat{S}^2. $$ The most general hamiltonians I have tried, inspired by the AKLT model, have the shape $$\hat{H_3} = \sum_{i,j} J_{ij}\hat{\vec{S}}_{i} \cdot \hat{\vec{S}}_{j} + \delta \hat{S}^2 + \sum_{i,j} Q_{ij}\left( \hat{\vec{S}}_{i} \cdot \hat{\vec{S}}_{j} \right)^2.$$
In some simple cases, ranging over the parameter space I can achieve the desired fitting or something sufficiently close. But in general I am at a lost and ranging over all the possible combinations of values of the parameters is usually unfeasible. Therefore, is there any clever way to, given correlations $C_{ij}$ and gap $\varepsilon$, find parameters $J_{ij}, Q_{ij}$ and $\delta$ such that the hamiltonians like $\hat{H_1}, \hat{H_2}$ or $\hat{H_3}$ reproduce the $\varepsilon$ and the $C_{ij}$ (for the ground state)?
