Given an action
$$S[\psi]=\int d^4x \mathcal{L}(\psi,\nabla_\mu \psi),$$
how do we find the equations of motion? The principle of stationary action says that the equations of motion can be found by extremizing the action, i.e., by setting its first-order change with respect to a change in its argument to zero. That is, we demand that
$$\left.\frac{dS[\psi+\alpha \delta \psi]}{d\alpha}\right|_{\alpha=0}=0.$$
Note that $\delta \psi$ need not be small, it is an arbitrary variation. The only thing we demand from $\delta\psi$ is that it must vanish at infinity. Let us look at the action you provided:
$$
S\left[\varphi, g_{\mu\nu}\right]=\int d^4x \sqrt{-g} e^{-2 \varphi}\left[R+4 g_{\mu \nu} \nabla^\mu \varphi \nabla^\nu \varphi\right]
$$
Set $$\left.\frac{d S\left[\varphi+\alpha \delta \varphi, g_{\mu \nu}\right]}{d \alpha}\right|_{\alpha=0}$$ equal to zero to get
\begin{equation}
\begin{aligned}
& 0=\frac{d}{d\alpha}|_{\alpha=0}\int d^4x \sqrt{-g} e^{-2 (\varphi+\alpha \delta \varphi)}\left[R+4 g_{\mu \nu} \nabla^\mu (\varphi+\alpha \delta \varphi) \nabla^\nu (\varphi+\alpha \delta \varphi)\right]\\
&=\int d^4 x \sqrt{-g} e^{-2 \varphi}\left\{\left[R+4 \nabla_\mu \varphi \nabla^\mu \varphi\right](-2 \delta \varphi)+8 \nabla^\mu \delta \varphi \nabla_\mu \varphi\right\} \\
&=\int d^4 x \sqrt{-g} e^{-2 \varphi}\left\{-2 R+8 \nabla^\mu \varphi \nabla_\mu \varphi-8 \square \phi\right\} \delta \varphi \\
& \Rightarrow R-4 \nabla^\mu \varphi \nabla_\mu \varphi+4 \square \varphi=0
\end{aligned}
\qquad(1)
\end{equation}
where $g_{\mu\nu}\nabla^{\mu}\nabla^{\nu}=\Box$. I have made use of
$$
\begin{aligned}
\nabla^\mu\left(\delta \varphi \nabla_\mu \varphi e^{-2 \varphi}\right)= & e^{-2 \varphi} \nabla^\mu \delta \varphi \nabla_\mu \varphi+\Box \phi \delta \varphi e^{-2 \varphi} -2 \nabla^\mu \varphi \nabla_\mu \varphi \delta \varphi e^{-2 \varphi} \\
\Rightarrow e^{-2 \varphi} \nabla^\mu \delta \varphi \nabla_\mu \varphi= & \delta \varphi e^{-2 \varphi}\left(2 \nabla^\mu \varphi \nabla_\mu \varphi-\Box \varphi\right)+\underbrace{\nabla^\mu\left(\delta \varphi \nabla_\mu \varphi e^{-2\varphi}\right)}_{\text{Integrates to an inconsequential boundary term}}
\end{aligned}
$$
Next, let us vary the action with respect to the metric tensor. As you probably know, we usually vary with respect to the inverse metric rather than the metric itself because it's usually easier to do it this way. Without further ado, let us write our action as
$$
S\left[\varphi, g^{\mu\nu}\right]=\int d^4x \sqrt{-g} e^{-2 \varphi}\left[R+4 g^{\mu \nu} \nabla_\mu \varphi \nabla_\nu \varphi\right]
$$
and set
$$\left.\frac{d S\left[\varphi, g^{\mu \nu}+\alpha \delta g^{\mu\nu}\right]}{d \alpha}\right|_{\alpha=0}$$ equal to zero:
$$
\begin{aligned}
& 0=\left.\frac{d}{d \alpha}\right|_{\alpha=0} \int d^4x \sqrt{-\operatorname{det}\left(g(\alpha))\right.}\left[R\left(g^{\mu \nu}+\alpha \delta g^{\mu\nu}\right)+4\left(g^{\mu \nu}+\alpha \delta g^{\mu\nu}\right) \nabla_\mu \varphi \nabla_\nu \varphi\right] e^{-2 \varphi} \\
& =\int d^4x\left[-\frac{1}{2} \sqrt{-g} g_{\mu \nu} \delta g^{\mu \nu}\left[R+4 g^{\alpha \beta} \nabla_\alpha \varphi \nabla_\beta \varphi \right]e^{-2 \varphi}+\sqrt{-g}\left[g^{\mu\nu}\frac{d R_{\mu\nu}\left(g^{\mu \nu}+\alpha \delta g^{\mu \nu}\right)}{d \alpha}\right|_{\alpha=0}\right. \\
& \left.+\delta g^{\mu\nu}R_{\mu\nu}+4 \delta g^{\mu \nu} \nabla_\mu \varphi \nabla_\nu \varphi\right] e^{-2 \varphi} \\
& =\int d^4 x \sqrt{-g} e^{-2 \varphi}\left[R_{\mu \nu}-\frac{1}{2} R g_{\mu \nu}+4 \nabla_\mu \varphi \nabla_\nu \varphi-2 g_{\mu\nu}g^{\alpha \beta} \nabla_\alpha \varphi \nabla_\beta \varphi \right]\delta g^{\mu\nu} \\
& \Rightarrow R_{\mu \nu}-\frac{1}{2} R g_{\mu \nu}=2 g_{\mu \nu} \nabla^\alpha \varphi \nabla_\alpha \varphi-4 \nabla_\mu \varphi \nabla_\nu \varphi. \\
&
\end{aligned}
$$
In going from the first equality to the second, I used
$$\frac{d}{d\alpha}\sqrt{-\operatorname{det}(g(\alpha))}|_{\alpha=0}=-\frac{1}{2}g_{\mu \nu}\delta g^{\mu\nu}\sqrt{-\operatorname{det}(g)}, $$
and
$$\left.\int d^4x \sqrt{-g}g^{\mu\nu}\frac{d R_{\mu\nu}\left(g^{\mu \nu}+\alpha \delta g^{\mu \nu}\right)}{d \alpha}\right|_{\alpha=0}=0.$$
Both of these are proven in every general relativity textbook. For example, you can take a look at chapter 4 of Carroll's book.
Now take the trace of both sides to find the Ricci scalar:
$$
\begin{aligned}
& R-2R=8 \nabla^\alpha \varphi \nabla_\alpha \varphi-4 \nabla^\alpha \varphi \nabla_\alpha \varphi=4 \nabla^\alpha \varphi \nabla_\alpha \varphi \\
& \Rightarrow R_{\mu \nu}=-4 \nabla_\mu \varphi \nabla_\nu \varphi.
\end{aligned}
$$
By the way, if we plug this back into the final equation in (1), we find that
\begin{equation}
\nabla^{\mu}\varphi\nabla_{\mu}\varphi=\frac{1}{2}\Box\varphi.
\qquad(2)
\end{equation}
To find the energy-momentum tensor, we simply write the equation of motion resulting from varying the action with respect to the metric in the form
$$R_{\mu \nu}-\frac{1}{2}Rg_{\mu\nu}=T_{\mu\nu}.$$
In our case,
$$T_{\mu\nu}=2 g_{\mu \nu} \nabla^\alpha \varphi \nabla_\alpha \varphi-4 \nabla_\mu \varphi \nabla_\nu \varphi=g_{\mu\nu}\Box\varphi-4 \nabla_\mu \varphi \nabla_\nu \varphi$$
where I have made use of (2).
I think the authors of the paper you quoted made a minor error. I think they took (2) to mean that
$$\nabla_{\mu}\varphi\nabla_{\nu}\varphi=\frac{1}{2}\nabla_{\mu}\nabla_{\nu}\varphi.$$
Go through everything and if you find a mistake, please let me know. I hope this helps.
