This is just an extended comment, giving reformulation of the problem and reducing it to just $p-1$ unknowns and $p-1$ quadratic equations over the Gaussian integers.

Consider the generating polynomials:
\begin{split}
A(x) &:= \sum_{i=0}^{p-1} \alpha_i x^i, \\
B(x) &:= \sum_{i=0}^{p-1} \beta_i x^i.
\end{split}

The linear equations $\sum_j \alpha_j = \sum_j \beta_j = 0$ are equivalent to $A(1)=B(1)=0$, i.e., both $A(x)=(x-1)\bar A(x)$ and $B(x)=(x-1)\bar B(x)$ are multiples of $x-1$.

Viewing indices modulo $p$ is equivalent to viewing the polynomials modulo $x^p - 1 = (x-1)\Phi_p(x)$, where $\Phi_p(x) := 1 + x + \dots + x^{p-1}$ is $p$-th cyclotomic polynomial.

For reciprocal polynomials (of fixed degree $p-1$) we have $A^\star(x):=x^{p-1}A(x^{-1})\equiv x^{p-1}A(x^{p-1})\pmod{x^p-1}$ and $B^\star(x):=x^{p-1}B(x^{-1})\equiv x^{p-1}B(x^{p-1})\pmod{x^p-1}$. Then the quadratic equations (under the condition $A(1)=B(1)=0$) translate into
$$\begin{cases}
A(x)B^\star(x) + A^\star(x)B(x) \equiv 0 \pmod{x^p-1},\\
-A(x)A^\star(x) + B(x)B^\star(x) \equiv -4x^{p-1} + 2x + 2x^{p-3} \equiv 2(x^2-1)^2x^{p-3} \pmod{x^p-1}
\end{cases}
$$
Dividing both congruences by $(x-1)x(\frac1x-1)=-(x-1)^2$, we get
$$\begin{cases}
\bar A(x)\bar B^\star(x) + \bar A^\star(x)\bar B(x) \equiv 0 \pmod{\Phi_p(x)},\\
-\bar A(x)\bar A^\star(x) + \bar B(x)\bar B^\star(x) \equiv -2(x+1)^2x^{p-3} \pmod{\Phi_p(x)}.
\end{cases}
$$

In terms of polynomials over Gaussian integers, we have
$$F(x)F^\star(x) \equiv -2(x+1)^2x^{p-3}\pmod{\Phi_p(x)},$$
where
$$F(x) := \bar B(x) + I\cdot \bar A(x)$$
is a polynomial of degree $p-1$ over Gaussian integers. (Notice that the polynomials at both sides of this congruence are [palindromic](https://en.wikipedia.org/wiki/Reciprocal_polynomial#Palindromic_and_antipalindromic_polynomials) if viewed as polynomials of the fixed degree $2p-4$.) 

The last congruence can be viewed as a system of $p-1$ quadratic equations over the the coefficients of $F(x)$.