%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  QVC-TEVC QUASIPARTICLE TUNNELING SIMULATION SUITE
%  ==================================================
%  Quantum Vacuum Catalyzation in Twisted Engineered Vacuum Crystals
%
%  Models:
%    1. Moire superlattice potential from twisted vdW bilayer
%    2. Catalyzon quasiparticle dispersion omega_cat(k, A_ZPF)
%    3. Transfer-matrix tunneling through TEVC barrier landscape
%    4. ZPF-enhanced tunneling (vacuum catalyzation of barrier)
%    5. Schwinger-analog pair production in reduced-gap metamaterial
%    6. Hopfion topological excitation nucleation
%
%  Material system: TMD heterobilayer (MoSe2/WSe2 or similar)
%  Framework:  Mineral Hoiland, QVC Theory (2024-2026)
%
%  Usage: Run this script. All figures are generated sequentially.
%         Adjust parameters in Section 1 below.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear; close all; clc;
fprintf('=== QVC-TEVC Quasiparticle Tunneling Simulation ===\n\n');

% --- Version check ---
v = ver('MATLAB');
if ~isempty(v)
    rel = v.Release;
    fprintf('MATLAB %s %s\n\n', v.Version, rel);
end

%% ========================================================================
%  SECTION 1: PHYSICAL PARAMETERS
%  ========================================================================

% --- Fundamental constants (SI) ---
const.hbar  = 1.0546e-34;   % reduced Planck constant [J.s]
const.e     = 1.6022e-19;   % elementary charge [C]
const.me    = 9.1094e-31;   % electron mass [kg]
const.kB    = 1.3806e-23;   % Boltzmann constant [J/K]
const.c     = 2.9979e8;     % speed of light [m/s]
const.eps0  = 8.8542e-12;   % vacuum permittivity [F/m]

% --- TMD bilayer material parameters ---
mat.a       = 0.328e-9;     % lattice constant [m] (MoSe2)
mat.mstar   = 0.45*const.me;% exciton effective mass [kg]
mat.epsilon = 4.5;          % effective dielectric constant
mat.d_layer = 0.65e-9;      % interlayer separation [m]

% --- TEVC lattice parameters ---
tevc.theta_deg = 1.1;       % twist angle [degrees]
tevc.theta  = tevc.theta_deg * pi/180;  % [radians]
tevc.lambda_M = mat.a / tevc.theta;     % moire period [m]
tevc.N_layers = 3;          % number of stacked layers

% --- Energy scales (in meV, then convert to Joules) ---
meV = const.e * 1e-3;       % 1 meV in Joules
tevc.J      = 15.0 * meV;   % interlayer tunneling [J]  (typ 5-50 meV)
tevc.V0     = 25.0 * meV;   % moire potential amplitude [J]
tevc.V1     =  5.0 * meV;   % second harmonic of moire potential [J]

% --- QVC drive parameters ---
qvc.A_ZPF_range = linspace(0, 0.58, 200);  % normalized ZPF amplitude
qvc.A_ZPF_default = 0.20;   % default A_ZPF (in units of 2J)
qvc.omega0_over_2J = 1.0;   % drive frequency ratio

% --- Boson coupling strengths (in units of J) ---
qvc.g_BJ  = 0.15;           % Josephson boson coupling
qvc.g_Bpl = 0.12;           % polariton boson coupling
qvc.g_BH  = 0.10;           % hopfion boson coupling

% --- Topology ---
topo.QH      = 1;           % Hopf charge
topo.aH_over_xi = 1.5;      % hopfion core size / coherence length
topo.alpha_R = 0.50;         % Rashba SOC (in units of J*a)
topo.Chern   = 1;            % Chern number |C|

% --- Temperature ---
thermo.T_range = linspace(0.01, 5.0, 200);  % temperature [K]
thermo.T_default = 0.08;    % default T (in units of J/kB)

% --- Coherence length ---
xi = const.hbar / sqrt(2 * mat.mstar * tevc.J);  % [m]
fprintf('Moire period:       lambda_M = %.1f nm\n', tevc.lambda_M*1e9);
fprintf('Coherence length:   xi       = %.2f nm\n', xi*1e9);
fprintf('Tunneling energy:   J        = %.1f meV\n', tevc.J/meV);
fprintf('Moire potential:    V0       = %.1f meV\n', tevc.V0/meV);
fprintf('\n');

%% ========================================================================
%  SECTION 2: MOIRE POTENTIAL LANDSCAPE
%  ========================================================================
fprintf('--- Computing moire potential landscape ---\n');

Nx = 300; Ny = 300;
x_range = linspace(-2*tevc.lambda_M, 2*tevc.lambda_M, Nx);
y_range = linspace(-2*tevc.lambda_M, 2*tevc.lambda_M, Ny);
[X, Y] = meshgrid(x_range, y_range);

[V_moire, V_type] = compute_moire_potential(X, Y, tevc.lambda_M, ...
    tevc.V0, tevc.V1, tevc.theta);

% --- Figure 1: Moire potential landscape ---
figure('Name','TEVC Moire Potential','Color','k','Position',[50 500 900 700]);
set(gcf, 'InvertHardcopy', 'off');

