How can I generate a sine wave with consistent "vibrato"

I am trying to create a .wav file which contains a 440Hz sine wave tone, with 10Hz vibrato that varies the pitch between 430Hz and 450Hz. Something must be wrong with my approach, because when I listen to the generated .wav file, it sounds like the "amplitude" of the vibrato (e.g. the highest/lowest pitch reached by the peaks and troughs of the waveform of the vibrato) just progressively increases over time, instead of staying between 430-450Hz. What is wrong with my approach here? Here is some minimal python code which illustrates the issue:

import math
import wave
import struct

SAMPLE_RATE = 44100

NOTE_PITCH_HZ = 440.0        # Note pitch, Hz
VIBRATO_HZ = 10.0             # Vibrato frequency, Hz
VIBRATO_VARIANCE_HZ = 10.0    # Vibrato +/- variance from note pitch, Hz

NOTE_LENGTH_SECS = 2.0      # Length of .wav file to generate, in seconds

NUM_SAMPLES = int(SAMPLE_RATE * NOTE_LENGTH_SECS)

# Generates a single point on a sine wave
def _sine_sample(freq: float, sine_index: int):
    return math.sin(2.0 * math.pi * float(freq) * (float(sine_index) / SAMPLE_RATE))

samples = []
for i in range(NUM_SAMPLES):
    # Generate sine point for vibrato, map to range -VIBRATO_VARIANCE_HZ:VIBRATO_VARIANCE_HZ
    vibrato_level = _sine_sample(VIBRATO_HZ, i)
    vibrato_change = vibrato_level * VIBRATO_VARIANCE_HZ

    # Mofidy note pitch based on vibrato state
    note_pitch = NOTE_PITCH_HZ + vibrato_change
    sample = _sine_sample(note_pitch, i) * 32767.0

    # Turn amplitude down to 80%
    samples.append(int(sample * 0.8))

# Create mono .wav file with a 2 second 440Hz tone, with 10Hz vibrato that varies the
# pitch by +/- 10Hz (between 430Hz and 450Hz)
with wave.open("vibrato.wav", "w") as wavfile:
    wavfile.setparams((1, 2, SAMPLE_RATE, NUM_SAMPLES, "NONE", "not compressed"))

    for sample in samples:
        wavfile.writeframes(struct.pack('h', sample))

A more straight forward approach that will accomplish what you want is to use a phasor (linear ramp that goes from 0 to 1 then shoots back down to 0) to look up the sin of that value. Then, you can control the amount the phasor increments (the frequency of vibrato).

Here is the code. I lowered the sampling rate to make it easier to look at:

import math
import matplotlib.pyplot as plt

SAMPLE_RATE = 10000

NOTE_PITCH_HZ = 100.0        # Note pitch, Hz
VIBRATO_HZ = 20.0             # Vibrato frequency, Hz
VIBRATO_VARIANCE_HZ = 20.0    # Vibrato +/- variance from note pitch, Hz

NOTE_LENGTH_SECS = 2.0      # Length of .wav file to generate, in seconds

NUM_SAMPLES = int(SAMPLE_RATE * NOTE_LENGTH_SECS)

# Generates a single point on a sine wave
def _sine_sample(freq: float, sine_index: int):
    return math.sin(2.0 * math.pi * float(freq) * (float(sine_index) / SAMPLE_RATE))

phasor_state = 0
phasored_samples = []
samples = []
unmodulated_samples = []
for i in range(NUM_SAMPLES):

    # Generate sine point for vibrato, map to range -VIBRATO_VARIANCE_HZ:VIBRATO_VARIANCE_HZ
    vibrato_level = _sine_sample(VIBRATO_HZ, i)
    vibrato_change = vibrato_level * VIBRATO_VARIANCE_HZ

    # Mofidy note pitch based on vibrato state
    note_pitch = NOTE_PITCH_HZ + vibrato_change
    samples.append(_sine_sample(note_pitch, i)+5)
    unmodulated_samples.append(_sine_sample(NOTE_PITCH_HZ, i))
    phasored_samples.append(math.sin(2*math.pi*phasor_state)+10)
    phasor_inc = note_pitch/SAMPLE_RATE
    phasor_state += phasor_inc
    if phasor_state>=1:
        phasor_state -=1
plt.plot(unmodulated_samples, label='unmodulated')
plt.plot(samples, label='not working')
plt.plot(phasored_samples, label='using phasor')
plt.legend()
plt.show()

A zoom in on the output shows you the difference between these approaches:

Keep in mind though, that this still isn't quite right. A violinist or vocalist will vibrate up and down in a more or less linear trajectory, not a sinusoidal one. To be more 'correct' (if that is what you are going for, that is) would be to compute the change in phase increment as a triangle wave, not a sinusoidal one.