Error in the diagonalization of the Hamiltonian of the Zeeman effect

Hi. I am calculating the energy shifts in the Cs-133 atom due to external magnetic fields by the Zeeman effect (my code is attached). For the ground state and the D1 line I have no problem, the code works fine, the plots give as the D. Steck paper. However, when I try to do it for the D2 line, it seems that the lines of different colors are mixed, which gives incorrect results, as seen in the attached image.

In the code, I construct the Hamiltonian for each magnetic field value and diagonalize the Hamiltonian. One way I got it to work is to diagonalize the Hamiltonian symbolically and then substitute the magnetic field values. However, this method takes too long, and when I tried to add the octupolar term, it does not work. I think it is a problem in the order of the eigenvalues, but I have not been able to solve it. If someone could help me, I would be eternally grateful.

This is the result I obtain.

This is the result of the reference document. In this case, the lines of different colors are not mixed.

The reference document is attached below.

Cesium D Line Data

clc
clear variables
close all

tic

%% Constants definition

line = 'D2'; % [Ground D1 D2]
I = 7/2; % Nuclear spin
g_I = -0.00039885395; % Nuclear g-factor
J = [1/2 3/2]; % Total angular momentum
g_J = [2.002540261 0.665900 1.33408749]; %  Fine structure g-factor
mu_B = (9.27400899E-24)*1E-4; % Bohr Magneton (J/G)
hbar = 1.054E-34; %  Reduced Planck's constant (Js)
Ahfs = 2*pi*hbar*[2.2981579425E9 291.9201E6 50.28827E6]; % Magnetic Dipole Constant
Bhfs = -2*pi*hbar*0.4934E6; % Electric Quadrupole Constant 5^2 P_3/2
Chfs = 2*pi*hbar*0.560E3; % Magnetic Octupole Constant 6^2 P_3/D2
Folder = pwd;

switch line   
    case 'Ground'
        A = Ahfs(1); 
        J = J(1); % Total angular momentum
        g_J = g_J(1);
        BF = linspace(0,15000,1000); % Magnetic field vector (G)   
        titleg = '$ ^{133}\mathrm{Cs} \, \, \, 6^2 \mathrm{S}_{1/2}$';
        yl = '$E/h$ (GHz)';
        Amp = 1E-9;
        yy = [-26 26];
        str = 'Ground state';
    case 'D1'
        A = Ahfs(2);
        J = J(1); % Total angular momentum
        g_J = g_J(2);
        BF = linspace(0,5000,1000); % Magnetic field vector (G)
        titleg = '$ ^{133}\mathrm{Cs} \, \, \, 6^2 \mathrm{P}_{1/2}$';
        yl = '$E/h$ (GHz)';
        Amp = 1E-9;
        yy = [-2.5 2.5];
        str = 'D1 line';
    case 'D2'
        A =  Ahfs(3);
        J = J(2); % Total angular momentum
        g_J =  g_J(3);
        BF = linspace(0,500,1000); % Magnetic field vector (G)
        titleg = '$ ^{133}\mathrm{Cs} \, \, \, 6^2 \mathrm{P}_{3/2}$';
        yl = '$E/h$ (MHz)';
        Amp = 1E-6;
        yy = [-1500 1500];
        str = 'D2 line';
    otherwise
        error('Unexpected value of state.');
end

%% 

% Get the spin operator matrices for I and J

[SxI, SyI, SzI] = espin(I);
[SxJ, SyJ, SzJ] = espin(J);

% Construct the spin operators, defined as the tensor product of the spin matrices S(I) and S(J) with 
% indetity matrices 1(I) and 1(J)

Ix = kron(SxI, id(J));
Iy = kron(SyI, id(J));
Iz = kron(SzI, id(J));
Jx = kron(id(I), SxJ);
Jy = kron(id(I), SyJ);
Jz = kron(id(I), SzJ);

F_values = abs(I-J):(I+J);
degeneracies = 2*F_values + 1; % Degeneracies for each F value

%% We diagonalize the Hamiltonian symbolically and then evaluate the eigenvalues in terms of the external magnetic field
energies = zeros(length(BF), (2*I+1)*(2*J+1)); % Preallocate energy array

% Check the value of J before entering the loop
if J == 1/2
    % Construct Hamiltonian for J = 1/2
    HZeemanBase = A*(Ix*Jx + Iy*Jy + Iz*Jz) + ...
        mu_B*BF(1)*(g_I*Iz + g_J*Jz);  % Use BF(1) just as a placeholder for the structure
