Recently, we have developed an unconventional subface-based Finite Volume (FV) method [Gallice et al., J. Comp. Phys., 2022] for solving the multidimensional Euler equations over general unstructured grids. The building block of this numerical algorithm consists of a parameterized approximate Riemann solver constructed in the normal direction to the subface. This Riemann solver is firstly written under Lagrangian representation which is more convenient to impose positivity preserving and entropy stability by monitoring its wave speeds. Then, its Eulerian counterpart is deduced applying the Lagrange-to-Euler mapping consistently with the systematic methodology introduced in [Gallice Numer. Math., 2003]. One of the main advantage of this approach lies in the fact that the good properties of the Lagrangian Riemann solver are automatically transfered to its Eulerian counterpart. Once the positivity preserving and entropy stability properties are ensured at the approximate Riemann solver level, it is quite easy to extend them at the global level of the FV scheme simply by expressing the updated cell-centered value of the conservative variables as a convex combination under an explicit time step constraint. It is well known that the robustness of numerical methods for hyperbolic conservation laws is greatly related to their ability of satisfying a so called entropy inequality which ensures their entropy stability. In this presentation, we aim at investigating thoroughly the entropy production embedded in a class of approximate Riemann solvers dedicated to gas dynamics. We shall provide an explicit expression of this entropy production and show how to control it by modifying slightly the approximate Riemann solvers under consideration. We shall also exhibit a class of approximate Riemann solvers which are able to conserve entropy rigorously. The use of such entropy conservative Riemann solvers results in a FV numerical method which remains able to run quite demanding test cases characterized by strong shock and rarefaction waves without any robustness problems. We infer that the entropy stability of the FV method in the classical sense is not really needed. Indeed, the main source of entropy production, which is not controlled, is due to computation of the kinetic energy by squaring the averaged velocity obtained from the intermediate states of the approximate Riemann solver.