$p\leq [(D-1)/2]$($[N]$ denotes the integral part of the number $N$), $c_k$ are arbitrary constants with dimension of $[\mathrm{Length}]^{2k-2}$, and ${\cal L}_k$ are the Euler densities $${\cal L}_k=\frac{1}{2^k}\sqrt{-g}\delta^{\mu_1\cdots
\mu_k\nu_1\cdots\nu_k}_{\lambda_1\cdots
\lambda_k\sigma_1\cdots\sigma_k}R^{\lambda_1\sigma_1}_{~~~~
\mu_1\nu_1}\cdots R^{\lambda_k\sigma_k}_{~~~~ \mu_k\nu_k}\, .
\label{Lk}$$ Here, the generalized delta is totally antisymmetric in both sets of indices. ${\cal L}_0=1$, so the constant $c_0$ is just the cosmological constant. ${\cal L}_1$ gives us the usual curvature scalar term, and for simplicity, we set $c_1=1$, while ${\cal L}_2$ is just the Gauss-Bonnet term. The Gauss-Bonnet term is argued to appear in the low-energy action of heterotic string theory with a positive coefficient [@Zwie]. The equations of motion following from the Lagrangian (\[1eq1\]) have the form $\mathcal{G}_{\mu\nu}=0$, where $$\mathcal{G}^{\mu}_{\nu}=\sum_{k=0}^{p}\frac{1}{2^{k+1}}
c_{k}\delta^{\mu\mu_1\cdots
\mu_k\nu_1\cdots\nu_k}_{\nu\lambda_1\cdots
\lambda_k\sigma_1\cdots\sigma_k}R^{\lambda_1\sigma_1}_{~~~~
\mu_1\nu_1}\cdots R^{\lambda_k\sigma_k}_{~~~~ \mu_k\nu_k} \,
.\label{em}$$ Since the action of Lovelock gravity is the sum of the dimensionally extended Euler densities, there are no more than second order derivatives with respect to metric in its equations of motion. Furthermore, the Lovelock gravity is shown to be free of ghost when expanded on a flat space, evading any problems with unitarity [@Zwie; @Des].

Finding exact analytic solutions of any gravity theory is an issue of long-standing interest. Indeed, there exist a lot of works 