elseif J == 3/2
    % Construct Hamiltonian for J = 3/2
    HZeemanBase = A*(Ix*Jx + Iy*Jy + Iz*Jz) + ...
        mu_B*BF(1)*(g_I*Iz + g_J*Jz) + ...
        (Bhfs/(2*I*(2*I-1)*J*(2*J-1)))*( 3*(Ix*Jx + Iy*Jy + Iz*Jz).^2 + (3/2)*(Ix*Jx + Iy*Jy + Iz*Jz) - (Ix.^2 + Iy.^2 + Iz.^2)*(Jx.^2 + Jy.^2 + Jz.^2) );
end

for j = 1:length(BF)
    % Modify the Hamiltonian for each BF(j)
    HZeeman = HZeemanBase;
    % Update the Zeeman term for the current magnetic field strength
    HZeeman = (Amp/(2*pi*hbar))*(HZeeman + mu_B*BF(j)*(g_I*Iz + g_J*Jz));
    
    % Diagonalize the Hamiltonian
    [eigenvectors, eigenvalues] = eig(HZeeman);

    % Store the eigenvalues (energies)
    energies(j, :) = sort(diag(real(eigenvalues))); % Sort energies for clarity
end

% Plot levels for each F group with different colors
colors = {'r', 'b', 'g', 'm', 'k'}; % Colors for each F group
level_start = 1; % Start index for each F group
legend_handles = []; % Store plot handles for the legend

Fig = figure;
hold on;
for f_idx = 1:length(F_values)
    num_levels_in_group = degeneracies(f_idx); % Number of levels for this F
    level_end = level_start + num_levels_in_group - 1; % End index for this F
    h = plot(BF, energies(:, level_start:level_end), 'LineWidth',2, 'Color', colors{f_idx});
    legend_handles = [legend_handles, h(1)]; % Store one handle per group
    level_start = level_end + 1; % Update start index for the next F group
end

set(gca, 'xminortick', 'on', 'TickLabelInterpreter','latex','fontsize',18);
set(gca, 'yminortick', 'on', 'TickLabelInterpreter','latex','fontsize',18);
xlabel('$|B|$ (G)', 'Interpreter', 'latex', 'Fontsize', 20);
ylabel(yl, 'Interpreter', 'latex','Fontsize', 20);
ylim(yy);
title(titleg,'Interpreter', 'latex', 'Fontsize', 22);
legend_labels = arrayfun(@(f) sprintf('$F = %.0f$', f), F_values, 'UniformOutput', false);
legend(legend_handles, legend_labels, 'Interpreter', 'latex', 'Location', 'northwest');
box on;
hold off;
Filename = sprintf('%s/Zeeman splitting %s.png', Folder, str);
print(Fig, Filename, '-dpng', '-r300');

toc

%%

function [Sx, Sy, Sz] = espin(S)
    
    m_values = -S:(S-1);  % Magnetic quantum numbers m from -S to S-1
    sp = sqrt(S*(S+1) - m_values.*(m_values+1));  % Eigenvalues of S_+
    
    Splus = diag(sp, 1);  % Construct the S_+ matrix

    % Construct the S_- matrix
    Sminus = Splus.';
    
    % Construct the Sx and Sy matrices
    Sx = (Splus + Sminus)/2;
    Sy = (Splus - Sminus)/(2*1i);

    spz = S:-1:-S; % Generate the vector [S, S-1, ..., -S] for Sz matrix
    
    % Create Sz as a diagonal matrix
    Sz = diag(spz);

end

function ident = id(S)
    ident = eye(2*S+1);
end

[Alan Stevens]

Where is your espin function?

[Sebastián Martínez]

Hi. I apologize, I've already updated the code in the publication, but I'm still attaching it here.

%%
function [Sx, Sy, Sz] = espin(S)
    
    m_values = -S:(S-1);  % Magnetic quantum numbers m from -S to S-1
    sp = sqrt(S*(S+1) - m_values.*(m_values+1));  % Eigenvalues of S_+
    
    Splus = diag(sp, 1);  % Construct the S_+ matrix
    % Construct the S_- matrix
    Sminus = Splus.';
    
    % Construct the Sx and Sy matrices
    Sx = (Splus + Sminus)/2;
    Sy = (Splus - Sminus)/(2*1i);
    spz = S:-1:-S; % Generate the vector [S, S-1, ..., -S] for Sz matrix
    
    % Create Sz as a diagonal matrix
    Sz = diag(spz);
