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’);
tocHi. 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 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 matrix, eigenvalues, diagonalization, atomic physics MATLAB Answers — New Questions
