Let $(M^3,g)$ be a complete, orientable and connected Riemannian $3$-manifold and let $\Sigma^2 \subset M$ be a compact orientable totally geodesic surface embedded in $M$. For $f \in C^{2,\alpha}(\Sigma)$ with small $C^0$ norm, let

$$\Sigma_f = \{ \operatorname{exp}_p(f(p)N(p)) : p \in \Sigma \}$$

be the normal graph of $f$ over $\Sigma$, where, of course, $N$ is a unit normal field for $\Sigma$. What equation does $f$ obey if $\Sigma_f$ is a minimal surface in $(M,g)$?

My thoughts: try to write the first variation formula for the area of $\Sigma_f$. For this, define

$$\varphi_t(p) = \operatorname{exp}_p((f(p) + t\eta(p))N(p)), \quad p \in \Sigma$$

where $\eta \in C^{\infty}(\Sigma)$ is any smooth function and $t$ is small. To compute the area of $\Sigma_f(t) := \varphi_t(\Sigma)$, we need to evaluate the following integral:

$$\operatorname{Area}(\Sigma_f(t)) = \int_\Sigma \sqrt{E_tG_t - F_t^2} dA$$

where $E_t, F_t$ and $G_t$ are the coefficients of the first fundamental form of $\varphi_t$. For instance, if $\{e_1, e_2\}$ is a local orthonormal referential for $\Sigma$, then

$$E_t = g\left( \mathrm{d} \varphi_t(e_1), \mathrm{d} \varphi_t(e_1) \right).$$

However, this involves the derivatives of the exponential map with respect to both base point and vector argument. How do I proceed? Do you have a faster argument?