I can recommend "Nonlinear Potential Theory of Degenerate Elliptic Equations" by Heinonen, Kilpeläinen, Martio; and "Function Spaces and Potential Theory" by Adams, Hedberg

The desired estimate cannot hold on arbitrary domains, you need some boundary regularity. However, your answer does not need boundary regularity and, consequently, it should be invalid.

My reason is already in the second comment. You just have to compute the derivative, which is a standard computation.

No, you fix $E_i$ and only vary $x_i$. You can think of the entire process as a two-stage problem: In the first stage, you fix the regions $E_i$ and in the second stage you just optimize over the $x_i$. In the second stage, the $E_i$ are fixed and the optimization w.r.t. $x_i$ can be solved in closed form. The remaining optimization over $E_i$ is delicate.

In your objective, you have the contribution $\int_{E_i} \| x_i^* - x\|^2 \, \mathrm d \ell(x)$. The derivative of this term w.r.t. $x_i^*$ is zero if and only if $x_i^*$ is the center of mass.

I think that the answer depends crucially on the definition of "subdifferential". If we take the standard definition $\partial W(x) = \{ u \mid \forall y \in K : W(y) \ge W(x) + \langle u, y - x \rangle\}$, then the answer is "no". Take $K = [0, \infty) \subset \mathbb R$ and $W(x) = \sqrt{x}$. Then, the subdifferential is empty everywhere (up to the origin) and, consequently, the subdifferential mapping is monotone.

No, but you can replace $\phi$ by $\phi^{cc}$.

Once you have the existence of a solution $\phi$ of (1), you can take $\phi^{cc}$. The objective value of $\phi^{cc}$ is at least that of $\phi$, hence, it is again a solution and further, the $L$-Lipschitz continuity of $c$ implies that $c$-conjugates are $L$-Lipschitz.