  We present details of a lattice QCD calculation of the $B_s\to D_s^*$ axial
form factor at zero recoil using the Highly Improved Staggered Quark (HISQ)
formalism on the second generation MILC gluon ensembles that include up, down,
strange and charm quarks in the sea. Using the HISQ action for all valence
quarks means that the lattice axial vector current that couples to the $W$ can
be renormalized fully non-perturbatively, giving a result free of the
perturbative matching errors that previous lattice QCD calculations have had.
We calculate correlation functions at three values of the lattice spacing, and
multiple `$b$'-quark masses, for physical $c$ and $s$. The functional
dependence on the $b$-quark mass can be determined and compared to Heavy Quark
Effective Theory expectations, and a result for the form factor obtained at the
physical value of the $b$-quark mass. We find $\mathcal{F}^{B_s\to D_s^*}(1) =
h^s_{A_1}(1) = 0.9020(96)_{\text{stat}}(90)_{\text{sys}}$. This is in agreement
with earlier lattice QCD results, which use NRQCD $b$ quarks, with a total
uncertainty reduced by more than a factor of two. We discuss implications of
this result for the $B\to D^*$ axial form factor at zero recoil and for
determinations of $V_{cb}$.
