Adding phase noise to a signal results in amplitude variation on a constellation - whats the mistake? *plots inside* I am attempting to add phase noise to a modulated signal.

I am doing this by taking a phase noise mask (single sided 0 Hz to 1
MHz), which shows dBc values away from central value or reference
value. Phase noise is decreasing as we move in frequency away from
the reference point which is 0 Hz. [Top left image in the 4-plot
image below]
I then make this double sided.
I convolve this phase noise plot with a spectrum of a clean complex
exponential at 1 MHz, this is my oscillator for converting frequency.
I think this has imparted the phase noise onto the complex
exponential.  [Top right image in the 4-plot image below]

 
The spectrums are convolved  to make a oscillator with phase noise:
$\mathcal{F}(e^{j \omega t}) * \mathcal{F}( e^{j \phi(t) }) = \mathcal{F}(e^{j \omega t})* \mathcal{F}(e^{j B}) = \mathcal{F}(e^{j( \omega t + \phi(t))}) $

where $\mathcal{F}$ is fourier transform, $ \omega$ = 1 MHz, $\mathcal{F}( e^{j \phi(t) }) $ is the double sided phase noise mask on 0 Hz, then also this can be interpreted as a bandwidth $B$ of the upper and lower sidebands of the phase noise centered on 0 Hz (the double sided spectrum in the top left plot below when not in logarithmic x-axis scale). The convolution result is the top right image, symmetric phase noise translated from 0 Hz onto my a oscillator or CW or sine wave at 1 MHz.


I take the Inverse Fourier Transform of this convolution to produce
an oscillator in the time domain, my thought is that the phase noise
is on it.
I then multiply my complex modulated signal (16-APSK) by this oscillator (complex exponential with phase noise) to produce an up conversion and then downconvert with a clean complex exponential
with no phase noise. This should leave the phase noise on the
modulated signal.

The math now being done in time domain:
Up-conversion using the inverse transform of the result above
$A(m) e^{ j( \phi (m))}  e^{j( \omega t + \phi(t))} $
Down-conversion of the upconversion
$ A(m) e^{ j( \phi (m))}  e^{j( \omega t + \phi(t))} e^{-i \omega t}  $
where $A(m) e^{ j \phi (m)} $is a modulated signal


However I am showing results that are not like classical phase noise,
which would be rotation acrcs around constellation points but it
appears as an amplitude modulation?  [Final constellation image]



I know this is wrong, but I stuck why. If I subtract in the time
domain the clean sinusoid without phase noise and the oscillator with
phase noise, I get this variation across the full time domain
simulation. [Central images are the real and imag values of the oscillators in time and the bottom images are the difference]

Why is this producing amplitude variation and not phase variation?
Is there another way to do this?

--Edit Updated constellation based on replies below--
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