I see, that's just my bad luck. I get $P_{21144099} −\sigma(21144097)=139^4$, but as luck would have it, $139^4-139^3+1$ is a prime :)

n = 21144097 seems to provide a counterexample

It is pretty rare for $\sigma(n)$ to not be even so, as Denis Shatrov points out the condition implies that $P_{n+2}-\sigma(n)$ is a power of an odd prime. If it is a cube of a prime, then $p^3-p^2+1$ is not going to be $3\mod 5$, so it has to be at least the fourth power of a prime. This is a rare event, which makes it hard to find by random search. Still, my guess is that the statement is unlikely to be true, as fourth powers of primes do occur among $P_{n+2}-\sigma(n)$.

I think that the problem can essentially be reduced to convex polytopes and piecewise linear functions. If there were a counterexample, one should be able to approximate both the function and the region, while keeping the centers of mass equality.

@MaxAlekseyev Thank you, interesting method!

I don't see it as straightforward. If one just looks at various cases of what you get modulo small primes, then the number of cases seems to be not feasible. In fact, if this could be done, then your original problem would also be manageable, since a reasonable proportion of $q$ would be a generator, hence quadratic non-residue.

What about just asking whether each of the first 1000 primes $q$ is a quadratic non-residue modulo $p$? It seems that this would happen roughly half the time, so you would get some 1000 bits of useful info. I don't know whether this is enough to reconstruct $p$, but it might be an easier problem.

Convexity is clearly important, even in dimension 3. Otherwise one can start with two pyramids from North and South poles, connected together by edges in the middle and then rotate one of the pyramids (while keeping the edge structure).

Thank you! I will first try Will Sawin's example, since I am comfortable with modular forms, but will come back to Will Chen's and your suggestions.

@WillChen Yes, definitely I should check that. I have seen some of your work, but not in detail, so if this turns out to be a $2$-torsion then it would definitely give me a reason to look further.

@WillChen Interesting reference, thank you! My $E$ is fairly ugly, with the $J$-invariant $(2661975021951 - 1086757026275 \sqrt{-7})/4598530048$, so computing Belyi maps could be computationally prohibitive. But it is definitely something to think about.

One can try to take the divisor $4(\hat p_1 + \hat p_2)$ where $\hat p_i$ are the preimages of $p_i$ which will likely give a map $X\to \mathbb P^2$, hopefully with image a nodal octic curve. The coefficients of the equation of the octic should be given by some algebraic functions in the two parameters of $M_{g=1,n=2}$, so if these can be understood to sufficiently high order near a specific point, then one should be able to find them.

I would deform to high order step by step and hopefully get something. Not a sure thing, but worth a shot. Thank you!

@WillSawin It's an interesting idea and I at first thought it might work, but now I think that the map $X(7)\to X(1)$ has other ramification points. So I am worried about additional ramification over $i$ and $\rho$. But it's easy to check!