Converting a discrete real time domain to a discrete complex time domain signal
My updated answer to OP's update

To do this specifically for the case of the FFT of a product of two sinusoids, you must choose a sampling rate or modulation frequency such that the upper sideband of the "positive" frequency component aliases to be on top of the lower sideband of the "negative" frequency component. For instance using a carrier frequency of 1.2MHz and a signal frequency of 50KHz would give you one peak in the center with two sidepeaks. (Although to note with 4095 points this wouldn't match the plot position since the center bin would be 2047 which is closer to where the lower sideband is in the plot).

Otherwise for an arbitrary signal the fft of the following would provide this result:

$$s[n] = Ae^{j\omega_c n}B\cos(\omega_m n)$$

Where

$\omega_c$ is the normalized radian frequency of the carrier ($0$ to $2\pi$)

$\omega_n$ is the normalized radian frequency of the modulation ($0$ to $2\pi$)

$A$, $B$ are arbitrary gain constants for each term (which simply scales the result by $AB$)

The opening part of the question was you are doing this to enhance your understanding of the frequency domain of AM modulation. So I am assuming you really want help understanding what AM modulation is supposed to look like:

When you multiply two sinusoids to get AM modulation, the result is DSB-SC (Double Side-band Suppressed Carrier): The carrier frequency is suppressed and there will only be two sidebands. This is what a cosine wave looks like in the FFT: two sidebands centered on the DC bin (use fftshift command if you want to center DC in the middle of the plot, then showing a positive and negative frequency axis which may help intuition), but since there is no DC in the cosine, that bin is supressed. So a cosine on its own is a DSB-SC signal with the carrier being DC, which is suppressed (we call that the baseband signal).  If you multiply this baseband signal with a exponential carrier $e^{j\omega_c n}$ as I did above, this causes the frequency to shift by $+\omega_c$   , moving what was centered at DC to now be centered at $\omega_c$. If you muiltiply the baseband signal with a sinusoid, you shift in two directions (moving the signal to the positive and negative frequency carrier) since a sinusoid contains two exponential frequencies as further explained below.  

What you show in the lower plot is Large Carrier AM, and specifically such a plot could be created from adding a DC offset prior to doing the upconversion:

$$s[n] = (C + A e^{j\omega_c n})B \cos(\omega_m n)$$

Converting the real signal from your product to a complex signal can of course be done, but the result of this would be the expected two sidebands only that you now have in either the positive frequency axis or negative frequency axis. This is the Hilbert Transform that will convert your real signal that has complex conjugate symmetric positive and negative frequencies to a complex signal that only has positive or negative frequencies depending on the sign of Q. 

My earlier answer regarding why the upper and lower sidebands are not equal magnitude in the FFT result

Either window your time domain data before taking the FFT or increase the length of your time domain data (and therefore your FFT). You are seeing the effects of “scalloping loss” which occurs when the actual frequency terms are not an bin centers.  Zero padding your time domain data before taking the FFT will also help to make the result look closer to what you would expect since that will interpolate more samples in the frequency domain resulting in balanced sum and difference frequency magnitudes.

Also the command fftshift() will move “DC” to the center of the axis to give you a plot with positive and negative frequencies which may be more intuitive. Keep in mind that a real signal will have positive and negative frequency components in the FFT result since each bin represents an exponential frequency term $e^{j\omega t}$, and as is clear in Euler's identity, a sinusoidal term is composed on two exponential frequency terms. For example:

$$2\cos(\omega t) = e^{j\omega t}+ e^{-j\omega t}$$

Your time domain repesentation would be a sine wave at 900 KHz and 1.1 MHz, so the result if you increase the time domain length by adding more samples or zero-padding will approximate the expected result showing two equally weighted frequency impulses in the positive frequency axis at 900 KHz and 1.1 MHz as well as two frequency impulses in the negative frequency axis at -900 KHz and -1.1 MHz.