# How to compute gauge variation of expression?

Suppose I have a symmetric tensor field $$h_{\mu\nu}$$ I want to implement somehow the following gauge variation of this tensor field as follows

$$\delta h_{\mu\nu} = \nabla_{\mu}\epsilon_{\nu} + \nabla_{\nu}\epsilon_{\mu}$$

and with this operator compute various gauge variations of different objects.

For example, suppose I have the following tensor in a flat background(which is gauge-invariant by definition)

$$k_{\mu\nu} = \Box h_{\mu\nu} - \nabla_\mu \nabla^{\lambda} h_{\nu\lambda} - \nabla_\nu \nabla^{\lambda} h_{\mu\lambda} + \nabla_{\mu\nu}h_{\lambda}^{\lambda}$$

where $$\Box = \nabla_{a}\nabla^{a}$$

The gauge variation of the $$k_{\mu\nu}$$ tensor should be zero

$${\delta k_{\mu\nu}} = 0$$

How can I implement the variation operator in Wolfram Mathematica using xAct package? Is it possilbe?