unique [e.g., @Wagner2010]. Following @Wagner2010, we assume that the data can be represented as $$\begin{aligned}
X &= X_{\rm G} + X_{\rm NL} \label{eqXgenerate},\end{aligned}$$ where $\XG$ is a Gaussian field and $\XNL$ a small nonlinear component. In the Fourier domain, the Gaussian field follows a complex Gaussian distribution, $$\begin{aligned}
X_{\rm G}(\bk) &= \mathcal{N}_{\bk} \sqrt{\frac{1}{h_k^2} P(\bk)}, \label{eqXGgenerate}\end{aligned}$$ where $\mathcal{N}_{\bk}$ are Gaussian random variables with independent real and imaginary parts, zero mean and unit variance. We set $\mathcal{N}_{-\bk} = \mathcal{N}_{\bk}^*$ to ensure that $X_{\rm G}$ is real, but otherwise the $\mathcal{N}_{\bk}$ are independent for different $\bk$.

We assume that the nonlinear component is obtained from the Gaussian one by an arbitrary quadratic nonlinearity, $$\begin{aligned}
X_{\rm NL}(\bmm) &= \sum_{\bk} A_{\bk, \bmm -\bk} X_{\rm G}(\bk) X_{\rm G}(\bmm -\bk) \label{eqXNLgenerate},\end{aligned}$$ where $A_{\bk,\bl}$ denote coupling coefficients of modes $\bk$ and $\bl=\bmm-\bk$. To guarantee that $\XNL$ is real, the coupling coefficients have to fulfill $A_{-\bk,-\bl} = A_{\bk,\bl}^*$.

Synthetic data with a given bispectrum $B_0(\bk,\bl)$ can then be generated by setting $$\begin{aligned}
A_{\bk,\bl} &= \frac{1}{6} \frac{h_k^2\, B_0(\bk,\bl) }{ P(\bk)P(\bl)}. \label{eqAkl}\end{aligned}$$ Appendix \[appModel\] outlines further details of this model.

The bispectrum of hexagons
--------------------------

For the patterns introduced in Section \[secPatternsIntro\], the bispectrum of the patterns is given by $B(\bk_0,\bl_0) 