This identity might be explicitly given in one of the following articles of Arthur Mattuck (Mattuck does prove many identities), but I could not find it.
MR0142553 (26 #122)
Mattuck, Arthur
Symmetric products and Jacobians.
Amer. J. Math. 83 1961 189–206.
14.20 (14.51)
MR0136608 (25 #76)
Mattuck, Arthur
On symmetric products of curves.
Proc. Amer. Math. Soc. 13 1962 82–87.
14.10 (14.20)
Just to repeat your notation, for every $d\in \mathbb{Z}$, define $u_d:C_d\to \text{Pic}^d_{C/k}$ to be the Abel map. Define $i_{d-1,d}:C_{d-1}\to C_d$ to be the morphism adding $\underline{p}_0$ to an effective divisor. Define $j_{d-1,d}:\text{Pic}^d_{C/k}$ to be the morphism that twists an invertible sheaf by $\mathcal{O}_C(\underline{p_0})$. The morphism $j_{d-1,d}$ is an isomorphism.
The image $\Theta = u_{g-1}(C_{g-1})\subset \text{Pic}^{g-1}_{C/k}$ is an effective, reduced Cartier divisor. For every $d$, define $\Theta_d\subset \text{Pic}^d_{C/k}$ to be the unique Cartier divisor such that $j_{d-1,d}^*\Theta_d$ equals $\Theta_{d-1}$. Define $\Delta_d \subset C_d$ to be the (reduced) Cartier divisor parameterizing effective divisors on $C$ that are nonreduced (i.e., at least one pair of points comes together). Finally, denote $H_d = i_{d-1}(C_{d-1})$ as an effective Cartier divisor on $C_d$.
The basic identity is that $i_{d+1,d}^*(\underline{\Delta_{d+1}})$ equals $\underline{\Delta_d} + 2\underline{H_d}$ as effective Cartier divisors on $C_d$. Set-theoretically this is clear. The fact that the coefficient of $\underline{H_d}$ equals $2$ follows from a local analysis: it is the fact that the quotient morphism $q:C^2\to C_2$ is simply branched along the diagonal (i.e., reduce to the case that $d$ equals $1$). By definition, $i_{d+1,d}^*u_{d+1}^* \Theta_{d+1}$ equals $u_d^*\Theta_d$.
Finally, the normal sheaf of the embedding $i_{d,d+1}$ equals $T_{C,p_0}\otimes_k \mathcal{O}_{C_d}(\underline{H}_d)$. To see this, consider the one-parameter family of deformations of $i_{d,d+1}$ obtained by varying $p_0\in C$. In particular, a first-order deformation in $T_{C,p_0}\setminus\{0\}$ defines a section of the normal sheaf of $i_{d,d+1}$. Away from $H_d$, it is straightforward to see that this section is nonzero. Since $H_d$ is irreducible, it follows that the normal sheaf equals $T_{C,p_0}\otimes_k \mathcal{O}_{C_d}(m\underline{H}_d)$ for some integer $m$. Working locally at a generic point of $H_d$, $m$ equals $1$ since $C^2 \to C_2$ is simply branched along the diagonal (the same reason for the coefficient $2$ in the previous paragraph).
That is all the necessary preparation. Because of the computations of $i_{d+1,d}^*$ on various divisors, on divisor classes we have the relation $$ i_{d+1,d}^* \mathcal{O}_{C_{d+1}}(-\underline{\Delta}_{d+1}-2\underline{\Theta}_{d+1} +2((d+1)+g-1)\underline{H}_{d+1})) \cong $$ $$ T_{C,p_0}^{\otimes 2(d+g)}\otimes_k \mathcal{O}_{C_d}(-\underline{\Delta}_d -2 \underline{H}_d - 2\underline{\Theta}_d + 2(d+g)\underline{H}_d). $$ Thus the claim for $d+1$ implies the claim for $d$.
Finally, when $d > 2g-2$, by Abel's Theorem / Jacobi Inversion and Riemann-Roch, the morphism $u_d$ is a smooth, projective morphism whose geometric generic fibers are projective spaces of relative dimension $d-g$. Moreover, $H_d$ gives a relative hyperplane class. Thus, for every invertible sheaf $\mathcal{L}$ on $C_d$ that has degree $0$ on a fiber, $(u_d)_*\mathcal{L}$ is an invertible sheaf on $\text{Pic}^d_{C/k}$, and the adjunction map $$u_d^*(u_d)_*\mathcal{L} \to \mathcal{L},$$ is an isomorphism. Now let $f:C\to \mathbb{P}^1$ be a tamely ramified, finite, flat morphism of degree $d> 2g-2$ (such morphisms exist by Riemann-Roch and Bertini type arguments). Let $\widetilde{f}:\mathbb{P}^1\to C_d$ be the morphism associated to the graph of $f$ in $C\times \mathbb{P}^1$ considered as a Cartier divisor of relative degree $d$ over $\mathbb{P}^1$. Then $\widetilde{f}^*H_d$ is the reduced divisor $f(p_0)$. By Riemann-Hurwitz, $\widetilde{f}^*\Delta_d$ is the branch divisor and has degree $2(d+g-1)$. Thus, the invertible sheaf from the previous paragraph has degree $0$ on fibers of $u_d$.
Thus, for every $C$, for every $p_0$, for every $d\geq 0$, there exists a unique invertible sheaf $\mathcal{L}_d$ on $\text{Pic}^d_{C/k}$ such that $j_{d,d+1}^*\mathcal{L}_{d+1}$ is isomorphic to $\mathcal{L}_d$, such that
$$u_d^*\mathcal{L}_d \cong \mathcal{O}_{C_d}(-\underline{\Delta}_d - 2\underline{\Theta}_d + 2(d+g-1)\underline{H}_d),
$$ for all $d$, and such that $$\mathcal{L}_d = (u_d)_* \mathcal{O}_{C_d}(-\underline{\Delta}_d - 2\underline{\Theta}_d + 2(d+g-1)\underline{H}_d),$$
for all $d>2g-2$. The claim is that $\mathcal{L}_d$ is an invertible sheaf on $\text{Pic}^d_{C/k}$ that is algebraically equivalent to zero.
By considering the case that $d=1$, it is already clear that $\mathcal{L}_d$ has degree $0$ on the Abel image of $C$. Moreover, since the construction works relatively in families, via properness and separatedness of the (components) of the relative Picard scheme, it suffices to prove for $C$ very general that every invertible sheaf on $\text{Pic}^1_{C/k}$ that has degree $0$ on the Abel image of $C$ is algebraically equivalent to zero. This follows, for instance, from the Franchetta conjecture (proved over the complex numbers by John Harer, and extended to positive characteristic by Stefan Scröer).
Edit. You do not need to use Franchetta's conjecture. If every correspondence on $C$ has valence, then you can use $d=2$ to prove that $u_2^*\mathcal{L}_2$ is algebraically equivalent to zero. By the theorem of the cube, already $u_2^*$ is injective on Picard groups. So it suffices to prove that there exists a curve $C$ such that every correspondence has valence. That is much easier than Franchetta's conjecture.
Second Edit. If you specialize to the case $d=g-1$, which again is sufficient to prove the claim when $g\geq 3$ so that $d\geq 2$, you can eliminate the dependence on $p_0$. Once you do that, the claim reduces to the following claim: for the universal divisor $\mathcal{D}\subset C\times C_{g-1}$, the following invertible sheaf on $C_d$ is algebraically equivalent to zero, $$\mathcal{B}_C := \text{det}(R^1\text{pr}_{1,*}\mathcal{O}(-\underline{\mathcal{D}}))^{\otimes 2}(-\underline{\Delta}).$$ Since this invertible sheaf is defined over the field of definition of $C$, in fact, it must be that $\mathcal{B}_C$ is isomorphic to the structure sheaf. There may be a direct proof that $\mathcal{B}_C$ is trivial using Grothendieck-Riemann-Roch (which would prove your result without specialization).
Final Edit. If you apply Grothendieck-Riemann-Roch to compute the class of $\text{det}(R^1\text{pr}_{1,*}\mathcal{O}(-\underline{\mathcal{D}}))$, it gives that the Cartier divisor class of the square of this line bundle equals the pushforward via the finite flat morphism, $$\text{pr}_{1,\mathcal{D}}:\mathcal{D} \to C_{d-1},$$ of the relative canonical divisor class. Via the pullback sequence of sheaves of differentials associated to $\text{pr}_{1,\mathcal{D}}$, there is an effective representative of the relative divisor class supported on the ramification locus of $\text{pr}_{1,\mathcal{D}}$. That maps to the diagonal $\Delta_{g-1}$ in $C_{g-1}$. By the same local analysis as above, the multiplicity equals $1$. This proves that $\mathcal{B}_C$ is isomorphic to the structure sheaf without using specialization. That proves that $u_{g-1}^*\mathcal{L}_{g-1}$ is algebraically equivalent to zero. When $g\geq 3$, that implies your result without specialization.
Final, Final Edit. Francesco Polizzi points out that the formula does indeed follow from Formulas (1) and (5) of Mattuck's article "On symmetric products of curves" listed above.