function [E_hopf, Gamma_hopf, rho_s_arr, T_BKT_arr, hopf_phase] = ...
    hopfion_nucleation(A_range, T_range, J, qvc, topo, const)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  HOPFION_NUCLEATION  Hopfion topological excitation energetics
%
%  A hopfion is a 3D topological soliton classified by the Hopf invariant
%  Q_H in pi_3(S^2) = Z. In the TEVC context, hopfions are spinor
%  textures in the excitonic condensate order parameter.
%
%  The hopfion energy is:
%    E_H = 4*pi^2 * rho_s * a_H * Q_H * J
%  where:
%    rho_s  = superfluid stiffness (vanishes at BKT transition)
%    a_H    = hopfion core size (in units of xi)
%    Q_H    = Hopf charge (integer topological invariant)
%
%  The nucleation rate follows an Arrhenius-type law with quantum
%  corrections from vacuum fluctuations:
%    Gamma_H = Gamma_0 * exp(-E_H / kBT) * [1 + f(A_ZPF)]
%
%  Phase boundaries:
%    1. Stable BEC + hopfion: T < T_BKT and E_H > kBT
%    2. Hopfion melting: E_H < kBT but T < T_BKT
%    3. BKT proliferation: T > T_BKT (vortex unbinding)
%    4. Normal: T >> T_BKT (no condensate)
%
%  Inputs:
%    A_range  - ZPF amplitude array
%    T_range  - temperature array (units of J/kB)
%    J        - tunneling energy [Joules]
%    qvc      - QVC parameters struct
%    topo     - topology parameters struct
%    const    - fundamental constants struct
%
%  Outputs:
%    E_hopf     - hopfion energy [J], size (NA, 3) for Q_H = 1,2,3
%    Gamma_hopf - nucleation rate, size (NT, NA)
%    rho_s_arr  - superfluid stiffness, size (NA,)
%    T_BKT_arr  - BKT temperature, size (NA,)
%    hopf_phase - phase diagram, size (NT, NA), values 0-3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

NA = length(A_range);
NT = length(T_range);

% Rabi splitting
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);

% Preallocate
E_hopf    = zeros(NA, 3);     % for Q_H = 1, 2, 3
rho_s_arr = zeros(1, NA);
T_BKT_arr = zeros(1, NA);

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

    % Excitonic gap
    Delta_ex = sqrt(max(0, 4*J^2 - 2*J*A*(2*J)));
    Delta_cat = max(0, Delta_ex - Omega_R/2);

    % Superfluid stiffness
    % rho_s ~ (1 - A/(2J))^2 * (condensate fraction)
    % Vanishes at the QCP when A -> A_cr = 2J
    A_cr = 1.0;  % critical A_ZPF (normalized) where gap closes
    rho_s_arr(ia) = max(0, (1 - A/A_cr)^2);

    % BKT temperature: T_BKT = (pi/2) * rho_s
    T_BKT_arr(ia) = (pi/2) * rho_s_arr(ia);

    % Hopfion energy for Q_H = 1, 2, 3
    for iq = 1:3
        QH = iq;
        % E_H = 4*pi^2 * rho_s * a_H * Q_H * J
        % The factor a_H/xi is the core size
        % Larger Q_H => more energy (more linked preimage circles)
        E_hopf(ia, iq) = 4*pi^2 * rho_s_arr(ia) * ...
            topo.aH_over_xi * QH * J;

        % Quantum correction: ZPF fluctuations reduce the topological
        % barrier, making nucleation easier near QCP
        zpf_correction = 1 - 0.3 * A;  % linear softening
        E_hopf(ia, iq) = E_hopf(ia, iq) * zpf_correction;
    end
end

% Nucleation rate Gamma(T, A_ZPF)
Gamma_hopf = zeros(NT, NA);
Gamma_0 = 1.0;  % attempt rate (normalized)

for ia = 1:NA
    for iT = 1:NT
        T = T_range(iT);
        kBT = T;  % T is already in units of J/kB

        % Arrhenius with quantum tunneling correction
        E_barrier = E_hopf(ia, 1) / J;  % for Q_H = 1, in units of J
        if kBT > 1e-6
            % Thermal nucleation + quantum assist
            thermal = exp(-E_barrier / kBT);
            % Quantum tunneling through the topological barrier
            % (WKB-like, significant when T << E_H)
            quantum = exp(-sqrt(max(0, E_barrier - kBT)) / max(kBT, 1e-6));
            Gamma_hopf(iT, ia) = Gamma_0 * (thermal + 0.1*quantum);
        else
            Gamma_hopf(iT, ia) = 0;
        end
    end
end

% Phase diagram: classify each (T, A) point
hopf_phase = zeros(NT, NA);
for ia = 1:NA
    for iT = 1:NT
        T = T_range(iT);
        T_bkt = T_BKT_arr(ia);
        E_h1 = E_hopf(ia, 1) / J;  % hopfion energy in units of J

        if T < T_bkt && E_h1 > T
            % Phase 0: Stable BEC with frozen hopfions
            hopf_phase(iT, ia) = 0;
        elseif T < T_bkt && E_h1 <= T
            % Phase 1: Hopfion melting (thermal excitation of hopfions)
            hopf_phase(iT, ia) = 1;
        elseif T >= T_bkt && T < 2*T_bkt
            % Phase 2: BKT proliferation (vortex-antivortex unbinding)
            hopf_phase(iT, ia) = 2;
        else
            % Phase 3: Normal state (no condensate)
            hopf_phase(iT, ia) = 3;
        end
    end
end

end