\begin{document}

\title{Two Groups in a Curie-Weiss Model with Heterogeneous Coupling}

\author{Werner Kirsch and Gabor Toth\\
 Fakult\"at f\"ur Mathematik und Informatik \\
 FernUniversit\"at Hagen, Germany}
\maketitle
\begin{abstract}
Following our paper ``Two Groups in a Curie-Weiss Model'' \cite{KT1},
where we showed bivariate laws of large numbers and central limit
theorems for two groups of spins with a homogeneous coupling matrix,
we now show similar results when the coupling matrix is heterogeneous
and positive definite, under the assumption that the coupling both
within groups and between groups is weak. 
\end{abstract}

\section{Introduction}

The Curie-Weiss model is probably the easiest model of magnetism which
shows a phase transition between a diamagnetic and a ferromagnetic
phase. In this model the spins can take values in $\{-1,1\}$ (or
up/down), each spin interacts with all the others in the same way.
More precisely, for finitely many spins $(X_{1},X_{2},\ldots,X_{N})\in\{-1,1\}$
the energy of the spins is given by 
\begin{align}
H~=~H(X_{1},\ldots,X_{N})~:=~-\frac{J}{2N}\,\big(\sum_{j=1}^{N}\,X_{j}\big)^{2}\,,
\end{align}
where $J$ is a positive real number.

Consequently, in the `canonical ensemble' with inverse temperature
$\beta\geq0$ the probability of a spin configuration is given by
\begin{align}
\mathbb{P}\big(X_{1}=x_{1},\ldots,X_{N}=x_{N}\big)~:=~Z^{-1}\;e^{-\beta H(x_{1},\ldots,x_{N})}\label{eq:Pbeta}
\end{align}
where $x_{i}\in\{-1,1\}$ and $Z$ is a normalisation constant which
depends on $N$, $J$ and $\beta$. Since only the product of $\beta$
and $J$ occurs in (\ref{eq:Pbeta}) we may set $J=1$ without loss
of generality.

The quantity 
\begin{align}
S_{N}~=~\sum_{j=1}^{N}X_{j}
\end{align}
is called the (total) magnetisation. It is well known (see e.\,g.
Ellis \cite{Ellis} or \cite{MM}) that the Curie-Weiss model has
a phase transition at $\beta=1$ in the following sense 
\begin{align}
\frac{1}{N}\,S_{N}~\Longrightarrow~\frac{1}{2}\,(\delta_{-m(\beta)}+\delta_{m(\beta)})\label{eq:lln}
\end{align}
where $\Rightarrow$ denotes convergence in distribution, $\delta_{x}$
the Dirac measure in $x$.

For $\beta\leq1$ we have $m(\beta)=0$ which is the unique solution
of 
\begin{align}
\tan(\beta x)=x\label{eq:mbeta}
\end{align}
for this case.

If $\beta>1$ equation (\ref{eq:mbeta}) has exactly three solutions
and $m(\beta)$ is the unique positive one.

Equation (\ref{eq:lln}) is a substitute for the law of large numbers
for i.i.d. random variables.

Moreover, for $\beta<1$ there is a central limit theorem, i.\,e.
\begin{align}
\frac{1}{\sqrt{N}}\,S_{N}~\Longrightarrow~\mathcal{N}(0,\frac{1}{1-\beta})
\end{align}

In this paper we form out of $N$ Curie-Weiss spins two disjoint groups
$X_{1},\ldots,X_{\tilde{N}_{1}}$ and $Y_{1},\ldots,Y_{\tilde{N}_{2}}$
with $\tilde{N}_{1}+\tilde{N}_{2}=N$. We let $\tilde{N}_{1}$ and
$\tilde{N}_{2}$ depend on $N$ in such a way that both $\tilde{N}_{1}$
and $\tilde{N}_{2}$ go to infinity as $N$ does. We may want to study
situations in which one or both groups of spins become very small
in relation to the entire set of $N$ spin. Therefore, we define $\{1,2,\ldots,N_{\nu}\}\subset\{1,2,\ldots,\tilde{N}_{\nu}\}$
for $\nu=1,2$. We set 
\begin{equation}
\alpha_{1}=\lim\frac{N_{1}}{N}\qquad\qquad\alpha_{2}=\lim_{N\to\infty}\frac{N_{2}}{N}
\end{equation}
and assume that these limits exist. Except for section \ref{sec:Sublin},
we will also assume that $\alpha_{\nu}>0$ for both groups.

The two groups $\{1,2,\ldots,\tilde{N}_{\nu}\}$ interact with each
other via a heterogeneous coupling matrix $J=\left[\begin{array}{cc}
J_{1} & \bar{J}\\
\bar{J} & J_{2}
\end{array}\right]$, which has all positive entries and is positive definite. In addition,
we shall assume unless otherwise noted that the interactions within
groups and between groups is weak as defined by

\begin{align}
J_{1} & <\frac{1}{\alpha_{1}},\label{eq:condJ1-1}\\
J_{2} & <\frac{1}{\alpha_{2}},\label{eq:condJ2-1}\\
\bar{J}^{2} & <(\frac{1}{\alpha_{1}}-J_{1})(\frac{1}{\alpha_{2}}-J_{2}).\label{eq:condJ3-1}
\end{align}
We shall refer to these conditions as the `high temperature regime'.
Note that if we use the symbol $\alpha$ for the diagonal $2\times2$
matrix with entries $\alpha_{1}$ and $\alpha_{2}$, we can formulate
these conditions equivalently in matrix form: the matrix

\[
J^{-1}-\alpha
\]
is positive definite if and only if we are in the high temperature
regime.

Let throughout this article $\Delta=\det J=J_{1}J_{2}-\bar{J}^{2}>0$.

Consequently, the energy of the spins is given by $H~:=~H(X,Y):=$

\begin{align}
-\frac{1}{2N}\left[J_{1}\big(\sum_{j=1}^{\tilde{N}_{1}}X_{j}\big)^{2}+J_{2}\big(\sum_{j=1}^{\tilde{N}_{2}}Y_{j}\big)^{2}+2\bar{J}\sum_{i=1}^{\tilde{N}_{1}}\sum_{j=1}^{\tilde{N}_{2}}X_{i}Y_{j}\right].
\end{align}

We consider the asymptotic behaviour of the two-dimensional random
variables 
\begin{equation}
\big(\sum_{i=1}^{N_{1}}\,X_{i}\,,\,\sum_{j=1}^{N_{2}}\,Y_{j}\,\big)
\end{equation}
as $N$ goes to infinity.

We prove 
\begin{thm}
\label{LLN} Under the above assumptions, we have

\[
\big(\frac{1}{N_{1}}\sum_{i=1}^{N_{1}}X_{i},\frac{1}{N_{2}}\sum_{j=1}^{N_{2}}Y_{j}\big)~\underset{N\to\infty}{\Longrightarrow}(0,0).
\]
\end{thm}

Above '$\Longrightarrow$' denotes convergence in distribution of
the $2$-dimensional random variable on the left hand side.

We also have a central limit theorem. 
\begin{thm}
\label{CLT} Under the above assumptions, we have 
\begin{equation}
(\frac{1}{\sqrt{N_{1}}}\sum_{i=1}^{N_{1}}X_{i},\frac{1}{\sqrt{N_{2}}}\sum_{j=1}^{N_{2}}Y_{j})~\underset{N\to\infty}{\Longrightarrow}~\mathcal{N}\big((0,0),C),\label{eq:CLT}
\end{equation}
where the covariance matrix $C$ is given by 
\begin{align}
C & =\frac{1}{(1-\alpha_{1}J_{1})(1-\alpha_{2}J_{2})-\alpha_{1}\alpha_{2}\bar{J}^{2}}\left[\begin{array}{cc}
1-\alpha_{2}J_{2} & \sqrt{\alpha_{1}\alpha_{2}}\bar{J}\\
\sqrt{\alpha_{1}\alpha_{2}}\bar{J} & 1-\alpha_{1}J_{1}
\end{array}\right]\label{eq:covariance}
\end{align}
\end{thm}

We prove the two theorems in section \ref{sec:ProofsLIT}. Thereafter,
we take a look at what these results imply when there is no interaction
between the two groups $(\bar{J}=0)$ in section \ref{sec:Indep},
and when there is sublinear population growth $(\alpha_{1}=0)$ in
section \ref{sec:Sublin}.

Finally, in section \ref{sec:Special-Case}, we discuss a special
case with low temperature. There is spontaneous magnetisation and
hence the central limit theorem does not hold, similarly to the one-dimensional
model when $\beta>1$. 
\begin{thm}
\label{thm:SpecialCase}For $J_{1}=J_{2}=J$, $\alpha_{1}=\alpha_{2}=\alpha>0$,
and $J+\bar{J}>1/\alpha$, we have

\[
\big(\frac{1}{N_{1}}\sum_{i=1}^{N_{1}}X_{i},\frac{1}{N_{2}}\sum_{j=1}^{N_{2}}Y_{j}\big)~\underset{N\to\infty}{\Longrightarrow}\frac{1}{2}\big(\delta_{(-m^{*},-m^{*})}+\delta_{(m^{*},m^{*})}\big),
\]
where $m^{*}$ is the positive solution to the equation in m

\[
\frac{J-\bar{J}}{\Delta}m-\sqrt{2}\alpha\tanh(\frac{m}{\sqrt{2}})=0.
\]
\end{thm}

We mention that the Curie-Weiss model is also used to model the behaviour
of voters who have the choice to vote `Yea' (spin=1, say) or `Nay'
(spin=-1) (see \cite{Penrose}).

In the proof of the results we employ the moment method (see e.\,g.
\cite{Breiman} or \cite{MM}). Thus, to show the convergence in distribution
of a sequence $(X_{n},Y_{n})$ of two-dimensional random variables
to a measure $\mu$ on $\mathbb{R}^{2}$ we prove that 
\begin{align}
\mathbb{E}\Big(X_{n}^{K}\cdot Y_{n}^{L}\Big)~\longrightarrow~\int x^{K}y^{L}\,\mu(dx,dy)\label{eq:momcon}
\end{align}
for all $K,L\in\mathbb{N}$.

Equation (\ref{eq:momcon}) implies convergence in distribution if
the moments of $\mu$ grow only moderately, namely if for some constant
$A$ and $C$ and all $K,L$ 
\begin{equation}
\int|x|^{K}|y|^{L}\,\mu(dx,dy)~\leq~A\,C^{K+L}\,(K+L)!
\end{equation}
holds.

While writing up this article, we learned that L\"owe and Schubert have
obtained closely related results by rather different methods \cite{LS}.
After uploading the first version of this article, we became aware
that this kind of question had already been considered by Micaela
Fedele and Pierluigi Contucci in \cite{FC,Fedele}. These papers deal
with an arbitrary number of groups and the methods used there are
quite different from ours. We are grateful to Francesca Collet for
bringing the articles to our attention.

\section{Preparation}

To use the method of moments we have to evaluate sums of the form
\begin{align}
 & \mathbb{E}\left[\left(\sum_{i=1}^{N_{1}}X_{i}\right)^{K}\left(\sum_{j=1}^{N_{2}}Y_{j}\right)^{L}\right]\nonumber \\
=~ & \sum_{i_{1},\ldots,i_{K}}\sum_{j_{1},\ldots,j_{L}}\,\mathbb{E}\Big(X_{i_{1}}\cdot X_{i_{2}}\cdot\ldots\cdot X_{i_{K}}\cdot Y_{j_{1}}\cdot Y_{j_{2}}\cdot\ldots\cdot Y_{j_{L}}\Big)\,.
\end{align}
To do the book-keeping for these huge sums we introduce a few combinatorial
concepts taken from \cite{MM}. 
\begin{defn}
We define a multiindex $\underline{i}=(i_{1},i_{2},\ldots,i_{L})\in\{1,2,\ldots,N\}^{L}$. 
\begin{enumerate}
\item For $j\in\{1,2,\ldots,N\}$ we set\\

\[
\nu_{j}(\underline{i})=|\{k\in\{1,2,\ldots,L\}|i_{k}=j\}|,
\]
where $|M|$ denotes the number of elements in the set $M$. 
\item For $l=0,1,\ldots,L$ we define\\

\[
\rho_{l}(\underline{i})=|\{j|\nu_{j}(\underline{i})=l\}|
\]
and\\

\[
\underline{\rho}(\underline{i})=(\rho_{1}(\underline{i}),\ldots,\rho_{L}(\underline{i})).
\]
\end{enumerate}
\end{defn}

The numbers $\nu_{j}(\underline{i})$ represent the multiplicity of
each index $j\in\{1,2,\ldots,N\}$ in the multiindex $\underline{i}$,
and $\rho_{l}(\underline{i})$ represents the number of indices in
$\underline{i}$ that occur exactly $l$ times. We shall call such
$\underline{\rho}(\underline{i})$ profile vectors. 
\begin{lem}[Lemma 3.8]
\label{lem:sum-l-rho}For all $\underline{i}=(i_{1},i_{2},\ldots,i_{L})\in\{1,2,\ldots,N\}^{L}$
we have\\
 $\sum_{l=1}^{L}l\rho_{l}(\underline{i})=L$. 
\end{lem}

\begin{defn}
Let $\underline{r}=(r_{1},\ldots,r_{L})$, $\sum_{l=1}^{L}lr_{l}=L$,
be a profile vector. We define

\[
w_{L}(\underline{r})=|\{\underline{i}\in\{1\ldots,N\}^{L}|\underline{\rho}(\underline{i})=\underline{r}\}|
\]
to represent the number of multiindices $\underline{i}$ that have
a given profile vector $\underline{r}$. 
\end{defn}

We now define the set of all profile vectors for a given $L\in\mathbb{N}$. 
\begin{defn}
Let $\Pi^{(L)}=\{\underline{r}\in\{0,1,\ldots,L\}^{L}|\sum_{l=1}^{L}lr_{l}=L\}$.
Some important subsets of $\Pi^{(L)}$ are $\Pi_{k}^{(L)}=\{\underline{r}\in\Pi^{(L)}|r_{1}=k\}$,
$\Pi^{0(L)}=\{\underline{r}\in\Pi^{(L)}|r_{l}=0\text{ for all }l\geq3\}$
and $\Pi^{+(L)}=\{\underline{r}\in\Pi^{(L)}|r_{l}>0\text{ for some }l\geq3\}$. 
\end{defn}

\begin{prop}
\label{thm:comb-coeff-multiindex}For $\underline{r}\in\Pi^{(L)}$
set $r_{0}=N-\sum_{l=1}^{L}r_{l}$, then\\

\[
w_{L}(\underline{r})=\frac{N!}{r_{1}!r_{2}!\ldots r_{L}!r_{0}!}\frac{L!}{1!^{r_{1}}2!^{r_{2}}\cdots L!^{r_{L}}}.
\]
\end{prop}


\section{\label{sec:ProofsLIT}Proofs of the High Temperature Results}

As a reminder, we are studying a Curie-Weiss model with two disjoint
groups of sizes $N_{1}$ and $N_{2}$ and a positive definite heterogeneous
coupling matrix $J=\left[\begin{array}{cc}
J_{1} & \bar{J}\\
\bar{J} & J_{2}
\end{array}\right]$, which has all positive entries.

We define the inverse matrix

\[
L=\left[\begin{array}{cc}
L_{1} & -\bar{L}\\
-\bar{L} & L_{2}
\end{array}\right]=\frac{1}{J_{1}J_{2}-\bar{J}^{2}}\left[\begin{array}{cc}
J_{2} & -\bar{J}\\
-\bar{J} & J_{1}
\end{array}\right]=J^{-1}.
\]

We are interested in estimating correlations $\mathbb{E}(X_{1}\cdots X_{K}Y_{1}\cdots Y_{L})$
and we use Laplace's Method. For any voter configuration $X$, we
define the function $h(S_{1},S_{2})$ $=\frac{1}{2N}(S_{1},S_{2})J\left(\begin{array}{c}
S_{1}\\
S_{2}
\end{array}\right)=\frac{1}{2N}s'Js$, where $S_{1}=\sum_{i=1}^{N_{1}}X_{i}$ and $S_{2}=\sum_{j=1}^{N_{2}}Y_{j}$,
and want to calculate the expression $e^{h(S_{1},S_{2})}$. We are
using a lower case $h$ to distinguish this function of $(S_{1},S_{2})$
from the Hamiltonian $H$ that is a function of the entire voter configuration
$X$.

For a $2\times2$ matrix $L$ and a point $x_{0}\in\mathbb{R}^{2}$
we can use the following equality to express a value of the exponential
function as an integral:

\[
e^{\frac{x_{0}'Lx_{0}}{2}}=\frac{1}{2\pi\sqrt{\det L}}\int_{\mathbb{R}^{2}}e^{-\frac{x'Lx}{2}}e^{-x'Lx_{0}}\mathrm{d}x.
\]

The above equality holds because of

\begin{align*}
2\pi\sqrt{\det L} & =\int_{\mathbb{R}^{2}}e^{-\frac{x'Lx}{2}}\mathrm{d}x\\
 & =\int_{\mathbb{R}^{2}}e^{-\frac{(x-x_{0})'L(x-x_{0})}{2}}\mathrm{d}x\\
 & =\int_{\mathbb{R}^{2}}e^{-\frac{x'Lx}{2}}e^{-\frac{x'Lx_{0}}{2}}e^{-\frac{x_{0}'Lx}{2}}e^{-\frac{x_{0}'Lx_{0}}{2}}\mathrm{d}x\\
 & =\int_{\mathbb{R}^{2}}e^{-\frac{x'Lx}{2}}e^{-x'Lx_{0}}e^{-\frac{x_{0}'Lx_{0}}{2}}\mathrm{d}x.
\end{align*}

According to this equality,

\[
e^{h(S_{1},S_{2})}=e^{\frac{s'Js}{2N}}=c\int_{\mathbb{R}^{2}}e^{-\frac{1}{2N}x'Jx}e^{-\frac{1}{N}x'Js}\mathrm{d}x,
\]
where $c=\frac{1}{2\pi\sqrt{\det J}}$.

We switch variables: $y=\frac{1}{N}Jx$ and sum over all possible
values of the variables $(X_{1},\ldots,X_{N_{1}},Y_{1},\ldots,Y_{N_{2}})$:

\[
\sum e^{h(S_{1},S_{2})}=c\sum\int_{\mathbb{R}^{2}}e^{-\frac{N}{2}y\cdot Ly}e^{s\cdot y}\mathrm{d}^{2}y,
\]
where $c$ is a term that depends on the matrix $L$ and $N$.

The expectation $\mathbb{E}(X_{1}\cdots X_{K}Y_{1}\cdots Y_{L})$
can be calculated as follows

\begin{align*}
 & =c\sum\int_{\mathbb{R}^{2}}e^{-\frac{N}{2}y\cdot Ly}e^{s\cdot y}\mathrm{d}^{2}yE_{y_{1},y_{2}}(X_{1}\cdots X_{K}Y_{1}\cdots Y_{L})\\
 & =c\sum\int_{\mathbb{R}^{2}}e^{-\frac{N}{2}y\cdot Ly}e^{s\cdot y}\mathrm{d}^{2}yE_{y_{1}}(X_{1})\cdots E_{y_{1}}(X_{K})E_{y_{2}}(Y_{1})\cdots E_{y_{2}}(Y_{L})\\
 & =c\sum\int_{\mathbb{R}^{2}}e^{-\frac{N}{2}y\cdot Ly}e^{s\cdot y}\mathrm{d}^{2}y\tanh(y_{1})\cdots\tanh(y_{1})\tanh(y_{2})\cdots\tanh(y_{2})\\
 & =c\sum\int_{\mathbb{R}^{2}}e^{-\frac{N}{2}y\cdot Ly}e^{s\cdot y}\mathrm{d}^{2}y\tanh^{K}(y_{1})\tanh^{L}(y_{2}),
\end{align*}
where the third equality follows from

\begin{align*}
E_{y_{1}}(X_{1}) & =+1\frac{1}{2}(1+\tanh(y_{1}))-1\frac{1}{2}(1-\tanh(y_{1}))\\
 & =\frac{1}{2}\tanh(y_{1})+\frac{1}{2}\tanh(y_{1})\\
 & =\tanh(y_{1}),
\end{align*}
and similarly for $E_{y_{2}}(Y_{1})$.

Note that summing over all values of $(X_{1},\ldots,X_{N_{1}},Y_{1},\ldots,Y_{N_{2}})$,
we obtain $\sum e^{s\cdot y}=(e^{+1y_{1}}+e^{-1y_{1}})\cdots(e^{+1y_{1}}+e^{-1y_{1}})(e^{+1y_{2}}+e^{-1y_{2}})\cdots(e^{+1y_{2}}+e^{-1y_{2}})=2^{N}\cosh^{N_{1}}y_{1}\cosh^{N_{2}}y_{2}$.
Including the factor $2^{N}$ in the constant $c$ we then have

\begin{align*}
 & =c\int_{\mathbb{R}^{2}}e^{-\frac{N}{2}y\cdot Ly}\mathrm{d}^{2}y\cosh^{N_{1}}y_{1}\cosh^{N_{2}}y_{2}\tanh^{K}(y_{1})\tanh^{L}(y_{2})\\
 & =c\int e^{-N(1/2y\cdot Ly-N_{1}/N\ln\cosh y_{1}-N_{2}/N\ln\cosh y_{2})}\tanh^{K}y_{1}\tanh^{L}y_{2}\mathrm{d}^{2}y.
\end{align*}

\subsection{Extrema of the function $F$}

In order to apply Laplace's Method, we need to determine the minima
of the function

\[
F(y_{1},y_{2})=\frac{1}{2}L_{1}y_{1}^{2}+\frac{1}{2}L_{2}y_{2}^{2}-\bar{L}y_{1}y_{2}-\alpha_{1}\ln\cosh y_{1}-\alpha_{2}\ln\cosh y_{2}.
\]
\begin{prop}
\label{prop:pos-def}If

\begin{align}
L_{1} & >\alpha_{1},\label{eq:posdef1}\\
(L_{1}-\alpha_{1})(L_{2}-\alpha_{2}) & >\bar{L}^{2},\label{eq:posdef2}
\end{align}

then the above function $F$ has a unique minimum at $(0,0)$. 
\end{prop}

\begin{proof}
We take derivatives with respect to both variables

\begin{align}
F_{1}(y_{1},y_{2}) & =L_{1}y_{1}-\bar{L}y_{2}-\alpha_{1}\tanh y_{1},\label{eq:FOC1}\\
F_{2}(y_{1},y_{2}) & =L_{2}y_{2}-\bar{L}y_{1}-\alpha_{2}\tanh y_{2},\label{eq:FOC2}\\
F_{11}(y_{1},y_{2}) & =L_{1}-\frac{\alpha_{1}}{\cosh^{2}y_{1}},\nonumber \\
F_{22}(y_{1},y_{2}) & =L_{2}-\frac{\alpha_{2}}{\cosh^{2}y_{2}}.\nonumber 
\end{align}

The Hessian matrix of $F$ is

\[
H=\left[\begin{array}{cc}
F_{11} & F_{12}\\
F_{21} & F_{22}
\end{array}\right]=\left[\begin{array}{cc}
L_{1}-\frac{\alpha_{1}}{\cosh^{2}y_{1}} & -\bar{L}\\
-\bar{L} & L_{2}-\frac{\alpha_{2}}{\cosh^{2}y_{2}}
\end{array}\right].
\]

One solution to the first order conditions (\ref{eq:FOC1}) and (\ref{eq:FOC2})
is $y_{1}=y_{2}=0$.

The matrix $H$ is positive definite at the origin if and only if
(\ref{eq:posdef1}) and (\ref{eq:posdef2}) hold. Hence there is a
local minimum at the origin. If the Hessian matrix at the origin is
positive definite, it is also positive definite at any other point
due to $\cosh|s|>\cosh|t|$ for all $|s|>|t|$. Thus $F$ is strictly
convex and it follows that the minimum is unique and global. 
\end{proof}
We now show that the conditions given in the above proposition are
equivalent to our definition of the high temperature regime. 
\begin{prop}
The conditions on the Hessian matrix H given in proposition \ref{prop:pos-def}
are equivalent to the following conditions on the coupling matrix
$J$:

\begin{align}
J_{1} & <\frac{1}{\alpha_{1}},\label{eq:condJ1}\\
J_{2} & <\frac{1}{\alpha_{2}},\label{eq:condJ2}\\
\bar{J}^{2} & <(\frac{1}{\alpha_{1}}-J_{1})(\frac{1}{\alpha_{2}}-J_{2}).\label{eq:condJ3}
\end{align}
\end{prop}

\begin{proof}
We start by showing the equivalence of (\ref{eq:posdef2}) and (\ref{eq:condJ3}):

\begin{align*}
(L_{1}-\alpha_{1})(L_{2}-\alpha_{2}) & >\bar{L}^{2}\iff\\
L_{1}L_{2}-\alpha_{1}L_{2}-\alpha_{2}L_{1}+\alpha_{1}\alpha_{2} & >\bar{L}^{2}\iff\\
\frac{J_{1}J_{2}}{\Delta^{2}}-\alpha_{1}\frac{J_{1}}{\Delta}-\alpha_{2}\frac{J_{2}}{\Delta}+\alpha_{1}\alpha_{2} & >\frac{\bar{J}^{2}}{\Delta^{2}}\iff\\
\frac{J_{1}J_{2}-\bar{J}^{2}}{\Delta^{2}}-\alpha_{1}\frac{J_{1}}{\Delta}-\alpha_{2}\frac{J_{2}}{\Delta}+\alpha_{1}\alpha_{2} & >0\iff\\
\frac{1}{\Delta}-\alpha_{1}\frac{J_{1}}{\Delta}-\alpha_{2}\frac{J_{2}}{\Delta}+\alpha_{1}\alpha_{2} & >0\iff\\
1-\alpha_{1}J_{1}-\alpha_{2}J_{2}+\alpha_{1}\alpha_{2}(J_{1}J_{2}-\bar{J}^{2}) & >0\iff\\
1-\alpha_{1}J_{1}-\alpha_{2}J_{2}+\alpha_{1}\alpha_{2}J_{1}J_{2} & >\alpha_{1}\alpha_{2}\bar{J}^{2}\iff\\
(1-\alpha_{1}J_{1})(1-\alpha_{2}J_{2}) & >\alpha_{1}\alpha_{2}\bar{J}^{2}\iff\\
(\frac{1}{\alpha_{1}}-J_{1})(\frac{1}{\alpha_{2}}-J_{2}) & >\bar{J}^{2}
\end{align*}

Next we show that (\ref{eq:posdef1}) and (\ref{eq:condJ3}) imply
both (\ref{eq:condJ1}) and (\ref{eq:condJ2}). Due to (\ref{eq:condJ3}),
which is implied by (\ref{eq:posdef2}), $\frac{1}{\alpha_{1}}-J_{1}$
and $\frac{1}{\alpha_{2}}-J_{2}$ must have the same sign. Now suppose
we had $\frac{1}{\alpha_{1}}-J_{1}\leq0$ and $\frac{1}{\alpha_{2}}-J_{2}\leq0$.
We calculate

\begin{align*}
L_{1} & >\alpha_{1}\iff\\
\frac{J_{2}}{J_{1}J_{2}-\bar{J}^{2}} & >\alpha_{1}\iff\\
J_{2} & >\alpha_{1}J_{1}J_{2}-\alpha_{1}\bar{J}^{2}\iff\\
1 & >\alpha_{1}J_{1}-\alpha_{1}\frac{\bar{J}^{2}}{J_{2}}\iff\\
1+\alpha_{1}\frac{\bar{J}^{2}}{J_{2}} & >\alpha_{1}J_{1}\iff\\
\frac{\bar{J}^{2}}{J_{2}} & >J_{1}-\frac{1}{\alpha_{1}}.
\end{align*}

Similarly, we have $\frac{\bar{J}^{2}}{J_{1}}>J_{2}-\frac{1}{\alpha_{2}}$.
Remember that the smaller sides of these two inequalities are assumed
non negative. Therefore, multiplying these two inequalities yields

\[
\frac{\bar{J}^{4}}{J_{1}J_{2}}>(J_{1}-\frac{1}{\alpha_{1}})(J_{2}-\frac{1}{\alpha_{2}}).
\]
Due to $\Delta>0$, the left hand side of this inequality is smaller
than $\bar{J}^{2}$. Hence we have

\[
\bar{J}^{2}>(J_{1}-\frac{1}{\alpha_{1}})(J_{2}-\frac{1}{\alpha_{2}})>\bar{J}^{2},
\]
a contradiction. It must be that $\frac{1}{\alpha_{1}}-J_{1}>0$ and
$\frac{1}{\alpha_{2}}-J_{2}>0$. We have therefore proved that the
conditions on $H$ imply the conditions on $J$.

Now we show the other direction: suppose the conditions on $J$ hold
and also assume $L_{1}-\alpha_{1}\le0$. Then

\begin{align*}
\frac{J_{2}}{J_{1}J_{2}-\bar{J}^{2}}-\alpha_{1} & \leq0\iff\\
J_{2} & \leq\alpha_{1}J_{1}J_{2}-\alpha_{1}\bar{J}^{2}\iff\\
1 & \leq\alpha_{1}J_{1}-\alpha_{1}\frac{\bar{J}^{2}}{J_{2}}\iff\\
\frac{1}{\alpha_{1}} & \leq J_{1}-\frac{\bar{J}^{2}}{J_{2}}\iff\\
\frac{\bar{J}^{2}}{J_{2}} & \leq J_{1}-\frac{1}{\alpha_{1}}.
\end{align*}
Note that $J_{1}-\frac{1}{\alpha_{1}}<0$ and we have a contradiction
because $\frac{\bar{J}^{2}}{J_{2}}>0$. Therefore, $L_{1}-\alpha_{1}>0$
must hold. Analogously, $L_{2}-\alpha_{2}>0$ can be shown. 
\end{proof}

\subsection{Correlations $\mathbb{E}(X_{1}\cdots X_{K}Y_{1}\cdots Y_{L})$}

To calculate correlations $\mathbb{E}(X_{1}\cdots X_{K}Y_{1}\cdots Y_{L})$,
we need to estimate the integral

\[
\int e^{-N(1/2y\cdot Ly-N_{1}/N\ln\cosh y_{1}-N_{2}/N\ln\cosh y_{2})}\tanh^{K}y_{1}\tanh^{L}y_{2}\mathrm{d}^{2}y.
\]

We switch variables and define $t_{k}=\tanh y_{k}$ for $k=1,2$.
Then the same results as in \cite{MM} hold: the Jacobian matrix for
$(t_{1},t_{2})=(\tanh y_{1},\tanh y_{1})$ is

\[
\left(\begin{array}{cc}
1-\tanh^{2}y_{1} & 0\\
0 & 1-\tanh^{2}y_{1}
\end{array}\right),
\]
hence we need to multiply by $((1-\tanh^{2}y_{1})(1-\tanh^{2}y_{1}))^{-1}=(1-t_{1}^{2})^{-1}(1-t_{2}^{2})^{-1}$:

\begin{align*}
\int_{\mathbb{R}^{2}}e^{-\frac{N}{2}y'Ly+N_{1}\ln\cosh y_{1}+N_{2}\ln\cosh y_{2}}\tanh^{K}y_{1}\tanh^{L}y_{2}\mathrm{d}y & =\\
=\int_{[-1,1]^{2}}e^{-\frac{N}{2}F(t_{1},t_{2})}t_{1}^{K}t_{2}^{L}(1-t_{1}^{2})^{-1}(1-t_{2}^{2})^{-1}\mathrm{d}t,
\end{align*}
where

\begin{equation}
F(t_{1},t_{2})=\frac{1}{4}(L_{1}l(t_{1})^{2}-2\bar{L}l(t_{1})l(t_{2})+L_{2}l(t_{2})^{2})+\alpha_{1}\ln(1-t_{1}^{2})+\alpha_{2}\ln(1-t_{2}^{2})\label{eq:F-fn}
\end{equation}
and $l(x)=\ln\frac{1+x}{1-x}$ for all $x\in(-1,1)$.

Therefore, the integral to estimate is:

\[
\int_{[-1,1]^{2}}e^{-\frac{N}{2}F(t_{1},t_{2})}t_{1}^{K}t_{2}^{L}(1-t_{1}^{2})^{-1}(1-t_{2}^{2})^{-1}\mathrm{d}t.
\]

As we have established, under certain assumptions, the function $F$
has a unique minimum at $(0,0)$, where its first partial derivatives
are both equal to 0 and its Hessian is positive definite. Additionally,
its second partial derivatives are positive everywhere. Unfortunately,
the function $F$ is not separable additively into functions $G_{1}$
and $G_{2}$ which depend only on $t_{1}$ and $t_{2}$, respectively.
If it were, we could apply the one-dimensional Laplace's Method to
estimate this integral as the product of two one-dimensional integrals.
Instead, we can apply a multivariate version of Taylor's theorem to
$F$:

\begin{align*}
F(t_{1},t_{2}) & =F(0,0)+F_{1}(0,0)t_{1}+F_{2}(0,0)t_{2}\\
 & +\frac{1}{2}(F_{11}(0,0)t_{1}^{2}+2F_{12}(0,0)t_{1}t_{2}+F_{22}(0,0)t_{2}^{2})+r(t_{1},t_{2}),
\end{align*}
where we used the notation $F_{k}$ to indicate a partial derivative
of $F$ with respect to $t_{k}$, $F_{kl}$, to indicate a second
partial derivative, and $r$ is a remainder function that satisfies

\begin{equation}
|r(t_{1},t_{2})|\leq M(|t_{1}|+|t_{2}|)^{3}\label{eq:RemIneq}
\end{equation}
for some $M>0$.

As we know $F(0,0)=F_{1}(0,0)=F_{2}(0,0)=0$. Hence,

\begin{equation}
F(t_{1},t_{2})=\frac{1}{2}(F_{11}(0,0)t_{1}^{2}+2F_{12}(0,0)t_{1}t_{2}+F_{22}(0,0)t_{2}^{2})+r(t_{1},t_{2}).\label{eq:Taylor}
\end{equation}

Let's choose a small $\delta>0$ and divide the integral 
\begin{equation}
\int_{[-1,1]^{2}}e^{-\frac{N}{2}F(t_{1},t_{2})}t_{1}^{K}t_{2}^{L}(1-t_{1}^{2})^{-1}(1-t_{2}^{2})^{-1}\mathrm{d}t\label{eq:LaplaceInt}
\end{equation}
into two parts:

\begin{align*}
\int_{|t_{1}|,|t_{2}|<\delta}e^{-\frac{N}{2}F(t_{1},t_{2})}t_{1}^{K}t_{2}^{L}(1-t_{1}^{2})^{-1}(1-t_{2}^{2})^{-1}\mathrm{d}t\\
+\int_{|t_{1}|,|t_{2}|\geq\delta}e^{-\frac{N}{2}F(t_{1},t_{2})}t_{1}^{K}t_{2}^{L}(1-t_{1}^{2})^{-1}(1-t_{2}^{2})^{-1}\mathrm{d}t.
\end{align*}

We shall refer to the first integral as $I_{1}$ and the second as
$I_{2}$. Since $F$ has a unique minimum at $(0,0)$, we have $\inf_{|t_{1}|,|t_{2}|\geq\delta}F(t_{1},t_{2})\geq B>F(0,0)$
for some $B>0$. Then we have the following upper bound for $I_{2}$:

\begin{align*}
\int_{|t_{1}|,|t_{2}|\geq\delta}e^{-\frac{N}{2}F(t_{1},t_{2})}t_{1}^{K}t_{2}^{L}(1-t_{1}^{2})^{-1}(1-t_{2}^{2})^{-1}\mathrm{d}t\\
\leq e^{-\frac{N-N_{0}}{2}B}\int_{|t_{1}|,|t_{2}|\geq\delta}e^{-\frac{N_{0}}{2}F(t_{1},t_{2})}t_{1}^{K}t_{2}^{L}(1-t_{1}^{2})^{-1}(1-t_{2}^{2})^{-1}\mathrm{d}t.
\end{align*}
The right hand side goes to 0 and therefore so does $I_{2}$ as $N\rightarrow\infty$.
As for $I_{1}$, we substitute the expression for $F(t_{1},t_{2})$
from (\ref{eq:Taylor}). Then $I_{1}$ can be written as

\[
\int_{|t_{1}|,|t_{2}|<\delta}e^{-\frac{N}{4}((F_{11}(0,0)t_{1}^{2}+2F_{12}(0,0)t_{1}t_{2}+F_{22}(0,0)t_{2}^{2})+2r(t_{1},t_{2}))}t_{1}^{K}t_{2}^{L}(1-t_{1}^{2})^{-1}(1-t_{2}^{2})^{-1}\mathrm{d}t.
\]

Now we use the symbol $H$ again to refer to the Hessian matrix

\[
H=\left(\begin{array}{cc}
F_{11}(0,0) & F_{12}(0,0)\\
F_{12}(0,0) & F_{22}(0,0)
\end{array}\right).
\]

Now we change variables $x=(\frac{N}{2})^{1/2}t$ and obtain

\begin{align*}
\frac{2}{N}\int_{|t_{1}|,|t_{2}|<\delta(\frac{N}{2})^{1/2}}e^{-\frac{1}{2}x'Hx-\frac{N}{2}r((\frac{N}{2})^{-1/2}x)}(\frac{2}{N})^{\frac{K+L}{2}}x_{1}^{K}x_{2}^{L}\\
\cdot(1-(\frac{2}{N}x_{1})^{2})^{-1}(1-(\frac{2}{N}x_{2})^{2})^{-1}\mathrm{d}x.
\end{align*}

The remainder function $r$ satisfies

\[
\frac{N}{2}r((\frac{N}{2})^{-1/2}x)\rightarrow0
\]
due to inequality (\ref{eq:RemIneq}).

The function $(1-(\frac{2}{N}x_{1})^{2})^{-1}(1-(\frac{2}{N}x_{2})^{2})^{-1}$
converges uniformly on\\
 $[-\delta(\frac{N}{2})^{1/2},\delta(\frac{N}{2})^{1/2}]$ as $N\rightarrow\infty$,
provided $\delta$ is small enough. The limiting function has the
value 1 at the origin.

Therefore, the integral asymptotically behaves like

\begin{align*}
(\frac{2}{N})^{\frac{K+L}{2}+1}\int_{|t_{1}|,|t_{2}|<\delta(\frac{N}{2})^{1/2}}e^{-\frac{1}{2}x'Hx}x_{1}^{K}x_{2}^{L}\mathrm{d}x.
\end{align*}

As $N\rightarrow\infty$, the limits of integration in the above integral
extend to $-\infty$ and $\infty$. If we define

\begin{equation}
\left(\begin{array}{cc}
\sigma_{11} & \sigma_{12}\\
\sigma_{12} & \sigma_{22}
\end{array}\right)=\Sigma=H^{-1}\label{eq:Sigma}
\end{equation}

and $m_{K,L}=m_{K,L}(0,\Sigma)$ as the moment of order $(K,L)$ of
a bivariate normal random vector $(Z_{1},Z_{2})\sim N(0,\Sigma)$,
then we see that the integral converges asymptotically to

\begin{align*}
(\frac{2}{N})^{\frac{K+L}{2}+1}2\pi(\det H)m_{K,L}.
\end{align*}

In order to calculate correlations $\mathbb{E}(X_{1}\cdots X_{K}Y_{1}\cdots Y_{L})$,
we take the above expression and divide by its value when $K=L=0$:

\[
\frac{(\frac{2}{N})^{\frac{K+L}{2}+1}2\pi(\det H)m_{K,L}}{\frac{2}{N}\cdot2\pi(\det H)m_{0,0}}=(\frac{2}{N})^{\frac{K+L}{2}}m_{K,L}.
\]

Since $m_{K,L}>0$ if and only if $K+L$ is even and 0 otherwise,
we have shown 
\begin{thm}
\label{thm:EXY-het}For all $K,L=0,1,2,\ldots$, the expected value
$\mathbb{E}(X_{1}\cdots X_{K}Y_{1}\cdots Y_{L})$ is asymptotically
zero if $K+L$ is odd and otherwise given by the expression 
\[
(\frac{2}{N})^{\frac{K+L}{2}}m_{K,L}(0,H^{-1}),
\]
where H is the Hessian matrix H at the origin of the function 
\[
F(t_{1},t_{2})=\frac{1}{4}(L_{1}l(t_{1})^{2}+2\bar{L}l(t_{1})l(t_{2})+L_{2}l(t_{2})^{2})+\alpha_{1}\ln(1-t_{1}^{2})+\alpha_{2}\ln(1-t_{2}^{2}),
\]
and $l(x)=\ln\frac{1+x}{1-x}$ for all $x\in(-1,1)$. 
\end{thm}


\subsection{Moments}

We calculate the moments $\mathbb{E}(\frac{1}{N_{1}^{K}}S_{1}^{K}\frac{1}{N_{2}^{L}}S_{2}^{L})$
and $M_{K,L}=\mathbb{E}(\frac{1}{N_{1}^{K/2}}S_{1}^{K}\frac{1}{N_{2}^{L/2}}S_{2}^{L})$
which we shall use to prove the law of large numbers and the central
limit theorem.

\subsubsection{Law of Large Numbers}

We start with $\mathbb{E}(\frac{1}{N_{1}^{K}}S_{1}^{K}\frac{1}{N_{2}^{L}}S_{2}^{L})$.
These moments $\mathbb{E}(\frac{1}{N_{1}^{K}}S_{1}^{K}\frac{1}{N_{2}^{L}}S_{2}^{L})$
can be calculated as the sum

\begin{equation}
\frac{1}{N_{1}^{K}N_{2}^{L}}\sum_{k=0}^{K}\sum_{\underline{i}\in\Pi_{k}^{(K)}}\sum_{l=0}^{L}\sum_{\underline{j}\in\Pi_{l}^{(L)}}w_{K}(\underline{i})w_{L}(\underline{j})\mathbb{E}(X_{\underline{i}}Y_{\underline{j}}).\label{eq:moments-LLN}
\end{equation}
By theorem \ref{thm:EXY-het}, the correlation $\mathbb{E}(X_{\underline{i}}Y_{\underline{j}})$
depends on the number of different indices that occur in $\underline{i}$
and $\underline{j}$. Since for any $n\in\mathbb{N}$ and any $X_{i}$
and $Y_{j}$ $X_{i}^{2n}=Y_{j}^{2n}=1$, only the indices which occur
an odd number of times contribute a power of $N$ asymptotically to
$\mathbb{E}(X_{\underline{i}}Y_{\underline{j}})$. That power is

\[
N^{-\frac{1}{2}(\sum_{m=0}^{\left\lfloor K/2\right\rfloor }\rho_{2m+1}(\underline{i})+\sum_{n=0}^{\left\lfloor L/2\right\rfloor }\rho_{2n+1}(\underline{j}))}.
\]

Each summand has thus the form $\frac{1}{N_{1}^{K}N_{2}^{L}}w_{K}(\underline{i})w_{L}(\underline{j})\mathbb{E}(X_{\underline{i}}Y_{\underline{j}})=$

\[
=\frac{1}{N_{1}^{K}N_{2}^{L}}bN_{1}^{\sum_{m=1}^{K}\rho_{m}(\underline{i})}N_{2}^{\sum_{n=1}^{L}\rho_{n}(\underline{j})}N^{-\frac{1}{2}(\sum_{m=0}^{\left\lfloor K/2\right\rfloor }\rho_{2m+1}(\underline{i})+\sum_{n=0}^{\left\lfloor L/2\right\rfloor }\rho_{2n+1}(\underline{j}))},
\]
where $b$ is a combinatorial constant that also includes all powers
of two from $\mathbb{E}(X_{\underline{i}}Y_{\underline{j}})$. We
disregard the constant, and, since $N_{\nu}\approx\alpha_{\nu}N$,
what remains are only powers of $N$:

\[
\frac{1}{N^{K}N^{L}}N^{\sum_{m=1}^{K}\rho_{m}(\underline{i})}N^{\sum_{n=1}^{L}\rho_{n}(\underline{j})}N^{-\frac{1}{2}(\sum_{m=0}^{\left\lfloor K/2\right\rfloor }\rho_{2m+1}(\underline{i})+\sum_{n=0}^{\left\lfloor L/2\right\rfloor }\rho_{2n+1}(\underline{j}))}.
\]
Hence the exponent of $N$ in each summand is given by

\begin{align*}
 & -K-L+\sum_{m=1}^{K}\rho_{m}(\underline{i})+\sum_{n=1}^{L}\rho_{n}(\underline{j})-\frac{1}{2}(\sum_{m=0}^{\left\lfloor K/2\right\rfloor }\rho_{2m+1}(\underline{i})+\sum_{n=0}^{\left\lfloor L/2\right\rfloor }\rho_{2n+1}(\underline{j}))\\
= & -K-L+\sum_{m=1}^{\left\lfloor K/2\right\rfloor }\rho_{2m}(\underline{i})+\sum_{n=1}^{\left\lfloor L/2\right\rfloor }\rho_{2n}(\underline{j})+\frac{1}{2}(\sum_{m=0}^{\left\lfloor K/2\right\rfloor }\rho_{2m+1}(\underline{i})+\sum_{n=0}^{\left\lfloor L/2\right\rfloor }\rho_{2n+1}(\underline{j})).
\end{align*}
\begin{lem}
For all $\underline{i}\in\Pi^{(K)}$ and all $\underline{j}\in\Pi^{(L)}$,
we have the inequality

\[
-K-L+\sum_{m=1}^{\left\lfloor K/2\right\rfloor }\rho_{2m}(\underline{i})+\sum_{n=1}^{\left\lfloor L/2\right\rfloor }\rho_{2n}(\underline{j})+\frac{1}{2}(\sum_{m=0}^{\left\lfloor K/2\right\rfloor }\rho_{2m+1}(\underline{i})+\sum_{n=0}^{\left\lfloor L/2\right\rfloor }\rho_{2n+1}(\underline{j}))<0.
\]
\end{lem}

\begin{proof}
By lemma 3.8 in \cite{MM}, we have $\sum_{m=1}^{K}m\rho_{m}(\underline{i})=K$.
If $\rho_{1}(\underline{i})=K$, then $\rho_{m}(\underline{i})=0$
for all $m>1$, and we have $-K+\sum_{m=1}^{\left\lfloor K/2\right\rfloor }\rho_{2m}(\underline{i})+\frac{1}{2}\sum_{m=0}^{\left\lfloor K/2\right\rfloor }\rho_{2m+1}(\underline{i})=-K+0+K/2<0$.
If $\rho_{1}(\underline{i})<K$, then there must be some $m>1$ with
$\rho_{m}(\underline{i})>0$ and therefore, $\sum_{m=1}^{\left\lfloor K/2\right\rfloor }\rho_{2m}(\underline{i})+\frac{1}{2}\sum_{m=0}^{\left\lfloor K/2\right\rfloor }\rho_{2m+1}(\underline{i})<\sum_{m=1}^{K}m\rho_{m}(\underline{i})=K$.
Similarly, $\sum_{n=1}^{\left\lfloor L/2\right\rfloor }\rho_{2n}(\underline{j})+\frac{1}{2}\sum_{n=0}^{\left\lfloor L/2\right\rfloor }\rho_{2n+1}(\underline{j})<L$.
By summing both inequalities, the result follows. 
\end{proof}
Since there are only finitely many summands in (\ref{eq:moments-LLN}),
each of which converges to 0, the above lemma proves the law of large
numbers (theorem \ref{LLN}).

\subsubsection{Central Limit Theorem}

Now we turn to $M_{K,L}$. These moments $M_{K,L}$ can be calculated
as the sum

\begin{equation}
\frac{1}{N_{1}^{K/2}N_{2}^{L/2}}\sum_{k=0}^{K}\sum_{\underline{i}\in\Pi_{k}^{(K)}}\sum_{l=0}^{L}\sum_{\underline{j}\in\Pi_{l}^{(L)}}w_{K}(\underline{i})w_{L}(\underline{j})\mathbb{E}(X_{\underline{i}}Y_{\underline{j}}).
\end{equation}

We separate the above sum into four summands:

\begin{align*}
A_{1} & =\frac{1}{N_{1}^{K/2}N_{2}^{L/2}}\sum_{k=0}^{K}\sum_{\underline{i}\in\Pi_{k}^{0(K)}}\sum_{l=0}^{L}\sum_{\underline{j}\in\Pi_{l}^{0(L)}}w_{K}(\underline{i})w_{L}(\underline{j})\mathbb{E}(X_{\underline{i}}Y_{\underline{j}}),\\
A_{2} & =\frac{1}{N_{1}^{K/2}N_{2}^{L/2}}\sum_{k=0}^{K}\sum_{\underline{i}\in\Pi_{k}^{0(K)}}\sum_{l=0}^{L}\sum_{\underline{j}\in\Pi_{l}^{+(L)}}w_{K}(\underline{i})w_{L}(\underline{j})\mathbb{E}(X_{\underline{i}}Y_{\underline{j}}),\\
A_{3} & =\frac{1}{N_{1}^{K/2}N_{2}^{L/2}}\sum_{k=0}^{K}\sum_{\underline{i}\in\Pi_{k}^{+(K)}}\sum_{l=0}^{L}\sum_{\underline{j}\in\Pi_{l}^{0(L)}}w_{K}(\underline{i})w_{L}(\underline{j})\mathbb{E}(X_{\underline{i}}Y_{\underline{j}}),\\
A_{4} & =\frac{1}{N_{1}^{K/2}N_{2}^{L/2}}\sum_{k=0}^{K}\sum_{\underline{i}\in\Pi_{k}^{+(K)}}\sum_{l=0}^{L}\sum_{\underline{j}\in\Pi_{l}^{+(L)}}w_{K}(\underline{i})w_{L}(\underline{j})\mathbb{E}(X_{\underline{i}}Y_{\underline{j}}).
\end{align*}

We claim that only those multiindices $\underline{i},\underline{j}$
contribute to this sum that only contain each index at most twice: 
\begin{prop}
The limits of $A_{2},A_{3}$ and $A_{4}$ as $N\rightarrow\infty$
are zero. 
\end{prop}

\begin{proof}
We only show $A_{2}\rightarrow\infty$; the other claims are proved
analogously. For fixed $k,l$ let $\underline{i}\in\Pi_{k}^{0(K)}$
and $\underline{j}\in\Pi_{l}^{+(L)}$ and let $\underline{r}$ and
$\underline{s}$ be their respective profile vectors. Then

\begin{align*}
\frac{1}{N_{1}^{K/2}N_{2}^{L/2}}w_{K}(\underline{i})w_{L}(\underline{j})\mathbb{E}(X_{\underline{i}}Y_{\underline{j}}) & \approx\\
\frac{1}{N_{1}^{K/2}N_{2}^{L/2}}\frac{N_{1}!}{r_{1}!r_{2}!\cdots r_{K}!r_{0}!}\frac{K!}{1!^{r_{1}}2!^{r_{2}}\cdots K!^{r_{K}}}\frac{N_{2}!}{s_{1}!s_{2}!\cdots s_{L}!s_{0}!} & \cdot\\
\cdot\frac{L!}{1!^{s_{1}}2!^{s_{2}}\cdots L!^{s_{L}}}(\frac{2}{N})^{\frac{k+l}{2}}m_{k,l} & =\\
\frac{1}{N_{1}^{K/2}N_{2}^{L/2}}\frac{N_{1}!}{k!(\frac{K-k}{2})!(N_{1}-\frac{K+k}{2})!}\frac{K!}{2!^{\frac{K-k}{2}}}\frac{N_{2}!}{l!s_{2}!\cdots s_{L}!s_{0}!} & \cdot\\
\cdot\frac{L!}{2!^{s_{2}}\cdots L!^{s_{L}}}(\frac{2}{N})^{\frac{k+l}{2}}m_{k,l} & =\\
\frac{1}{N_{2}^{L/2}}\frac{K!}{k!(\frac{K-k}{2})!2!^{\frac{K-k}{2}}}\frac{N_{2}!}{l!s_{2}!\cdots s_{L}!s_{0}!} & \cdot\\
\cdot\frac{L!}{2!^{s_{2}}\cdots L!^{s_{L}}}2^{\frac{k+l}{2}}\frac{1}{N^{\frac{l}{2}}}m_{k,l} & \approx\\
\frac{c}{N_{2}^{L/2}}N_{2}^{\sum_{j=1}^{L}s_{j}}N^{-\frac{l}{2}} & \leq\\
\frac{c}{N_{2}^{L/2}}N_{2}^{\frac{L+l-1}{2}}N^{-\frac{l}{2}} & \approx\\
c\alpha_{1}^{k/2}\alpha_{2}^{l/2}N_{2}^{-\frac{1}{2}} & \rightarrow0.
\end{align*}

The constant $c$ in the fourth line above represents the product
of all those factors that do not depend on $N_{2}$ and $N$. In the
inequality above we used that

\[
\sum_{j=1}^{L}s_{j}\leq l+\frac{1}{2}\sum_{j=2}^{L}js_{j}-\frac{1}{2}=\frac{l}{2}+\frac{1}{2}\sum_{j=1}^{L}js_{j}-\frac{1}{2}=\frac{L+l-1}{2}.
\]

Since each summand goes to zero, and there are only finitely many
summands, $A_{2}$ goes to zero as $N_{2}$ goes to infinity. 
\end{proof}
It follows from the last proposition that $M_{K,L}$ is asymptotically
equal to $A_{1}$. We note that for $K+L$ odd, we have $M_{K,L}=0$
due to $\mathbb{E}(X_{\underline{i}}Y_{\underline{j}})=0$ in that
case. For $K+L$ even, we have to distinguish the case where $K,L$
are even and the case where they are both odd. We proceed with the
case where they are both even. Then $M_{K,L}=$

\begin{align*}
 & =\frac{1}{N_{1}^{K/2}N_{2}^{L/2}}\sum_{k=0}^{K/2}\sum_{l=0}^{L/2}\frac{N_{1}!}{(2k)!(\frac{K}{2}-k)!(N_{1}-\frac{K}{2}-k)!}\frac{K!}{2^{\frac{K}{2}-k}}\cdot\\
 & \cdot\frac{N_{2}!}{(2l)!(\frac{L}{2}-l)!(N_{2}-\frac{L}{2}-l)!}\frac{L!}{2^{\frac{L}{2}-l}}(\frac{2}{N})^{k+l}m_{2k,2l}\\
 & \approx\frac{1}{N_{1}^{K/2}N_{2}^{L/2}}\sum_{k=0}^{K/2}\sum_{l=0}^{L/2}\frac{N_{1}^{\frac{K}{2}+k}}{(2k)!(\frac{K}{2}-k)!}\frac{K!}{2^{\frac{K}{2}-k}}\cdot\\
 & \cdot\frac{N_{2}^{\frac{L}{2}+l}}{(2l)!(\frac{L}{2}-l)!}\frac{L!}{2^{\frac{L}{2}-l}}(\frac{2}{N})^{k+l}m_{2k,2l}\\
 & \approx\sum_{k=0}^{K/2}\sum_{l=0}^{L/2}\frac{\alpha_{1}^{k}}{(2k)!(\frac{K}{2}-k)!}\frac{K!}{2^{\frac{K}{2}}}\frac{\alpha_{2}^{l}}{(2l)!(\frac{L}{2}-l)!}\frac{L!}{2^{\frac{L}{2}}}2^{2(k+l)}m_{2k,2l}.
\end{align*}

We now turn to the normal moment $m_{2k,2l}$. We shall use Isserlis's
theorem to express this moment as the sum of products of covariances.

Let $(Z_{1},Z_{2})$ be a bivariate normal random vector with zero
mean and covariance matrix $\Sigma$ defined previously.

Isserlis's Theorem states that higher moments of the bivariate normal
distribution can be calculated

\[
\mathbb{E}(Z_{1}^{K}Z_{2}^{L})=\sum_{\pi\in\mathcal{P}_{2}}\prod_{i=1}^{\frac{K+L}{2}}\mathbb{E}(\xi_{i1}\xi_{i2}),
\]

where $\mathcal{P}_{2}$ is the set of all pair partitions on the
set $\{1,2,\ldots,K+L\}$ and $\xi_{i1}\xi_{i2}$ are two of the $K+L$
variables grouped together by a particular pair partition $\pi$.

We start by using Isserlis's Theorem to express higher moments of
the bivariate normal distribution recursively. 
\begin{lem}
\label{lem:Isserlis}For all $K,L\in\mathbb{N}_{0}$, the moments
$m_{K,L}=\mathbb{E}(Z_{1}^{K}Z_{2}^{L})$ satisfy the following equalities: 
\begin{enumerate}
\item $m_{K,L+2}=Km_{1,1}m_{K-1,L+1}+(L+1)m_{0,2}m_{K,L},$ 
\item $m_{K+2,L}=(K+1)m_{2,0}m_{K,L}+Lm_{1,1}m_{K+1,L-1}.$ 
\end{enumerate}
\end{lem}

If $K=0$ or $L=0$, then the formulas still hold, setting any moments
with negative indices equal to 0. Note that these two formulas suffice
to calculate any higher moment as a function of $K,L$ and $m_{2,0},m_{1,1},m_{0,2}$.
For a proof of the above lemma see \cite{KT1}. 
\begin{lem}
The moment $m_{K,L}$ can be expressed as

\[
\sigma_{11}^{\frac{K}{2}}\sigma_{22}^{\frac{L}{2}}\sum_{r=0}^{\min\{\frac{K}{2},\frac{L}{2}\}}\frac{K!L!}{(2r)!(\frac{K}{2}-r)!(\frac{L}{2}-r)!2^{\frac{K+L}{2}}}\left(\frac{2^{2}\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{r},
\]
when K,L are even, and as

\[
\sigma_{11}^{\frac{K-1}{2}}\sigma_{22}^{\frac{L-1}{2}}\sigma_{12}\sum_{r=0}^{\min\{\frac{K-1}{2},\frac{L-1}{2}\}}\frac{K!L!}{(2r+1)!(\frac{K-1}{2}-r)!(\frac{L-1}{2}-r)!2^{\frac{K+L}{2}-1}}\left(\frac{2^{2}\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{r},
\]
when K,L are odd. 
\end{lem}

\begin{proof}
We shall only prove the expression for even $K,L$. The odd case is
very similar. When we group the $K+L$ variables into pairs, there
are three types of pairs: both variables are $Z_{1}$, both are $Z_{2}$
and one of each variable. Depending on this, $\mathbb{E}(\xi_{i1}\xi_{i2})$
is $\sigma_{11},\sigma_{22}$ and $\sigma_{12}$, respectively. For
even $K,L$ the number of mixed pairs must be even. Let $2r$ stand
for this number. From the $K$ copies of $Z_{1}$, we pick $K-2r$
to be paired up among themselves: there are $\left(\begin{array}{c}
K\\
K-2r
\end{array}\right)$ ways to achieve this. Similarly, we pick from the $L$ copies of
$Z_{2}$ $L-2r$ to be paired up among themselves: $\left(\begin{array}{c}
L\\
L-2r
\end{array}\right)$. Then we proceed to pair the $K-2r$ $Z_{1}$: $(K-2r-1)!!$ and
the $L-2r$ $Z_{2}$: $(L-2r-1)!!$. Lastly, we pair the $2r$ $Z_{1}$
and the $2r$ $Z_{2}$: $(2r)!$. Multiplying all these factors and
summing over all possible values of mixed pairs, we obtain

\begin{align*}
m_{K,L} & =\sum_{r=0}^{\min\{\frac{K}{2},\frac{L}{2}\}}\frac{K!}{(2r)!(K-2r)!}\frac{L!}{(2r)!(L-2r)!}\frac{(K-2r)!}{(\frac{K}{2}-r)!2^{\frac{K}{2}-r}}\frac{(L-2r)!}{(\frac{L}{2}-r)!2^{\frac{L}{2}-r}}\cdot\\
 & \cdot(2r)!\sigma_{11}^{\frac{K}{2}-r}\sigma_{12}^{2r}\sigma_{22}^{\frac{L}{2}-r},
\end{align*}
and the result follows after some simplification. 
\end{proof}
Now we return to the problem of calculating $M_{K,L}$:

\begin{align*}
 & \approx\sum_{k=0}^{K/2}\sum_{l=0}^{L/2}\frac{\alpha_{1}^{k}}{(2k)!(\frac{K}{2}-k)!}\frac{K!}{2^{\frac{K}{2}}}\frac{\alpha_{2}^{l}}{(2l)!(\frac{L}{2}-l)!}\frac{L!}{2^{\frac{L}{2}}}2^{2(k+l)}m_{2k,2l}\\
 & =\sum_{k=0}^{K/2}\sum_{l=0}^{L/2}\frac{\alpha_{1}^{k}}{(2k)!(\frac{K}{2}-k)!}\frac{K!}{2^{\frac{K}{2}}}\frac{\alpha_{2}^{l}}{(2l)!(\frac{L}{2}-l)!}\frac{L!}{2^{\frac{L}{2}}}2^{2(k+l)}\cdot\\
 & \cdot\sigma_{11}^{k}\sigma_{22}^{l}\sum_{r=0}^{\min\{k,l\}}\frac{(2k)!(2l)!}{(2r)!(k-r)!(l-r)!2^{k+l}}\left(\frac{2^{2}\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{r}\\
 & =\frac{K!L!}{2^{\frac{K+L}{2}}}\sum_{k=0}^{K/2}\sum_{l=0}^{L/2}\frac{(2\alpha_{1}\sigma_{11})^{k}}{(\frac{K}{2}-k)!}\frac{(2\alpha_{2}\sigma_{22})^{l}}{(\frac{L}{2}-l)!}\sum_{r=0}^{\min\{k,l\}}\frac{1}{(2r)!(k-r)!(l-r)!}\left(\frac{2^{2}\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{r}.
\end{align*}
We have thus shown the first part of the following 
\begin{prop}
Asymptotically, the moments $\mathbb{E}(\frac{1}{N_{1}^{K/2}}S_{1}^{K}\frac{1}{N_{2}^{L/2}}S_{2}^{L})$
are given by 
\[
\frac{K!L!}{2^{\frac{K+L}{2}}}\sum_{k=0}^{K/2}\sum_{l=0}^{L/2}\frac{(2\alpha_{1}\sigma_{11})^{k}}{(\frac{K}{2}-k)!}\frac{(2\alpha_{2}\sigma_{22})^{l}}{(\frac{L}{2}-l)!}\sum_{r=0}^{\min\{k,l\}}\frac{1}{(2r)!(k-r)!(l-r)!}\left(\frac{2^{2}\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{r},
\]
when $K,L$ and $K+L$ are even, 
\begin{align*}
\frac{K!L!}{2^{\frac{K+L}{2}}}\cdot4\sqrt{\alpha_{1}\alpha_{2}}\sigma_{12}\sum_{k=0}^{\frac{K-1}{2}}\sum_{l=0}^{\frac{L-1}{2}}\frac{(2\alpha_{1}\sigma_{11})^{k}}{(\frac{K-1}{2}-k)!}\frac{(2\alpha_{2}\sigma_{22})^{l}}{(\frac{L-1}{2}-l)!}\cdot\\
\cdot\sum_{r=0}^{\min\{k,l\}}\frac{1}{(2r+1)!(k-r)!(l-r)!}\left(\frac{2^{2}\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{r},
\end{align*}
when $K,L$ are odd and $K+L$ even, and zero when $K+L$ is odd. 
\end{prop}

The second part is shown in a very similar fashion. We shall need
both formulas for the proof of the central limit theorem, which we
shall prove using lemma \ref{lem:Isserlis}. Using the above proposition,
we calculate all moments $M_{K,L}$ of order 2:

\begin{align*}
M_{2,0} & =\frac{2!0!}{2^{\frac{2+0}{2}}}\sum_{k=0}^{1}\sum_{l=0}^{0}\frac{(2\alpha_{1}\sigma_{11})^{k}}{(1-k)!}\frac{(2\alpha_{2}\sigma_{22})^{l}}{(0-l)!}\sum_{r=0}^{\min\{k,l\}}\frac{1}{(2r)!(k-r)!(l-r)!}\left(\frac{2^{2}\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{r}\\
 & =\sum_{k=0}^{1}\frac{(2\alpha_{1}\sigma_{11})^{k}}{(1-k)!}\frac{1}{0!(k-0)!(0-0)!}\left(\frac{2^{2}\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{0}\\
 & =1+2\alpha_{1}\sigma_{11},\\
M_{0,2} & =1+2\alpha_{2}\sigma_{22},\\
M_{1,1} & =\frac{1!1!}{2^{\frac{1+1}{2}}}\cdot4\sqrt{\alpha_{1}\alpha_{2}}\sigma_{12}\sum_{k=0}^{0}\sum_{l=0}^{0}\frac{(2\alpha_{1}\sigma_{11})^{k}}{(\frac{K-1}{2}-k)!}\frac{(2\alpha_{2}\sigma_{22})^{l}}{(\frac{L-1}{2}-l)!}\cdot\\
 & \cdot\sum_{r=0}^{\min\{k,l\}}\frac{1}{(2r+1)!(k-r)!(l-r)!}\left(\frac{2^{2}\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{r}\\
 & =2\sqrt{\alpha_{1}\alpha_{2}}\sigma_{12}.
\end{align*}

Now let $(Z_{1},Z_{2})$ be a bivariate normal random vector with
zero mean and covariance matrix

\[
\left(\begin{array}{cc}
M_{2,0} & M_{1,1}\\
M_{1,1} & M_{0,2}
\end{array}\right).
\]

If we can show that the moments $M_{K,L}$ are asymptotically equal
to the moments $\mathbb{E}(Z_{1}^{K}Z_{2}^{L})$, then the central
limit theorem is proved. Since $(Z_{1},Z_{2})$ is bivariate normal,
we can do this by showing that the $M_{K,L}$ satisfy the following
recursive relations: 
\begin{enumerate}
\item $M_{K,L+2}=KM_{1,1}M_{K-1,L+1}+(L+1)M_{0,2}M_{K,L},$ 
\item $M_{K+2,L}=(K+1)M_{2,0}M_{K,L}+LM_{1,1}M_{K+1,L-1}.$ 
\end{enumerate}
Since the formulas for $M_{K,L}$ are completely symmetric with respect
to the two variables, we shall limit ourselves to the second of these
equations.

The left hand side $(LHS)$ is given by

\[
\frac{(K+2)!L!}{2^{\frac{K+L}{2}+1}}\sum_{k=0}^{\frac{K}{2}+1}\sum_{l=0}^{\frac{L}{2}}\frac{(2\alpha_{1}\sigma_{11})^{k}}{(\frac{K}{2}-k)!}\frac{(2\alpha_{2}\sigma_{22})^{l}}{(\frac{L}{2}-l)!}\sum_{r=0}^{\min\{k,l\}}\frac{1}{(2r)!(k-r)!(l-r)!}\left(\frac{2^{2}\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{r},
\]
whereas the right hand side $(RHS)$ is composed of the sum

\begin{align*}
(K+1)M_{2,0}M_{K,L}+LM_{1,1}M_{K+1,L-1} & =\\
(K+1)(1+2\alpha_{1}\sigma_{11})M_{K,L}+L\cdot2\sqrt{\alpha_{1}\alpha_{2}}\sigma_{12}M_{K+1,L-1}.
\end{align*}

Hence, $RHS$ can be separated into three sums:

\begin{align*}
A & =\frac{(K+1)!L!}{2^{\frac{K+L}{2}}}\sum_{k=0}^{K/2}\sum_{l=0}^{L/2}\frac{(2\alpha_{1}\sigma_{11})^{k}}{(\frac{K}{2}-k)!}\frac{(2\alpha_{2}\sigma_{22})^{l}}{(\frac{L}{2}-l)!}\cdot\\
 & \cdot\sum_{r=0}^{\min\{k,l\}}\frac{1}{(2r)!(k-r)!(l-r)!}\left(\frac{2^{2}\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{r},\\
B & =\frac{(K+1)!L!}{2^{\frac{K+L}{2}}}\sum_{k=0}^{K/2}\sum_{l=0}^{L/2}\frac{(2\alpha_{1}\sigma_{11})^{k+1}}{(\frac{K}{2}-k)!}\frac{(2\alpha_{2}\sigma_{22})^{l}}{(\frac{L}{2}-l)!}\cdot\\
 & \cdot\sum_{r=0}^{\min\{k,l\}}\frac{1}{(2r)!(k-r)!(l-r)!}\left(\frac{2^{2}\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{r},\\
C & =2\frac{(K+1)!L!}{2^{\frac{K+L}{2}}}\sum_{k=0}^{\frac{K}{2}}\sum_{l=0}^{\frac{L}{2}-1}\frac{(2\alpha_{1})^{k+1}\sigma_{11}^{k}}{(\frac{K-1}{2}-k)!}\frac{(2\alpha_{2})^{l+1}\sigma_{22}^{l}}{(\frac{L-1}{2}-l)!}\cdot\\
 & \cdot\sum_{r=0}^{\min\{k,l\}}\frac{1}{(2r+1)!(k-r)!(l-r)!}\left(\frac{2^{2}}{\sigma_{11}\sigma_{22}}\right)^{r}\sigma_{12}^{2(r+1)}.
\end{align*}
We can simplify the equation by dividing both sides by the constant
$\frac{(K+1)!L!}{2^{\frac{K+L}{2}}}$. Then $LHS$ equals

\[
\frac{K+2}{2}\sum_{k=0}^{\frac{K}{2}+1}\sum_{l=0}^{\frac{L}{2}}\frac{(2\alpha_{1}\sigma_{11})^{k}}{(\frac{K}{2}-k)!}\frac{(2\alpha_{2}\sigma_{22})^{l}}{(\frac{L}{2}-l)!}\sum_{r=0}^{\min\{k,l\}}\frac{1}{(2r)!(k-r)!(l-r)!}\left(\frac{2^{2}\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{r},
\]

\begin{align*}
A & =\sum_{k=0}^{K/2}\sum_{l=0}^{L/2}\frac{(2\alpha_{1}\sigma_{11})^{k}}{(\frac{K}{2}-k)!}\frac{(2\alpha_{2}\sigma_{22})^{l}}{(\frac{L}{2}-l)!}\cdot\\
 & \cdot\sum_{r=0}^{\min\{k,l\}}\frac{1}{(2r)!(k-r)!(l-r)!}\left(\frac{2^{2}\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{r},\\
B & =\sum_{k=0}^{K/2}\sum_{l=0}^{L/2}\frac{(2\alpha_{1}\sigma_{11})^{k+1}}{(\frac{K}{2}-k)!}\frac{(2\alpha_{2}\sigma_{22})^{l}}{(\frac{L}{2}-l)!}\cdot\\
 & \cdot\sum_{r=0}^{\min\{k,l\}}\frac{1}{(2r)!(k-r)!(l-r)!}\left(\frac{2^{2}\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{r},\\
C & =2\sum_{k=0}^{\frac{K}{2}}\sum_{l=0}^{\frac{L}{2}-1}\frac{(2\alpha_{1})^{k+1}\sigma_{11}^{k}}{(\frac{K-1}{2}-k)!}\frac{(2\alpha_{2})^{l+1}\sigma_{22}^{l}}{(\frac{L-1}{2}-l)!}\cdot\\
 & \cdot\sum_{r=0}^{\min\{k,l\}}\frac{1}{(2r+1)!(k-r)!(l-r)!}\left(\frac{2^{2}}{\sigma_{11}\sigma_{22}}\right)^{r}\sigma_{12}^{2(r+1)}.
\end{align*}
Note that the powers $a_{1},a_{2},s_{11},s_{22},s_{12}$ of the parameters
$\alpha_{1},\alpha_{2},\sigma_{11},\sigma_{22},\sigma_{12}$, respectively
in each summand are as follows:

\selectlanguage{english}
$LHS$:

\selectlanguage{british}
\[
a_{1}=k,a_{2}=l,s_{11}=k-r,s_{22}=l-r,s_{12}=2r
\]

$A$:

\[
a_{1}=k,a_{2}=l,s_{11}=k-r,s_{22}=l-r,s_{12}=2r
\]

$B$:

\[
a_{1}=k+1,a_{2}=l,s_{11}=k+1-r,s_{22}=l-r,s_{12}=2r
\]

$C$:

\[
a_{1}=k+1,a_{2}=l+1,s_{11}=k-r,s_{22}=l-r,s_{12}=2(r+1)
\]

In order to investigate systematically whether $LHS=RHS$, we propose
the following scheme: 
\begin{enumerate}
\item Pick $a_{1}\in\{0,1,\ldots,\frac{K}{2}+1\}$ and $a_{2}\in\{0,1,\ldots,\frac{L}{2}\}$.
These choices fix the values of $k,l$ in each of the four sums. 
\item Pick $s_{12}\in\{0,2,\ldots,2\min\{a_{1},a_{2}\}\}$. This fixes the
value of $r$ in each of the four sums, as well as the values of $s_{11}=a_{1}-\frac{s_{12}}{2}$
and $s_{22}=a_{2}-\frac{s_{12}}{2}$. 
\end{enumerate}
Following this, we distinguish six main cases, as well as some subcases
if applicable: 
\begin{enumerate}
\item If $a_{1}=a_{2}=0$, then we obtain $s_{11}=s_{22}=s_{12}=0$. 
\item If $0<a_{1}<\frac{K}{2}+1$ and $a_{2}=0$, then we obtain $s_{11}=a_{1}$
and $s_{22}=s_{12}=0$. 
\item If $a_{1}=\frac{K}{2}+1$ and $a_{2}=0$, then we obtain $s_{11}=\frac{K}{2}+1$
and $s_{22}=s_{12}=0$. 
\item If $a_{1}=0$ and $a_{2}>0$, then we obtain $s_{11}=0$, $s_{22}=a_{2}$
and $s_{12}=0$. 
\item If $0<a_{1}<\frac{K}{2}+1$ and $a_{2}>0$ and $s_{12}\in\{0,2,\ldots,2\min\{a_{1},a_{2}\}\}$,
then we obtain $s_{11}=a_{1}-\frac{s_{12}}{2}$ and $s_{22}=a_{2}-\frac{s_{12}}{2}$. 
\item If $a_{1}=\frac{K}{2}+1$ and $a_{2}>0$ and $s_{12}\in\{0,2,\ldots,2\min\{a_{1},a_{2}\}\}$,
then we obtain $s_{11}=a_{1}-\frac{s_{12}}{2}$ and $s_{22}=a_{2}-\frac{s_{12}}{2}$. 
\end{enumerate}
In the cases 5 and 6, we further have to distinguish whether $s_{12}<2\min\{a_{1},a_{2}\}$
or $s_{12}=2\min\{a_{1},a_{2}\}$ and in the second of these cases
whether $a_{1}\leq a_{2}$ or $a_{1}>a_{2}$. This makes cases 5 and
6 the most extensive to prove. Since the technique is the same in
each case, we shall limit ourselves to case 5 with $s_{12}<2\min\{a_{1},a_{2}\}$
which incidentally covers the majority of all possible values of $a_{1},a_{2},s_{11},s_{22},s_{12}$. 
\begin{proof}
Due to the chosen values of $a_{1},a_{2},s_{11},s_{22},s_{12}$, we
have

$LHS$:

\[
k=a_{1},l=a_{2},r=\frac{s_{12}}{2}
\]

Hence the summand reads

\[
(\frac{K}{2}+1)\frac{(2\alpha_{1}\sigma_{11})^{a_{1}}}{(\frac{K}{2}+1-a_{1})!}\frac{(2\alpha_{2}\sigma_{22})^{a_{2}}}{(\frac{L}{2}-a_{2})!}\frac{1}{s_{12}!(a_{1}-\frac{s_{12}}{2})!(a_{2}-\frac{s_{12}}{2})!}\left(\frac{2^{2}\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{\frac{s_{12}}{2}}.
\]

$A$:

\[
k=a_{1},l=a_{2},r=\frac{s_{12}}{2}
\]

Hence the summand reads

\[
\frac{(2\alpha_{1}\sigma_{11})^{a_{1}}}{(\frac{K}{2}-a_{1})!}\frac{(2\alpha_{2}\sigma_{22})^{a_{2}}}{(\frac{L}{2}-a_{2})!}\frac{1}{s_{12}!(a_{1}-\frac{s_{12}}{2})!(a_{2}-\frac{s_{12}}{2})!}\left(\frac{2^{2}\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{\frac{s_{12}}{2}}.
\]

$B$:

\[
k=a_{1}-1,l=a_{2},r=\frac{s_{12}}{2}
\]

Hence the summand reads

\[
\frac{(2\alpha_{1}\sigma_{11})^{a_{1}}}{(\frac{K}{2}+1-a_{1})!}\frac{(2\alpha_{2}\sigma_{22})^{a_{2}}}{(\frac{L}{2}-a_{2})!}\frac{1}{s_{12}!(a_{1}-1-\frac{s_{12}}{2})!(a_{2}-\frac{s_{12}}{2})!}\left(\frac{2^{2}\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{\frac{s_{12}}{2}}.
\]

$C$:

\[
k=a_{1}-1,l=a_{2}-1,r=\frac{s_{12}}{2}-1
\]

Hence the summand reads

\[
2\frac{(2\alpha_{1})^{a_{1}}\sigma_{11}^{a_{1}-1}}{(\frac{K}{2}+1-a_{1})!}\frac{(2\alpha_{2})^{a_{2}}\sigma_{22}^{a_{2}-1}}{(\frac{L}{2}-a_{2})!}\frac{1}{(s_{12}-1)!(a_{1}-\frac{s_{12}}{2})!(a_{2}-\frac{s_{12}}{2})!}\left(\frac{2^{2}}{\sigma_{11}\sigma_{22}}\right)^{\frac{s_{12}}{2}-1}\sigma_{12}^{s_{12}}.
\]

Now we divide all four terms by

\[
\frac{(2\alpha_{1}\sigma_{11})^{a_{1}}}{(\frac{K}{2}-a_{1})!}\frac{(2\alpha_{2}\sigma_{22})^{a_{2}}}{(\frac{L}{2}-a_{2})!}\frac{1}{(s_{12}-1)!(a_{1}-1-\frac{s_{12}}{2})!(a_{2}-\frac{s_{12}}{2})!}\left(\frac{\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{\frac{s_{12}}{2}}2^{s_{12}-1}.
\]

After this simplification, we have

$LHS$:

\[
(\frac{K}{2}+1)\frac{1}{\frac{K}{2}+1-a_{1}}\frac{1}{s_{12}(a_{1}-\frac{s_{12}}{2})}\cdot2.
\]

$A$:

\[
\frac{1}{s_{12}(a_{1}-\frac{s_{12}}{2})}\cdot2.
\]

$B$:

\[
\frac{1}{\frac{K}{2}+1-a_{1}}\frac{1}{s_{12}}\cdot2.
\]

$C$:

\[
\frac{1}{\frac{K}{2}+1-a_{1}}\frac{1}{a_{1}-\frac{s_{12}}{2}}.
\]

To see that indeed $LHS=RHS=A+B+C$, we multiply all four of these
expressions by

\[
(\frac{K}{2}+1-a_{1})s_{12}(a_{1}-\frac{s_{12}}{2}).
\]
Then we obtain

$LHS$:

\[
2(\frac{K}{2}+1).
\]

$A$:

\[
2(\frac{K}{2}+1-a_{1}).
\]

$B$:

\[
2(a_{1}-\frac{s_{12}}{2}).
\]

$C$:

\[
s_{12}.
\]

Adding up the expressions obtained from $A,B,C,$ we see that what
remains of $RHS$ is $K+2$, same as of $LHS$. 
\end{proof}
This concludes the proof of the central limit theorem (theorem \ref{CLT}).

\section{\label{sec:Indep}Independent Case}

In this section we shall analyse a particular type of coupling matrix:

\[
J=\left[\begin{array}{cc}
J_{1} & 0\\
0 & J_{2}
\end{array}\right],
\]
where $J_{1},J_{2}>0$ and, additionally, $J_{\nu}<\frac{1}{\alpha_{\nu}}$.
In this case, there is no interaction between spins in different groups,
hence we expect stochastic independence of the random variables $\frac{1}{\sqrt{N_{1}}}\sum_{i=1}^{N_{1}}X_{i},\frac{1}{\sqrt{N_{2}}}\sum_{j=1}^{N_{2}}Y_{j}$.
The inverse matrix $L=J^{-1}$ is given by

\[
L=\left[\begin{array}{cc}
\frac{1}{J_{1}} & 0\\
0 & \frac{1}{J_{2}}
\end{array}\right].
\]

The conditions on the coupling matrix reduce to

\begin{align*}
J_{1} & <\frac{1}{\alpha_{1}},\\
J_{2} & <\frac{1}{\alpha_{2}}.
\end{align*}

Under these assumptions, we calculate the correlations $\mathbb{E}(X_{1}\cdots X_{K}Y_{1}\cdots Y_{L})$
and the moments $M_{K,L}$.

\subsection{Correlations $\mathbb{E}(X_{1}\cdots X_{K}Y_{1}\cdots Y_{L})$ }

We show the following theorem: 
\begin{thm}
\label{thm:correlations}For all $K,L=\mathbb{N}_{0}$, the expected
value $\mathbb{E}(X_{1}\cdots X_{K}Y_{1}\cdots Y_{L})$ is asymptotically
given by the expression

\[
(K-1)!!\left(\frac{J_{1}}{1-\alpha_{1}J_{1}}\right)^{K/2}(L-1)!!\left(\frac{J_{2}}{1-\alpha_{2}J_{2}}\right)^{L/2}N^{-\frac{K+L}{2}}
\]
if both $K$ and $L$ are even and zero otherwise. 
\end{thm}

Note that this result is consistent with the correlations given for
a single group of spins, hence $\alpha=1$, in a Curie-Weiss model
in Method moments, with $J_{1}$ and $J_{2}$ taking the place of
the inverse temperature $\beta$. As expected, we have $\mathbb{E}(X_{1}\cdots X_{K}Y_{1}\cdots Y_{L})=\mathbb{E}(X_{1}\cdots X_{K})\mathbb{E}(Y_{1}\cdots Y_{L})$. 
\begin{proof}
We calculate the second derivatives of function $F$ defined in (\ref{eq:F-fn})

\begin{align*}
F_{11} & =2(L_{1}-\alpha_{1}),\\
F_{22} & =2(L_{2}-\alpha_{2}),\\
F_{12} & =2\bar{L},
\end{align*}
which we will need in this proof.

According to theorem \ref{thm:EXY-het}, the expected value $\mathbb{E}(X_{1}\cdots X_{K}Y_{1}\cdots Y_{L})$
is asymptotically equal to

\[
(\frac{2}{N})^{\frac{K+L}{2}}m_{K,L},
\]
where $m_{K,L}$ is the moment of a bivariate normal distribution
with zero mean and covariance matrix

\begin{align*}
H^{-1} & =\frac{1}{\det H}\left[\begin{array}{cc}
F_{22} & -F_{12}\\
-F_{12} & F_{11}
\end{array}\right]\\
 & =\frac{1}{F_{11}F_{22}}\left[\begin{array}{cc}
F_{22} & 0\\
0 & F_{11}
\end{array}\right]\\
 & =\left[\begin{array}{cc}
\frac{1}{F_{11}} & 0\\
0 & \frac{1}{F_{22}}
\end{array}\right],
\end{align*}
where we used $F_{12}=-\bar{L}=-\frac{\bar{J}}{J_{1}J_{2}-\bar{J}^{2}}=0$.
Since $H^{-1}$ is diagonal, the moment $m_{K,L}$ is zero, unless
both $K$ and $L$ are even. Then

\begin{align*}
m_{K,L} & =(K-1)!!\left(\frac{1}{F_{11}}\right)^{\frac{K}{2}}(L-1)!!\left(\frac{1}{F_{22}}\right)^{\frac{L}{2}}\\
 & =(K-1)!!\left(\frac{1}{2(L_{1}-\alpha_{1})}\right)^{\frac{K}{2}}(L-1)!!\left(\frac{1}{2(L_{2}-\alpha_{2})}\right)^{\frac{L}{2}}\\
 & =(K-1)!!\left(\frac{1}{\frac{1}{J_{1}}-\alpha_{1}}\right)^{\frac{K}{2}}2^{-\frac{K}{2}}(L-1)!!\left(\frac{1}{\frac{1}{J_{2}}-\alpha_{2}}\right)^{\frac{L}{2}}2^{-\frac{L}{2}}\\
 & =2^{-\frac{K+L}{2}}(K-1)!!\left(\frac{J_{1}}{1-\alpha_{1}J_{1}}\right)^{\frac{K}{2}}(L-1)!!\left(\frac{J_{2}}{1-\alpha_{2}J_{2}}\right)^{\frac{L}{2}}.
\end{align*}
Multiplying by $(\frac{2}{N})^{\frac{K+L}{2}}$ yields the desired
result. 
\end{proof}

\subsection{Moments $M_{K,L}$ and the Central Limit Theorem}

We use the correlations given by theorem \ref{thm:correlations} to
calculate the moments

\[
\mathbb{E}(\frac{1}{N_{1}^{K/2}}(\sum_{i=1}^{N_{1}}X_{i})^{K}\frac{1}{N_{2}^{L/2}}(\sum_{j=1}^{N_{2}}Y_{j})^{L})
\]
and show 
\begin{thm}
The moments for $K$ and $L$ even are given asymptotically by

\begin{align*}
\mathbb{E}(\frac{1}{N_{1}^{K/2}}(\sum_{i=1}^{N_{1}}X_{i})^{K}\frac{1}{N_{2}^{L/2}}(\sum_{j=1}^{N_{2}}Y_{j})^{L}) & \approx\\
(K-1)!!\left(\frac{1}{1-\alpha_{1}J_{1}}\right)^{\frac{K}{2}}(L-1)!!\left(\frac{1}{1-\alpha_{2}J_{2}}\right)^{\frac{L}{2}}.
\end{align*}
Therefore, the random vectors

\[
\left(\frac{1}{N_{1}^{1/2}}\sum_{i=1}^{N_{1}}X_{i},\frac{1}{N_{2}^{1/2}}\sum_{j=1}^{N_{2}}Y_{j}\right){}_{N}
\]
converge in distribution to a bivariate normal random vector with
zero mean and covariance matrix

\[
\left[\begin{array}{cc}
\frac{1}{1-\alpha_{1}J_{1}} & 0\\
0 & \frac{1}{1-\alpha_{2}J_{2}}
\end{array}\right].
\]
\end{thm}

\begin{proof}
We calculate $M_{K,L}$ when $K,L$ are both even.

\begin{align*}
M_{K,L} & \approx\frac{K!L!}{2^{\frac{K+L}{2}}}\sum_{k=0}^{K/2}\sum_{l=0}^{L/2}\frac{(2\alpha_{1}\sigma_{11})^{k}}{(\frac{K}{2}-k)!}\frac{(2\alpha_{2}\sigma_{22})^{l}}{(\frac{L}{2}-l)!}\cdot\\
 & \cdot\sum_{r=0}^{\min\{k,l\}}\frac{1}{(2r)!(k-r)!(l-r)!}\left(\frac{2^{2}\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{r}\\
 & =\frac{K!L!}{2^{\frac{K+L}{2}}}\sum_{k=0}^{K/2}\sum_{l=0}^{L/2}\frac{(2\alpha_{1}\sigma_{11})^{k}}{(\frac{K}{2}-k)!}\frac{(2\alpha_{2}\sigma_{22})^{l}}{(\frac{L}{2}-l)!}\sum_{r=0}^{0}\frac{1}{(2r)!(k-r)!(l-r)!}\left(\frac{2^{2}\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{r}\\
 & =\frac{K!L!}{2^{\frac{K+L}{2}}}\sum_{k=0}^{K/2}\sum_{l=0}^{L/2}\frac{(2\alpha_{1}\sigma_{11})^{k}}{(\frac{K}{2}-k)!}\frac{(2\alpha_{2}\sigma_{22})^{l}}{(\frac{L}{2}-l)!}\frac{1}{k!l!}\\
 & =\frac{K!}{(\frac{K}{2})!2^{\frac{K}{2}}}\frac{L!}{(\frac{L}{2})!2^{\frac{L}{2}}}\sum_{k=0}^{K/2}\sum_{l=0}^{L/2}\frac{(\frac{K}{2})!(2\alpha_{1}\sigma_{11})^{k}}{(\frac{K}{2}-k)!k!}\frac{(\frac{L}{2})!(2\alpha_{2}\sigma_{22})^{l}}{(\frac{L}{2}-l)!l!}\\
 & =(K-1)!!(L-1)!!\sum_{k=0}^{K/2}\left(\begin{array}{c}
\frac{K}{2}\\
k
\end{array}\right)(2\alpha_{1}\sigma_{11})^{k}\sum_{l=0}^{L/2}\left(\begin{array}{c}
\frac{K}{2}\\
k
\end{array}\right)(2\alpha_{2}\sigma_{22})^{l}\\
 & =(K-1)!!(1+2\alpha_{1}\sigma_{11})^{\frac{K}{2}}(L-1)!!(1+2\alpha_{2}\sigma_{22})^{\frac{L}{2}}\\
 & =(K-1)!!(1+2\alpha_{1}\frac{F_{22}}{F_{11}F_{22}-F_{12}^{2}})^{\frac{K}{2}}(L-1)!!(1+2\alpha_{2}\frac{F_{11}}{F_{11}F_{22}-F_{12}^{2}})^{\frac{L}{2}}\\
 & =(K-1)!!(1+2\alpha_{1}\frac{1}{F_{11}})^{\frac{K}{2}}(L-1)!!(1+2\alpha_{2}\frac{1}{F_{22}})^{\frac{L}{2}}\\
 & =(K-1)!!(1+\alpha_{1}\frac{1}{L_{1}-\alpha_{1}})^{\frac{K}{2}}(L-1)!!(1+\alpha_{2}\frac{1}{L_{2}-\alpha_{2}})^{\frac{L}{2}}\\
 & =(K-1)!!(1+\frac{\alpha_{1}}{\frac{1}{J_{1}}-\alpha_{1}})^{\frac{K}{2}}(L-1)!!(1+\frac{\alpha_{2}}{\frac{1}{J_{2}}-\alpha_{2}})^{\frac{L}{2}}\\
 & =(K-1)!!(1+\frac{\alpha_{1}J_{1}}{1-\alpha_{1}J_{1}})^{\frac{K}{2}}(L-1)!!(1+\frac{\alpha_{2}J_{2}}{1-\alpha_{2}J_{2}})^{\frac{L}{2}}\\
 & =(K-1)!!(\frac{1}{1-\alpha_{1}J_{1}})^{\frac{K}{2}}(L-1)!!(\frac{1}{1-\alpha_{2}J_{2}})^{\frac{L}{2}}
\end{align*}
The second equality above is due to the definition of $\Sigma$ in
equation (\ref{eq:Sigma}): $\sigma_{12}=\frac{-F_{12}}{\det H}=0$.
All summands with $r>0$ are zero. 
\end{proof}
Since the independent case is essentially the product of two one-dimensional
Curie-Weiss models, all the other results from that model hold as
well: 
\begin{thm}
\label{LLN-Indep} If $\bar{J}=0$, we have

\[
\big(\frac{1}{N_{1}}\sum_{i=1}^{N_{1}}X_{i},\frac{1}{N_{2}}\sum_{j=1}^{N_{2}}Y_{j}\big)~\underset{N\to\infty}{\Longrightarrow}\frac{1}{4}\sum_{t_{1},t_{2}\in\{-1,1\}}\delta_{t_{1}m(\alpha_{1}J_{1}),t_{2}m(\alpha_{2}J_{2})},
\]
where $m(\alpha_{\nu}J_{\nu})$ is the largest solution to the equation
in m

\[
\tanh(\alpha_{\nu}J_{\nu}m)=m.
\]
\end{thm}

Note that contrary to the model with homogeneous coupling matrices,
there is the possibility of the normalized sums converging to up to
four different points with probability $1/4$ each, depending on the
magnitude of $\alpha_{\nu}J_{\nu}$ (cf. theorem 1 in \cite{KT1}). 
\begin{thm}
If $\bar{J}=0$ and $\alpha_{\nu}J_{\nu}=1$, let $S_{\nu}$ be the
sum of spins in group $\{1,\ldots,N_{\nu}\}$. Then the normalized
sums

\[
\big(\frac{S_{\nu}}{N_{\nu}^{3/4}}\big)
\]
converge in distribution to a random variable with density function

\[
f(x)=2(12\alpha_{\nu})^{-1/4}\Gamma(1/4)^{-1}e^{-\frac{1}{12}x^{4}}.
\]
\end{thm}

This is a consequence of theorem 5.26 in \cite{MM}.

\section{Sublinear Population Growth\label{sec:Sublin}}

In this section we shall analyse the limiting distribution of the
sums

\[
(S_{1},S_{2})_{N}=(\frac{1}{\sqrt{N_{1}}}\sum_{i=1}^{N_{1}}X_{i},\frac{1}{\sqrt{N_{2}}}\sum_{j=1}^{N_{2}}Y_{j})_{N}
\]

if one or both groups grow at lower rates than the overall population
$N$. As in the previous sections, we allow for the presence of a
remainder population, i.e. $N_{1}+N_{2}\leq N$, where equality need
not hold.

Assume again that $\alpha_{1}=\lim\frac{N_{1}}{N}$ and $\alpha_{2}=\lim\frac{N_{2}}{N}$.
Now we allow one or both of these limits to be 0: let $\alpha_{1}=0$
and $0\leq\alpha_{2}\leq1$.

We shall show the following results: $S_{1}$ is asymptotically standard
normal, even though for finite $N$ the variance is of course strictly
greater than 1. In the large $N$ limit, $S_{1}$ and $S_{2}$ become
independent. Note that it suffices for this that one of the two groups
grow more slowly than at linear speed. Hence, contrary to the last
section, where we had positive moments for $K,L$ odd (such as the
covariance, e.g.), here only moments for $K,L$ even are positive.

The marginal moments $\mathbb{E}(S_{1}^{K})$ are given by

\begin{align*}
\approx & \frac{K!}{2^{\frac{K}{2}}}\sum_{k=0}^{K/2}\frac{(2\alpha_{1}\sigma_{11})^{k}}{(\frac{K}{2}-k)!}\sum_{r=0}^{\min\{k,0\}}\frac{1}{(2r)!(k-r)!(l-r)!}\left(\frac{2^{2}\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{r}\\
= & \frac{K!}{2^{\frac{K}{2}}}\sum_{k=0}^{K/2}\frac{(2\alpha_{1}\sigma_{11})^{k}}{(\frac{K}{2}-k)!}\frac{1}{k!}\\
= & \frac{K!}{2^{\frac{K}{2}}}\sum_{k=0}^{0}\frac{(2\alpha_{1}\sigma_{11})^{k}}{(\frac{K}{2}-k)!}\frac{1}{k!}\\
= & \frac{K!}{2^{\frac{K}{2}}}\frac{1}{(\frac{K}{2})!}\\
= & (K-1)!!
\end{align*}

The third equality is due to $\alpha_{1}=0$; only the summand with
$k=0$ contributes asymptotically to the value of the moment $\mathbb{E}(S_{1}^{K})$
. This shows that $S_{1}$ tends to a standard normal distribution.

The bivariate moment $M_{K,L}=\mathbb{E}(S_{1}^{K}S_{2}^{L})$ is
similar:

\begin{align*}
 & \frac{K!L!}{2^{\frac{K+L}{2}}}\sum_{k=0}^{K/2}\sum_{l=0}^{L/2}\frac{(2\alpha_{1}\sigma_{11})^{k}}{(\frac{K}{2}-k)!}\frac{(2\alpha_{2}\sigma_{22})^{l}}{(\frac{L}{2}-l)!}\sum_{r=0}^{\min\{k,l\}}\frac{1}{(2r)!(k-r)!(l-r)!}\left(\frac{2^{2}\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{r}\\
= & \frac{K!L!}{2^{\frac{K+L}{2}}}\sum_{k=0}^{0}\sum_{l=0}^{L/2}\frac{(2\alpha_{1}\sigma_{11})^{k}}{(\frac{K}{2}-k)!}\frac{(2\alpha_{2}\sigma_{22})^{l}}{(\frac{L}{2}-l)!}\sum_{r=0}^{\min\{0,l\}}\frac{1}{(2r)!(k-r)!(l-r)!}\left(\frac{2^{2}\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{r}\\
= & (K-1)!!\frac{L!}{2^{\frac{L}{2}}}\sum_{l=0}^{L/2}\frac{(2\alpha_{2}\sigma_{22})^{l}}{(\frac{L}{2}-l)!l!}\\
= & (K-1)!!\frac{L!}{(\frac{L}{2})!2^{\frac{L}{2}}}\sum_{l=0}^{L/2}\frac{(\frac{L}{2})!(2\alpha_{2}\sigma_{22})^{l}}{(\frac{L}{2}-l)!l!}\\
= & (K-1)!!(L-1)!!(1+2\alpha_{2}\sigma_{22})^{\frac{L}{2}},
\end{align*}

which implies a limiting bivariate normal distribution with zero mean
and covariance matrix

\[
\left[\begin{array}{cc}
1 & 0\\
0 & 1+2\alpha_{2}\sigma_{22}
\end{array}\right]=\left[\begin{array}{cc}
1 & 0\\
0 & \frac{1}{1-\alpha_{2}J_{2}}
\end{array}\right]
\]

As mentioned previously, it is enough if one of the two groups grows
more slowly for asymptotic independence to occur.

If we inspect the formula for odd $K,L$, $M_{K,L}=$

\begin{align*}
\frac{K!L!}{2^{\frac{K+L}{2}}}\cdot4\sqrt{\alpha_{1}\alpha_{2}}\sigma_{12}\sum_{k=0}^{\frac{K-1}{2}}\sum_{l=0}^{\frac{L-1}{2}}\frac{(2\alpha_{1}\sigma_{11})^{k}}{(\frac{K-1}{2}-k)!}\frac{(2\alpha_{2}\sigma_{22})^{l}}{(\frac{L-1}{2}-l)!}\cdot\\
\cdot\sum_{r=0}^{\min\{k,l\}}\frac{1}{(2r+1)!(k-r)!(l-r)!}\left(\frac{2^{2}\sigma_{12}^{2}}{\sigma_{11}\sigma_{22}}\right)^{r},
\end{align*}

we notice that owing to $\alpha_{1}=0$ the moment disappears.

\section{\label{sec:Special-Case}Special Case}

We first determine the extrema of the function $F$ defined in (\ref{eq:F-fn})
under the assumptions made in theorem \ref{thm:SpecialCase}: 
\begin{prop}
\label{prop:SC-Fminima}For $J_{1}=J_{2}=J$, $\alpha_{1}=\alpha_{2}=\alpha>0$,
and $J+\bar{J}>1/\alpha$, we set $L=\frac{J}{\Delta}$ and $\bar{L}=\frac{\bar{J}}{\Delta}$.
Then F has a local maximum at the origin and two global minima at
$(-m^{*},-m^{*})$ and $(m^{*},m^{*})$, where $m^{*}$ is the positive
solution to the equation

\[
(L-\bar{L})m-\sqrt{2}\alpha\tanh(\frac{m}{\sqrt{2}})=0.
\]
\end{prop}

In the proof of this proposition, we shall use the following lemma: 
\begin{lem}
These properties hold: 
\begin{enumerate}
\item $\sup_{t}\sin t\cos t=\frac{1}{2}$ 
\item For all $x\geq0$, we have $x\leq\sinh x\cosh x$, with equality if
and only if $x=0$. 
\item The function $g:[0,\infty)\rightarrow\mathbb{R},g(x)=\begin{cases}
\frac{1}{x}\tanh x, & x>0\\
1, & x=0
\end{cases}$, is strictly decreasing. 
\item For a fixed $r>0$, $\sup_{t}\cosh(r\sin t)\cosh(r\cos t)=\cosh(\frac{r}{\sqrt{2}})\cosh(\frac{r}{\sqrt{2}})$. 
\end{enumerate}
\end{lem}

\begin{proof}
\begin{enumerate}
\item Set $\varphi(t)=\sin t\cos t$. Then $\varphi'(t)=\cos^{2}t-\sin^{2}t=0$
implies $\cos t=\pm\sin t$. If $\cos t$ and $\sin t$ are of opposite
signs, then the product $\sin t\cos t$ will be negative, hence that
is a minimum of $\varphi$ and $\cos t=\sin t=\frac{1}{\sqrt{2}}$
is a maximum. 
\item We set $\varphi(x)=x$ and $\psi(x)=\sinh x\cosh x$. At the origin,
both functions are 0, hence the claim holds at the origin. The derivatives
are given by $\varphi'(x)=1$ and $\psi'(x)=\cosh^{2}x+\sinh^{2}x=(\cosh^{2}x-\sinh^{2}x)+2\sinh^{2}x=1+2\sinh^{2}x>1=\varphi'(x)$
for all $x>0$. Therefore, $\psi'(x)\geq\varphi'(x)$ for all $x\geq0$
and the claim holds. 
\item The function $g$ is differentiable at all $x>0$ and $g'(x)=\frac{1}{x}\frac{1}{\cosh^{2}x}-\frac{1}{x^{2}}\tanh x=\frac{1}{x^{2}}\frac{1}{\cosh^{2}x}(x-\sinh x\cosh x)<0$
due to statement 2. 
\item We derive the function $\varphi(t)=\cosh(r\sin t)\cosh(r\cos t)$:
\begin{align*}
\varphi'(t) & =r\cos t\cosh(r\cos t)\sinh(r\sin t)-r\sin t\sinh(r\cos)\cosh(r\sin t)\\
 & =r\cos t\sin t\cosh(r\cos t)\cosh(r\sin t)\cdot\\
 & \cdot(\frac{1}{\sin t}\tanh(r\sin t)-\frac{1}{\cos t}\tanh(r\cos t))\\
 & =r^{2}\cos t\sin t\cosh(r\cos t)\cosh(r\sin t)(g(r\sin t)-g(r\cos t)).
\end{align*}
The inequality $\sin t\leq\cos t$ holds on the interval $[0,\pi/4]$.
It follows from statement 3 that $g(r\sin t)\geq g(r\cos t)$. On
the interval $[\pi/4,\pi/2]$, we have $\sin t\geq\cos t$ and hence
$g(r\sin t)\leq g(r\cos t)$. The function $\varphi$ is increasing
on $[0,\pi/4]$ and decreasing on $[\pi/4,\pi/2]$. It thus reaches
a maximum at $t=\pi/4$. Since $\cosh$ is an even function, the sign
of the arguments $r\sin t$ and $r\cos t$ does not affect the value
$\cosh(r\sin t)\cosh(r\cos t)$. Hence this pattern repeats on each
interval $[\frac{k\pi}{2},\frac{(k+1)\pi}{2}]$ and claim 4 follows. 
\end{enumerate}
\end{proof}
Now we turn to the proof of the proposition. 
\begin{proof}[Proof of Proposition \ref{prop:SC-Fminima}]
We change to polar coordinates: $y_{1}=r\cos t,y_{2}=r\sin t$, $F(r,t)=\frac{1}{2}(L-2\bar{L}\sin t\cos t)r^{2}-\alpha\ln(\cosh(r\sin t)\cosh(r\cos t))$.
Since the value $r$ is non negative, and taking into account the
statements 1 and 4 from the previous lemma, $F$ reaches a minimum
whenever the function

\[
f(r)=\frac{1}{2}(L-\bar{L})r^{2}-2\alpha\ln(\cosh(\frac{r}{\sqrt{2}}))
\]
does. We calculate its derivative

\begin{equation}
f'(r)=(L-\bar{L})r-\sqrt{2}\alpha\tanh(\frac{r}{\sqrt{2}}).\label{eq:foc}
\end{equation}

We set $\varphi(r)=(L-\bar{L})r$ and $\psi(r)=\sqrt{2}\alpha\tanh(\frac{r}{\sqrt{2}})$.
$\varphi(0)=\psi(0)$ so the equation holds at the origin. We calculate

\begin{align*}
\varphi'(r) & =L-\bar{L},\\
\psi'(r) & =\frac{\alpha}{\cosh^{2}(\frac{r}{\sqrt{2}})}.
\end{align*}

If $L-\bar{L}\geq\alpha$, then $\varphi'(r)>\psi'(r)$ for all $r>0$
and 0 is the only solution to (\ref{eq:foc}) and $f$ reaches a global
minimum at the origin. If on the other hand $L-\bar{L}<\alpha$, then
due to the strict concavity of the tanh function on $[0,\infty)$
and its strict convexity on $(-\infty,0]$ there are two other solutions
to (\ref{eq:foc}). Call the strictly positive solution $m$. The
other solution is then $-m$ and $f$ reaches its global minima at
$-m$ and $m$ and it reaches a local maximum at the origin. Since
the condition $L-\bar{L}<\alpha$ is equivalent to $\alpha>\frac{J}{\det J}-\frac{\bar{J}}{\det J}=\frac{J-\bar{J}}{J^{2}-\bar{J}^{2}}=\frac{1}{J+\bar{J}}$,
the proposition follows. 
\end{proof}
{"config": "arxiv", "file": "1806.06708.tex"}