The rate of energy loss and orbital period decay of quasi- stable compact binary systems
are derived in f(R) theory of gravity using the method of a single vertex graviton emission
process from a classical source. After linearising the f(R) action written in an equivalent
scalar-tensor format in the Einstein frame, we identify the appropriate interaction terms between
the massless spin-2 tensor mode, massive scalar mode, and the energy momentum tensor. The
definition of the scalar field is related to the f(R) models. Then using the interaction vertex
we compute the rate of energy loss due to spin-2 quadrupole radiation, which comes out to be the
same as the Peter-Mathews formula with a multiplication factor, and also the energy loss due to
the scalar dipole radiation. The total energy loss is the sum of these two contributions. Our
derivation is most general as it is applicable for both arbitrary eccentricity of the binary
orbits and arbitrary mass of the scalar field. Using the derived theoretical formula for the
period decay of the binary systems, we compare the predictions of f(R) gravity and general
relativity for the observations of four binary systems, i.e. Hulse-Taylor Binary, PSR J1141-6545,
PSR J1738+0333, and PSR J0348+0432. Thus we put bound on three well-known f(R) dark energy
models, namely the Hu-Sawicki, the Starobinsky, and the Tsujikawa model. We get the best
constraint on f'(R
0)-1 (where R
0 is the scalar curvature of the Universe at the present
epoch) from the Tsujikawa model, i.e |f'(R
0)-1| < 2.09 × 10-4. This bound is
stronger than those from most of the astrophysical observations and even some
cosmological observations.