The usual approach to cosmological perturbations relies on the so-called “slicing and threading” [1]}, which ruins the covariance of the picture. This raises the issue of gauge-dependence of cosmological perturbations, which is tackled either by choosing a particular gauge or by cleverly defining gauge-invariant perturbation variables [2]}. An alternative covariant approach to cosmological perturbations was first suggested by Hawking [3]} and later developed by others [4]}, [5]}, [6]}, [7]}. This approach is based on the so-called 1+3 covariant decomposition approach, in which one writes down a generic form of the field equations in terms of a set of covariant quantities defined locally with respect to an arbitrary timelike observer [8]}. Homogeneous and isotropic FLRW background corresponds to the vanishing of some of these covariant quantities, which, when assumed to be up to first order of smallness, correspond to linear perturbations. Because all the quantities are by definition covariant, including the perturbation quantities, they are also manifestly gauge-invariant. This approach has been successfully implemented to \(f(R)\)  cosmology [9]}, [10]}, [11]}, [12]}. There are five scalar perturbation variables
\(\Delta _m \equiv \frac{a^2}{\rho _m}\tilde{\nabla }^2\rho _m, \quad Z \equiv 3a^2 \tilde{\nabla }^2 H, \quad C \equiv a^4 \tilde{\nabla }^2 \tilde{R}, \quad \mathcal {R} \equiv a^2 \tilde{\nabla }^2 R, \quad \Re \equiv a^2 \tilde{\nabla }\dot{R}.\) 
