Work in progress:
So, my question is: why can't do we use any value (even and odd, spaced by 1 fd), to reduce even more the bandwidth of our signal? Has it relation to interference between symbols in frequency spectrum?
Let's first consider binary or $2$-FSK. We have a carrier signal at fixed frequency $f_c$. Modulation by the binary data changes the frequency of the transmitted signal to $f_c + f_d$ or to $f_c-f_d$. Thus the change in the frequency of the transmitted signal which occurs when the transmitter transmits a $1$ followed by a $0$ is $\pm 2f_d$; the two possible FSK tones at frequencies $f_c + f_d$ and $f_c-f_d$ are spaced $2f_d$ apart in frequency; not $f_d$ apart as the OP apparently believes. The transmitted signal consists of a succession of RF pulses of duration $D$. e assume that $f_c \gg \dfrac 1D$ so that each RF pulse has several periods of the RF sinusoid.
What we would like to arrange is for the two possible tones to be orthogonal over each signaling epoch, that is, the RF pulses $\cos(2\pi(f_c + f_d)t+\theta_1)$ and $\cos(2\pi(f_c - f_d)t+\theta_0)$ should satisfy $$\int_{kD}^{(k+1)D}\cos(2\pi(f_c + f_d)t+\theta_1)\cos(2\pi(f_c - f_d)t+\theta_0) \,\mathrm dt = 0$$ for all integers $k$. Why the insistence on orthogonality? Well, if the receiver consists of two matched filters for the two RF pulses and a comparator to determine which filter has larger response at the sampling instant, then the output of each matched filter is the cross-correlation function of the input RF pulse with the RF pulse to which the filter is matched, and at the sampling instant, this cross-correlation has value $0$ if the input is the other pulse, and large value (the peak value of the autocorrelation function, in fact) if the input is the one to which the filter is matched (cf. this answer for details.)
Now, assuming without loss of generality that $f_1 > f_0$, we have that \begin{align} & ~~~~~~~~\int_{0}^{D}2\cos(2\pi f_1t+\theta_1)\cos(2\pi f_0t+\theta_0) \,\mathrm dt\\ &= \int_{0}^{D}\cos(2\pi (f_1-f_0)t+\theta_1-\theta_0) +\cos(2\pi (f_1+f_0)t+\theta_1+\theta_0) \,\mathrm dt\\ &= \left.\frac{\sin(2\pi (f_1-f_0)t+\theta_1-\theta_0)}{2\pi (f_1-f_0)} + \frac{\sin(2\pi (f_1+f_0)t+\theta_1+\theta_0)}{2\pi (f_1+f_0)}\right\vert_{0}^{)D}\\ &= ~~~~~\frac{\sin(2\pi (f_1-f_0)D+\theta_1-\theta_0)-\sin(\theta_1-\theta_0)}{2\pi (f_1-f_0)}\\ &~~~~+\frac{\sin(2\pi (f_1+f_0)D+\theta_1+\theta_0)-\sin(\theta_1+\theta_0)}{2\pi (f_1+f_0)} \end{align} The numerators of both fractions in the last RHS above are bounded by $2$ while $2\pi (f_1+f_0)\gg 2\pi (f_1-f_0)$ and so the second fraction is considerably smaller than the first. But, if $f_1$ and $f_0$ are such that $(f_1-f_0)D$ is an integer, then the first fraction has value $0$ while if $(f_1+f_0)D$ is an integer, the second fraction has value $0$. Note that both conditions hold whenever $f_0D$ and $f_1D$ both are integers ($n$ and $m$, say) or both are half-integers ($\frac{2m-1}{2}$ and $\frac{2n-1}{2}$, say), and give us the desired orthogonality. We will restrict attention to the integer case only.
Applying the above notion to $$\int_{0}^{D}\cos(2\pi(f_c + f_d)t+\theta_1)\cos(2\pi(f_c - f_d)t+\theta_0) \,\mathrm dt = 0$$ we see that