  [Abridged]We present high-resolution three-dimensional direct numerical
simulations of breaking waves solving the two-phase Navier-Stokes equations. We
investigate the role of the Reynolds and Bond numbers on the energy, bubble and
droplet statistics of strong plunging breakers, and explore the asymptotic
regimes at high Reynolds and Bond numbers to be compared with laboratory
breaking waves. Energetically, we show that the breaking wave transitions from
laminar to three-dimensional turbulent flow on a timescale that depends on the
turbulent Reynolds number up to a limiting value of $Re_\lambda \sim 100$,
consistent with the mixing transition observed in other canonical turbulent
flows. We characterize the role of capillary effects on the impacting jet and
ingested main cavity shape and subsequent fragmentation process. We confirm two
distinct power-law regimes in the bubble size distribution, separated by the
Hinze scale $r_H$. We also confirm the interplay between the generation of
bubbles with the turbulent dissipation rate and extend the buoyant-energetic
scaling of Deike et al. (2016) to account for further Bond number effects. We
show resolved bubbles up to one order of magnitude below the Hinze scale and
observe a good collapse of the numerical data, both above and below the Hinze
scale, compared to laboratory breaking wave Deane and Stokes (2002). We resolve
droplet statistics and our data show good agreement with recent experiments
(Erinin et al., 2009) in various statistics. We present velocity distributions
for the droplets, finding ejection velocities up to four times the phase speed
of the wave, which are produced during the most intense splashing events of the
breaking process. This study paves the way for future integrated studies of the
statistics of energetics, bubbles and droplets and their interdependences in
breaking ocean waves and multiphase turbulent flows.
