I already posted this question on MSE.
Using theorem $IV$ from this article, it is possible to prove that when $p$ is a prime $p ≡ 3\bmod4$, $x ≢ y\bmod{p}$ and $\gcd(x,y) = 1$, then the equation $ \dfrac{x^p - y^p}{x - y} = a^2+pb^2 $ ever admits solutions for some integers $\{a, b\}$ (see proof below).
I want to know how $\{a, b\}$ can be expressed as a function of $\{x, y\}$, the case $p=3$ is easy
$ a = \dfrac{x-y}{2}$ , $ b = \dfrac{x+y}{2}$ (of course not all the solutions have this form). With the great help of Will Jagy on MSE, the closed form for the following cases are
$p = 7$
$ a = \frac{ (x-y) (2x^2 + 3xy + 2y^2 )}{ 2}$ , $ b = \frac{xy(x+y)}{2}$
$p = 11$
$ a = \frac{( x-y) (2x^4 + 3x^3 y + x^2 y^2 + 3 x y^3 + 2 y^4)}{ 2}$ , $ b = \frac{ xy(x+y)(x^2 - xy+y^2)}{ 2}$
$p = 19$
$ a = \frac{ (x-y) (2x^8 + 3x^7y - x^6y^2 + 2x^5y^3 + 7x^4y^4 + 2x^3y^5 - x^2y^6 + 3xy^7 + 2y^8 )}{ 2}$, $ b = \frac{ x y (x +y) ( x^2 -xy +y^2)(x^4 - x^2 y^2 + y^4 )}{ 2} $
What can be said in general?
I guess it is possible to restate the problem in terms of cyclotomic binary forms: if we define $ f(k) = y^{\phi(k)} \Phi_k(x/y) $ where $ \phi(k)$ is Euler's phi function and $ \Phi_k(x/y) $ is the $k$-th cyclotomic polynomial, so $ f(p) = \dfrac{x^p - y^p}{x - y} $ and
$ f(7) = a^2 + 7\left (\dfrac{xyf(2)}{2} \right )^2 $
$ f(11) = a^2 + 11\left (\dfrac{xyf(2)f(6)}{2} \right )^2 $
$ f(19) = a^2 + 19\left (\dfrac{xyf(2)f(6)f(12)}{2} \right )^2 $
Proof:
When $(x - y) ≢ {0}\bmod{p}$
$ \dfrac{x^p - y^p}{x - y} = P(p) $
where $P(p)$ is the arithmetic primitive factor of $ (x^p - y^p) $ defined in that article (p. $175$). With other considerations we can say that each one of those primitive prime factors is of the form $(2k_i p+1)$ where $k_i$ $∈ N$.
Now $q\mid{a^2+pb^2} $, with $q\nmid p$, iff $ \left(\frac{-p}{q}\right) = 1$ and equivalently $q ≡ \{q_1, q_2, q_3,..\}\bmod{4p}$ (see corollary $1.19$ in David Cox's book Primes of the Form $x^2+ny^2$), so being $q=2k_i p+1$, we have just $q_1=1$ and $q_2 = 2p+1$. Using that $p ≡ 3\bmod4$, we have $ \left(\frac{-p}{2p+1}\right) = 1$ and so $ \dfrac{x^p - y^p}{x - y} = a^2+pb^2 $ ever admits solutions for some integers $\{a, b\}$.