  We present a fully relativistic calculation of the matter bispectrum at
second order in cosmological perturbation theory assuming a Gaussian primordial
curvature perturbation. For the first time we perform a full numerical
integration of the bispectrum for both baryons and cold dark matter using the
second-order Einstein-Boltzmann code, SONG. We review previous analytical
results and provide an improved analytic approximation for the second-order
kernel in Poisson gauge which incorporates Newtonian nonlinear evolution,
relativistic initial conditions, the effect of radiation at early times and the
cosmological constant at late times. Our improved kernel provides a percent
level fit to the full numerical result at late times for most configurations,
including both equilateral shapes and the squeezed limit. We show that baryon
acoustic oscillations leave an imprint in the matter bispectrum, making a
significant impact on squeezed shapes.
