Introduction
I will start from the end, namely answering why the two plots look so different. I point out first that the resulting sampled signals (one as a "nicer plot" sampled directly at 210 Hz, and the other as a "bad plot" first sampled at 2.1 Hz and then interpolated linearly to 210 Hz) are both sampled at the same rate of 210 Hz, and not just a plotting artifact as may be incorrectly concluded in haste. The "nicer plot" is also consistent with using better interpolation techniques to recreate the higher sampled signal, while the "bad plot" has been done with a First Order Hold (FOH) interpolation only.
One common approach of interpolating a signal to a higher sampling rate when the interpolation is an integer amount ($I$ times the original sampling rate with $I$ as any positive integer) is to first upsample through the insertion of $I-1$ zeros between every sample, followed by filtering to remove the resulting spectral images. Combining the two is referred to as "interpolation". This approach will lend itself nicely to a straight forward approach for quantifying the error, that once understood can easily be extended to non-integer cases.
In the OP's case, the filter is being done with a FOH filter. In general, although the Nyquist criteria is met, the closer the sampling is to the Nyquist criteria, the more complex is the filtering to remove those spectral images, and an FOH filter is not a very good one.
As to the side note question about the appearance of a superposition of sine waves, I refer to the OP to this post which answers that specifically.
Details
The OP is interested in the limitations of reconstruction using a First-Order-Hold (FOH) only. To understand this, I will first detail the zero insert and filter method of resampling to a higher sampling rate. An FOH interpolation creates new samples between the existing samples that land on a linearly connected line between the samples. This is in contrast to a Zero-Order-Hold (ZOH) which adds to new samples between the existing samples that are at the same value as the last previous sample at the input rate (similar to a "sample and hold"). For further details on the digital interpolation process through zero insert followed by filtering, please see this post which also quantifies the ZOH and FOH filters.
Once the resampling process is clear (if not, read post just linked above!), it is helpful to understand first the "perfect and unrealizable" case, and then from that we can best understand the the trade space between distortion, occupied bandwidth, filter complexity and filter delay with realizable implementations (what we can implement on actual signals). For this I will bring up the term "First Nyquist Zone", which refers to the spectrum that extends from $\pm f_s/2$, where $f_s$ is the sampling rate, and "Nyquist Frequency" which is $f_s/2$. The Nyquist Criteria states that perfect reconstruction is feasible if the spectrum of an analog waveform is completely limited to the First Nyquist Zone (no energy beyond $\pm f_s/2$) prior to digital sampling. For a real waveform, all of the spectral single-sided bandwidth $B$ must be strictly less than the Nyquist frequency: $|B| < f_s/2$. That's what we usually first learn. The practical details that we may often not really get until much later in our education, often when we try to build real hardware, are summarized below:
It is impossible to have no spectral energy beyond the Nyquist frequency. There will always be some amount of aliasing and folding. To minimize this distortion, we enlist the use of anti-alias filters prior to any analog filtering. (Anti-alias does not mean no aliasing, we just reduce it to levels we no longer care about). (To understand anti-alias filter considerations in more detail please see this post).
Filter complexity and delay are proportional to the inverse of transition bandwidth. The steeper a filter cut-off is, the more time delay will be introduced and the more complex (costly) the filter will be. For going from analog to digital, this effects the order or number of poles in the analog anti-alias filter mentioned above. More pertinent to the OP's question for going from a lower sampling rate to a higher sampling rate, this effects the complexity in terms of number of filter coefficients needed in the digital interpolation filter used (more details on this specifically below). Also similarly for going from digital to analog, this effects the order of the analog reconstruction filter used after the digital to analog converter. This relationship is key to the OP's question on distortion, as ultimately for all of the above filter designs, we are trading that filter complexity (and delay) with an allowable distortion to meet a given specification. (To see common "rules of thumb" for estimating the number of coefficients in a digital FIR filter, see this post.)
An approach to perfect (and unrealizable) reconstruction of a waveform with spectral content anywhere within the first Nyquist band ($|B| < f_s/2$) to an integer higher sampling rate $I f_s$ is inserting $I-1$ additional zeros in between the samples at the original sampling rate $f_s$ to get to the higher sampling rate $I f_s$. This increases the sampling rate, but introduces spectral images at higher frequencies which must be filtered out to complete the interpolation process. This filter, if correctly done, essentially "grows" the zeros that were inserted to the proper values to represent the original waveform. Digitally, such filtering is done by convolving those resulting samples with a sampled Sinc function as the inverse Fourier Transform of a brick-wall filter with single-sided bandwidth of $f_s/2$. This equivalently means digital filtering using a time-domain filter with an impulse response that is samples of a normalized Sinc function given as:
$$h[n] = \frac{\sin(2\pi n/I)}{\pi n}, \space -\infty < n < \infty$$
Such a filter would perfectly reconstruct or interpolate the digital samples to any higher rate with no distortion. Unfortunately, since $h[n]$ has infinite time support, and is non-causal, it is not physically realizable as an actual filter. What is instead implemented very effectively is to do the following to modify $h[n]$, completing the digital filter design:
Reduce the filter frequency cutoff to not extend all the way to Nyquist (this allows for an increased transition band, simplifying the filter. 70 to 90% is typical).
Delay and truncate $h[n]$ for causality and finite time support.
Multiply the truncated and delayed Sinc with a window (a Kaiser window is a good choice) to reduce edge effects from truncation.
The closer the filter cutoff is to Nyquist, the more challenging the actual implementation of the filter will be. Typical implementations use a sampling rate such that the filter cutoff is at 70% to 80% of Nyquist, the actual value is traded with filter complexity.
(UPDATE IN PROCESS, REMAINS OF ORIGINAL ANSWER BELOW)
One good way to have more intuition with this is to take a windowed FFT of the sampled result to see the resulting waveform in frequency, before and after inserting the zeros as described above. The ideal reconstruction filter will pass the original signal without distortion and reject the images that the zero insert will create. A ZOH or FOH interpolation as a reconstruction filter is not very effective in both not distorting the original signal and rejecting the images. It is easy to see this from the frequency response of a ZOH (which approximately has a Sinc shaped frequency response as the "Dirichlet Kernel") or FOH (which approximately has a Sinc squared frequency response and specifically the Dirichlet Kernel convolved with itself). I demonstrate this below:
First is a plot of the FFT of the OP's signal, where I both extended the total time duration from 8 seconds to 80 seconds and then windowed the result with a Kaiser window prior to taking the FFT. Extending the time results in better frequency resolution (the width of the peaks shown is proportional to the inverse of the time duration), and windowing, which degrades the frequency resolution, is necessary to reduce the sidelobes (spectral leakage). The lower plot is the result after inserting 9 zeros in between every zero (upsampling). We see the desired increase in sampling rate (10x) but the introduction of many images which at this point is distortion.
Note this important detail, which is a motivation to do "zero-insert" for interpolation: By inserting zeros, we minimize the distortion of our primary signal within the spectrum that it occupies. The distortion manifests itself as images only (I say minimize as in practical implementation where we don't have infinite time, the resulting sidelobes, which in this case are > -60 dB down will be a shared spectrum distortion). See the plot below where I super-imposed the zoomed-in view of the interpolated spectrum with the original spectrum (we could do more detailed tests which I don't show here that would confirm the distortion level is as I described). This demonstrates how with zero-insert, our goal then is to come up with a filter to pass that original signal, and reject the closest image shown, as well as all higher images.
At this point, let's first see what a ZOH or FOH would do in this case in meeting this objective. For this interpolate by 10 case, the ZOH is done with an FIR filter with coefficients as ten ones; we convolve our zero-inserted signal with [1,1,1,1,1,1,1,1,1,1] which has the same effect as holding each value for 10 samples. Thus the resulting distortion is easily viewed in the frequency domain, showing the zero-insert signal together with our ZOH filter:
What we see here is how miserably the ZOH filter fails in both not distorting the original spectrum (instead we get "passband droop") and eliminating the images. With a single tone, the passband droop itself isn't a problem, but clearly we see how poorly the images will be reduced. Not a good result. I don't think I need to elaborate further to actually show the FOH case (we could get it easily by passing it through the ZOH filter twice, that would match the FOH result). Instead I will focus on a realizable solution using a least squares filter designed to work for all signals from DC up to the frequency of the OP's cosine wave at 1 Hz, when the sampling rate is only 2.1 Hz. To do such a feat, down to the distortion level I was willing to accept (>-70 dB), required 1021 taps for a least squares designed filter (a Kaiser windowed Sinc as the coefficients with 1021 taps would have a similar result), with the comparative results shown in the plot below:
Zooming in at the transition band shows the effectiveness of this approach, with complete predictability of the resulting distortion with this frequency domain analysis technique:
(Note for purposes of visibility in the above plot I had increased the total time duration of the cosine wave for the FFT observation; reducing the resolution bandwidth - this has no effect on the filtering approach, it is just for viewing the FFT).
So with that after confirming filtering performance, I filter the zero-inserted signal which used only the OP's original signal, with the following result below, showing the reconstruction of the 1 Hz cosine wave, and the resulting delay and initial settling time as the only deviations from representing our original signal (along with distortions that are > 70 dB down in this case).
If we zoom in on the reconstructed signal after 40 seconds, we see how well the original 1 Hz cosine wave has been recovered from the OP's samples, which in this case is the cosine wave sampled 10x higher (and connecting those samples as an unimportant artifact of the plotting method I used).
We note again that this is all done from the original samples only, of the 1 Hz cosine sampled at 2.1 Hz. Those samples describe perfectly the 1 Hz cosine as the unique solution within the "First Nyquist Zone" extending from $\pm f_s/2$ where $f_s$ is the sampling rate. Beyond that there are infinitely many other solutions that could also be represented by these samples.
What would have happened if we did actually proceed and use that ZOH as the filter, what would the final result look like instead? I am glad you asked! Basically it's an unfinished job. We can of course proceed from here with additional filtering, along with compensation for the passband droop previously mentioned (important if our actual waveform occupies more of that first Nyquist zone than just a single tone). Here is the result:
Python code used in addition to OP's (using cool code format descriptor OLGD introduced me to):
import scipy.fft as fft
tf = 80 # duration increased to achieve sufficient freq res in FFT plots
N=len(y1)
pad = 50*N # zero pad for visual interpolation of more freq samples
freqout = fft.fftfreq(pad, 1/ws1) # freq axis
fout = fft.fft(y1 * sig.kaiser(N, 6), pad)/ N # FFT with window
# upsampling by 10 by taking original samples only and from that recreating
# higher sampled version of same cosine signal
# insert 9 zeros between each original sample
interp = 10 # interpolation rate
y2 = np.zeros(interp*N)
y2[::interp] = y1
freqout2 = fft.fftfreq(pad, 1/(interp*ws1))
fout2 = fft.fft(y2 * sig.kaiser(len(y2), 6), pad)/ N
plt.figure(figsize=(6,4))
plt.subplot(2,1,1)
plt.plot(fft.fftshift(freqout),20*np.log10(fft.fftshift(np.abs(fout))))
plt.axis([-1.2, 1.2, -80, 0])
plt.grid()
plt.title("FFT of Cosine Sampled at 2.1 Hz")
plt.ylabel('dB')
plt.subplot(2,1,2)
plt.plot(fft.fftshift(freqout2),20*np.log10(fft.fftshift(np.abs(fout2))))
plt.axis([-12, 12, -80, 0])
plt.grid()
plt.title("FFT of Samples after Zero Insert- Upsample by 10")
plt.xlabel("Freq (Hz)")
plt.ylabel('dB')
plt.tight_layout()
# reconstruction filter: pass 1 Hz, block 1.2 Hz
# ZOH filter
coeff = np.ones(interp)
# filter orginal signal
output = sig.lfilter(coeff, 1, y2)
# determine and plot frequency spectrum
w, h = sig.freqz(coeff, interp, whole=True, worN=2**16)
plt.figure(figsize= (6,4))
plt.plot(fft.fftshift(freqout2), 20*np.log10(fft.fftshift(np.abs(fout2))))
plt.plot((w-np.pi)/(2*np.pi) * interp * ws1 , 20*np.log10(np.abs(fft.fftshift(h))))
plt.axis([-12, 12, -80, 0])
plt.grid()
plt.title("Spectrum after Zero-Insert and Freq Response of ZOH")
plt.xlabel("Freq (Hz)")
plt.ylabel('dB')
plt.tight_layout()
# least squares filter
coeff = interp*sig.firls(1021, [0, f0, ws1-f0, interp*ws1/2], [1,1,0,0], fs=interp*ws1)
# (repeat above with these coefficients to get high performance result)
