We present details of a lattice QCD calculation of the $B_sto 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_sto 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 $Bto D^*$ axial form factor at zero recoil and for determinations of $V_{cb}$.