The frequency spectrum of cosine function is as below.
While adopting a square window function to the cosine function, it is same as implementing convolution of their frequency spectrum as below. The thick dash line is the sum of them. Is it correct to say the points on this dash line are the FFT result of corresponding frequency? I feel it is not true, but don't know where the problem is.
There is no problem in your plot if by fft you mean DTFT (or Discrete Time Fourier transform) and it is true to say "the points on this dash line are the DTFT result of corresponding frequency":
If your function is a cosine:
$$ x[n] = \cos(\omega_0n) \Longrightarrow X(e^{j\omega}) = \pi\delta(\omega-\omega_0)+\pi\delta(\omega+\omega_0)\ ; \qquad|\omega|<\pi $$
and your window function is defined as follow:
$$ w[n] = \left\{ \begin{array}{ll} 1 & \quad |n| \leq N \\ 0 & \quad |n| > N \end{array} \right. \Longrightarrow W(e^{j\omega}) = \frac{\sin\left(\omega(\frac{N}{2}+1)\right)}{\sin\left(\frac{\omega}{2}\right)} $$ as you mentioned if we multiply two functions in time domain, it's equal to periodic convolution in frequency domain: $$ y[n] = x[n]w[n]\\ \Longrightarrow Y(e^{j\omega})=X(e^{j\omega})\circledast W(e^{j\omega})=\frac{1}{2\pi}\int_{2\pi}^{}X(e^{j\omega})W(e^{j(\omega-\theta)}) d\theta $$ in this case (in one period): $$ Y(e^{j\omega})=\frac{1}{2}\frac{\sin\left((\omega-\omega_0)\left(\frac{N}{2}+1\right)\right)}{\sin\left(\frac{\omega-\omega_0}{2}\right)}+\frac{1}{2}\frac{\sin\left((\omega+\omega_0)\left(\frac{N}{2}+1\right)\right)}{\sin\left(\frac{\omega+\omega_0}{2}\right)} $$
if you plot this function in one period for $\omega_0=\frac{\pi}{5}$ and $N=2$ it look like this:
but if you change the window length to $N=4$ the result changes to:
and the cause for this, is that two main lobes overlap with one another.
That's one way to look at it conceptually. The DFT definition is only concerned with the bin values at integer $k$'s. As far as it is concerned there is nothing in between the bins.
However, there is nothing that mathematically prevents you from using a fractional k value in the DFT definition. This called the "trigonometric interpolation function" and is equivalent to a Fourier Series evaluation.
When you do, it is as if you set
$$ \omega = \frac{k}{N}2\pi $$
in the DTFT and applied a window function (implicitly) by using the samples 0 to N-1. Very similar to this with the FT
Amplitude after Fourier transform
The points on the dashed curve can also be found by summing two Dirichlet kernel functions (discrete version of the sinc) or you can use my more efficient combined forms found here:
The latter is for frequencies very close to a bin, or a different perspective.
They calculate the DFT at the bin values and will give a different answer in between for fractional $k$ values.
This is what is called "leakage" if you are off bin (not an integer number of cycles per frame) in the DFT. (The DFT is the proper name of the operation, FFT is a faster way to do it.)
