JFET mismatch between mathematical model and SPICE

Question

As a sequel of this question, I built a mathematical model for the Id current of the following (same) JFET circuit:

For a simplified model, I took LAMBDA=0 and set the internal resistances Rd=0 and Rs=0. I then proceeded to calculate the roots for Vg_limit (the value of the input voltage that takes the JFET to the limit between saturation and triode regions) and Id_sat (the drain current in the saturation region). For simplification, I didn't consider the triode region in this analysis.

After doing some math, the roots for the Vg and the equation for Id_sat are obtained as:

$$Vg_{limit} = \dfrac{2 \beta R_{ds} V_{dd} \pm \sqrt{4 \beta R_{ds} V_{dd} + 1} + 1}{2 \beta R_{ds}^2} $$

$$Id_{sat} = \dfrac{2 \beta R_{s} V_{gt} \pm \sqrt{2 \beta R_{S} V_{gt} + 1} + 1}{2 \beta R_{s}^2} $$

where $$ R_{ds} = R_d + R_s \quad and \quad V_{gt} = V_g - V_{to}$$

Out of the two roots for each equation above, only one makes sense physically.

I tested this for many NJF models from LTspice and, while it worked for quite a few, I faced some discrepancies for some other devices, in which the BETA was particularly low. As an example, if I take a simplified version of the 2N4118A, you can see from the plots below that Id_sat from the model is totally different than the one from SPICE. Moreover, the Vds - Vgst curve doesn't touch the origin, which should happen, as the circuit is fed with a sine wave having amplitude equal to Vg_limit.

Can you maybe spot what I'm doing wrong here? As I mentioned, this works for other devices, in which the BETA is larger. As a matter of fact, if I replace the value in this model with a higher one (0.4m or higher, instead of 0.095m), the plots start to match.

Another funny fact that I noticed was that, when increasing the temperatures to absurdly high values, for instance 250ºC, the plots also start to match. Strange indeed, as I forced BETATCE=0 and VTOTC=0.

Here is the PySpice code:

import numpy as np
import matplotlib.pyplot as plt
from PySpice.Spice.Netlist import Circuit
from prefixed import Float


def main():
    F0 = 20  # [Hz]
    FS = 48000  # [Hz]
    temp_c = 25  # [celsius]
    T = 0.1  # [seconds]

    # Ground truth: spice model
    circuit = Circuit('JFET')
    my2N4118A = circuit.model("my2N4118A", "NJF",
                              Beta="0.095m",
                              Betatce="0",
                              Vto="-1.299",
                              Vtotc="0",
                              Lambda="0",
                              Rd="0",
                              Rs="0",
                              )
    # component values
    VDD = 9.0
    RG = 1E6
    RD = 4.4E3
    RS = 1E3
    RDS = RD+RS

    jfet_spice_params = {
        p.upper(): my2N4118A._parameters[p]
        for p in my2N4118A._parameters
    }

    BETA = Float(jfet_spice_params["BETA"])
    VTO = Float(jfet_spice_params["VTO"])

    # Roots for ID at the limit region, where
    # ID - BETA * (VDS)**2 = ID - BETA * (VDD - ID * RDS)**2 = 0
    id_limit = np.array([
        (2*BETA*RDS*VDD - np.sqrt(4*BETA*RDS*VDD + 1) + 1)/(2*BETA*(RDS**2)),
        (2*BETA*RDS*VDD + np.sqrt(4*BETA*RDS*VDD + 1) + 1)/(2*BETA*(RDS**2))
    ])
    vd = VDD - id_limit * RD
    vg = vd + VTO
    vs = id_limit * RS
    vds = vd - vs
    if (id_limit[0] > 0 and vds[0] > 0) and (id_limit[1] <= 0 or vds[1] <= 0):
        vg_limit = vg[0]
    elif (id_limit[1] > 0 and vds[1] > 0) and (id_limit[0] <= 0 or vds[0] <= 0):
        vg_limit = vg[1]
    else:
        vg_limit = None
        print("Error! VG limit sanity check failed")

    # Netlist
    circuit.V('Vdd', 'vdd', circuit.gnd, VDD)
    circuit.SinusoidalVoltageSource(
        'in', 'gate', circuit.gnd, amplitude=vg_limit, frequency=F0)
    circuit.R('Rg', 'gate', circuit.gnd, RG)
    circuit.R('Rd', 'vdd', 'drain', RD)
    circuit.R('Rs', 'source', circuit.gnd, RS)
    circuit.J('JFET', 'drain', 'gate', 'source', model="my2N4118A")

    plt.figure()
    plt.title('Id x t')

    # Spice model
    simulator = circuit.simulator(
        temperature=temp_c, nominal_temperature=temp_c)
    # analysis = simulator.dc(Vin=Vsl)
    analysis = simulator.transient(step_time=1/FS, end_time=T)
    vin_spice = np.array(analysis["gate"])
    vgs_spice = np.array(analysis["gate"]) - np.array(analysis["source"])
    id_spice = (np.array(analysis["vdd"]) - np.array(analysis["drain"])) / RD
    vgst_spice = vgs_spice - VTO
    vds_spice = np.array(analysis["drain"]) - np.array(analysis["source"])
    # if > 0 then saturation, else triode
    v_limit_spice = vds_spice - vgst_spice
    t = np.array(analysis.time)

    # Roots for IDsat, when LAMBDA = 0
    vgt = vin_spice - VTO
    BETA_2 = 2 * BETA
    BETA_2_RS_VGT = BETA_2 * RS * vgt
    # Out of the 2 existent roots, only the one with "- np.sqrt(...)" makes sense physically
    id_sat = (BETA_2_RS_VGT - np.sqrt(2 * BETA_2_RS_VGT + 1) + 1) / \
        (BETA_2*(RS**2))
    id_sat[vgt <= 0] = 0

    plt.plot(t, id_spice, label="Id_spice")
    plt.plot(t, id_sat, label="Id_model")

    plt.xlabel('T [s]')
    plt.ylabel('Id [A]')
    plt.margins(0, 0.1)
    plt.grid(which='both', axis='both')
    plt.legend()

    plt.figure()
    plt.title('V x t')
    plt.plot(t, vgst_spice, label=f"Vgst")
    plt.plot(t, v_limit_spice, label=f"Vds-Vgst")
    plt.xlabel('T [s]')
    plt.ylabel('V [V]')
    plt.margins(0, 0.1)
    plt.grid(which='both', axis='both')
    plt.legend()

    plt.show()


if __name__ == "__main__":
    main()

Answer 1

It turns out that I didn't consider the Id behavior when Vgs > 0. In this region, the gate-source PN junction (as hinted by Ste Kulov) starts to conduct and the Igs starts to play a major role in the circuit, effectively clamping Id as Vin further increases.

For my needs, the use-case when Vgs > 0 is not a valid one. Hence, I just need to limit my analysis for an input voltage below a threshold Vin_limit_pn, where Vgs=0. Hence, as the JFET is in the saturation region:

$$ I_{d\_limit\_pn} = \beta \times (0+V_{TO}^2) $$ $$ V_{in\_limit\_pn} = V_{gs} + V_s = 0 + R_s \times I_{d\_limit\_pn} = R_S \times \beta \times V_{TO}^2 $$