  Given a cusp form $f$ which is supersingular at a fixed prime $p$ away from
the level, and a Coleman family $F$ through one of its $p$-stabilisations, we
construct a $2$-variable meromorphic $p$-adic $L$-function for the symmetric
square of $F$, denoted $L^{\mathrm{imp}}_p(\mathrm{Sym}^2 F)$. We prove that
this new $p$-adic $L$-function interpolates values of complex imprimitive
symmetric square $L$-functions, for the various specialisations of the family
$F$. It is in fact uniquely determined by its interpolation properties. We also
prove that the function $L^{\mathrm{imp}}_p(\mathrm{Sym}^2 F)$ satisfies a
functional equation. We use this $p$-adic $L$-function to prove a $p$-adic
factorisation formula, expressing the geometric $p$-adic $L$-function attached
to the self-convolution of $F$, as the product of
$L^{\mathrm{imp}}_p(\mathrm{Sym}^2 F)$ and a Kubota-Leopoldt $L$-function. This
extends a result of Dasgupta in the ordinary case.
  Using Beilinson-Flach classes constructed by Kings, Zerbes and the second
author we construct motivic cohomology classes $b_f$, and prove that, under
some hypotheses, they differ by a scalar factor from the higher cyclotomic
classes constructed by Beilinson. Using this relation, we prove the
interpolation formulae for $L^{\mathrm{imp}}_p(\mathrm{Sym}^2 F)$ and the
factorisation formula.
