Given a partition function of the form $Z\sim \Sigma_{\sigma} \exp(J_{ij}\sigma_{i}\sigma_{j})$

and a bimodal distribution of the couplings,

\begin{equation}P(J_{ij})=p\,\delta(J_{ij}+J)+(1-p)\,\delta(J_{ij}-J)\end{equation}

I'm attempting to calculate $\langle Z^{n}\rangle$ using the replica method, however it seems that the replicas under this distribution (as opposed to a Gaussian distribution, for example) don't become coupled together, leaving me with the result that $\langle Z^{n}\rangle =\langle Z\rangle ^{n}$ which doesn't seem right. (Maybe it is, but it strikes me as too simple of a solution)

**My attempt:**

$H^{a}=\sum_{ij} J^{a}_{ij}\sigma^{a}_{i}\sigma^{a}_{j}$ ($a$ is the replica index).

Then after averaging over the coupling disorder, I have the thermal average of something like

$\prod_{a=1}^{n}\prod_{\langle ij\rangle}\left[p\exp(-J^{a}\sigma^{a}_{i}\sigma^{a}_{j})+(1-p)\exp(J^{a}\sigma^{a}_{i}\sigma^{a}_{j})\right]$

Do I need to perform an additional average for each replica? As it stands, the replicas are independent so $\langle Z^{n}\rangle$ seems wrong (it indicates the system is self-averaging which it shouldn't be). For other distributions, you typically get a term of the form $H\sim J\sigma^{a}_{i}\sigma^{b}_{i}$ or something similar.