Adding phase noise to a signal results in amplitude variation on a constellation - whats the mistake? *plots inside*
As the OP has described, the noise as shown is AM noise only as the constellation points are not changing in phase at all. A double sided spectrum with symmetric sidebands in phase with the carrier is amplitude noise:
Consider a reference carrier as 1 at angle 0 degrees for which we do the following to create sinusoidal sideband noise by adding to it a small sinuoid with with $\alpha$:
$$c_{n(am)} = 1 + \alpha cos(\omega_nt)$$
This is a real waveform which only changes in amplitude and all samples will be of zero phase (AM noise only). We can gain further insight by decomposing the cosine using Euler's Formula (of which the magnitude and starting phase of each are the two sidebands in the Fourier spectrum):
$$c_{n(am)} = 1 + \alpha/2 (e^{j\omega_nt} +e^{-j\omega_nt})  $$
However if the two sidebands were in quadrature relationship to the carrier, then for small $\alpha$ where the small angle approximation holds, the noise will be dominantly phase noise:
$$c_{n(pm)} = 1 + j\alpha cos(\omega_nt)$$
$$= 1 + j\alpha/2 (e^{j\omega_nt} +e^{-j\omega_nt})$$
For the case of noise density where the noise is relatively small so that the small angle approximation holds, the OP can simply multiply the sidebands representing the noise pattern by $j$, mapping the result to be PM instead of AM. Note from the graphics below that negating half the spectrum (for example multiplying all components in the negative half spectrum by $e^{j\pi}$ or 180°) will also accomplish converting AM to PM (to do this with an FFT use an odd number of bins, multiply all the $(N-1)/2$ upper bins by $-1$ and multiply the DC bin by j to have this come out perfectly):

For adding phase noise from a power spectral density in MATLAB or Octave, please see this great function with all that worked out from Alex Bur-Guy:
https://www.mathworks.com/matlabcentral/fileexchange/8844-phase-noise