function zpf = zpf_vacuum_catalyzation(A_range, k_range, E_range, ...
    tevc, qvc, topo, mat, const, meV)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  ZPF_VACUUM_CATALYZATION  Unified view of vacuum catalyzation effects
%
%  Computes all key QVC observables as functions of A_ZPF (the zero-point
%  field amplitude that drives the system toward the quantum critical point).
%
%  This function aggregates:
%    - Catalyzon gap Delta_cat(A)
%    - Tunneling probability T(A) at fixed energy
%    - Schwinger pair production rate Gamma(A)
%    - Hopfion energy E_H(A)
%    - Superfluid stiffness rho_s(A) and BKT temperature
%    - Vacuum entanglement entropy S(A)
%    - Effective coherence length xi(A)
%
%  These collectively define the "catalyzation landscape" — the parameter
%  space over which vacuum effects become experimentally accessible.
%
%  Inputs:
%    A_range  - ZPF amplitude array [normalized]
%    k_range  - wavevector array [1/xi units]
%    E_range  - energy array for tunneling [J]
%    tevc     - TEVC lattice parameters struct
%    qvc      - QVC drive parameters struct
%    topo     - topology parameters struct
%    mat      - material parameters struct
%    const    - fundamental constants struct
%    meV      - conversion factor (1 meV in Joules)
%
%  Output:
%    zpf      - struct with all computed observables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

NA = length(A_range);
J = tevc.J;
xi_0 = const.hbar / sqrt(2 * mat.mstar * J);

% Boson coupling
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);

% Interaction parameter
gn = 1.3;

% Preallocate outputs
zpf.Delta_ex         = zeros(1, NA);
zpf.Delta_cat        = zeros(1, NA);
zpf.T_tunnel         = zeros(1, NA);
zpf.Gamma_schwinger  = zeros(1, NA);
zpf.E_hopfion        = zeros(1, NA);
zpf.rho_s            = zeros(1, NA);
zpf.T_BKT            = zeros(1, NA);
zpf.S_entanglement   = zeros(1, NA);
zpf.xi_eff           = zeros(1, NA);

for ia = 1:NA
    A = A_range(ia);

    % --- Gap structure ---
    zpf.Delta_ex(ia)  = sqrt(max(0, 4*J^2 - 2*J*A*(2*J)));
    zpf.Delta_cat(ia) = max(0, zpf.Delta_ex(ia) - Omega_R/2);

    % --- Tunneling probability (at E = V0/2 through 3 periods) ---
    V_eff = tevc.V0 * (1 - A);
    V1_eff = tevc.V1 * (1 - A);
    E_probe = tevc.V0 * 0.5;
    zpf.T_tunnel(ia) = transfer_matrix_1d(E_probe, V_eff, V1_eff, ...
        tevc.lambda_M, 3, mat.mstar, const.hbar);

    % --- Schwinger rate (at E/E_cr = 0.5) ---
    E_cr = zpf.Delta_cat(ia)^2 / (const.e * xi_0);
    if zpf.Delta_cat(ia) > 1e-30
        prefactor = zpf.Delta_cat(ia)^2 / (4*pi^3 * const.hbar * xi_0^2);
        zpf.Gamma_schwinger(ia) = prefactor * 0.5 * exp(-pi / 0.5);
    else
        % At QCP: copious production (gap collapsed)
        zpf.Gamma_schwinger(ia) = J^2 / (4*pi^3 * const.hbar * xi_0^2);
    end

    % --- Hopfion energy (Q_H = 1) ---
    zpf.rho_s(ia) = max(0, (1 - A)^2);
    zpf.E_hopfion(ia) = 4*pi^2 * zpf.rho_s(ia) * ...
        topo.aH_over_xi * topo.QH * J * (1 - 0.3*A);

    % --- BKT temperature ---
    zpf.T_BKT(ia) = (pi/2) * zpf.rho_s(ia);

    % --- Vacuum entanglement entropy ---
    % S = sum_k S_k, where S_k = (n_k+1)ln(n_k+1) - n_k ln(n_k)
    % and n_k = |v_k|^2 is the Bogoliubov occupation
    S_total = 0;
    Dg = zpf.Delta_ex(ia) / J;
    for ik = 1:min(length(k_range), 100)  % sample subset for speed
        k = k_range(ik);
        eps_k = k^2 / 4;
        Ek = sqrt(max(0, eps_k*(eps_k + 2*gn) + Dg^2/4));
        if Ek > 1e-10
            nk = 0.5 * ((eps_k + gn + Dg/2) / Ek - 1);
            nk = max(nk, 1e-15);
            S_k = (nk+1)*log(nk+1) - nk*log(max(nk, 1e-15));
            S_total = S_total + S_k;
        end
    end
    zpf.S_entanglement(ia) = S_total;

    % --- Effective coherence length ---
    % xi_eff diverges at QCP as Delta_cat -> 0
    % xi_eff = hbar / sqrt(2 m* Delta_cat)
    if zpf.Delta_cat(ia) > 1e-30
        zpf.xi_eff(ia) = const.hbar / sqrt(2 * mat.mstar * zpf.Delta_cat(ia));
    else
        zpf.xi_eff(ia) = 1e-6;  % cap at 1 micron (QCP divergence)
    end
end

% Normalize Schwinger rate for plotting
Gamma_max = max(zpf.Gamma_schwinger);
if Gamma_max > 0
    zpf.Gamma_schwinger = zpf.Gamma_schwinger / Gamma_max;
end

end