Something interesting about the quintic $x^5 + x^4 - 4 x^3 - 3 x^2 + 3 x + 1=0$ and its cousins

I think the depressed quintic in the question is what Emma Lehmer called the reduced quintic in her paper, The quintic character of 2 and 3, Duke Math. J. Volume 18, Number 1 (1951), 11-18, MR0040338. In the proof of Theorem 4, she writes that the reduced quintic is $$F(z)=z^5-10pz^3-5pxz^2-5p\left({x^2-125w^2\over4}-p\right)z+p^2x-{x^3+625(u^2-v^2)w\over8}$$ where $16p=x^2+50u^2+50v^2+125w^2$, $xw=v^2-u^2-4uv$ determines $x$ uniquely. You might plug these into your conjectured diophantine relations, to see whether they hold.

EDIT:

A slightly different formula is given by Berndt and Evans, The determination of Gauss sums, Bull Amer Math Soc 5, no. 2 (1981) 107-129, on page 119 (formula 5.2). It goes $$ F(z)=z^5-10pz^3-5pxz^2-5p\left({x^2-125w^2\over4}-p\right)z+p^2x-p{x^3+625(u^2-v^2)w\over8} $$ the difference being the factor of $p$ in the last term. Berndt and Evans note that Lehmer inadvertently omitted this factor.

Update: (by OP)

Using the edited coefficients of $F(z)$ by Berndt and Evans, if we express it as, $$z^5+az^3+bz^2+cz+d=0$$ and taking into account $16p=x^2+50u^2+50v^2+125w^2,\,x=(v^2-u^2-4uv)/w$, then $$a^3 + 10 b^2 - 20a c =2\times125^2p^2w^2=2z_1^2$$ $$5(a^2 - 4c)^2 + 32a b^2 = 500^2 p^2(u^2 - u v - v^2)^2=z_2^2$$ $$(a^3 + 10 b^2 - 20 a c)\,\big(5 (a^2 - 4 c)^2 + 32 a b^2\big) = 2(a^2 b + 20 b c - 100 a d)^2$$ which confirms all three Diophantine relations by the OP.

Question 1: Yes; in fact $$ \sum_{n\bmod 5} y_n y_{n+1} = \! \sum_{n\bmod 5} y_n y_{n+2} = -p/5 $$ for the "depressed" $y_n$, using an order consistent with the action of the cyclic Galois group. Likewise "for other values of 5".

This follows from the formula $|\tau|^2 = p$ for the absolute value of a nontrivial Gauss sum $\tau$, together with the observation that the quintic Gauss sums are the values at $m \neq 0$ of the discrete Fourier transform $\hat y$ of $\vec y$: $$ \tau_m = \sum_{n \bmod 5} e^{2\pi i mn/5} y_n $$ (and $\tau_0 = 0$ thanks to the "depression"). Now the $y_n$ are real, so $\bar\tau$ is the Fourier transform of the map $n \mapsto y_{-n}$; thus the autocorrelation function $$ c_j = \sum_{n \bmod 5} y_n y_{n+j} $$ is the convolution of $\vec y$ with $n \mapsto y_{-n}$, whence the Fourier transform $\hat c$ is $\tau \bar \tau = |\tau|^2$. But that's $0$ at $m=0$, and $p$ otherwise; therefore $c_0 = 4p/5$ and $c_j = -p/5$ for other $j$. In particular $c_1 = c_2$, QED.

P.S. I was using $x = y - \frac15$, not your $x = \frac{y-1}{5}$, so most formulas above will have to be scaled by suitable powers of $5$ to apply to your $y$'s.