ソニーのマイクの改造をしていた時に気が付いたのですが、増幅器の入力容量が変化すると、歪みが出ることがあるようです。
これについては過去に技術論文を見たことがないので、ひょっとしたら私が最初に声を上げているのかもしれません。

［2022.8.13 追記］
すみません。　その後の調査で私が最初ではないことが分かりました。　==> "Operational Amplifier Distortion, page 13" で Samuel Groner 氏が述べている第４の歪みがこれに該当します。

もう少し正確に言いますとJFET などの入力デバイスでは信号の電位によって入力容量が異なるため、 高いインピーダンスの信号源から受けると入力端でひずみが発生するということです。　 この歪みは入力端で既に発生してしまうので、増幅器内部で NF などによる改善はできません。
また、入力容量だけでなく、入力微分抵抗が変化するときにも同様のひずみが発生します。　 これについては過去に記事にしています。

この記事では、FET の入力容量変化時について、具体的な数字を入れてどの程度その影響が出るか検証してみたいと思います。

入力容量の変化がわかっていないと検証できませんので、素性が分かっている 2SK30A を先頭に置いたモデルを作ります。
差動にして、ソース・ゲート間電圧が一定になるように構成しました。

この回路はコンデンサ・マイクを信号源にしましたが、単に高い抵抗値の信号源の場合でも似たような結果になるでしょう。　 そちらの方がもっと計算が単純になるでしょうが、マイクの改造の時に気が付いたのでこの回路で解析することにします。

入力付近だけをモデル化したものです。　計算に使用するバイアスは K30A の Y を仮定し、1mA ほどのドレイン電流の位置とします。
ソース電圧はゲート電圧に従属して上下しますので、C4(6.5pF) には充放電は行われず、実質的に 0pF になります（ブートストラップ効果）。
C3 は、ドレインに負荷抵抗がある場合はミラー効果によって影響が大きくなりますが、今回はカスコードにして 4V に固定していますので単純に仕様書から読み取った値にします。
配線などの浮遊容量 C2 は 3pF を仮定しました。
ドレインからの漏れ電流がありますが、今回の趣旨にはあまり関係がないので無視することにします。

このグラフは東芝の 2SK30A の仕様書から抜き出したものです。
ゲート・ドレイン間の容量はグラフの Crss を使いますが、これは変化しますから、このグラフから -3V ないし -5V 辺りの近似式（ C3=f(e2) ）を導き出しておきます。
Crss が変化する原因は、ゲート電圧が下がると空乏層が深くなり、必然的に接合容量が少なくなるためです。　 FET である限り避けられません。　またバイポーラでも同様の現象があります。

C3 の近似式を作ります。　空乏層ができる原理から、双曲線に近いだろうと思って描いてみました。
k30A の仕様書の絵を Excel に張り付け、その上に半透明のグラフを乗せ、双曲線のカーブをフィットさせました。　 　3V から 5V の範囲だけですので、まあこれで十分でしょう。
x がドレインとゲートの電位差ですので極性をひっくり返し、4V のバイアスを加えて式は：
　　C3=6.9/(5.1-e2)+1.63　　　[pF]
になります。

歪みを調べるので入力 e1 をサインカーブとします。　 リーク抵抗への漏れ電流 i1 によって C 群の電荷が放電され、時々の電圧が変化しますから、周波数によって歪みが異なることが予想されます。　 なので、1Hz と 1kHz の２つで調べることにします。

１サイクルを 512 点の離散数値で表し、キャパシタへの充電電流を調べて順次端子電圧 e2(k+1) を計算していきます。　 512 という値は、後でフーリエ変換するときに２の累乗でないといけないので採用しました。　 この数字は多いほどフーリエ変換時の解像度が上がりますが、256 もあれば十分でしょう。　 ８波ほど作成しようと思います（全部で 4096 ポイント）。　波数が多いほど窓関数の影響が少なくなります。

パッと見て予想すると、リーク抵抗が無視できる高い周波数なら e2 の上端が +0.625V で下端が -0.633V ぐらいではないでしょうか。　 あ、直流分が出てくるから、どちらもわずかに下へシフトするかも。
これらの値は計算した後の検算に使います。

計算で気を付けないといけないのは、e2 が C3 の関数なのに、C3 が e2 の関数であるという、鶏と卵の関係にあるということです。　 i1、i2 も同様。　移項して e2 でまとめると長い２次式になりそうです。　 でも、この循環は収束が速いので２分木を使えば何とかなるでしょう。
２分木法を使うなら、 e2(k+1) を仮に etemp と置いて C3 や i1～i3 を求め、それらを使って e2 を試算してみます。　 そして etemp が その e2 と異なったら修正して一致するように追いかけるのです。　 別の方法もいくつかあります。

こうやって、e2 の電圧を時間軸で順次計算して積み上げていきます。
結果のオレンジの線のピーク値は予想値に近いものでした。

一つの点での計算結果には誤差があって、それを次点の計算に使いますから誤差が拡大する懸念があります。　 しかし、拡大することはなく、誤差は自動的に修正されて行きますから心配ありません。

もう一つ問題になることがあります。　初期値とそれ以前（過去）の電圧を e1 も e2 も０としているのですが、本来は違います。　 このため、e2 は頭の方でゆがんだ形になってしまいます。　 それを少しでも避けるため２サイクル目からの 512 点を捕まえることにします。
もう一つの手段として、いったん計算した後、513 番目か 1025 番目かの値を 1 番目の値（初期値）に放り込み、再計算する方法があります。　 これは効果的です。
両方やってみましたが、最終結果はどちらも同じでした。　 右図は 1Hz の頭の部分を拡大したもので、 Method-1 が前者の、Method-2 が後者の方法です。　 50 ステップもすれば両者は全く同じになります。　上に書いた「誤差は自動的に修正」が効いています。
どちらの方法かで捕まえた 512 点を８波コピーし、それにレイズドコサインの窓をかけてフーリエ変換します。

フーリエ変換した結果です。
見てわかる通り、R よりも C の効果の方が優勢となる高い周波数では、２次のひずみが 0.3% ほど入ることが分かります。

また、非直線であるため波形が上下で非対称になりますから、リーク抵抗を通して C 群に充電する電流と放電する電流が異なって、直流分が現れます。　 これは混変調歪みの原因になります。

1kHz のグラフをよく見ると、第５と第７の調波が第３よりも大きくなっています。　面白いですね。　 理由は追及していませんが、もしかしたら、C3=f(e2) が直線ではなく反っていることが原因かもしれません。

対策は：

1. 入力容量の小さいデバイスを使うこと。　これは高周波用のデバイスで実現できます。
2. 入力は大きい振幅にしないこと（といっても無理ですよね）。
3. ドレイン電圧を高くすること。　だだし、漏れ電流が増えてフリッカが増える可能性があります。
4. この例に挙げたようにカスコード(*)にしてドレイン電圧を振らないこと。　 ドレインから増幅結果を取り出すと C3 が［増幅率＋１］倍のミラー容量になって歪みが増大します。
5. 反転入力を使って増幅する。　 非反転入力を接地し、フィードバックを反転入力につなぐとゲート電圧が振れないのでこの問題は発生しません。　 可能ならば検討する価値はあります。

ついでの情報ですが、TI の OPA1641、1642、1644 という OPAMP は Cin が変動しない構造になっているようです。