subplot(2,2,1);
imagesc(x_range*1e9, y_range*1e9, V_moire/meV);
axis equal tight; colormap(gca, parula);
cb = colorbar; ylabel(cb, 'V_{moire} (meV)');
xlabel('x (nm)'); ylabel('y (nm)');
title('Moire Superlattice Potential', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');

subplot(2,2,2);
imagesc(x_range*1e9, y_range*1e9, V_type);
axis equal tight; colormap(gca, hot);
cb = colorbar; ylabel(cb, 'Stacking type');
xlabel('x (nm)'); ylabel('y (nm)');
title('AA/AB/BA Stacking Regions', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');

subplot(2,2,3);
% 1D slice through moire potential
ix_mid = round(Ny/2);
V_slice = V_moire(ix_mid, :);
plot(x_range*1e9, V_slice/meV, 'c', 'LineWidth', 1.5);
hold on;
% Mark AA sites (minima) — manual peak detection (no toolbox needed)
aa_locs = [];
for ip = 2:length(V_slice)-1
    if V_slice(ip) < V_slice(ip-1) && V_slice(ip) < V_slice(ip+1)
        aa_locs(end+1) = ip; %#ok<SAGROW>
    end
end
plot(x_range(aa_locs)*1e9, V_slice(aa_locs)/meV, 'rv', 'MarkerSize',8, 'MarkerFaceColor','r');
xlabel('x (nm)'); ylabel('V (meV)');
title('1D Slice: Tunneling Path', 'Color','w');
legend('V_{moire}(x)', 'AA nodes', 'TextColor','w', 'Color',[0.1 0.1 0.1]);
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

subplot(2,2,4);
% 3D surface of moire potential
surf(X*1e9, Y*1e9, V_moire/meV, 'EdgeColor','none', 'FaceAlpha',0.85);
colormap(gca, parula);
xlabel('x (nm)'); ylabel('y (nm)'); zlabel('V (meV)');
title('3D Moire Potential Surface', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
view(35, 30); lighting gouraud; camlight headlight;

sgtitle('TEVC Moire Superlattice  (\theta = 1.1°, \lambda_M = 17.1 nm)', ...
    'Color','w', 'FontSize', 14);

%% ========================================================================
%  SECTION 3: CATALYZON DISPERSION
%  ========================================================================
fprintf('--- Computing catalyzon dispersion ---\n');

k_range = linspace(0, 5, 500);   % k in units of 1/xi
A_range = qvc.A_ZPF_range;       % A_ZPF normalized

% Compute dispersion surfaces
[omega_cat, omega_plus, omega_minus, Delta_cat, Delta_ex, Omega_R] = ...
    catalyzon_dispersion(k_range, A_range, tevc.J, qvc);

% --- Figure 2: Catalyzon dispersion ---
figure('Name','Catalyzon Dispersion','Color','k','Position',[100 450 1100 750]);
set(gcf, 'InvertHardcopy', 'off');

subplot(2,3,1);
% omega_cat(k) for several A_ZPF values
A_vals = [0.0, 0.10, 0.20, 0.35, 0.50];
colors_A = [0.3 0.5 1.0; 0.4 0.7 0.9; 0.6 0.5 0.9; 0.9 0.4 0.5; 1.0 0.2 0.2];
hold on;
for ia = 1:length(A_vals)
    [~,idx] = min(abs(A_range - A_vals(ia)));
    plot(k_range, omega_cat(:,idx)/tevc.J, ...
        'Color', colors_A(ia,:), 'LineWidth', 1.5);
end
xlabel('k\xi'); ylabel('\omega_{cat} / J');
title('Catalyzon Dispersion \omega_{cat}(k)', 'Color','w');
legend(arrayfun(@(a) sprintf('A_{ZPF}=%.2f',a), A_vals, 'Uni',0), ...
    'TextColor','w', 'Color',[0.1 0.1 0.1], 'Location','northwest');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);
ylim([0, max(omega_cat(:,1)/tevc.J)*1.1]);

subplot(2,3,2);
% Gap Delta_cat vs A_ZPF
plot(A_range, Delta_cat/meV, 'm', 'LineWidth', 2);
hold on;
plot(A_range, Delta_ex/meV, 'c--', 'LineWidth', 1.2);
yline(Omega_R/(2*meV), 'r:', 'LineWidth', 1.2);
xlabel('A_{ZPF}'); ylabel('Energy (meV)');
title('Gap Collapse: \Delta_{cat} vs A_{ZPF}', 'Color','w');
legend('\Delta_{cat}', '\Delta_{ex}^{QVC}', '\Omega_R/2', ...
    'TextColor','w', 'Color',[0.1 0.1 0.1]);
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

subplot(2,3,3);
% Effective mass ratio
mstar_cat = zeros(size(A_range));
for ia = 1:length(A_range)
    % m*_cat from curvature at k=0: 1/m* = d^2 omega / dk^2
    if length(k_range) > 2
        dk = k_range(2) - k_range(1);
        d2w = (omega_cat(3,ia) - 2*omega_cat(2,ia) + omega_cat(1,ia)) / dk^2;
        if d2w > 0
            mstar_cat(ia) = tevc.J / (2*d2w);
        else
            mstar_cat(ia) = NaN;
        end
    end
end
plot(A_range, mstar_cat, 'Color', [0.5 1.0 0.4], 'LineWidth', 1.5);
xlabel('A_{ZPF}'); ylabel('m^*_{cat} / m^*');
title('Catalyzon Effective Mass', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

subplot(2,3,4);
% 2D dispersion surface omega_cat(k, A_ZPF)
[K_grid, A_grid] = meshgrid(k_range, A_range);
surf(K_grid, A_grid, omega_cat'/tevc.J, 'EdgeColor','none', 'FaceAlpha',0.85);
colormap(gca, hot);
xlabel('k\xi'); ylabel('A_{ZPF}'); zlabel('\omega_{cat}/J');
title('Dispersion Surface \omega_{cat}(k, A_{ZPF})', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
view(35, 25);

subplot(2,3,5);
% Bogoliubov comparison: omega_+ and omega_- at A_ZPF = 0.20
[~, idx20] = min(abs(A_range - 0.20));
plot(k_range, omega_plus(:,idx20)/tevc.J, 'Color',[0.1 0.8 0.6], 'LineWidth',1.5);
hold on;
plot(k_range, omega_minus(:,idx20)/tevc.J, 'Color',[0.3 0.4 1.0], 'LineWidth',1.5);
plot(k_range, omega_cat(:,idx20)/tevc.J, 'm', 'LineWidth',2);
xlabel('k\xi'); ylabel('\omega / J');
title('Bogoliubov Branches at A_{ZPF}=0.20', 'Color','w');
legend('\omega_+(k)', '\omega_-(k)', '\omega_{cat}(k)', ...
    'TextColor','w', 'Color',[0.1 0.1 0.1]);
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

subplot(2,3,6);
% Spectral weight / vacuum entanglement |v_k|^2
gn = 1.3;  % interaction strength (units of J)
v_k_sq = zeros(length(k_range), 1);
for ik = 1:length(k_range)
    eps_k = k_range(ik)^2 / 4;
    Ek = sqrt(eps_k * (eps_k + 2*gn));
    if Ek > 1e-12
        v_k_sq(ik) = 0.5 * ((eps_k + gn) / Ek - 1);
    end
end
plot(k_range, v_k_sq, 'Color', [0.9 0.4 0.6], 'LineWidth', 1.5);
xlabel('k\xi'); ylabel('|v_k|^2');
title('Vacuum Entanglement |v_k^-|^2', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

sgtitle('Catalyzon Quasiparticle Dispersion & Gap Structure', ...
    'Color','w', 'FontSize', 14);

%% ========================================================================
%  SECTION 4: TRANSFER-MATRIX TUNNELING THROUGH MOIRE BARRIERS
%  ========================================================================
fprintf('--- Computing transfer-matrix tunneling ---\n');

% Energy range for tunneling (in units of V0)
E_range = linspace(0.01, 2.5, 600) * tevc.V0;  % [J]

% --- Case A: Bare moire barrier (no ZPF) ---
N_barriers = 3;  % number of moire periods to tunnel through
T_bare = transfer_matrix_1d(E_range, tevc.V0, 0, tevc.lambda_M, ...
    N_barriers, mat.mstar, const.hbar);

% --- Case B: ZPF-enhanced tunneling at several A_ZPF ---
A_zpf_vals = [0.0, 0.10, 0.20, 0.35, 0.50];
T_zpf = zeros(length(E_range), length(A_zpf_vals));
for ia = 1:length(A_zpf_vals)
    % ZPF reduces effective barrier: V_eff = V0 * (1 - A_ZPF * f(x))
    % At AA nodes, f(x) -> max => barrier is lowered
    V_eff = tevc.V0 * (1 - A_zpf_vals(ia));
    T_zpf(:,ia) = transfer_matrix_1d(E_range, V_eff, tevc.V1*(1-A_zpf_vals(ia)), ...
        tevc.lambda_M, N_barriers, mat.mstar, const.hbar);
end

% --- Case C: Tunneling vs number of moire periods ---
N_range = 1:8;
T_vs_N = zeros(length(E_range), length(N_range));
for iN = 1:length(N_range)
    T_vs_N(:,iN) = transfer_matrix_1d(E_range, tevc.V0, tevc.V1, ...
        tevc.lambda_M, N_range(iN), mat.mstar, const.hbar);
end

% --- Figure 3: Tunneling Transmission ---
figure('Name','Transfer Matrix Tunneling','Color','k','Position',[150 400 1100 750]);
set(gcf, 'InvertHardcopy', 'off');

subplot(2,3,1);
semilogy(E_range/tevc.V0, T_bare, 'c', 'LineWidth', 1.5);
xlabel('E / V_0'); ylabel('T(E)');
title('Bare Moire Tunneling (N=3)', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);
ylim([1e-15, 1.5]);

subplot(2,3,2);
hold on;
for ia = 1:length(A_zpf_vals)
    semilogy(E_range/tevc.V0, T_zpf(:,ia), 'Color', colors_A(ia,:), 'LineWidth', 1.3);
end
xlabel('E / V_0'); ylabel('T(E)');
title('ZPF-Enhanced Tunneling', 'Color','w');
legend(arrayfun(@(a) sprintf('A_{ZPF}=%.2f',a), A_zpf_vals, 'Uni',0), ...
    'TextColor','w', 'Color',[0.1 0.1 0.1], 'Location','southeast');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);
ylim([1e-15, 1.5]);

subplot(2,3,3);
% Enhancement ratio T(A_ZPF) / T(0)
hold on;
for ia = 2:length(A_zpf_vals)
    ratio = T_zpf(:,ia) ./ max(T_zpf(:,1), 1e-30);
    ratio(ratio > 1e15) = NaN;
    semilogy(E_range/tevc.V0, ratio, 'Color', colors_A(ia,:), 'LineWidth', 1.3);
end
xlabel('E / V_0'); ylabel('T(A_{ZPF}) / T(0)');
title('Tunneling Enhancement Factor', 'Color','w');
legend(arrayfun(@(a) sprintf('A_{ZPF}=%.2f',a), A_zpf_vals(2:end), 'Uni',0), ...
    'TextColor','w', 'Color',[0.1 0.1 0.1], 'Location','northeast');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

subplot(2,3,4);
% Tunneling vs N (barrier count)
hold on;
N_colors = cool(length(N_range));
for iN = 1:length(N_range)
    semilogy(E_range/tevc.V0, T_vs_N(:,iN), 'Color', N_colors(iN,:), 'LineWidth', 1.2);
end
xlabel('E / V_0'); ylabel('T(E)');
title('Tunneling vs Barrier Count N', 'Color','w');
legend(arrayfun(@(n) sprintf('N=%d',n), N_range, 'Uni',0), ...
    'TextColor','w', 'Color',[0.1 0.1 0.1], 'Location','southeast');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);
ylim([1e-30, 1.5]);

subplot(2,3,5);
% Tunneling current J(V) = integral T(E) [f_L(E) - f_R(E)] dE
V_bias = linspace(0, 50, 200) * meV;  % bias voltage [J]
T_eff = 0.5;  % temperature in K
J_current = zeros(length(V_bias), length(A_zpf_vals));
for iv = 1:length(V_bias)
    for ia = 1:length(A_zpf_vals)
        % Landauer-type integral
        f_L = 1 ./ (1 + exp((E_range - V_bias(iv)/2) / (const.kB*T_eff)));
        f_R = 1 ./ (1 + exp((E_range + V_bias(iv)/2) / (const.kB*T_eff)));
        dE = E_range(2) - E_range(1);
        J_current(iv,ia) = sum(T_zpf(:,ia)' .* (f_L - f_R)) * dE;
    end
end
J_current = J_current / max(abs(J_current(:)));  % normalize
hold on;
for ia = 1:length(A_zpf_vals)
    plot(V_bias/meV, J_current(:,ia), 'Color', colors_A(ia,:), 'LineWidth', 1.3);
end
xlabel('V_{bias} (meV)'); ylabel('J / J_{max}');
title('Tunneling Current J(V)', 'Color','w');
legend(arrayfun(@(a) sprintf('A_{ZPF}=%.2f',a), A_zpf_vals, 'Uni',0), ...
    'TextColor','w', 'Color',[0.1 0.1 0.1], 'Location','northwest');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

subplot(2,3,6);
% Tunneling probability at E = V0/2 as function of A_ZPF (continuous)
A_fine = linspace(0, 0.58, 500);
T_at_half = zeros(size(A_fine));
E_probe = tevc.V0 * 0.5;  % probe at half-barrier
for ia = 1:length(A_fine)
    V_eff = tevc.V0 * (1 - A_fine(ia));
    T_temp = transfer_matrix_1d(E_probe, V_eff, tevc.V1*(1-A_fine(ia)), ...
        tevc.lambda_M, 3, mat.mstar, const.hbar);
    T_at_half(ia) = T_temp;
end
semilogy(A_fine, T_at_half, 'Color', [0.8 0.5 1.0], 'LineWidth', 2);
xlabel('A_{ZPF}'); ylabel('T(E = V_0/2)');
title('Vacuum Catalyzation of Tunneling', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

sgtitle('Transfer-Matrix Tunneling Through TEVC Moire Barriers', ...
    'Color','w', 'FontSize', 14);

%% ========================================================================
%  SECTION 5: SCHWINGER-ANALOG PAIR PRODUCTION
%  ========================================================================
fprintf('--- Computing Schwinger pair production rates ---\n');

% In a narrow-gap system, the Schwinger critical field maps to:
%   E_cr^eff = Delta_cat^2 / (e * xi)
% where Delta_cat plays the role of 2mc^2 and xi is the coherence length.
%
% The production rate per unit volume per unit time:
%   Gamma = (Delta_cat^2 / 4 pi^3 xi^2 hbar) * (E/E_cr) * exp(-pi E_cr / E)
%
% The "applied field" E_eff comes from the interlayer electric field
% at AA-stacking nodes, enhanced by dielectric mismatch.

E_field_range = linspace(0.01, 5.0, 400);  % E / E_cr

[Gamma_sch, E_cr_eff, Gamma_vs_A] = schwinger_rate(A_range, E_field_range, ...
    tevc.J, qvc, xi, const);

% --- Figure 4: Schwinger pair production ---
figure('Name','Schwinger Pair Production','Color','k','Position',[200 350 1100 700]);
set(gcf, 'InvertHardcopy', 'off');

subplot(2,3,1);
% Schwinger rate vs E/E_cr at fixed A_ZPF
A_sch_vals = [0.05, 0.15, 0.25, 0.40, 0.52];
colors_sch = [0.3 0.4 1.0; 0.3 0.7 0.9; 0.7 0.5 0.9; 1.0 0.5 0.3; 1.0 0.2 0.2];
hold on;
for ia = 1:length(A_sch_vals)
    [~,idx] = min(abs(A_range - A_sch_vals(ia)));
    semilogy(E_field_range, Gamma_sch(:,idx), 'Color', colors_sch(ia,:), 'LineWidth',1.3);
end
xlabel('E / E_{cr}^{eff}'); ylabel('\Gamma (arb. units)');
title('Schwinger Rate vs Field', 'Color','w');
legend(arrayfun(@(a) sprintf('A_{ZPF}=%.2f',a), A_sch_vals, 'Uni',0), ...
    'TextColor','w', 'Color',[0.1 0.1 0.1], 'Location','northwest');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

subplot(2,3,2);
% Critical field E_cr vs A_ZPF
semilogy(A_range, E_cr_eff/1e6, 'Color', [1.0 0.4 0.3], 'LineWidth', 2);
xlabel('A_{ZPF}'); ylabel('E_{cr}^{eff} (MV/m)');
title('Effective Schwinger Threshold', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

subplot(2,3,3);
% Schwinger enhancement at FIXED absolute field E_abs
% As A_ZPF increases, E_cr drops, so E_abs/E_cr grows => rate explodes
E_abs_fixed = E_cr_eff(1) * 0.3;  % 30% of threshold at A=0
enhancement = zeros(size(A_range));
for ia = 1:length(A_range)
    Ecr_local = max(E_cr_eff(ia), 1e-10);
    ratio = E_abs_fixed / Ecr_local;
    if ratio > 0.01
        enhancement(ia) = ratio * exp(-pi / ratio);
    else
        enhancement(ia) = 1e-50;
    end
end
enhancement = enhancement / max(enhancement(1), 1e-100);
semilogy(A_range, enhancement, 'Color', [0.9 0.6 0.2], 'LineWidth', 2);
xlabel('A_{ZPF}'); ylabel('\Gamma(A) / \Gamma(0)');
title('Schwinger Enhancement (fixed |E|)', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

subplot(2,3,4);
% Phase diagram: log10(Gamma) in (E/E_cr, A_ZPF) space
[E_grid, A_grid2] = meshgrid(E_field_range, A_range);
log_Gamma = log10(max(Gamma_vs_A', 1e-50));
imagesc(E_field_range, A_range, log_Gamma);
axis xy; colormap(gca, hot);
cb = colorbar; ylabel(cb, 'log_{10} \Gamma');
xlabel('E / E_{cr}^{eff}'); ylabel('A_{ZPF}');
title('Pair Production Phase Diagram', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
% Mark the QCP line
hold on;
A_qcp = 0.50;  % approximate QCP
yline(A_qcp, 'w--', 'LineWidth', 1.5);
text(0.5, A_qcp+0.02, 'QCP', 'Color','w', 'FontSize',10);

subplot(2,3,5);
% Compare QED Schwinger vs TEVC-analog Schwinger
% QED: E_cr = m_e^2 c^3 / (e hbar) ~ 1.3e18 V/m
E_cr_QED = const.me^2 * const.c^3 / (const.e * const.hbar);
% TEVC analog with various gaps
gap_meV = linspace(0.1, 30, 200);
gap_J = gap_meV * meV;
E_cr_analog = gap_J.^2 ./ (const.e * xi);
semilogy(gap_meV, E_cr_analog/1e6, 'Color', [0.5 0.8 1.0], 'LineWidth', 2);
hold on;
semilogy(gap_meV, ones(size(gap_meV))*E_cr_QED/1e6, 'r--', 'LineWidth', 1.2);
xlabel('\Delta_{cat} (meV)'); ylabel('E_{cr} (MV/m)');
title('TEVC vs QED Schwinger Threshold', 'Color','w');
legend('TEVC E_{cr}^{eff}(\Delta)', 'QED E_{cr} = 1.3\times10^{18} V/m', ...
    'TextColor','w', 'Color',[0.1 0.1 0.1], 'Location','southeast');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);
ylim([1e-2, 1e14]);

subplot(2,3,6);
% Pair production rate vs temperature (thermal-assisted Schwinger)
T_K = linspace(0.01, 5, 300);
E_over_Ecr = 0.3;  % sub-threshold field
Delta_fixed = Delta_cat(100);  % pick a representative gap
Gamma_thermal = zeros(size(T_K));
for iT = 1:length(T_K)
    kBT = const.kB * T_K(iT);
    % Thermal enhancement: exp(-pi (E_cr/E - kBT/Delta))
    exponent = pi * (1/E_over_Ecr - kBT/max(Delta_fixed, 1e-25));
    Gamma_thermal(iT) = exp(-max(exponent, 0));
end
semilogy(T_K, Gamma_thermal, 'Color', [1.0 0.5 0.7], 'LineWidth', 2);
xlabel('Temperature (K)'); ylabel('\Gamma_{thermal}');
title('Thermal-Assisted Schwinger (E/E_{cr}=0.3)', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

sgtitle('Schwinger-Analog Pair Production in TEVC', ...
    'Color','w', 'FontSize', 14);

%% ========================================================================
%  SECTION 6: HOPFION NUCLEATION & TOPOLOGY
%  ========================================================================
fprintf('--- Computing hopfion nucleation energetics ---\n');

[E_hopf, Gamma_hopf, rho_s_arr, T_BKT_arr, hopf_phase] = ...
    hopfion_nucleation(A_range, thermo.T_range, tevc.J, qvc, topo, const);

% --- Figure 5: Hopfion topology ---
figure('Name','Hopfion Nucleation','Color','k','Position',[250 300 1100 750]);
set(gcf, 'InvertHardcopy', 'off');

subplot(2,3,1);
% Hopfion energy vs A_ZPF for different Q_H
QH_vals = [1, 2, 3];
colors_QH = [0.6 0.5 0.9; 0.9 0.4 0.6; 1.0 0.7 0.2];
hold on;
for iq = 1:length(QH_vals)
    E_h = E_hopf(:, iq);
    plot(A_range, E_h/meV, 'Color', colors_QH(iq,:), 'LineWidth', 1.8);
end
xlabel('A_{ZPF}'); ylabel('E_{hopfion} (meV)');
title('Hopfion Energy vs ZPF Amplitude', 'Color','w');
legend(arrayfun(@(q) sprintf('Q_H = %d',q), QH_vals, 'Uni',0), ...
    'TextColor','w', 'Color',[0.1 0.1 0.1]);
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

subplot(2,3,2);
% Superfluid stiffness rho_s vs A_ZPF
plot(A_range, rho_s_arr, 'Color', [0.1 0.8 0.6], 'LineWidth', 2);
xlabel('A_{ZPF}'); ylabel('\rho_s / \rho_0');
title('Superfluid Stiffness (BEC Order)', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);
ylim([0, 1.1]);

subplot(2,3,3);
% BKT transition temperature vs A_ZPF
plot(A_range, T_BKT_arr, 'Color', [0.3 0.9 0.4], 'LineWidth', 2);
xlabel('A_{ZPF}'); ylabel('T_{BKT} (arb.)');
title('BKT Transition Temperature', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

subplot(2,3,4);
% Hopfion nucleation rate (thermal) vs T at A_ZPF = 0.20
[~, idx_a20] = min(abs(A_range - 0.20));
semilogy(thermo.T_range, Gamma_hopf(:, idx_a20), 'Color', [0.8 0.5 1.0], 'LineWidth', 2);
xlabel('T (arb.)'); ylabel('\Gamma_{nucleation}');
title('Hopfion Nucleation Rate (A_{ZPF}=0.20)', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

subplot(2,3,5);
% Phase diagram: hopfion stability in (A_ZPF, T) plane
imagesc(A_range, thermo.T_range, hopf_phase);
axis xy;
colormap(gca, [0.05 0.05 0.15; 0.2 0.15 0.5; 0.5 0.3 0.8; 1.0 0.2 0.2]);
caxis([-0.5 3.5]);
cb = colorbar; ylabel(cb, 'Phase');
set(cb, 'Ticks', [0, 1, 2, 3], ...
    'TickLabels', {'Stable BEC+Hopfion', 'Hopfion melting', 'BKT proliferation', 'Normal'});
xlabel('A_{ZPF}'); ylabel('T (arb.)');
title('Hopfion-BEC Phase Diagram', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');

subplot(2,3,6);
% Hopf fibration visualization: linking number
% Show Q_H = 1 preimage circles on S^2
phi_range = linspace(0, 2*pi, 200);
theta_fib = linspace(0, pi, 100);
[PHI, THETA] = meshgrid(phi_range, theta_fib);
% Stereographic projection of Hopf fiber
R_fib = topo.aH_over_xi * 1.0;
r_fib = topo.aH_over_xi * 0.3;
X_torus = (R_fib + r_fib*cos(THETA)) .* cos(PHI);
Y_torus = (R_fib + r_fib*cos(THETA)) .* sin(PHI);
Z_torus = r_fib * sin(THETA);
% Color by Hopf phase
C_hopf = mod(PHI + THETA, 2*pi) / (2*pi);
surf(X_torus, Y_torus, Z_torus, C_hopf, 'EdgeColor','none', 'FaceAlpha', 0.7);
colormap(gca, hsv);
axis equal; view(30, 25);
xlabel('x/\xi'); ylabel('y/\xi'); zlabel('z/\xi');
title(sprintf('Hopfion Texture (Q_H = %d)', topo.QH), 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
lighting gouraud; camlight headlight;

sgtitle('Hopfion Topological Excitations in TEVC', ...
    'Color','w', 'FontSize', 14);

%% ========================================================================
%  SECTION 7: ZPF VACUUM CATALYZATION — UNIFIED VIEW
%  ========================================================================
fprintf('--- Computing unified ZPF catalyzation effects ---\n');

[zpf_data] = zpf_vacuum_catalyzation(A_range, k_range, E_range, ...
    tevc, qvc, topo, mat, const, meV);

% --- Figure 6: Unified ZPF Effects ---
figure('Name','ZPF Vacuum Catalyzation','Color','k','Position',[300 250 1200 800]);
set(gcf, 'InvertHardcopy', 'off');

subplot(2,4,1);
plot(A_range, zpf_data.Delta_cat/meV, 'Color', [0.6 0.4 0.9], 'LineWidth', 2);
hold on;
plot(A_range, zpf_data.Delta_ex/meV, 'c--', 'LineWidth', 1.2);
xlabel('A_{ZPF}'); ylabel('Gap (meV)');
title('\Delta_{cat} Gap', 'Color','w');
legend('\Delta_{cat}','\Delta_{ex}','TextColor','w','Color',[0.1 0.1 0.1]);
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

subplot(2,4,2);
semilogy(A_range, zpf_data.T_tunnel, 'Color', [0.4 0.8 1.0], 'LineWidth', 2);
xlabel('A_{ZPF}'); ylabel('T(E=V_0/2)');
title('Tunneling Probability', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

subplot(2,4,3);
semilogy(A_range, zpf_data.Gamma_schwinger, 'Color', [1.0 0.4 0.3], 'LineWidth', 2);
xlabel('A_{ZPF}'); ylabel('\Gamma_{Schwinger}');
title('Pair Production Rate', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

subplot(2,4,4);
plot(A_range, zpf_data.E_hopfion/meV, 'Color', [0.8 0.5 1.0], 'LineWidth', 2);
xlabel('A_{ZPF}'); ylabel('E_{hopfion} (meV)');
title('Hopfion Energy', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

subplot(2,4,5);
plot(A_range, zpf_data.rho_s, 'Color', [0.1 0.8 0.6], 'LineWidth', 2);
hold on;
plot(A_range, zpf_data.T_BKT, 'Color', [0.3 0.9 0.4], 'LineWidth', 1.5);
xlabel('A_{ZPF}'); ylabel('Order parameter');
title('\rho_s and T_{BKT}', 'Color','w');
legend('\rho_s','T_{BKT}','TextColor','w','Color',[0.1 0.1 0.1]);
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

subplot(2,4,6);
% Negentropy: S_vac = -Tr(rho_red ln rho_red) for Bogoliubov vacuum
% For each k mode: S_k = -|v_k|^2 ln|v_k|^2 - (1+|v_k|^2)ln(1+|v_k|^2) + ...
% Simplified: total entanglement entropy
plot(A_range, zpf_data.S_entanglement, 'Color', [0.9 0.6 0.2], 'LineWidth', 2);
xlabel('A_{ZPF}'); ylabel('S_{entanglement}');
title('Vacuum Entanglement Entropy', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

subplot(2,4,7);
% Coherence length vs A_ZPF
plot(A_range, zpf_data.xi_eff*1e9, 'Color', [0.5 1.0 0.8], 'LineWidth', 2);
xlabel('A_{ZPF}'); ylabel('\xi_{eff} (nm)');
title('Effective Coherence Length', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

subplot(2,4,8);
% Combined "catalyzation figure of merit"
% FOM = T_tunnel * Gamma_schwinger / Delta_cat (higher = more catalyzed)
FOM = zpf_data.T_tunnel .* zpf_data.Gamma_schwinger ./ max(zpf_data.Delta_cat, 1e-30);
FOM = FOM / max(FOM);
semilogy(A_range, FOM, 'Color', [1.0 0.8 0.2], 'LineWidth', 2);
xlabel('A_{ZPF}'); ylabel('FOM (normalized)');
title('Catalyzation Figure of Merit', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

sgtitle('Unified ZPF Vacuum Catalyzation Effects in TEVC', ...
    'Color','w', 'FontSize', 14);

%% ========================================================================
%  SECTION 8: TWIST ANGLE DEPENDENCE
%  ========================================================================
fprintf('--- Computing twist angle dependence ---\n');

theta_range_deg = linspace(0.3, 5.0, 300);
theta_range_rad = theta_range_deg * pi/180;
lambda_M_arr = mat.a ./ theta_range_rad;

% For each theta, compute key observables at fixed A_ZPF = 0.25
A_fixed = 0.25;
Delta_vs_theta = zeros(size(theta_range_deg));
T_tunnel_vs_theta = zeros(size(theta_range_deg));
E_hopf_vs_theta = zeros(size(theta_range_deg));

for it = 1:length(theta_range_deg)
    lam = lambda_M_arr(it);
    % Tunneling J depends on theta (flat band enhancement near magic angle)
    J_eff = tevc.J * (1 + 2*exp(-(theta_range_deg(it) - 1.1)^2 / 0.1));

    % Gap
    D_ex = sqrt(max(0, 4*J_eff^2 - 2*J_eff*A_fixed*2*J_eff));
    OR2 = compute_rabi(J_eff, qvc);
    Delta_vs_theta(it) = max(0, D_ex - OR2) / meV;

    % Tunneling through N=3 periods
    V_eff = tevc.V0 * (1 - A_fixed);
    E_probe = tevc.V0 * 0.5;
    T_temp = transfer_matrix_1d(E_probe, V_eff, tevc.V1*(1-A_fixed), ...
        lam, 3, mat.mstar, const.hbar);
    T_tunnel_vs_theta(it) = T_temp;

    % Hopfion energy: E_H = 4*pi^2 * rho_s * (a_H/xi) * Q_H * J
    rhos_t = max(0, (1 - A_fixed)^2);
    E_hopf_vs_theta(it) = 4*pi^2 * rhos_t * J_eff * topo.aH_over_xi * topo.QH / meV;
end

% --- Figure 7: Twist angle dependence ---
figure('Name','Twist Angle Dependence','Color','k','Position',[350 200 1000 650]);
set(gcf, 'InvertHardcopy', 'off');

subplot(2,2,1);
plot(theta_range_deg, lambda_M_arr*1e9, 'c', 'LineWidth', 1.5);
xlabel('\theta (degrees)'); ylabel('\lambda_M (nm)');
title('Moire Period vs Twist Angle', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

subplot(2,2,2);
plot(theta_range_deg, Delta_vs_theta, 'Color', [0.6 0.4 0.9], 'LineWidth', 2);
xlabel('\theta (degrees)'); ylabel('\Delta_{cat} (meV)');
title('Catalyzon Gap vs Twist Angle', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);
% Mark magic angle
xline(1.1, 'r--', 'Magic', 'Color','r', 'LineWidth', 1.2, 'LabelColor','r');

subplot(2,2,3);
semilogy(theta_range_deg, T_tunnel_vs_theta, 'Color', [0.4 0.8 1.0], 'LineWidth', 2);
xlabel('\theta (degrees)'); ylabel('T(E=V_0/2)');
title('Tunneling vs Twist Angle', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);
xline(1.1, 'r--', 'Magic', 'Color','r', 'LineWidth', 1.2, 'LabelColor','r');

subplot(2,2,4);
plot(theta_range_deg, E_hopf_vs_theta, 'Color', [0.8 0.5 1.0], 'LineWidth', 2);
xlabel('\theta (degrees)'); ylabel('E_{hopfion} (meV)');
title('Hopfion Energy vs Twist Angle', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);
xline(1.1, 'r--', 'Magic', 'Color','r', 'LineWidth', 1.2, 'LabelColor','r');

sgtitle('Twist Angle Dependence (A_{ZPF} = 0.25)', ...
    'Color','w', 'FontSize', 14);

%% ========================================================================
%  SECTION 9: RASHBA SOC + CHIRAL TRANSPORT
%  ========================================================================
fprintf('--- Computing Rashba spin-orbit and chiral effects ---\n');

% Rashba splits the Bogoliubov dispersion into spin branches:
%   omega_pm(k) = omega_-(k) +/- alpha_R * k
% This creates a spin-momentum locked texture enabling topological
% transport and chiral anomaly analogs.

k_fine = linspace(0, 4, 400);
alpha_R_vals = [0, 0.25, 0.50, 1.0, 1.5];
A_fixed2 = 0.20;

% --- Figure 8: Rashba + Chirality ---
figure('Name','Rashba SOC & Chirality','Color','k','Position',[400 150 1000 650]);
set(gcf, 'InvertHardcopy', 'off');

subplot(2,2,1);
hold on;
gn = 1.3;
D_bg = sqrt(max(0, 4 - 2*A_fixed2));  % gap in units of J
for ia = 1:length(alpha_R_vals)
    aR = alpha_R_vals(ia);
    om_base = zeros(size(k_fine));
    om_up = zeros(size(k_fine));
    om_dn = zeros(size(k_fine));
    for ik = 1:length(k_fine)
        eps_k = k_fine(ik)^2 / 4;
        om_base(ik) = sqrt(max(0, eps_k*(eps_k + 2*gn) + D_bg^2/4));
        om_up(ik) = om_base(ik) + aR * k_fine(ik) * 0.2;
        om_dn(ik) = om_base(ik) - aR * k_fine(ik) * 0.2;
    end
    col = colors_A(ia,:);
    plot(k_fine, om_up, '-', 'Color', col, 'LineWidth', 1.3);
    plot(k_fine, om_dn, '--', 'Color', col, 'LineWidth', 1.0);
end
xlabel('k\xi'); ylabel('\omega / J');
title('Rashba-Split Dispersion \omega_{\uparrow\downarrow}(k)', 'Color','w');
legend('\alpha_R=0','',...
    '\alpha_R=0.25','',...
    '\alpha_R=0.50','',...
    '\alpha_R=1.0','',...
    '\alpha_R=1.5','', ...
    'TextColor','w','Color',[0.1 0.1 0.1],'Location','northwest');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

subplot(2,2,2);
% Berry curvature Omega_z(k) for Rashba-coupled Bogoliubov system
% Omega_z = -2 Im(<du/dk_x | du/dk_y>)
% For our model: Omega_z ~ alpha_R * Delta / (2*(epsilon_k^2 + Delta^2)^{3/2})
Omega_z = zeros(length(k_fine), length(alpha_R_vals));
for ia = 1:length(alpha_R_vals)
    aR = alpha_R_vals(ia);
    for ik = 1:length(k_fine)
        eps_k = k_fine(ik)^2/4;
        denom = (eps_k^2 + D_bg^2)^1.5 + 1e-10;
        Omega_z(ik,ia) = aR * D_bg / (2*denom);
    end
end
hold on;
for ia = 1:length(alpha_R_vals)
    plot(k_fine, Omega_z(:,ia), 'Color', colors_A(ia,:), 'LineWidth', 1.3);
end
xlabel('k\xi'); ylabel('\Omega_z(k)');
title('Berry Curvature (Anomalous Hall)', 'Color','w');
legend(arrayfun(@(a) sprintf('\\alpha_R=%.2f',a), alpha_R_vals, 'Uni',0), ...
    'TextColor','w','Color',[0.1 0.1 0.1],'Location','northeast');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
grid on; set(gca, 'GridColor', [0.3 0.3 0.3]);

subplot(2,2,3);
% Chern number integration: C = (1/2pi) int Omega_z dk_x dk_y
% For 2D: C = (1/2pi) * 2pi * int_0^inf Omega_z(k) * k dk
Chern_computed = zeros(1, length(alpha_R_vals));
for ia = 1:length(alpha_R_vals)
    dk = k_fine(2) - k_fine(1);
    integrand = Omega_z(:,ia) .* k_fine';
    Chern_computed(ia) = sum(integrand) * dk;  % factor of 2pi cancels
end
bar(alpha_R_vals, Chern_computed, 0.5, 'FaceColor', [0.6 0.4 0.9], 'EdgeColor', [0.8 0.7 1.0]);
xlabel('\alpha_R'); ylabel('Chern number C');
title('Integrated Chern Number', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');

subplot(2,2,4);
% Spin texture: <sigma_x>, <sigma_y> around a constant-energy contour
% For Rashba: <sigma> perpendicular to k (spin-momentum locking)
n_pts = 200;
k_circ = 1.5;  % constant energy contour
phi = linspace(0, 2*pi, n_pts);
kx = k_circ * cos(phi);
ky = k_circ * sin(phi);
% Spin expectation: sigma_x ~ -sin(phi), sigma_y ~ cos(phi) (Rashba)
sx = -sin(phi) * topo.alpha_R;
sy = cos(phi) * topo.alpha_R;
quiver(kx(1:5:end), ky(1:5:end), sx(1:5:end), sy(1:5:end), 0.4, ...
    'Color', [0.9 0.5 0.9], 'LineWidth', 1.2);
hold on;
plot(kx, ky, 'c--', 'LineWidth', 0.8);
axis equal;
xlabel('k_x\xi'); ylabel('k_y\xi');
title('Spin Texture (Rashba Locking)', 'Color','w');
set(gca, 'Color','k', 'XColor','w', 'YColor','w');
xlim([-2.5 2.5]); ylim([-2.5 2.5]);

sgtitle('Rashba SOC, Berry Curvature & Chiral Transport', ...
    'Color','w', 'FontSize', 14);

%% ========================================================================
%  DONE
%  ========================================================================
fprintf('\n=== Simulation complete. %d figures generated. ===\n', 8);
fprintf('    Moire period:        %.1f nm\n', tevc.lambda_M*1e9);
fprintf('    Coherence length:    %.2f nm\n', xi*1e9);
fprintf('    Catalyzon gap at A=0.20: %.3f meV\n', Delta_cat(100)/meV);
fprintf('    Schwinger E_cr(eff): %.2e V/m (vs QED %.2e V/m)\n', ...
    E_cr_eff(100), E_cr_QED);
fprintf('    Enhancement ratio:   %.0f orders of magnitude\n', ...
    log10(E_cr_QED / max(E_cr_eff(100), 1)));


%% ========================================================================
%  LOCAL HELPER: compute Rabi splitting (used in Section 8)
%  ========================================================================
function OR2 = compute_rabi(J, qvc)
    g_BJ  = qvc.g_BJ * J;
    g_Bpl = qvc.g_Bpl * J;
    g_BH  = qvc.g_BH * J;
    OR = 2*sqrt(g_BJ^2 + g_Bpl^2 + g_BH^2 + g_BJ*g_Bpl + g_BJ*g_BH + g_Bpl*g_BH);
    OR2 = OR / 2;
end