In string theory, the Polyakov string action is given by \begin{align} S = -\frac{T}{2}\int d^2\sigma\:\sqrt{-\gamma}\gamma^{ab}\eta_{\mu\nu}\frac{\partial X^{\mu}}{\partial \sigma^a}\frac{\partial X^{\nu}}{\partial \sigma^b}, \end{align} where $\gamma^{ab}$ is an independent worldsheet metric that has a residual conformal symmetry, and $\gamma = \text{det}\gamma_{ab}$. By performing an infinitesimal worldsheet reparameterisation of the form $\tilde{\sigma}^a = \sigma^a - \epsilon^a(\sigma)$, I have been told that when the variation of the action vanishes, I should find \begin{align} \partial^a(T_{ab}\epsilon^b) = 0. \end{align} $T_{ab}$ is the stress-energy tensor given by \begin{align} T_{ab} = \eta_{\mu\nu}\frac{\partial X^{\mu}}{\partial \sigma^a}\frac{\partial X^{\nu}}{\partial \sigma^b} - \frac{1}{2}\gamma_{ab}\gamma^{cd}\eta_{\mu\nu}\frac{\partial X^{\mu}}{\partial \sigma^c}\frac{\partial X^{\nu}}{\partial \sigma^d}. \end{align}
I've not been able to prove this so far, and at this point I'm stuck. I was wondering if anyone would be able to show me how it is done?
My attempt:
First calculate \begin{align} \tilde{\gamma}^{ab}(\tilde{\sigma}) &= \frac{\partial\tilde{\sigma}^a}{\partial \sigma^c} \frac{\partial\tilde{\sigma}^b}{\partial \sigma^d}\gamma^{cd}(\sigma) \nonumber\\ &= (\delta^a_c - \partial_c\epsilon^a(\sigma))(\delta^b_d - \partial_d\epsilon^b(\sigma))\gamma^{cd}(\sigma) \nonumber\\ &= \gamma^{ab}(\sigma) - \partial_c\epsilon^a(\sigma)\gamma^{cb}(\sigma) - \partial_d\epsilon^b(\sigma)\gamma^{ad}(\sigma)\nonumber\\ &= - \partial^b\epsilon^a(\sigma) - \partial^a\epsilon^b(\sigma), \end{align} and so \begin{align} \delta\gamma^{ab} = - \partial^b\epsilon^a(\sigma) - \partial^a\epsilon^b(\sigma). \end{align} Use the result \begin{align} \delta \sqrt{-\gamma} = -\frac{1}{2}\sqrt{-\gamma} \gamma_{ab}(\sigma)\delta\gamma^{ab}. \end{align} Vary the action with respect to the reparameterisation: \begin{align} \delta S =& \delta\left( -\frac{T}{2}\int d\sigma^2\:\sqrt{-\gamma}\gamma^{ab}(\sigma)\eta_{\mu\nu}\frac{\partial X^{\mu}}{\partial \sigma^a}\frac{\partial X^{\nu}}{\partial \sigma^b}\right) \nonumber\\ =& -\frac{T}{2}\int d\sigma^2\:\left[\delta\sqrt{-\gamma}\gamma^{ab}(\sigma) + \sqrt{-\gamma}\delta\gamma^{ab}\right]\eta_{\mu\nu}\frac{\partial X^{\mu}}{\partial \sigma^a}\frac{\partial X^{\nu}}{\partial \sigma^b} \nonumber\\ =& -\frac{T}{2}\int d\sigma^2\:\left[-\frac{1}{2}\sqrt{-\gamma} \gamma_{cd}(\sigma)\delta\gamma^{cd}\gamma^{ab}(\sigma) + \sqrt{-\gamma}\delta\gamma^{ab}\right]\eta_{\mu\nu}\frac{\partial X^{\mu}}{\partial \sigma^a}\frac{\partial X^{\nu}}{\partial \sigma^b} \nonumber\\ =& -\frac{T}{2}\int d\sigma^2\:\sqrt{-\gamma} \underbrace{\left[ -\frac{1}{2}\gamma_{ab}(\sigma)\gamma^{cd}(\sigma)\eta_{\mu\nu}\frac{\partial X^{\mu}}{\partial \sigma^c}\frac{\partial X^{\nu}}{\partial \sigma^d} + \eta_{\mu\nu}\frac{\partial X^{\mu}}{\partial \sigma^a}\frac{\partial X^{\nu}}{\partial \sigma^b}\right]}_{T_{ab}(\sigma)} \delta\gamma^{ab} \nonumber\\ = & -\frac{T}{2}\int d\sigma^2\:\sqrt{-\gamma}T_{ab}(\sigma)(- \partial^b\epsilon^a(\sigma) - \partial^a\epsilon^b(\sigma)) \nonumber\\ =& T\int d\sigma^2\:\sqrt{-\gamma}T_{ab}(\sigma)\partial^a\epsilon^b(\sigma), \end{align} where in the last line I used the symmetry of $T_{ab}(\sigma)$.
Now using the continuity equation $\partial^aT_{ab}(\sigma) = 0$ for the stress-energy tensor, \begin{align} \delta S_P =& T\int d\sigma^2\:\sqrt{-\gamma}T_{ab}(\sigma)\partial^a\epsilon^b(\sigma)\nonumber\\ =& T\int d\sigma^2\:\sqrt{-\gamma}\left[\partial^a(T_{ab}(\sigma)\epsilon^b(\sigma)) - \underbrace{\partial^aT_{ab}(\sigma)}_{0}\epsilon^b(\sigma)\right]\nonumber\\ =& T\int d\sigma^2\:\sqrt{-\gamma}\partial^a(T_{ab}(\sigma)\epsilon^b(\sigma)). \end{align}