The OP's opening statement is incorrect:
$f_s > f_{max}/2$ prevents frequency aliasing for a bandlimited
signal, but not amplitude aliasing
$f_s > 2 f_{max}$ prevents aliasing. It's as simple as that. There is no such distinction as "amplitude aliasing". Since the OP has stated the signal is band-limited; as long as we can assume that means all spectrum is confined to be less than $f_s/2$, then there is no aliasing involved in the resulting samples. Those samples represent the sampled sinusoidal tone exactly with no additional distortion (beyond quantization based on whatever precision was used in the sampling).
Since the OP has brought up aliasing specifically and not reconstruction, the focus of my answer here is to show how aliasing can result in the same digital samples given for the cases when the sampled waveform is not band-limited, and specifically and significantly why aliasing can never occur in a band-limited waveform with all spectral content below $f_s/2$.
"Aliasing" or "Imaging"? Aliasing occurs in the process of going from continuous-time to discrete-time (sampling an analog waveform). When doing the opposite, going from the discrete samples back to an analog waveform, then reconstruction and interpolation explanations where imaging is a factor would be applicable and of interest: see Fat32's good answer and RBJ's comments under the original question which is from this perspective. This post also has a nice graphic to further help distinguish the two: https://electronics.stackexchange.com/questions/267408/aliasing-vs-imaging-what-is-imaging#:~:text=Aliasing%20is%20when%20a%20higher,zero%20padding%22%20and%20so%20forth.
Aliasing is when a waveform at another frequency location appears in (aliases into) our unique digital spectrum that extends from DC to half the sampling rate (for real signals). Aliasing can indeed distort signal amplitude when one waveform that does exist below $f_s/2$ in frequency interacts with an aliased waveform with spectral content above $f_s/2$ that wasn't filtered out prior to sampling. However, that is not what is occurring here as the OP appears to suspect: The amplitude of a pure sinusoid is NOT being distorted. Each of the samples given (to the extent of rounding precision) are the exact mathematical results for the samples of a sine wave at each of those points in time.
The interesting points in this question that I will further detail are as follows:
I show that the resulting samples can certainly be explained by aliasing in that the resulting samples would result from either a tone with a frequency < $f_s/2$ (direct sampling in the first Nyquist zone), OR from a more interesting under-sampled waveform (aliased) that has an envelope similar to what the OP has plotted (AM modulated waveform). This does not imply the signal the OP has plotted is distorted from aliasing; I am merely showing how the resultant digital samples, when we have no other information as to their origin, can be equally created from a waveform that was aliased (which alone would require the analog signal to not be band-limited).
The direct representation of the digital samples as an analog waveform is a stream of weighted impulses in time so would be better represented in the OP's plot using a stem plot rather than a line plot as was done. A stream of impulses in time is represented in frequency also as a stream of impulses. Selecting and averaging two tones from this stream results in the amplitude-modulated waveform cases I have plotted here; as we add more and more of the frequency tones the result would approach what would be shown in a stem plot as the quantity of tones approaches infinity. I demonstrate this as a final plot as well.
I also quantify the error between the closest maximum sample and the actual peak of a sampled tone, as a function of sampling rate and tone frequency (as I first thought this was of actual mathematical interest to the OP). It must be emphasized with regards to a sampled sinusoid that the result isn't an actual "error" since according to Nyquist theory we can recreate the noise free signal exactly for the cases when we know the sampling rate is greater than twice the maximum bandwidth of the bandlimited signal.
How Aliasing Can Occur
Aliasing occurs in the the process of going from an analog waveform to the discrete samples presented, and any explanation with regards to aliasing is with that process in mind.
A further intuitive view is further detailed below explaining how one might view the resulting appearance of an envelope in the OP's plot as an "alias" when there is no bandwidth limiting restriction on the signal that was sampled, here specifically it is a alias of an amplitude modulated waveform that "aliases" to a single tone.
What is visible here from sampling a real sinusoidal tone as that tone located at $f_s/2-f$ approaches $f_s/2$, is the interaction of that tone with a image located at $f_s/2 + f$, resulting in what would be identical to an amplitude modulated signal with a suppressed carrier at $f_s/2$ and two sidebands each at $f_s/2 \pm f$ (referred to as Double Sideband Suppressed Carrier or DBS-SC), where similarly the envelope would be sinusoidal with frequency equal to $(f_s/2-f)$, where $f_s$ is the sampling rate and $f$ is the frequency of the real tone. [And as RBJ points out and may be easier to visualize for some is the simple additive beating of those two tones].
This is demonstrated in the graphic below showing two cases of a sinusoidal tone sampled at $f_s=20 Hz$. In each case the resulting digital spectrum is the convolution of the input analog spectrum with the spectrum of the sampling process (since the sampled signal is the result of multiply in time a stream of impulses with the time domain sine wave waveform, and multiplication in time is convolution in frequency). The first case is with the tone at $3$ Hz while the second case is with the tone at $8$ Hz. Those familiar with DSB-SC would recognize the spectrums, in that within the unique digital sampling span of $f=0$ to $f_s$ the resulting waveforms are identical to a DSB-SC AM modulated waveform with the carrier at $f_s/2$. This is also apparent in the time domain waveform that the OP presented; that it is identical to such a modulation, which has an envelope whose frequency matches the modulation signal.
The approach to completely prevent this artifact is to use a complex signal and not a real signal, meaning the tone would be represented as $e^{j\omega t}$ instead of $cos(\omega t)$ and the sampled result given as real and imaginary components would have a constant magnitude for all samples, and there will be only one tone in the unique digital span from $f=0$ to $f_s$, as depicted in the plot below for this case.
To do this with a real signal implies that on every cycle, the peak of the sine wave is hit exactly, which can only occur with sampling rates that are an integer multiple of the tone itself (and aligned as such in phase).
Note the additional plots below in case they offer further insight into the mathematical equivalence of a sampled AM modulation (DBS-SC) with the envelope predicted by $f_s/2-f$ as the envelope line in these plots show.
These plots show the actual signal at frequency $f$ in blue, ranging from $9.7$ Hz to $1.7$ Hz, and the resulting samples when sampled with a $f_s = 20$ Hz clock (to match the figures above). In beige I also plot the envelope in each case with the frequency as predicted by $(f_s/2-f)$. The second plot in each case is the equivalent DSB-SC modulation at the carrier of $fs/2$ and sidebands $\pm(f_s/2-f)$ also sampled by the same clock (same samples) along with beige envelope waveform superimposed.
Notice in all plots the equivalence to a sampled DSB-SC AM modulated signal with a carrier at $f_s/2$ Hz and sidebands offset from the carrier by $\pm(f_s/2-f)$ Hz, which is the frequency of the envelope shown. Also observe when we sampled at an integer multiple of the tone, (2 Hz case) we could land on the peaks consistently on every cycle.
In Summary
These plots were done by selecting just two of the frequencies from the spectrum plots above showing the sampled spectrum as a stream of impulses. We could select ANY number of these (and only these) and average to create all possible analog waveforms that would result in the same digital sequence once sampled. The plots above were done with two of the tones, and another similar plot below shows what would occur in the first case with a 9.7 Hz sampled tone by selecting 10 of the tones in the extended digital spectrum ("Extended Digital Spectrum" meaning the digital spectrum we would get if we extended the frequency axis beyond the unique span out to $f_s/2$ as I had shown in the spectrums of the sampled signal). This is to the perspective that Fat32 makes in his answer and RBJ makes in the comments under the original question with reference to interpolation filters and reconstruction: this would apply when going from the discrete samples in time to the equivalent analog signal, and the ideal interpolation filter (given by Sinc) eliminates all but the first of these tones. The digital samples alone do not specify what is in between (which is the reason we do have viable implementations under-sampling higher frequency waveforms- here the OP said "band-limiting" so it can only be a single tone less than $f_s/2$ but in general nothing dictates that the originating analog waveform must be as such.) When we begin to interpret the plot given by discrete samples by adding lines between the points is when any of these effects and explanations would come into play, depending on how we choose to connect those dots. When we only have the dots, we cannot generally say it was a single tone that was sampled or any of these other combinations of alias waveforms, but if we also say that the signal is bandlimited less than $f_s/2$ as the OP has done, then the only possible answer in this case would be a single tone.
Derivation of Magnitude Difference Signal:
This is what I first thought the OP was primarily interested in, was quantifying the magnitude of the peak difference between the actual peak magnitude of a sampled sinusoid and the magnitude of the closest sample. The OP had suspected that there would be some frequency where this difference would go to zero. That is not true as will be shown, but we can make that difference arbitrarily small through over-sampling. (and through proper reconstruction, this difference is no form of a distortion since we can perfectly reconstruct the signal according to the Nyquist-Shannon sampling theorem). I don't yet see a practical utility for what follows, but it would perhaps be of interest if anyone needed a mathematical description and further properties of this difference signal.
The difference in magnitude is given by the following equation:
$$\epsilon_p = 1 - \cos\bigg(\pi\frac{f_T}{f_s}\bigg)$$
$$ = 2\sin^2\bigg(\pi\frac{f_T}{2f_s}\bigg) \tag{1} \label{1}$$
From the small angle approximation ( $\sin(\theta) \approx \theta$ ) this becomes:
$$\epsilon_p \approx \bigg(\frac{\pi f_T}{\sqrt{2}f_s}\bigg)^2 \tag{2} \label{2}$$
Where:
$\epsilon_p$: peak error relative to the peak magnitude of the sine wave
$f_T$: frequency of the sinewave (in any units) $f_T \in [0, f_s/2] $ for \ref{2}, $f_T \ll f_s/2$
$f_s$: frequency of the sampling clock (in same units as $f_T$)
In units of normalized radian frequency this becomes:
$$\epsilon_p \approx \bigg(\frac{\omega_T}{2\sqrt{2}}\bigg)^2 \tag{3} \label{3}$$
Where:
$\omega_T =2\pi f_T/f_s$: normalized frequency in units of radians/sample, $\omega \ll \pi$.
Equations \ref{2} and \ref{3} are reasonably accurate for small $\omega$ as demonstrated in the error plot below, while \ref{1} is accurate for all frequencies in the first Nyquist zone ($f_T \in [0, f_s/2]$).
This relationship is derived from the graphic below showing two samples (the red dots) of the sample sinusoidal waveform (the blue curve) at the peak difference condition, illustrating a bound for all conditions; a sample will always be at most within $\delta/2$ radians of the peak location of the tone, where $\delta$ is the sample period in time converted to phase in radians for a tone of a given frequency according to:
$$\delta = T_s \text{ sec} \times \frac{2\pi}{T} \text{ rad/sec} = \frac{2\pi f_T}{f_s} \text{ rad}$$
This is plotted below on a dB and log frequency scale, where we see the slope is approximately 40 dB/decade consistent with the $(f_T/f_s)^2$ relationship given by \ref{2}, and the close match of this approximation over most of this range when viewed on a log frequency axis. The right end of this plot is $fs/2$ and this shows the maximum possible difference relative to the actual full scale of a tone for the closest sample on any given cycle of the tone that is sampled.
This plot is the envelope of the maximum possible peak difference, where the actual curve would be less than this under certain rational relationships between the frequency of the tone and the frequency of the sampling clock as further detailed below (for instance we can easily see how the difference can be zero when the sampling clock is an integer multiple of the tone frequency with zero time offset of the sampling clock and the location of the peak of the tone). The difference will also be reduced below this bound for other rational sampling relationships such that the difference forms a repeating pattern with time (avoiding the max difference conditions where the samples are symmetrical around the peak such as depicted in the diagram above showing the samples and the peak location on the sinusoid).
We will now consider an error waveform as being the resulting difference computed for each positive and negative peak of the sinusoid (at each peak we update the error waveform with a new result representing the magnitude of this difference such that we get a sequence of errors versus time), we would see that this error waveform is represented by samples on a parabolic function which itself is given by an equation of identical form to \ref{1}, specifically as depicted in the diagram below.
Where in this diagram $\Delta f$ is the difference in frequency between the sinusoidal tone of frequency f, and the closest sub-integer of the sampling clock of frequency $f_s$ (when the sampling clock is an integer multiple of the sinusoidal tone, $\Delta f = 0$), and $T_r$ is the repetition period of the error waveform.
Thus, $\epsilon_p$ is the peak magnitude of the underlying continuous-time error waveform that itself would repeat with time depending on the offset in frequency between the tone and the closest integer fraction of the sampling rate: When the sampling frequency is an integer multiple of the sampling rate, this error will be a constant value, (since the sample closest to the peak of the tone will land in the same position relative to that peak on every cycle), and will be zero error specifically only when the sampling clock is aligned it time to land exactly on the peak; shift the time offset by half a sample in this condition and the constant error will be $\epsilon_p$. Thus for the integer sampling case, the actual error will be constant on a point given by the red curve in the plot above, where $t$ would be the time offset of the sampling clock. If we changed the frequency slightly (of either the tone or the sampling clock), this time offset will continuously increase, and thus the resulting samples of the error versus time will be on this error function depicted here (cyclically repeating at rate $\Delta f$).
Python script for DSB-SC graphics:
import numpy as np
import matplotlib.pyplot as plt
fs = 20 # sampling rate
interp = 10 # interpolation emulating "continuous time" waveform
T = 2 # total time in seconds
f = 0.5 # frequency of tone
N = T*fs # number of samples
N2 = N * interp # number of high precision samples
t = np.arange(N)/fs # time increment for samples
t2 = np.arange(N2)/(fs*interp) # high precision time increment
sig = np.cos(2*np.pi*f*t2) # underlying cont-time tone
sig1 = np.cos(2*np.pi*f*t) # sampled tone
sig2 =np.cos(2*np.pi*(fs/2 - f)*t2) # envelope
dsbsc = np.cos(2*np.pi*fs/2*t2)*sig2 # DBS-SC signal
plt.figure()
plt.subplot(2,1,1)
plt.title(f"$f_s = 20 $ Hz, f = {f} Hz ")
plt.plot(t2, sig, label = "Actual Sig")
plt.plot(t2, sig2, linewidth =2, color= 'bisque',label = "Envelope")
plt.plot(t2, -sig2, linewidth =2, color = 'bisque', label = "Envelope")
plt.plot(t, sig1, 'ro', markersize = 4, label= "Samples")
plt.xlabel('Time [s]')
plt.ylabel('Magnitude')
plt.legend()
plt.subplot(2,1,2)
plt.title(f"DSB-SC 10 Hz Carrier, fs/2 - f = {fs/2 - f:.1f} Hz Mod ")
plt.plot(t2, sig2, linewidth =2, color= 'bisque',label = "Envelope")
plt.plot(t2, -sig2, linewidth =2, color = 'bisque', label = "Envelope")
plt.plot(t, sig1, 'ro', markersize = 4, label= "Samples")
plt.plot(t2, dsbsc)
plt.xlabel('Time [s]')
plt.ylabel('Magnitude')
plt.tight_layout()
