function [omega_cat, omega_plus, omega_minus, Delta_cat, Delta_ex, Omega_R] = ...
    catalyzon_dispersion(k_range, A_range, J, qvc)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  CATALYZON_DISPERSION  Compute the catalyzon quasiparticle dispersion
%
%  The catalyzon is the emergent quasiparticle in the QVC framework,
%  arising from the hybridization of three boson-mediated channels
%  (Josephson, polariton, hopfion) with the excitonic BEC vacuum.
%
%  The excitonic gap is modified by ZPF coupling:
%    Delta_ex^QVC(A) = sqrt(max(0, 4J^2 - 2*J*A*2J))
%
%  Three boson channels create a Rabi splitting:
%    Omega_R = 2*sqrt(g_J^2 + g_pl^2 + g_H^2 + cross-terms)
%
%  The catalyzon gap:
%    Delta_cat = max(0, Delta_ex - Omega_R/2)
%
%  The full dispersion hybridizes Bogoliubov, Josephson-plasma,
%  and polariton modes:
%    omega_cat(k) = effective dressed quasiparticle energy
%
%  Inputs:
%    k_range  - wavevector array [1/xi units]
%    A_range  - ZPF amplitude array [normalized]
%    J        - tunneling energy [Joules]
%    qvc      - struct with boson coupling parameters
%
%  Outputs:
%    omega_cat   - catalyzon dispersion [J units], size (Nk, NA)
%    omega_plus  - in-phase Bogoliubov mode
%    omega_minus - out-of-phase Bogoliubov mode
%    Delta_cat   - catalyzon gap [J], size (NA,)
%    Delta_ex    - excitonic gap [J], size (NA,)
%    Omega_R     - Rabi frequency [J], scalar
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Nk = length(k_range);
NA = length(A_range);

% Boson coupling strengths (in Joules)
g_BJ  = qvc.g_BJ * J;
g_Bpl = qvc.g_Bpl * J;
g_BH  = qvc.g_BH * J;

% Rabi splitting from three-channel hybridization
% Including interference cross-terms (bosons are not orthogonal channels)
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 parameters (in units of J)
gn_minus = 1.3;   % repulsive interaction (out-of-phase)
gn_plus  = 0.7;   % attractive interaction (in-phase)

% Preallocate
omega_cat   = zeros(Nk, NA);
omega_plus  = zeros(Nk, NA);
omega_minus = zeros(Nk, NA);
Delta_ex    = zeros(1, NA);
Delta_cat   = zeros(1, NA);

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

    % --- Excitonic gap (modified by ZPF) ---
    % Delta_ex = sqrt(4J^2 - 2J * A_ZPF * (2J))
    % A_ZPF couples to the interlayer field, reducing the excitonic binding
    Delta_ex(ia) = sqrt(max(0, 4*J^2 - 2*J * A * (2*J)));

    % --- Catalyzon gap ---
    Delta_cat(ia) = max(0, Delta_ex(ia) - Omega_R/2);

    Dg = Delta_ex(ia) / J;  % dimensionless gap

    for ik = 1:Nk
        k = k_range(ik);
        eps_k = k^2 / 4;  % free-particle dispersion (units of J)

        % --- Bogoliubov branches ---
        % omega_+(k): in-phase (amplitude) mode
        omega_plus(ik,ia) = J * sqrt(max(0, eps_k * (eps_k + 2*gn_plus)));

        % omega_-(k): out-of-phase (phase) mode, gapped by Delta
        omega_minus(ik,ia) = J * sqrt(max(0, eps_k*(eps_k + 2*gn_minus) + Dg^2/4));

        % --- Josephson plasma mode ---
        % omega_J(k) = sqrt(omega_J0^2 + v_J^2 k^2)
        omega_J0 = Dg;
        v_J = 0.4;
        omega_J = J * sqrt(omega_J0^2 + v_J^2 * k^2);

        % --- Polariton mode ---
        % omega_pl(k) = omega_0 * sqrt(k/k_0 + epsilon)
        omega_pl = J * max(Dg, 0.01) * sqrt(k/1.5 + 0.001);

        % --- Catalyzon: three-mode hybridization ---
        % The catalyzon emerges as the dressed quasiparticle from
        % avoided crossings of omega_-, omega_J, and omega_pl
        % with Rabi couplings g_BJ, g_Bpl, g_BH
        %
        % Simplified: weighted average minus Rabi shift
        w_m = omega_minus(ik,ia);
        w_J = omega_J;
        w_pl = omega_pl;
        OR_half = Omega_R / 2;

        % Level repulsion formula (two-level avoided crossing extended to 3)
        omega_avg = (w_m + w_J + w_pl) / 3;
        omega_cat(ik,ia) = max(0.001*J, omega_avg - OR_half + ...
            Delta_cat(ia)/3 * (1 - exp(-k)));
    end
end

end