結論

周波数 $f$ が時刻変化するとき、sin関数等に与えるべき位相は「 $2 \pi$ ×周波数」の積算値（積分）であって、 $2 \pi f t$ ではない。

実装例（C#）

時刻 $t$ における周波数 $f(t)$ は、 $440 \times 2^{\frac{2}{12}\sin(2 \pi \times 1 \times t)}$
時刻 $t$ における位相 $\phi(t)$ は、 $2 \pi \left(\int_{0}^{t} f(t) dt \right)$
よって、 $y = A \sin(\phi(t))$

const int SAMPLE_PER_SEC = 48000;
float length = 2;
float sample_size = SAMPLE_PER_SEC * length;
float[] waveform = new float[sample_size] ;
float phase = 0;
float base_freq = 440; // A4(hiA) 「栄光の架橋」の最高音らしい
float base_amp = 1;
float LFO_freq = 1;
float LFO_amp = 2 / 12; // 上下に全音一個分ずつ揺らす
for (int sample = 0; sample < sample_size; sample++)
{
    float t = sample / SAMPLE_PER_SEC;
   
    float LFO_value = LFO_amp * MathF.Sin(2 * MathF.PI * LFO_freq * t);
    float frequency = base_freq * MathF.Pow(2, LFO_value);
    float dt = 1.0F / SAMPLE_PER_SEC;


    phase += 2 * MathF.PI * frequency * dt;
    // 周波数を 2 * PI に収める処理
    if (phase > 2 * MathF.PI){
        phase -= 2 * MathF.PI;
    }

    float val = base_amp * MathF.Sin(phase) ;

    waveform[sample] = val;

}

備考

「周波数を $2 \pi$ に収める処理」について
これをしないと、２～３秒経ってから数値計算をミスっていそうな音が鳴るようになります。たぶん丸め誤差です。

失敗例

時刻 $t$ における周波数 $f(t)$ は、 $440 \times 2^{\frac{2}{12}\sin(2 \pi \times 1 \times t)}$ （定義）
波の方程式、 $y = A \sin(2\pi f t)$ （ $f$ が一定値なら正しい）
よって、 $y = A \sin(2\pi f(t) t)$ （間違い）

const int SAMPLE_PER_SEC = 48000;
float length = 2;
float sample_size = SAMPLE_PER_SEC * length;
float[] waveform = new float[sample_size] ;
float phase = 0;
float base_freq = 440;  
float base_amp = 1;
float LFO_freq = 1;
float LFO_amp = 2 / 12; 
for (int sample = 0; sample < sample_size; sample++)
{
    float t = sample / SAMPLE_PER_SEC;
   
    float LFO_value = LFO_amp * MathF.Sin(2 * MathF.PI * LFO_freq * t);
    float frequency = base_freq * MathF.Pow(2, LFO_value);

    phase = 2 * MathF.PI * frequency * t; // y = A sin(2 * pi * f * t) <-- 間違い

    float val = base_amp * MathF.Sin(phase) ;

    waveform[sample] = val;

}

これをやると、ピッチが宇宙の彼方へと飛んでいきます。それもそのはず、位相が $\phi (t) = 2\pi f(t) t$ のとき、周波数はこれを微分した値 $(2 \pi)(f(t) + f'(t) t)$ で、tが大きくなるとどんどん大きくなっていきます。

最初LFOを $f(t) = f_0 + \sin(2 \pi f_1 t)$ で実装してしまい、これが下手に解析的に積分できてしまうせいで、正しく $f(t) = f_0 \times 2^{\sin(2 \pi f_1 t)}$ （解析的に積分できないらしい）で実装したりADSRエンベロープを実装したりしようとしたときに、数値積分という発想に至らず、しばらく泥沼にはまってしまったので、書き残してみます。