function [Gamma, E_cr_eff, Gamma_vs_A] = schwinger_rate(A_range, E_field_range, J, qvc, xi, const)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  SCHWINGER_RATE  Schwinger-analog pair production in TEVC metamaterial
%
%  In QED, the Schwinger mechanism produces e+/e- pairs from vacuum
%  when the electric field exceeds the critical threshold:
%    E_cr^QED = m_e^2 c^3 / (e * hbar) ~ 1.3e18 V/m
%
%  In a narrow-gap metamaterial (TEVC), the excitonic gap Delta plays
%  the role of 2mc^2, and the coherence length xi plays the role of
%  the Compton wavelength. The effective critical field becomes:
%    E_cr^eff = Delta_cat^2 / (e * xi)
%
%  For Delta_cat ~ 1 meV and xi ~ 10 nm, this gives E_cr^eff ~ 10^4 V/m,
%  which is 14 orders of magnitude below the QED threshold!
%
%  The pair production rate (per unit volume per unit time):
%    Gamma = (Delta_cat^2 / (4*pi^3 * hbar * xi^2)) *
%            (E/E_cr) * exp(-pi * E_cr / E)
%
%  QVC enhancement: as A_ZPF increases, Delta_cat decreases toward zero
%  (the QCP), causing E_cr^eff to drop and Gamma to diverge.
%
%  Inputs:
%    A_range       - ZPF amplitude array [normalized]
%    E_field_range - applied field array [in units of E_cr]
%    J             - tunneling energy [J]
%    qvc           - struct with coupling parameters
%    xi            - coherence length [m]
%    const         - fundamental constants struct
%
%  Outputs:
%    Gamma       - pair production rate [arb.], size (NE, NA)
%    E_cr_eff    - effective critical field [V/m], size (NA,)
%    Gamma_vs_A  - rate in (E, A) parameter space (for phase diagram)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

NE = length(E_field_range);
NA = length(A_range);

% Boson coupling strengths
g_BJ  = qvc.g_BJ * J;
g_Bpl = qvc.g_Bpl * J;
g_BH  = qvc.g_BH * J;
Omega_R = 2*sqrt(g_BJ^2 + g_Bpl^2 + g_BH^2 + ...
    g_BJ*g_Bpl + g_BJ*g_BH + g_Bpl*g_BH);

% Compute gap and critical field for each A_ZPF
Delta_ex  = zeros(1, NA);
Delta_cat = zeros(1, NA);
E_cr_eff  = zeros(1, NA);

for ia = 1:NA
    A = A_range(ia);
    Delta_ex(ia)  = sqrt(max(0, 4*J^2 - 2*J*A*(2*J)));
    Delta_cat(ia) = max(0, Delta_ex(ia) - Omega_R/2);

    % Effective Schwinger critical field
    % E_cr = Delta^2 / (e * xi)
    % As Delta -> 0, E_cr -> 0 (pair production threshold vanishes at QCP)
    E_cr_eff(ia) = Delta_cat(ia)^2 / (const.e * xi);
end

% Pair production rate: Gamma(E, A_ZPF)
Gamma = zeros(NE, NA);
Gamma_vs_A = zeros(NE, NA);

% Reference prefactor scale (use J as gap scale when Delta_cat -> 0)
prefactor_max = J^2 / (4*pi^3 * const.hbar * xi^2);

for ia = 1:NA
    Dcl = Delta_cat(ia);

    for iE = 1:NE
        E_over_Ecr = E_field_range(iE);

        if Dcl > 1e-30 && E_over_Ecr > 1e-6
            % Standard Schwinger exponential
            prefactor = Dcl^2 / (4*pi^3 * const.hbar * xi^2);
            exponent = -pi / E_over_Ecr;
            rate = prefactor * E_over_Ecr * exp(exponent);
        elseif Dcl <= 1e-30
            % At QCP, gap = 0 => threshold vanishes => copious production
            rate = prefactor_max * E_over_Ecr;
        else
            rate = 0;
        end

        Gamma(iE, ia) = rate;
        Gamma_vs_A(iE, ia) = rate;
    end
end

% Normalize for plotting
Gamma_max = max(Gamma(:));
if Gamma_max > 0
    Gamma = Gamma / Gamma_max;
    Gamma_vs_A = Gamma_vs_A / Gamma_max;
end

end