What is the total square on the dual Steenrod algebra?
We bothered to write it down in our paper. Look at pg 6, we give some of references that we know of.
I did not find a formula for the left action of the $Sq$ on $\zeta_i$s in the literature. But from the formula for left action of $Sq$ on $\xi_i$ and formulas relating $\xi_i$s and $\zeta_i$s one can do an extensive combinatorial argument to see that $$ Sq(\zeta_i) = \zeta_i + \zeta_{i-1} + \dots + \zeta_1 + 1$$.
(In my experience the combinatorial inductive argument was tedious but straightforward!)
[ For example, let's consider the first nontrivial case, ie calculate $Sq(\zeta_2)$. Keep in mind that $\zeta_2 = \xi_2 + \xi_1^3$ and $\zeta_1 = \xi_1$. Then $$Sq(\zeta_2) = Sq(\xi_2 + \xi_1^3) = (\xi_2 + \xi_1^2) + (\xi_1 +1)^3 = \zeta_2 + \zeta_1 + 1.$$ Keep going inductively to get the formulas for $Sq(\zeta_i)$... ]
I don't have a reference for you, but here is a comment on how to prove these formulas using the Kronecker pairing that you alluded to.
The Steenrod operation $Sq^m$ is dual to the element $\xi_1^m$ in the monomial basis of the dual Steenrod algebra; the left and right actions of the Steenrod algebra on $\mathcal{A}_*$ are composites of the coproduct in the dual Steenrod algebra and the action on the right or left side. If the coproduct satisfies $\Delta x = \sum x' \otimes x''$, we then get $$ \begin{align*} x \cdot Sq^m &= \sum (\xi_1^m)^*(x') x'',\\ Sq^m \cdot x &= \sum x' (\xi_1^m)^* (x''). \end{align*} $$ (The apparent order reversal is necessary to make this into a left/right action.) We'd like to apply this to the comultiplication formulas $\Delta \xi_n = \sum_{i+j=n} \xi_i^{2^j} \xi_j$ and $\Delta \zeta_n = \sum_{i+j=n} \zeta_i \zeta_j^{2^i}$. Here by convention $\xi_0 = \zeta_0 = 1$.
To apply this to the $\xi_n$, we first remark that $$ \sum_m (\xi_1^m)^*(\xi_i^{2^j}) = \begin{cases}1 &\text{if }i=0,1,\\0&\text{otherwise.}\end{cases} $$ Therefore: $$ \begin{align*} \xi_n \cdot Sq &= \sum (\xi_1^m)^* (\xi_i^{2^j}) \xi_j = \xi_n + \xi_{n-1}\\ Sq \cdot \xi_n &= \sum \xi_i^{2^j} (\xi_1^m)^* (\xi_j) = \xi_n + \xi_{n-1}^2. \end{align*} $$ To figure out the corresponding result for the $\zeta_n$, we have to figure out what the coefficient of $\xi_1^{2^n-1}$ is in the formula for $\zeta_n$. The $\zeta_i$ are defined inductively, for $n > 0$, using the formula $$ \sum_{i+j=n} \xi_i^{2^j} \zeta_j = 0. $$ If we take the quotient by the ideal generated by $\xi_2, \xi_3, \dots$ we find that this formula reduces to $$ \zeta_n + \xi_1^{2^{n-1}} \zeta_{n-1} \equiv 0 $$ and so inductively $\zeta_n \equiv \xi_1^{2^n - 1}$ mod the higher $\xi_i$. This means $$ \sum_m (\xi_1^m)^*(\zeta_j^{2^i}) = 1 $$ for any $i$ and $j$.
Therefore: $$ \begin{align*} \zeta_n \cdot Sq &= \sum (\xi_1^m)^* (\zeta_i) \zeta_j^{2^i} = \zeta_n + \zeta_{n-1}^2 + \dots + \zeta_1^{2^{n-1}} + 1\\ Sq \cdot \zeta_n &= \sum \zeta_i (\xi_1^m)^* (\zeta_j^{2^i}) = \zeta_n + \zeta_{n-1} + \dots + \zeta_1 + 1. \end{align*} $$