Let' s calculate one-loop counter terms of squared vector field operator using methods described in Calculating one-loop counterterms. The operator is: \begin{eqnarray*} &&D^2{}_\alpha{}^\beta = \delta_\alpha{}^\beta \Box^2 - \lambda \nabla_\alpha \nabla^\beta \Box + 2 P_\alpha{}^\beta \Box - \lambda \Box \nabla_\alpha \nabla^\beta + \lambda^2 \nabla_\alpha \Box \nabla^\beta \vphantom{\frac{1}{2}} -\nonumber\\ && - \lambda P_\alpha{}^\mu \nabla_\mu \nabla^\beta - \lambda P_\mu{}^\beta \nabla_\alpha \nabla^\mu + (\Box P_\alpha{}^\beta) + 2 (\nabla_\mu P_\alpha{}^\beta) \nabla^\mu - \lambda (\nabla_\alpha \nabla_\mu P^{\mu\beta}) \vphantom{\frac{1}{2}} -\nonumber\\ && - \lambda (\nabla_\alpha P_\mu{}^\beta) \nabla^\mu - \lambda (\nabla_\mu P^{\mu\beta}) \nabla^\alpha + P_\alpha{}^\mu P_\mu{}^\beta \vphantom{\frac{1}{2}} \end{eqnarray*}
By commuting covariant derivatives, it is easy to obtain the input quantities needed for the algorithm described in Calculating one-loop counterterms: \begin{multline} K^{\mu\nu\gamma\delta}{}_{\alpha}{}^{\beta} = \frac{1}{12} (\lambda^{2}-2 \lambda) (\delta_{\alpha}{}^{\gamma} g^{\mu\delta} g^{\beta\nu}+\delta_{\alpha}{}^{\delta} g^{\nu\gamma} g^{\beta\mu}+\delta_{\alpha}{}^{\gamma} g^{\beta\delta} g^{\mu\nu}+\delta_{\alpha}{}^{\nu} g^{\gamma\delta} g^{\beta\mu} +\\+ \delta_{\alpha}{}^{\nu} g^{\beta\gamma} g^{\mu\delta}+\delta_{\alpha}{}^{\gamma} g^{\beta\mu} g^{\nu\delta}+\delta_{\alpha}{}^{\mu} g^{\beta\delta} g^{\nu\gamma}+\delta_{\alpha}{}^{\mu} g^{\beta\gamma} g^{\nu\delta}+\delta_{\alpha}{}^{\nu} g^{\beta\delta} g^{\mu\gamma} +\\+ \delta_{\alpha}{}^{\delta} g^{\mu\gamma} g^{\beta\nu}+\delta_{\alpha}{}^{\mu} g^{\gamma\delta} g^{\beta\nu}+\delta_{\alpha}{}^{\delta} g^{\mu\nu} g^{\beta\gamma})+\frac{1}{3} \delta_{\alpha}{}^{\beta} (g^{\gamma\delta} g^{\mu\nu}+g^{\mu\gamma} g^{\nu\delta}+g^{\nu\gamma} g^{\mu\delta}), \end{multline}
\[ S^{\mu\nu\rho\alpha\beta} = 0, \]
\begin{multline} W^{\mu\nu}{}_{\alpha}{}^{\beta} = -\frac{1}{2} \lambda P_{\alpha}{}^{\mu} g^{\nu\beta}-\frac{1}{2} \lambda P_{\alpha}{}^{\nu} g^{\mu\beta}-\frac{1}{2} \lambda P^{\beta\mu} d^{\nu}{}_{\alpha}-\frac{1}{2} \lambda P^{\beta\nu} d^{\mu}{}_{\alpha}-\frac{2}{3} R^{\mu\nu} \delta_{\alpha}{}^{\beta} +\\+ 2 P_{\alpha}{}^{\beta} g^{\mu\nu}+\frac{1}{2} (2 \lambda-\lambda^{2}) R_{\alpha}{}^{\beta} g^{\mu\nu}+\frac{1}{6} (-2 \lambda^{2}+\lambda) (R^{\beta\mu} \delta^{\nu}{}_{\alpha}+R_{\alpha}{}^{\nu} g^{\mu\beta} +\\+ R_{\alpha}{}^{\mu} g^{\nu\beta}+R^{\beta\nu} \delta^{\mu}{}_{\alpha})+\frac{1}{6} (2 \lambda-\lambda^{2}) (R_{\alpha}{}^{\mu\beta\nu}+R_{\alpha}{}^{\nu\beta\mu}), \end{multline}
\[ N^{\rho\alpha\beta} = 0, \]
\begin{multline} M_{\alpha}{}^{\beta} = \frac{1}{4} (2 \lambda-\lambda^{2}) R^{\gamma\mu\nu\beta} R_{\alpha\mu\nu\gamma}+\frac{1}{12} (7 \lambda^{2}+4 \lambda) R^{\mu\nu} R_{\mu\alpha\nu}{}^{\beta}+\frac{1}{2} \lambda R^{\mu\beta} P_{\alpha\mu} +\\+ P_{\alpha\mu} P^{\mu\beta}-\frac{1}{2} R_{\mu\nu\gamma\alpha} R^{\mu\nu\gamma\beta}+\frac{1}{6} (-2 \lambda^{2}+\lambda) R^{\mu\beta} R_{\alpha\mu}+\frac{1}{2} \lambda R^{\mu}{}_{\alpha}{}^{\nu\beta} P_{\mu\nu}, \end{multline}
\[ (Kn)^{-1}{}_{\alpha}{}^{\beta} = (2 \gamma+\gamma^{2}) n_{\alpha} n^{\beta}+\delta_{\alpha}{}^{\beta}, \]
\[ F_{\mu\nu\alpha\beta} = R_{\mu\nu\alpha\beta}, \]
where $\gamma = \frac{\lambda}{1-\lambda}$.
The following code produces the result:
//setup symmetries of Riemann tensor
addSymmetries 'R_abcd', -[[0, 1]].p, [[0, 2], [1, 3]].p
setSymmetric 'R_ab', 'P_ab'
def iK = 'iK_a^b=d_a^b+(2*g+g**2)*n_a*n^b'.t
def K = ('K^{lmcd}_a^{b}=' +
'd_a^b*1/3*(g^{lm}*g^{cd} ' +
'+ g^{lc}*g^{md} + g^{ld}*g^{mc})' +
'+ 1/12*(-2*l+l**2)*(' +
'g^{lm}*d_a^c*g^{bd} + g^{lm}*d_a^d*g^{bc} +' +
'g^{lc}*d_a^m*g^{bd} + g^{lc}*d_a^d*g^{bm} +' +
'g^{ld}*d_a^m*g^{bc} + g^{ld}*d_a^c*g^{bm} +' +
'g^{mc}*d_a^l*g^{bd} + g^{mc}*d_a^d*g^{bl} +' +
'g^{md}*d_a^l*g^{bc} + g^{md}*d_a^c*g^{bl} +' +
'g^{cd}*d_a^l*g^{bm} + g^{cd}*d_a^m*g^{bl})').t
def S = 'S^lmpab=0'.t
def W = ('W^{lm}_a^b=' +
'2*P_{a}^{b}*g^{lm}-2/3*R^lm*d_a^b' +
'-l/2*P_a^l*g^mb -l/2*P_a^m*g^lb' +
'-l/2*P^bl*d^m_a -l/2*P^bm*d^l_a' +
'+1/6*(l-2*l**2)*(' +
'R_a^l*g^mb + R_a^m*g^lb +' +
'R^bl*d^m_a + R^bm*d^l_a)' +
'+1/6*(2*l-l**2)*(R_a^lbm+R_a^mbl)' +
'+1/2*(2*l-l**2)*g^lm*R_a^b').t
def N = 'N^pab=0'.t
def M = ('M_a^b = P_al*P^lb' +
'-1/2*R_lmca*R^lmcb + l/2*P_al*R^lb' +
'+l/2*P_lm*R^l_a^mb' +
'+1/6*(l-2*l**2)*R_al*R^lb' +
'+1/12*(4*l+7*l**2)*R_lam^b*R^lm' +
'+1/4*(2*l-l**2)*R_almc*R^clmb').t
def F = 'F_lmab=R_lmab'.t
def lambda = 'l = g/(1+g)'.t
(K, W, M) = [K, W, M].collect { lambda >> it }
def div = oneloopdiv4(iK, K, S, W, N, M, F)
def counterterms = div.counterterms
counterterms = (Factor[[FactorScalars: false]] & 'P_l^l = P'.t) >> counterterms
println counterterms
> counterterms = (1/12)*R*P*(4+g**2+2*g)+(1/12)*(12+g**2+6*g)*P^{m}_{l}*P^{l}_{m}
+(1/120)*(28+5*g**2+20*g)*R**2+(1/6)*g*(4+g)*P_{lm}*R^{lm}
+(1/60)*(-32+5*g**2+10*g)*R^{m}_{l}*R^{l}_{m}+(1/24)*g**2*P**2
Multiplying the produced result by $1\left/16\pi(d-4)\right.$ and integrating over the space-time volume gives: \[ \Gamma^{(1)}_{\infty} = \frac{1}{16\pi(d-4)} \int d^4 x \sqrt{-g} \left( \frac{1}{60}(-32+5 \gamma^2+10 \gamma) R_{\epsilon\mu} R^{\epsilon\mu} +\frac{1}{24}\gamma^2 P^2 +\right.\\ +\frac{1}{120} R^2 (28+5 \gamma^2+20 \gamma)+\frac{1}{12} (\gamma^2+12+6 \gamma) P_{\beta\alpha} P^{\alpha\beta} +\\ \left.+\frac{1}{6}\gamma (4+\gamma) R_{\nu\epsilon} P^{\nu\epsilon} +\frac{1}{12} R (\gamma^2+4+2 \gamma) P\right) \] where $R_{\mu\nu}$ is a Ricci tensor, $R$ is a Riemann scalar curvature, $P = P^{\mu}{}_\mu$ and notation $\lambda = \gamma/(1+\gamma)$ is used.
The result is two times larger then in the case of Vector field, because of the following obvious relation \begin{equation*} \mbox{Tr} \ln D^2 = 2 \mbox{Tr} \ln D. \end{equation*}