The Poisson-kernel in the plane and polynomials Let
\begin{align*}
 p(z) & = \ \prod_{j = 1}^{n} (z - z_j) \, ; \ |z_j| \ = \ 1  \\
      & = \ \prod_{j = 1}^{l+1} (z - z_j)^{M_j}
\end{align*}
be a non-constant complex polynomial with $l+1$ distinct zeros of multiplicity $M_j$, which all lie on the unit circle. Then there are $l$ not necessarily distinct zeros $z'_k$ of the derivative $p'(z)$ in the interior of the unit disc:
\begin{align*}
 p'(z) & = \ n \, \prod_{k = 1}^{n-1} (z - z'_k) \, ; \qquad |z'_k| \ < \ 1 \ : \ 1 \ \leq \ k \ \leq \ l  \, ,\\
\end{align*}
where $l \, = \, n-1$, if $p(z)$ has only simple zeros. The $z'_k$ in the interior of the unit disc are known as critical points of the second kind, as p(z) itself is not zero there.

The Poisson-kernel $k(a,b)$ of two complex numbers, where $|a| \ = \ 1$ and $|b| \ < \ 1$,
$$ 
k(a,b) \ = \ \frac{1 - |b|^2}{|a - b|^2} \, ,
$$
is well-known from the solution of the Dirichlet-Problem of the Laplace equation for the unit disc.

Question, equation (1):
Is the equation
\begin{equation}
 \sum^l_{k = 1} \frac{1 - |z'_k|^2}{|z_j - z'_k|^2} \ = \ \frac{n - M_j}{M_j}
\end{equation}
true for any zero $z_j$ of any polynomial $p(z)$ with at least two distinct zeros, all of whose zeros lie on the unit circle?

Some background:

