  Context. The spatial power spectrum of supergranulation does not fully
characterize the underlying physics of turbulent convection. For example, it
does not describe the non-Gaussianity in the horizontal flow divergence.
  Aims. Our aim is to statistically characterize the spatial pattern of solar
supergranulation beyond the power spectrum. The next-order statistic is the
bispectrum. It measures correlations of three Fourier components and is related
to the nonlinearities in the underlying physics.
  Methods. We estimated the bispectrum of supergranular horizontal surface
divergence maps that were obtained using local correlation tracking (LCT) and
time-distance helioseismology (TD) from one year of data from the Helioseismic
and Magnetic Imager on-board the Solar Dynamics Observatory starting in May
2010.
  Results. We find significantly nonzero and consistent estimates for the
bispectrum. The strongest nonlinearity is present when the three coupling wave
vectors are at the supergranular scale. These are the same wave vectors that
are present in regular hexagons, which were used in analytical studies of solar
convection. At these Fourier components, the bispectrum is positive, consistent
with the positive skewness in the data and with supergranules preferentially
consisting of outflows surrounded by a network of inflows. We use the
bispectrum to generate synthetic divergence maps that are very similar to the
data by a model that consists of a Gaussian term and a weaker quadratic
nonlinear component. Thereby, we estimate the fraction of the variance in the
divergence maps from the nonlinear component to be of the order of 4-6%.
  Conclusions. We propose that bispectral analysis is useful for understanding
solar turbulent convection, for example for comparing observations and
numerical models of supergranular flows. This analysis may also be useful to
generate synthetic flow fields.
