Artifacts signal conversion from float to integer
The quantization error is dependent on the input signal. In your case the input signal is a sine wave and the quantization error is a combination of its 8 odd harmonics that fit into the available bandwidth. You can reduce the dependency of the quantization error to the input signal by adding triangular dither before quantizing: Add to each sample a uniformly distributed random number between 0 and the quantization step, and subtract from the result another uniformly distributed random number between 0 and the quantization step. Then quantize. This way you'll get a flat spectral noise floor (no distortion of the signal) and no modulation of the quantization noise by the signal.

For further reading on dither: Robert Alexander Wannamaker's PhD thesis "The Theory of Dithered Quantization", 2003.

Modifying your script to include dithering fixes the problem:

import numpy as np
import random
import matplotlib.pyplot as plt
from scipy.signal import get_window
from scipy.fftpack import fft

l = 128
f = 4
x = np.arange(l) / float(l)
max_int = np.iinfo(np.int16).max
sine = np.sin(2 * np.pi * f * x)
dither = [random.triangular(-1, 1, 0) for _ in xrange(l)]
sine_96 = (sine / sine.max() * (max_int - 1) + dither).astype(np.int16)

w = get_window('hann', l)
S = np.abs(fft(sine * w))
S_96 = np.abs(fft(sine_96 * w))

fig, ax = plt.subplots(2, 2, figsize=(8, 7))
ax[0, 0].plot(sine * w)
ax[0, 1].plot(sine_96 * w)

ax[1, 0].plot(20 * np.log10(S[:l // 2]))
ax[1, 1].plot(20 * np.log10(S_96[:l // 2]))
fig.show()


Figure 1. Top left: windowed input signal, top right: dithered, quantized, and windowed input signal, bottom left: frequency spectrum of the windowed input signal (vertical axis in dB scale), bottom right: frequency spectrum of the dithered, quantized, and windowed input signal (vertical axis in dB scale).

Let's plot the quantization residual (error) of the undithered case before windowing, and the frequency spectrum of the residual after windowing it:

sine_96 = (sine / sine.max() * max_int).astype(np.int16)
residual = sine_96 - (sine / sine.max() * max_int)
fig, ax = plt.subplots(1, 2, figsize=(8, 3))
ax[0].plot(residual)
R = np.abs(fft(residual * w))
ax[1].plot(20 * np.log10(R[:l // 2]))
fig.show()


Figure 2. Quantization residual and its spectrum (window function applied before transforming).

Quantization is equivalent to a memoryless waveshaper shaped like a staircase. For a periodic input signal a memoryless waveshaper can only produce harmonics of the input signal, because the resulting output signal will be identically periodic:

$$f(x + s) = f(x) \text{ for all } x \in \mathbf{R}\\
\Rightarrow g\big(f(x + s)\big) = g\big(f(x)\big) \text{ for all } x \in \mathbf{R}$$

Here the input signal is a sine wave with a period of 128/4 = 32 samples. If the period would not be an integer number of samples, then also other frequencies than the harmonics would be generated, because quantization operates in discrete time and non-integer periodicities don't exist in discrete time.