In the Ising model, the two-spin correlation function is $$ C(\vec{r}) = \langle \sigma_{\vec{r}_0+\vec{r}}\sigma_{\vec{r}_0}\rangle - \langle \sigma_{\vec{r}_0+\vec{r}}\rangle \langle \sigma_{\vec{r}_0} \rangle. $$ This quantity doesn't depend on $\vec{r}_0$ due to the translational invariance. When $r = |\vec{r}|$ is large compared to the lattice spacing, we expect the following approximate form $$ C(\vec{r}) \sim \exp(-r/\xi), $$ where $\xi$ is the correlation length.
Different directions on the lattice are not equivalent. For example, in the Ising model on the square lattice, there are two directions, say vertical and horizontal, along which neighboring spins interact. I see no reasons to think that other directions are equivalent to these two. In the anisotropic Ising model, vertical and horizontal directions are also not equivalent.
Then the correlation length $\xi$ should depend on the direction of $\vec{r}$. Is the analytical form of this dependence known at least for the square lattice? The Ising model is probably the most studied model of statistical physics, but I was not able to find corresponding formulas. So any references would be appreciated.
P.S. I know that in the scaling limit the Ising model becomes isotropic. The question above is for systems far enough from the critical point.