  We theoretically study the thermodynamic properties of a strongly interacting
Fermi gas at the crossover from a Bardeen-Cooper-Schrieffer (BCS) superfluid to
a Bose-Einstein condensate (BEC), by applying a recently outlined
strong-coupling theory that includes pair fluctuations beyond the commonly-used
many-body $T$-matrix or ladder approximation at the Gaussian level. The beyond
Gaussian pair fluctuation (GPF) theory always respects the exact thermodynamic
relations and recovers the Bogoliubov theory of molecules in the BEC limit with
a nearly correct molecule-molecule scattering length. We show that the
beyond-GPF theory predicts quantitatively accurate ground-state properties at
the BEC-BCS crossover, in good agreement with the recent measurement by
Horikoshi \textit{et al.} in Phys. Rev. X \textbf{7}, 041004 (2017). In the
unitary limit with infinitely large $s$-wave scattering length, the beyond-GPF
theory predicts a reliable universal energy equation of state up to 0.6$T_c$,
where $T_c$ is the superfluid transition temperature at unitarity. The theory
predicts a Bertsch parameter $\xi \simeq 0.365$ at zero temperature, in good
agreement with the latest quantum Monte Carlo result $\xi = 0.367(7)$ and the
latest experimental measurement $\xi = 0.367(9)$. We attribute the excellent
and wide applicability of the beyond-GPF theory in the broken-symmetry phase to
the reasonable re-summation of Feynman diagrams following a dimensional
$\epsilon$-expansion analysis near four dimensions ($d=4-\epsilon$), which
gives rise to accurate predictions at the second order
$\mathcal{O}(\epsilon^2)$. Our work indicates the possibility of further
improving the strong-coupling theory of strongly interacting fermions based on
the systematic inclusion of large-loop Feynman diagrams at higher orders
$\mathcal{O}(\epsilon^n)$ with $n\ge 3$.
