The following can explain the shape of the expression and also why you always get pairs (why the minimal polynomials for the two elements in a pair differ by few coefficients, I don't know). Let $p$ be a prime congruent to $1$ mod $6$, and consider the rational function $f(X):=g(X^p)$, where $g(X):=\frac{X^3-3X+1}{X-X^2}$ (note that setting $t:=g(X)$ and clearing denominators gives exactly your generic cubic, except that I like $t$ better than $n$ as a parameter). Then $X\mapsto f(X)$ yields a cover $\mathbb{P}^1\to \mathbb{P}^1$, and we can try to compute its monodromy group (i.e., the Galois group of $f(X)-t=0$ over $\mathbb{C}(t)$ - and, as a matter of fact, still over $\mathbb{Q}(\zeta_p)(t)$). Since $Gal(g(X)-t/\mathbb{C}(t)) = C_3$ and $Gal(X^p-t/\mathbb{C}(t))=C_p$, this monodromy group will be a subgroup of $C_p\wr C_3 (=(C_p\times C_p\times C_p)\rtimes C_3) \le S_{3p}$. We can be more explicit: IF $a,b,c$ denote the three roots of $g(X)-t=0$, and $\alpha:=a^{1/p}$, $\beta:=b^{1/p}$, $\gamma:=c^{1/p}$, then the roots of $f(X)-t=0$ fall into three blocks $\Delta_1, \Delta_2, \Delta_3$, with $\Delta_1:=\{\alpha, \zeta\alpha, \dots, \zeta^{p-1}\alpha\}$ (and analogously for $\Delta_2,\Delta_3$). The inertia groups at the branch points of the cover $X\mapsto f(X)$ will now give more precise information about the Galois group. It's easy to see that $X\mapsto g(X)$ has exactly two (non-rational) branch points (namely the roots of your discriminant $D$), whereas $X\mapsto X^p$ has of course exactly the branch points $0$ and $\infty$, and these two get mapped to the same point $\infty$ under $X\mapsto g(X)$. So in total, the composition has exactly three branch points: the two irrational ones (whose inertia group generators are of order $3$ (permuting $a$, $b$, $c$ cyclically and thus permuting the three blocks $\Delta_i$ in the same way), and $t=\infty$ at which the inertia group is a double-$p$-cycle inside the block kernel (since unramified under $X\mapsto g(X)$, but with two totally ramified preimages on the "upper" level). This permutation group can then be identified as $G:=K\rtimes C_3$, where $K:=\{(x_1,x_2,x_3)\mid x_1,x_2,x_3\in \mathbb{F}_p; \sum x_i=0\}$ is the ``augmentation ideal". This group $G$ has two normal subgroups $N_1$ and $N_2$ isomorphic to $C_p$, both with quotient group $G/N_i\cong C_p\rtimes C_3$: namely, let $y_1,y_2$ be the two primitive third roots of unity in $\mathbb{F}_p$ and let $N_i:=\{(x,y_i x, y_i^2 x): x \in \mathbb{F}_p\}\trianglelefteq K$. We want to find an expression in $\alpha,\beta,\gamma$ which is fixed by $N_i$, and claim that this can be done in the form $\alpha^j\beta^k$ for $j,k \in \{1,\dots, p-1\}$ suitable, so that the minimal polynomial of $\alpha^j\beta^k (=a^jb^k)^{1/p}$ over $\mathbb{C}(t)$ will have Galois group $G/N_i\cong C_p\rtimes C_3$ (the remaining symmetrization $\alpha^j\beta^k + \beta^j\gamma^k + \gamma^j\alpha^k$ the only drops the degree by a factor $3$, reaching a degree-$p$ polynomial, but still with the same splitting field, so this is not so important anymore). So let $\sigma_i:=(1,y_i,y_i^2)\in N_i$ be a generator of $N_i$. This acts quite explicitly on the three blocks $\Delta_1,\Delta_2,\Delta_3$ by $\sigma_i(\alpha) = \zeta_p \alpha$, $\sigma_i(\beta) = \zeta_p^{y_i}\beta$, $\sigma_i(\gamma) = \zeta_p^{y_i^2}\gamma$. Hence, $\sigma_i(\alpha^j\beta^k) = \zeta_p^{j+k\cdot y_i} \alpha^j\beta_k$, and we're laughing if we can get $j+k\cdot y_i$ divisible by $p$. There's obviously some freedom here in choosing $j$. E.g., for $p=7$, the third roots of unity in $\mathbb{F}_7$ are $y_1=2$ and $y_2=4$; now if you choose $j=3$ with $y_2=4$, you get $k=1$, which after the symmetrization yields your first solution $x_7$; whereas $j=1$ with $y_1=2$ gives $k=3$, and that, after symmetrization is your second solution $x_7'$. Just to check one more case, for $p=37$, the third roots of unity are $10$ and $26$. Choosing $j=7$ with $y_1=10$ gives $k=3$. You could of course choose other exponents $j$ (whether choosing to increase $j$ by $1$ for each new prime is actually a natural choice, I don't know; it happens to give small $k$ for the first few terms, but if I made no mistake, choosing $j=8$ for $p=43$ should give you $k=13$, whereas $j=1$ would have allowed $k=7$), although these would of course not generate new field extensions.