Heisenberg ferromagnet in continuum limit I consider the case of the simple, say 2D, Heisenberg ferromagnet with exchange interaction between the nearest neighbors. The Hamiltonian is:
$$H = -J \sum_{<ij>} \mathbf S_i \mathbf S_j,$$
where $\mathbf S_{i,j}$ are classical spin vectors. I often encounter the following chain of arguments while transiting to the continuum limit. First, we express the scalar product of spins in terms of the small angle between them: 
$$\mathbf S_i \mathbf S_j  \approx S_0^2  \cos \theta_{ij} \approx S_0^2 -\frac{1}{2}S_0^2 \theta_{ij}^2,$$ 
where $S_0=|\mathbf S|$ and $\theta_{ij}$ is the angle between $i$th and $j$th spins. The Hamiltonian then turns into:
$$H =const +\frac{J}{2}S_0^2  \sum_{<ij>}\theta_{ij}^2 \rightarrow_{cont. \ lim.} \rightarrow \frac{J}{2}S_0^2  \int d^2 r (\nabla \theta(\mathbf r))^2.$$
The question is about the last transition: what is the correct way to perform the transition to the continuum limit rigorously?

Here I provide my attempt to obtain the continuum-case expression for a square lattice, which obviously has a flaw and leads to an incorrect answer. However, I don't see where the problem is.
First, I substitute an integral instead of the sum: $\sum_{<ij>} \rightarrow \int \frac{dx dy}{a^2}$, where $a$ is the distance between the neighbors. The expression for the angle may approximately be rewritten as 
$$\theta_{ij} \rightarrow \theta(\mathbf r+\Delta \mathbf r) - \theta(\mathbf r) \approx (\Delta\mathbf r, \nabla )\theta(\mathbf r).$$
Next, I need to write the square of it in the Hamiltonian:
$$H \rightarrow  \frac{J}{2} \int  \frac{d^2 r}{a^2} \left [  (\Delta \mathbf r, \nabla )\theta( \mathbf r)\right ]^2 =  \frac{J}{2} \int  \frac{d^2 r}{a^2} \left [  a_x \frac{\partial \theta }{\partial x} + a_y \frac{\partial \theta }{\partial x}\right ]^2 = \\ =\frac{J}{2}\int  \frac{d^2 r}{a^2} \left (  a_x^2 \left(\frac{\partial \theta }{\partial x} \right )^2+ a_y^2 \left(\frac{\partial \theta }{\partial y} \right )^2 + 2a_x a_y \frac{\partial \theta }{\partial x} \frac{\partial \theta }{\partial y}\right ),$$
where  $a_{x,y}=\Delta  r_{x,y}$. The last term in the round brackets clearly is absent in the correct expression. Where am I wrong in my calculations?
 A: Remember that the Hamiltonian involves a sum over all pairs of neighbouring sites. Assume that the sites are located on a square lattice so that their positions are given by $\vec r=a(n_x\vec e_x+n_y\vec e_y)$ where the n's are integer. The Hamiltonian reads
    $$\eqalign{
H&=-J\sum_{\vec r}  \big(\vec S_{\vec r+a\vec e_x}.\vec S_{\vec r}
+\vec S_{\vec r+a\vec e_y}.\vec S_{\vec r}\big)\cr
&=-JS^2\sum_{\vec r} \Big[\cos\big(\theta_{\vec r+a\vec e_x}
-\theta_{\vec r}\big)+\cos\big(\theta_{\vec r+a\vec e_y}
-\theta_{\vec r}\big)\Big]\cr
&\simeq -JS^2\sum_{\vec r} \Big[\cos\big(a\partial_x\theta\big)
+\cos\big(a\partial_y\theta\big)\Big]\cr
&\simeq -{JS^2\over 2}\sum_{\vec r} \Big[1-a^2
\big(\partial_x\theta\big)^2-a^2
\big(\partial_y\theta\big)^2\Big]\cr
&={\rm Cst}+{JS^2a^2\over 2}\sum_{\vec r} ||\vec\nabla\theta||^2
}$$
On an hexagonal lattice, one needs to distinguish two sublattices. On one of them, the three neighbours are at $\vec a=\vec e_x$, $\vec b=-{1\over 2}\vec e_x+{\sqrt 3\over 2}\vec e_y$ and $\vec c=-{1\over 2}\vec e_x-{\sqrt 3\over 2}\vec e_y$. The energy on site $\vec r$ then reads
    $$\eqalign{
&-J\big(\vec S_{\vec r+a\vec a}.\vec S_{\vec r}
+\vec S_{\vec r+a\vec b}.\vec S_{\vec r}
+\vec S_{\vec r+a\vec c}.\vec S_{\vec r}\big)\cr
&=-JS^2\Big[\cos\big(\theta_{\vec r+a\vec a}
-\theta_{\vec r}\big)+\cos\big(\theta_{\vec r+a\vec b}
-\theta_{\vec r}\big)+\cos\big(\theta_{\vec r+a\vec c}
-\theta_{\vec r}\big)\Big]\cr
&\simeq -JS^2\Big[\cos\big(a\partial_x\theta\big)
+\cos\Big(-{a\over 2}\partial_x\theta+a{\sqrt 3\over 2}
\partial_y\theta\Big)+\cos\Big(-{a\over 2}\partial_x\theta
-a{\sqrt 3\over 2}\partial_y\theta\Big)\Big]\cr
&\simeq -{JS^2\over 2}\sum_{\vec r} \Big[1-a^2
\big(a\partial_x\theta\big)^2-a^2\Big(-{1\over 2}\partial_x\theta
+{\sqrt 3\over 2}\partial_y\theta\Big)^2-a^2\Big(-{1\over 2}
\partial_x\theta-{\sqrt 3\over 2}\partial_y\theta\Big)^2\Big]\cr
&={\rm Cst}+{JS^2a^2\over 2}\times  {3\over 2}||\vec\nabla\theta||^2
}$$
The only difference with the square lattice is a geometric factor $3/2$.