This question is related to the matrix
\begin{equation}
\mathbf{D_{pu}}  \ := \ 
\begin{pmatrix} 
\frac{M_1}{n}           &        \ldots        &        \frac{M_{l+1}}{n}                 \\
\frac{M_1}{n}  \,   \frac{H(z'_1)}{|z_1 - z'_1|^2}                             &       \ldots         &         
\frac{M_{l+1}}{n}  \,   \frac{H(z'_1)}{|z_{l+1} - z'_1|^2}         \\
\vdots & & \vdots \\
\frac{M_1}{n}  \,   \frac{H(z'_l)}{|z_1 - z'_l|^2}         &       \ldots         &         
\frac{M_{l+1}}{n}  \,   \frac{H(z'_l)}{|z_{l+1} - z'_l|^2}         \\
\end{pmatrix} \, ,
\end{equation}
where $H(z'_k)$ is the harmonic mean  of the squares of the distances of the $z_j$ to a critical point of the second kind $z'_k$, defined as follows:
$$
H(z'_k) \ := \ \frac{n}{\sum^{l+1}_{j=1} \frac{M_j}{|z_j - z'_k|^2}} \, .
$$
$D_{pu}$ appears, with an immaterial difference in the definition, in the paper
"Borcea’s variance conjectures on the critical points of polynomials"
by D. Khavinson, R.Pereira, M. Putinar, E. Saff and S. Shimorin in Notions of Positivity and the Geometry of Polynomials, Series: Trends in Mathematics,
Birkhäuser, Basel, 2011, pp. 283 – 309, see https://arxiv.org/abs/1010.5167.
It is called "Augmented Gauss-Lucas matrix" there, as its rows 2 to l+1 give the coefficients in a convex combination of the zeros expressing a critical point of the second kind:
$D_{pu}(z_1, \ldots, z_{l+1})^t \ = \ (S, z'_1, \ldots, z'_{l})^t$ ,
where $S$ is the centroid of the zeros. These convex combinations are used in an early and well-known proof of the Gauss-Lucas-Theorem, see
"Analytic Theory of Polynomials" by Q. I. Rahman and G. Schmeisser, Oxford University Press, 2002, section 2.6, notes to section 2.1,
hence the name of the matrix given in the paper above. The authors of the paper anticipated me in proving, that $D_{pu}$ is orthostochastic if all zeros lie on a line and unistochastic if there are at most three different zeros. I arrived at these results independently with a slightly different approach.
The connection to the Poisson-kernel for unimodular zeros arises by the equation:
$$ (2) \qquad H(z'_k) \ = \ 1 - |z'_k|^2 \, ,$$
which is valid, if all zeros are unimodular, see equation $(\ast)$ in my new question Question about an exact expression for the root-mean-square of the distances of the critical points to a given zero of a polynomial. I edited this one, because it linked to an old MathStackexchange question, which was deleted, because it drew no attention.
Equation (2) implies that the stochastic matrix $D_{pu}$ is doubly stochastic for any polynomial $p(z)$ with unimodular zeros if and only if equation (1) is true. As $D_{pu}$ is stochastic, equation (1) is true "in the mean". Double stochasticity would prove a (slight) refinement of Lemma 7.3.3 in Rahman/Schmeisser, showing that there is a critical point of the second kind in a certain circle associated with a zero, implying the correctness of Sendov's conjecture for the polynomial $p(z)$, Sendov's conjecture. According to Profs. Khavinson and Pereira, this connection was not recognized before.

All my (rather few) investigations of examples of low degree showed, that $D_{pu}$ is doubly stochastic if and only if the zeros of $p(z)$ lie on a line or a circle, $p(z) = z^4 - z$ is an instructive counterexample, so a generalization of my problem above is the following

Further Conjecture/Question (*):
$D_{pu}$ is doubly stochastic if and only if all zeros of $p(z)$ lie on a line or a circle.
If this wrong: When is $D_{pu}$ doubly stochastic?

The paper by Khavinson et al. also gives an example of a $D_{pu}$ which is not doubly stochastic and where the zeros do not lie on a line or circle. The authors mention questions about the column-sums of $D_{pu}$, which are of the same nature as my own, without the connection to the Poisson-kernel and conjecture (*), though. According to Profs. Khavinson and Pereira, the problems have not been resolved up to now.

The following question is now discussed more throroughly in my new question cited above, I leave it here, because I refered to it there:
Is the inequality
\begin{equation}
 \sum^l_{k = 1} \frac{H(z'_k)}{|z_j - z'_k|^2} \ \leq \ \frac{n}{M_j}
\end{equation}
true for any complex polynomial $p(z)$ with at least two distinct zeros, without restrictions on its zeros?
 A: I would suggest the following strategy to attack your question. We fix some zero of $p$ say $z_1$ and we define following functions: $C(z):=p'(z)/p(z)$; 
$$u(z):=\log |C(z)|^2-\log|C(\frac 1{\bar z})|^2+\log|z|^2; \quad v(z):=\frac{1-|z|^2}{|z-z_1|^2}. $$
The function $u$ has the following properties: its logarithmic singularities in the unit disk are 
points  $z'_k$; it vanishes at the unit circle except points $z_j$; its Laplacian  satisfies 
$$\Delta u(z)= \sum_k \delta_{z_k'}(z)$$
inside the unit disk in the sense of distributions (subject to right normalization of the Laplacian). Now, we apply Green's formula 
$$\iint_D \left(\Delta u\cdot v- u\cdot \Delta v\right)\, dA = 
\int_{\partial D}\left(\frac{\partial u}{\partial n}\cdot v - u\cdot \frac {\partial v}{\partial n}\right) \, ds $$
where the domain $D$ is the unit disk with deleted small disks around points $z_j$. The left hand side of this Green's formula gives  the left hand side of your equation (1). The right hand side vanishes along arcs of the unit circle and it remains to analyse its asymptotic behaviour at small arcs around points $z_j$. The non-zero limit will be only for the integral along the arc around $z_1$ which is subject to a piece of hard analysis. 
A: Here is an easy algebraic proof of your identity. I present it in the case of a simple zero, a general case requires some obvious modifications. 
So, let $z_1$ be a simple zero of $p$. Your equation (1)  is a special case of more general equation: 
$$\sum_{l=1}^{n-1}\frac{|z_1|^2-|z_l'|^2}{|z_1-z_l'|^2}= n-1 + 2\sum_{k=2}^n\frac{|z_1|^2-|z_k|^2}{|z_1-z_k|^2}. \tag 3$$
To prove it, we write first 
$$\frac{p''(z_1)}{p'(z_1)}=\sum_{l=1}^{n-1}\frac 1{z_1-z_l'}$$
which is true by partial fraction decomposition of $\frac{p''(z)}{p'(z)}$. 
On the other hand, 
$$\frac{p''(z_1)}{p'(z_1)}=2\lim_{z\to z_1}\left(\frac{p'(z)}{p(z)}-\frac 1{z-z_1}\right) = 2\sum_{k=2}^n\frac 1{z_1-z_k}$$
by partial fraction decomposition of $\frac{p'(z)}{p(z)}$. 
This gives the identity 
$$\sum_{l=1}^{n-1}\frac 1{z_1-z_l'}=2\sum_{k=2}^n\frac 1{z_1-z_k}.$$
It remains to multiply it by $z_1$ and to take the real part using identity $2\text{Re}\left(\frac {z_1}{z_1-a}\right) =1+\frac{|z_1|^2-|a|^2}{|z_1-a|^2}$. 
A: This answer fills in some details not explicitly given in the accepted answer by user Karl Marx and treats arbitrary $M_j$ for my reference. Any upvotes should go to the accepted answer first, as the relevant idea is there already.
For ease of notation it is convenient to consider $z_1$ of multiplicity $M_1$ without loss of generality. The formular to be derived applies for arbitrary polynomials with at least two distinct zeros. We have
$$ p(z) \ = \ (z - z_1)^{M_1} q(z) $$ and
$$ p'(z) \ = \ (z - z_1)^{M_1-1} r(z) \, , $$
where $q(z)$ and $r(z)$ are polynomials of degree $n-M_1$ and
$$ \frac{q'(z)}{q(z)} \ = \ \sum_{j=M_1 + 1}^{n} \frac{1}{z - z_j} $$ and 
$$ \frac{r'(z)}{r(z)} \ = \ \sum_{k=1}^{n-M_1} \frac{1}{z - z'_k} \, . $$
The idea is to compute 
$$ \frac{p^{(M_1+1)}(z_1)}{p^{(M_1)}(z_1)} $$
in two different ways by using the given formulars for $p(z)$ and $p'(z)$. We get
\begin{align} p^{(M_1)}(z) & = \ M_1! \, q(z) + M_1 \, M_1! \, (z-z_1) \, q'(z) + O((z-z_1)^2) \\
& = \ (M_1-1)! \,  r(z) + (M_1 - 1) \, (M_1-1)! \,  (z-z_1) \, r'(z) + O((z-z_1)^2)
\end{align}
and
\begin{align} p^{(M_1+1)}(z) & = \ (M_1 + 1)! \, q'(z) + O((z-z_1)) \\
& = \ M_1! \,  r'(z) + O((z-z_1)) \, .
\end{align}
Inserting $z_1$ into those expressions and dividing, we obtain
$$ \left( \frac{p^{(M_1+1)}(z_1)}{p^{(M_1)}(z_1)} \ = \ \right)  M_1 \sum_{k=1}^{n-M_1} \frac{1}{z_1 - z'_k}  \ = \ (M_1 + 1) \sum_{j=M_1 + 1}^{n} \frac{1}{z_1 - z_j} \, . $$
Then we proceed as in the accepted answer by multiplying $z_1$ and looking at the real part to obtain
$$  \sum_{k=1}^{n-M_1} \frac{|z_1|^2 - |z'_k|^2}{|z_1 - z'_k|^2} \ = \ \frac{n-M_1}{M_1} + \frac{M_1 + 1}{M_1} \sum_{j=M_1 + 1}^{n} \frac{|z_1|^2 - |z_j|^2}{|z_1 - z_j|^2} \, , $$ 
which is valid for arbitrary polynomials of degree $n$ with at least two distinct zeros and we are done for the case of solely unimodular zero.
