We investigate the first and second order cosmological perturbation equations in f(R)
modified gravity theory and provide the equation of motion of second order scalar induced
gravitational waves. We find that the effects of modified gravity not only change the form of the
equation of motion of second order scalar induced gravitational waves but also contribute an
additional anisotropic stress tensor, composed of first order scalar perturbations, to the source
term of the gravitational waves. We calculate the energy density spectrum of second order scalar
induced gravitational waves in the HS model. Utilizing current pulsar timing array observational
data, we perform a rigorous Bayesian analysis of the parameter space of the HS model.