  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 $\vert f'(R_0)-1\vert < 2.09\times 10^{-4}$. This bound is
stronger than those from most of the astrophysical observations and even some
cosmological observations.
