Setup: Consider a space (or spacetime, the signature doesn't matter) $(M,g_{ab})$ of arbitrary dimension, and let $\Sigma_\lambda$ be a one-parameter (with parameter $\lambda$) family of surfaces of arbitrary dimension embedded in $(M,g_{ab})$. Let $h_{ab}$ denote the induced metric on the $\Sigma_\lambda$ (I assume this metric is nondegenerate, i.e. that the $\Sigma_\lambda$ are not null); then the extrinsic curvature of the $\Sigma_\lambda$ is given by $$ K^a_{bc} = {h_b}^d {h_c}^e \nabla_d {h_e}^a, $$ where $\nabla_a$ is the covariant derivative compatible with $g_{ab}$. The mean curvature is $K^a = h^{bc} K^a_{bc}$.
Define $\eta^a \equiv (\partial_\lambda)^a$ (so that $\eta^a$ is the "deviation vector field between infinitesimally nearby surfaces"); we may without loss of generality take $\eta^a$ to be normal to the $\Sigma_\lambda$ since any component of $\eta^a$ tangential to the surfaces just generates some diffeomorphism, which is uninteresting.
Goal: I want an elliptic equation for $\eta^a$ on the $\Sigma_\lambda$. Ultimately I'm interested in the case where the $\Sigma_\lambda$ are all minimal (more generally, extremal) surfaces.
Issue: I have been able to obtain such an equation, which I'm pretty sure is correct (see Method 1 below). However, this equation disagrees by an annoying factor of two with the equation I expect I should get from the second variation formula for minimal surfaces (Method 2 below). I'm asking if someone understand the reason for this discrepancy (did I screw something up? Is a step I've made incorrect?)
Method 1: By my own work, I can show that $\eta^a$ should satisfy the equation $$ {P^a}_b \Delta_\Sigma \eta^b + 2{S^a}_b \eta^b + R_{bcde} h^{bd} \eta^c P^{ea} - {P^a}_c \eta^b \nabla_b K^c = 0, $$ where $P_{ab} \equiv g_{ab} - h_{ab}$ is the projector onto the subspace orthogonal to the $\Sigma_\lambda$, ${P^a}_b \Delta_\Sigma \eta^b \equiv {P^a}_b h^{cd} \nabla_c({h_d}^e \nabla_e \eta^b)$ is the Laplacian on the normal bundle, $S^{ab} \equiv K^a_{cd} K^b_{ef} g^{ce} g^{df}$ is Simons' operator, and $R_{abcd}$ is the Riemann tensor of $g_{ab}$ (defined via $(\nabla_a \nabla_b - \nabla_b \nabla_a)u_c = {R_{abc}}^du_d$). It took me about two pages of work to derive this equation, so instead of writing it all up here let me instead provide a check of this equation (though I can provide the details if someone wishes). Consider the case where the $\Sigma_\lambda$ are codimension-one and normal to a geodesic congruence with unit normal vector field $\eta^a$. Then $h_{ab} = g_{ab} - \epsilon \eta_a \eta_b$, where $\epsilon = \eta^2 = \pm 1$ (depending on the signature of $\eta^a$). One can show that $\nabla_a \eta_b = -\eta_c K^c_{ab}$, and thus the expansion and shear of this congruence are $\theta = -\eta_a K^a$ and $\sigma_{ab} = -\eta_c K^c_{ab} - \theta h_{ab}/(d-1)$ (with $d$ the dimension of $M$). Contracting the above equation with $\eta_d$, it is then easy enough to show that one obtains $$ \frac{d\theta}{d\lambda} = -\frac{1}{d-1} \theta^2 - \sigma_{ab} \sigma^{ab} - R_{ab} \eta^a \eta^b, $$ where $R_{ab}$ is the Ricci tensor of $g_{ab}$. This is precisely the Raychaudhuri equation (with vanishing twist $\omega_{ab}$, as must be the case since by construction the geodesic congruence here is hypersurface orthogonal), which verifies that at least in this case the equation I derived is correct. (Even more explicitly, you can check the equation for e.g. concentric spheres in flat Euclidean space).
Note that in the particular case that the $\Sigma_\lambda$ are all taken to be minimal (or extremal), we have that $K^a = 0$, and thus $$ {P^a}_b \Delta_\Sigma \eta^b + 2{S^a}_b \eta^b + R_{bcde} h^{bd} \eta^c P^{ea} = 0. $$
Method 2: Denote the areas of the $\Sigma_\lambda$ by $A(\lambda) = \int_{\Sigma_\lambda} \epsilon_\lambda$, where $\epsilon_\lambda$ is the natural volume form on $\Sigma_\lambda$ (also assume the $\Sigma_\lambda$ are compact). From e.g. Colding and Minicozzi's book "A Course in Minimal Surfaces", chapter 1 section 8, if $\Sigma_0$ is a minimal surface, we have that $$ \frac{d^2 A(0)}{d\lambda^2} = -\int_{\Sigma_0} \eta_a L\eta^a |_{\lambda = 0}, $$ where the Jacobi operator $L$ has the form $$ L \eta^a = {P^a}_b \Delta_\Sigma \eta^b + {S^a}_b \eta^b + R_{bcde} h^{bd} \eta^c P^{ea} $$ (I had to convert from math notation to my physics notation here, but I think this expression is the same as Colding and Minicozzi's). Now, if all the $\Sigma_\lambda$ are minimal, then we must have $d^2 A(\lambda)/d\lambda^2 = 0$, and thus we should have that on every $\Sigma_\lambda$, $$ L\eta^a = {P^a}_b \Delta_\Sigma \eta^b + {S^a}_b \eta^b + R_{bcde} h^{bd} \eta^c P^{ea} = 0. $$ Note the issue here: this formula looks almost identical to the one I derived above, except for a factor of two discrepancy in the coefficient of Simons' operator. My formula must have that factor of two, otherwise it wouldn't reproduce the Raychaudhuri equation. However, the second variation formula from Colding and Minicozzi is well-established, and it seems to me like it ought to reproduce the same equation for $\eta^a$ on each $\Sigma_\lambda$. What's going on?