end
function ident = id(S)
    ident = eye(2*S+1);
end

Answer 1

[David]

Yes, the issue is ordering the eigenvalues. I am working on a simlar analysis for Rubidium 87 and solved this issue by tracking the eigenvectors and taking the dot product with the previous step's eigenvectors.

[eig_vec, eigvals] = eig(HZeeman);
eigvals = real(diag(HZeeman));  % Extract real eigenvalues

if j == 1
    % Store initial eigenvalues and eigenvectors
    energies(:, j) = eigvals;
    eigenvectors(:, :, j) = eig_vec;
else
    % Compute overlap matrix with previous step's eigenvectors
    overlap = abs(eigenvectors(:, :, j-1)' * eig_vec);
    
    % Find best-matching eigenvectors using maximum overlap
    [~, indices] = max(overlap, [], 2);
    
    % Reorder eigenvalues and eigenvectors to maintain continuity
    energy2(:, j) = eigvals(indices);
    eigenvectors(:, :, j) = eig_vec(:, indices);

Answer 2

[David Goodmanson]

Hi Sebastian,

The best way would be to find four subspaces that are not coupled by the Hamiltonian when BF is included. That may not be that easy since F and Jz are no longer good quantum numbers. Flipping the colors when lines cross seems to be the next best thing. The code below does that for 34 crossings of different colors and produces this:

oo Before plotting I used your code and just stuck with the 'D2' option with J = 3/2, making no attempt to generalize to other values of J and I, although it would not be very hard to do. There were three changes: The B field is

delB = .02;    BF = 0:delB:500;          % Magnetic field vector (G)

The sorting is descending since that seemed easier to think about,

energies(j, :) = sort(diag(real(eigenvalues)),'descend');

and in the statement

HZeemanBase = A*(Ix*Jx + Iy*Jy + Iz*Jz) + ...
    mu_B*BF(1)*(g_I*Iz + g_J*Jz) + ...

I replaced BF(1) with 0. The BF(1) placeholder idea is very dangerous because the intended value is zero and the code only works when the BF starting value is 0. If someone wants to check out the situation using, say, BF = 100:.01:200, the answer goes bad in a nonobvious way.

oo In the Hamiltonian, the code has a factor

( 3*(Ix*Jx + Iy*Jy + Iz*Jz).^2 

should this be ^2 without the dot? Either way seems to make no difference in the answer.

oo The plotting code uses the fact that when energies cross and are sorted, there is a discontinuity in the slope, so there should be a large value for the second derivative. There is a some crossing below BF = 50 but since this involves curves of the same intended color, no correction has to be done.

d = diff(energies);
d2 = diff(d);
d2(1:round(50/delB),:) = 0;   % crossings below BF = 50 don't count
% figure(8)
% plot(BF,(energies))
% figure(9)
% plot(BF(2:end),d)
figure(10)
plot(BF(2:end-1),d2,'-o')  % crossing point detection 
grid on
xlim([50 250])
[a1 b1] = find(d2> .001);  % crossing point indices
[a2 b2] = find(d2<-.001);  % find the 'opposite side' of the crossing point,
                           % not strictly necessary but good as a check
inds = [a1 b1 a2 b2];
inds = sortrows(inds);
inds = inds(1:2:end,:);      % lower index of crossing points
% denote colors by F value even though
% it is not a good quantum number when BF~=0
C = zeros(size(energies));   
C(:,1:11) = 5;
C(:,12:20) = 4;
C(:,21:27) = 3;
C(:,28:32) = 2;
for k = 1:size(inds,1)
    BFind = inds(k,1);
    levlo = inds(k,2);
    levhi = inds(k,4);       % levhi = levlo +1, good as a check
    Clo = C(BFind,levlo);
    Chi = C(BFind,levhi);
    C(BFind+1:end,levlo) = Chi;  % swap colors from here on out
    C(BFind+1:end,levhi) = Clo;
end
figure(11)
str = 'xrbgm'
for c = 5:-1:2
    energiestemp = energies;
    energiestemp(C~=c) = nan;
    energiestemp(:,all(isnan(energiestemp))) = [];
    plot(BF,energiestemp,str(c))    % nans are not plotted
    hold on
end
grid on
hold off

crossing point detection, fig 10: