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.