Constraints on the Einstein-Hilbert action The Einstein-Hilbert action is given by
$$S_{EH}=\frac{1}{2\kappa}\int\text{d}^4x\sqrt{-g}R$$
Let's say we want to replace it with an action that doesn't reproduce singularities for example $R\rightarrow \frac{1}{\alpha}\tanh\left(\alpha R\right)$ with $\alpha$ a dimensionless parameter. For small $\alpha$ this reproduces the regular EH action. I don't actually know if this removes singularities so ignore this bold claim. My question is what are the constraints that we can put on these replacements? Can you put any function there that reduces to the EH action in some regime? This is called $f(R)$ gravity I think. And could you, in theory, remove singularities this way?
 A: There are several constraints, it depends on what you want to do.
For example there is the Dolgov-Kawasaki stability criterion which states that:
$$f_{RR}(R)>0$$
where $f_R = \cfrac{df(R)}{dR}$, for the $f(R)$ model to be tachyonically stable. See for example the papers:
(2+1)-dimensional solutions in F(R) gravity
Energy conditions and stability in
f
(
R
)
theories of gravity with nonminimal coupling to matter
f(R) gravity: successes and challenges
Matter instability in modified gravity
There are also cosmological constraints:
The pattern of growth in viable f(R) cosmologies

and are discussed in great detail at section B. The action in this last paper contains the pure Einstein-Hilbert term so you may need to modify the constraints the authors give.
Moreover if ones dealing with black holes, since the entropy in $f(R)$ gravity is given by (see: Black hole entropy in scalar-tensor and f(R) gravity: an overview
 and Horizon thermodynamics in f(R) theory):
$$S(r_h) = \frac{A f_R(r_h)}{4G}$$
one may impose that $f_R(r_h) >0$ for a reasonable entropy which will lead to a relation for the parameters of your theory.
For the singularity part in the paper Non-trivial black hole solutions in f(R) gravitational theory
,the authors consider an arbitrary $f(R)$ theory and found that the leading order of the Kretschmann scalar is $O(r^{-2})$ which is softer than the Schwarzchild Kretschmann scalar $O(r^{-6})$.
A: Yes, you can modify the EH action by writing
$$\mathcal{S}_{EH} = \displaystyle \int \mathrm{d}^4 x\sqrt{-g} f(R)$$
where $f(R)$ is an arbitrary function of the Ricci scalar $R$. The appealing feature of this action is that it combines mathematical simplicity with a fair amount of generality. For instance if we take a series expansion of $f$
$$f(R) = \ldots + \frac{a_2}{R^2}+\frac{a_1}{R} − 2\Lambda + R + b_2 R^2 + b_3 R^3 + \ldots $$
Where the $a_i$ and $b_j$ coefficients have the appropriate dimensions. You ask about placing some function of the Ricci scalar in there, well, I would add caution here, for example the Starobinsky gravity formalism has the following form:
$$f(R)=R+\frac {R^{2}}{6M^{2}}$$
where $M$ has the dimensions of mass. This action corresponds to the potential
$$V(\phi) = \Lambda^4 \left(1 - e^{-\sqrt{2/3} \phi/M_p} \right)^2$$
If there is a good reason for writing a hyperbolic tangent then do so, but observational test seem to go along with the Taylor expansion form I wrote above.
