We use the functional renormalization group approach with partial bosonization in the particleparticle channel to study the effect of order parameter fluctuations on the BCS-BEC crossover of superfluid fermions in three dimensions. Our approach is based on a new truncation of the vertex expansion where the renormalization group flow of bosonic two-point functions is closed by means of Dyson-Schwinger equations and the superfluid order parameter is related to the single particle gap via a Ward identity. We explicitly calculate the chemical potential, the single-particle gap, and the superfluid order parameter at the unitary point and compare our results with experiments and previous calculations.PACS numbers: 03.75.Ss, 03.75.Hh, 74.20.Fg The BCS-BEC crossover in a two-component Fermi gas has attracted the attention of theorists for several decades [1,2,3,4,5]. It is generally accepted that the nature of the superfluid state exhibits a smooth crossover as a function of the dimensionless parameter 1/k F a s , where k F is the Fermi momentum and a s is the s-wave scattering length in vacuum. While for a small negative scattering length, i.e., 1/k F a s ≪ −1, the paired state generated by the attractive interaction is a collection of spatially extended Cooper pairs (BCS limit), in the opposite limit 1/k F a s ≫ 1 the superfluid state can be viewed as a Bose-Einstein condensate of tightly bound fermion pairs (BEC limit). Of particular interest is the unitary point 1/k F a s = 0 where the scattering length diverges and k F sets the only length scale of the system. In this regime quantitative calculations are difficult because there is no small parameter to justify approximations.In the past few years several observables such as the chemical potential and the quasi-particle gap have been determined experimentally at the unitary point [6,7,8,9,10], but there is still some uncertainty in the precise numerical values of these quantities. The unitary point has also been studied theoretically using Monte Carlo simulations [11,12] and various analytical methods based on field theoretical techniques [13,14,15,16,17] or the functional renormalization group (FRG) [18,19,20,21], but also in this case the theoretical results have not converged yet. In such a situation it is desirable to study this problem using new approximation strategies which are complementary to previous calculations. In this work we shall therefore develop a novel FRG approach for the BCS-BEC crossover which is based on a suitable truncation of the vertex expansion using skeleton equations and Ward identities. For fixed density we explicitly calculate the chemical potential, the single-particle gap, and the superfluid order-parameter at the unitary point in three dimensions, and compare our results with experiments and with previous calculations.We consider a system of neutral fermions with energy dispersion ǫ k = k 2 /(2m) and a short-range attractive two-body interaction g p depending on the total momentum p of a fermion pair. After decoupling the intera...