karthik/Calculations
1
0
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity

Latest - added new analysis and plotting routines.

Browse Source
This commit is contained in:
Karthik 2025-01-22 15:58:05 +01:00
parent 0889d1178b
commit 1170f6a8a6
4 changed files with 2047 additions and 67 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

355
Estimations/DipolarDispersionAndRotonInstabilityBoundary/AnalyzeResults.m
View File

@ -1,4 +1,8 @@
load('ExtractingKRoton_Result.mat')
%% Across range of a_s, n
% load('.\Results\ExtractingKRoton_Result_Below1000.mat')
% load('.\Results\ExtractingKRoton_Result_Above1000.mat')
load('.\Results\ExtractingKRoton_Result_Above10000.mat')
PlanckConstantReduced = 6.62607015E-34/(2*pi);
AtomicMassUnit = 1.660539066E-27;
@ -10,6 +14,293 @@ DyMagneticMoment = 9.93*BohrMagneton;
add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
% Create a tiled layout with tighter spacing
figure(17)
clf
set(gcf,'Position',[50 50 1800 500])
t = tiledlayout(1, 3, 'TileSpacing', 'compact', 'Padding', 'compact'); % 2x2 grid
% First subplot
nexttile; % Equivalent to subplot(2, 2, 1)
for idx = 1:length(data_struct)
theta_values = data_struct(idx).theta_values;
eps_dd_values = data_struct(idx).eps_dd_values;
plot(theta_values, eps_dd_values, '-o', 'LineWidth', 2.0, 'DisplayName', ['$w_z = 2 \pi \times $', num2str(data_struct(idx).wz_value), ' Hz']);
hold on
end
xlabel('$\theta$', 'fontsize', 16, 'interpreter', 'latex');
ylabel('$\epsilon_{dd}$', 'fontsize', 16, 'interpreter', 'latex');
grid on
legend('location', 'northeast', 'fontsize', 10, 'Interpreter', 'latex'); % Reduced font size
% Second subplot
nexttile; % Equivalent to subplot(2, 2, 2)
for idx = 1:length(data_struct)
theta_values = data_struct(idx).theta_values;
n_values = data_struct(idx).n_values;
plot(theta_values, n_values * 1E-15, '-o', 'LineWidth', 2.0, 'DisplayName', ['$w_z = 2 \pi \times $', num2str(data_struct(idx).wz_value), ' Hz']);
hold on
end
xlabel('$\theta$', 'fontsize', 16, 'interpreter', 'latex');
ylabel('$n (\times 10^{3} \mu m^{-2})$', 'fontsize', 16, 'interpreter', 'latex');
grid on
legend('location', 'northeast', 'fontsize', 10, 'Interpreter', 'latex'); % Reduced font size
% Third subplot
nexttile; % Equivalent to subplot(2, 2, 3)
for idx = 1:length(data_struct)
theta_values = data_struct(idx).theta_values;
k_roton_values = data_struct(idx).k_roton_values;
plot(theta_values, k_roton_values * 1E-6, '-o', LineWidth=2.0, DisplayName=['$w_z = 2 \pi \times $', num2str(data_struct(idx).wz_value), ' Hz']);
hold on
end
xlabel('$\theta$','fontsize',16,'interpreter','latex');
ylabel('$k_{roton} (\mu m^{-1})$','fontsize',16,'interpreter','latex');
grid on
legend('location', 'northeast','fontsize', 10, 'Interpreter','latex')
% Adjust layout to minimize space
t.TileSpacing = 'compact'; % Minimize space between tiles
t.Padding = 'compact'; % Minimize padding around the layout
% Convert to units relevant to experiment
% Create a tiled layout with tighter spacing
figure(18)
clf
set(gcf,'Position',[50 50 1800 500])
t = tiledlayout(1, 3, 'TileSpacing', 'compact', 'Padding', 'compact'); % 2x2 grid
% First subplot
nexttile; % Equivalent to subplot(2, 2, 1)
for idx = 1:length(data_struct)
theta_values = data_struct(idx).theta_values;
eps_dd_values = data_struct(idx).eps_dd_values;
plot(theta_values, (1 ./ eps_dd_values) * (add / BohrRadius), '-o', 'LineWidth', 2.0, 'DisplayName', ['$w_z = 2 \pi \times $', num2str(data_struct(idx).wz_value), ' Hz']);
hold on
end
xlabel('$\theta$', 'fontsize', 16, 'interpreter', 'latex');
ylabel('$a_s (\times a_o)$', 'fontsize', 16, 'interpreter', 'latex');
grid on
legend('location', 'southeast', 'fontsize', 10, 'Interpreter', 'latex'); % Reduced font size
% Second subplot
nexttile; % Equivalent to subplot(2, 2, 2)
for idx = 1:length(data_struct)
theta_values = data_struct(idx).theta_values;
n_values = data_struct(idx).n_values;
Lx = 10e-6;
Ly = 10e-6;
AtomNumber = n_values .* Lx * Ly;
plot(theta_values, AtomNumber * 1e-5, '-o', 'LineWidth', 2.0, 'DisplayName', ['$w_z = 2 \pi \times $', num2str(data_struct(idx).wz_value), ' Hz']);
hold on
end
xlabel('$\theta$', 'fontsize', 16, 'interpreter', 'latex');
ylabel('Atom number in a trap of area 100 $\mu m^2 (\times 10^{5})$', 'fontsize', 16, 'interpreter', 'latex');
grid on
legend('location', 'northeast', 'fontsize', 10, 'Interpreter', 'latex'); % Reduced font size
% Third subplot
nexttile; % Equivalent to subplot(2, 2, 3)
for idx = 1:length(data_struct)
theta_values = data_struct(idx).theta_values;
lambda_roton_values = (2 * pi) ./ data_struct(idx).k_roton_values;
plot(theta_values, lambda_roton_values * 1E6, '-o', LineWidth=2.0, DisplayName=['$w_z = 2 \pi \times $', num2str(data_struct(idx).wz_value), ' Hz']);
hold on
end
xlabel('$\theta$','fontsize',16,'interpreter','latex');
ylabel('$\lambda_{roton} (\mu m)$','fontsize',16,'interpreter','latex');
grid on
legend('location', 'northeast','fontsize', 10, 'Interpreter','latex')
% Adjust layout to minimize space
t.TileSpacing = 'compact'; % Minimize space between tiles
t.Padding = 'compact'; % Minimize padding around the layout
%% Fixed Density results
load('.\Results\ExtractingKRoton_Result_FixedDensity_phi0.mat')
PlanckConstantReduced = 6.62607015E-34/(2*pi);
AtomicMassUnit = 1.660539066E-27;
Dy164Mass = 163.929174751*AtomicMassUnit;
VacuumPermeability = 1.25663706212E-6;
BohrMagneton = 9.274009994E-24;
BohrRadius = 5.2917721067E-11;
DyMagneticMoment = 9.93*BohrMagneton;
add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
% Create a tiled layout with tighter spacing
figure(19)
clf
set(gcf,'Position',[50 50 1200 500])
t = tiledlayout(1, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 2x2 grid
% First subplot
nexttile;
for idx = 1:length(data_struct)
theta_values = data_struct(idx).theta_values;
eps_dd_values = data_struct(idx).eps_dd_values;
plot(theta_values, eps_dd_values, '-o', 'LineWidth', 2.0, 'DisplayName', ['$w_z = 2 \pi \times $', num2str(data_struct(idx).wz_value), ' Hz']);
hold on
end
xlabel('$\theta$', 'fontsize', 16, 'interpreter', 'latex');
ylabel('$\epsilon_{dd}$', 'fontsize', 16, 'interpreter', 'latex');
grid on
legend('location', 'northeast', 'fontsize', 10, 'Interpreter', 'latex'); % Reduced font size
% Second subplot
nexttile;
for idx = 1:length(data_struct)
theta_values = data_struct(idx).theta_values;
k_roton_values = data_struct(idx).k_roton_values;
plot(theta_values, k_roton_values * 1E-6, '-o', LineWidth=2.0, DisplayName=['$w_z = 2 \pi \times $', num2str(data_struct(idx).wz_value), ' Hz']);
hold on
end
xlabel('$\theta$','fontsize',16,'interpreter','latex');
ylabel('$k_{roton} (\mu m^{-1})$','fontsize',16,'interpreter','latex');
grid on
legend('location', 'northeast','fontsize', 10, 'Interpreter','latex')
% Adjust layout to minimize space
t.TileSpacing = 'compact'; % Minimize space between tiles
t.Padding = 'compact'; % Minimize padding around the layout
% Create a tiled layout with tighter spacing
figure(20)
clf
set(gcf,'Position',[50 50 1200 500])
t = tiledlayout(1, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 2x2 grid
% First subplot
nexttile;
for idx = 1:length(data_struct)
theta_values = data_struct(idx).theta_values;
eps_dd_values = data_struct(idx).eps_dd_values;
plot(theta_values, (1 ./ eps_dd_values) * (add / BohrRadius), '-o', 'LineWidth', 2.0, 'DisplayName', ['$w_z = 2 \pi \times $', num2str(data_struct(idx).wz_value), ' Hz']);
hold on
end
xlabel('$\theta$', 'fontsize', 16, 'interpreter', 'latex');
ylabel('$a_s (\times a_o)$', 'fontsize', 16, 'interpreter', 'latex');
grid on
legend('location', 'northwest', 'fontsize', 10, 'Interpreter', 'latex'); % Reduced font size
% Second subplot
nexttile;
for idx = 1:length(data_struct)
theta_values = data_struct(idx).theta_values;
lambda_roton_values = (2 * pi) ./ data_struct(idx).k_roton_values;
semilogy(theta_values, lambda_roton_values * 1E6, '-o', LineWidth=2.0, DisplayName=['$w_z = 2 \pi \times $', num2str(data_struct(idx).wz_value), ' Hz']);
hold on
end
% ylim([0 2])
xlabel('$\theta$','fontsize',16,'interpreter','latex');
ylabel('$\lambda_{roton} (\mu m)$','fontsize',16,'interpreter','latex');
grid on
legend('location', 'southeast','fontsize', 10, 'Interpreter','latex')
% Adjust layout to minimize space
t.TileSpacing = 'compact'; % Minimize space between tiles
t.Padding = 'compact'; % Minimize padding around the layout
%% Fixed Density results - compare two orthogonal directions
data0 = load('.\Results\ExtractingKRoton_Result_FixedDensity_phi0.mat');
data90 = load('.\Results\ExtractingKRoton_Result_FixedDensity_phi90.mat');
PlanckConstantReduced = 6.62607015E-34/(2*pi);
AtomicMassUnit = 1.660539066E-27;
Dy164Mass = 163.929174751*AtomicMassUnit;
VacuumPermeability = 1.25663706212E-6;
BohrMagneton = 9.274009994E-24;
BohrRadius = 5.2917721067E-11;
DyMagneticMoment = 9.93*BohrMagneton;
add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
% Create a tiled layout with tighter spacing
figure(21)
clf
set(gcf,'Position',[50 50 1200 500])
t = tiledlayout(1, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 2x2 grid
idx = 4;
% First subplot
nexttile;
theta_values = data0.data_struct(idx).theta_values;
eps_dd_values = data0.data_struct(idx).eps_dd_values;
plot(theta_values, eps_dd_values, '-o', 'LineWidth', 2.0, 'DisplayName', ['$w_z = 2 \pi \times $', num2str(data0.data_struct(idx).wz_value), ' Hz; $\phi = 0^\circ$']);
hold on
theta_values = data90.data_struct(idx).theta_values;
eps_dd_values = data90.data_struct(idx).eps_dd_values;
plot(theta_values, eps_dd_values, '-o', 'LineWidth', 2.0, 'DisplayName', ['$w_z = 2 \pi \times $', num2str(data90.data_struct(idx).wz_value), ' Hz; $\phi = 90^\circ$']);
xlabel('$\theta$', 'fontsize', 16, 'interpreter', 'latex');
ylabel('$\epsilon_{dd}$', 'fontsize', 16, 'interpreter', 'latex');
grid on
legend('location', 'northeast', 'fontsize', 10, 'Interpreter', 'latex'); % Reduced font size
% Second subplot
nexttile;
theta_values = data0.data_struct(idx).theta_values;
k_roton_values = data0.data_struct(idx).k_roton_values;
plot(theta_values, k_roton_values * 1E-6, '-o', LineWidth=2.0, DisplayName=['$w_z = 2 \pi \times $', num2str(data0.data_struct(idx).wz_value), ' Hz; $\phi = 0^\circ$']);
hold on
theta_values = data90.data_struct(idx).theta_values;
k_roton_values = data90.data_struct(idx).k_roton_values;
plot(theta_values, k_roton_values * 1E-6, '-o', LineWidth=2.0, DisplayName=['$w_z = 2 \pi \times $', num2str(data90.data_struct(idx).wz_value), ' Hz; $\phi = 90^\circ$']);
xlabel('$\theta$','fontsize',16,'interpreter','latex');
ylabel('$k_{roton} (\mu m^{-1})$','fontsize',16,'interpreter','latex');
grid on
legend('location', 'northeast','fontsize', 10, 'Interpreter','latex')
% Adjust layout to minimize space
t.TileSpacing = 'compact'; % Minimize space between tiles
t.Padding = 'compact'; % Minimize padding around the layout
% Create a tiled layout with tighter spacing
figure(22)
clf
set(gcf,'Position',[50 50 1200 500])
t = tiledlayout(1, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 2x2 grid
% First subplot
nexttile;
theta_values = data0.data_struct(idx).theta_values;
eps_dd_values = data0.data_struct(idx).eps_dd_values;
plot(theta_values, (1 ./ eps_dd_values) * (add / BohrRadius), '-o', 'LineWidth', 2.0, 'DisplayName', ['$w_z = 2 \pi \times $', num2str(data0.data_struct(idx).wz_value), ' Hz; $\phi = 0^\circ$']);
hold on
theta_values = data90.data_struct(idx).theta_values;
eps_dd_values = data90.data_struct(idx).eps_dd_values;
plot(theta_values, (1 ./ eps_dd_values) * (add / BohrRadius), '-o', 'LineWidth', 2.0, 'DisplayName', ['$w_z = 2 \pi \times $', num2str(data90.data_struct(idx).wz_value), ' Hz; $\phi = 90^\circ$']);
xlabel('$\theta$', 'fontsize', 16, 'interpreter', 'latex');
ylabel('$a_s (\times a_o)$', 'fontsize', 16, 'interpreter', 'latex');
grid on
legend('location', 'northwest', 'fontsize', 10, 'Interpreter', 'latex'); % Reduced font size
% Second subplot
nexttile;
theta_values = data0.data_struct(idx).theta_values;
k_roton_values = data0.data_struct(idx).k_roton_values;
lambda_roton_values = (2 * pi) ./ k_roton_values;
semilogy(theta_values, lambda_roton_values * 1E6, '-o', LineWidth=2.0, DisplayName=['$w_z = 2 \pi \times $', num2str(data0.data_struct(idx).wz_value), ' Hz; $\phi = 0^\circ$']);
hold on
theta_values = data90.data_struct(idx).theta_values;
k_roton_values = data90.data_struct(idx).k_roton_values;
lambda_roton_values = (2 * pi) ./ k_roton_values;
semilogy(theta_values, lambda_roton_values * 1E6, '-o', LineWidth=2.0, DisplayName=['$w_z = 2 \pi \times $', num2str(data90.data_struct(idx).wz_value), ' Hz; $\phi = 90^\circ$']);
xlabel('$\theta$','fontsize',16,'interpreter','latex');
ylabel('$\lambda_{roton} (\mu m)$','fontsize',16,'interpreter','latex');
grid on
legend('location', 'northwest','fontsize', 10, 'Interpreter','latex')
% Adjust layout to minimize space
t.TileSpacing = 'compact'; % Minimize space between tiles
t.Padding = 'compact'; % Minimize padding around the layout
%%
%{
figure(13)
clf
@ -71,65 +362,3 @@ ylabel('$k_{roton} (\mu m^{-1})$','fontsize',16,'interpreter','latex');
grid on
legend('location', 'northeast','fontsize', 16, 'Interpreter','latex')
%}
% Create a tiled layout with tighter spacing
figure(17)
clf
set(gcf,'Position',[50 50 1200 900])
t = tiledlayout(2, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 2x2 grid
% First subplot
nexttile; % Equivalent to subplot(2, 2, 1)
for idx = 1:length(data_struct)
theta_values = data_struct(idx).theta_values;
eps_dd_values = data_struct(idx).eps_dd_values;
plot(theta_values, eps_dd_values, '-o', 'LineWidth', 2.0, 'DisplayName', ['$w_z = 2 \pi \times $', num2str(data_struct(idx).wz_value), ' Hz']);
hold on
end
xlabel('$\theta$', 'fontsize', 16, 'interpreter', 'latex');
ylabel('$\epsilon_{dd}$', 'fontsize', 16, 'interpreter', 'latex');
grid on
legend('location', 'northeast', 'fontsize', 10, 'Interpreter', 'latex'); % Reduced font size
% Second subplot
nexttile; % Equivalent to subplot(2, 2, 2)
for idx = 1:length(data_struct)
theta_values = data_struct(idx).theta_values;
eps_dd_values = data_struct(idx).eps_dd_values;
plot(theta_values, (1 ./ eps_dd_values) * (add / BohrRadius), '-o', 'LineWidth', 2.0, 'DisplayName', ['$w_z = 2 \pi \times $', num2str(data_struct(idx).wz_value), ' Hz']);
hold on
end
xlabel('$\theta$', 'fontsize', 16, 'interpreter', 'latex');
ylabel('$a_s (\times a_o)$', 'fontsize', 16, 'interpreter', 'latex');
grid on
legend('location', 'southeast', 'fontsize', 10, 'Interpreter', 'latex'); % Reduced font size
% Third subplot
nexttile; % Equivalent to subplot(2, 2, 3)
for idx = 1:length(data_struct)
theta_values = data_struct(idx).theta_values;
n_values = data_struct(idx).n_values;
plot(theta_values, n_values * 1E-15, '-o', 'LineWidth', 2.0, 'DisplayName', ['$w_z = 2 \pi \times $', num2str(data_struct(idx).wz_value), ' Hz']);
hold on
end
xlabel('$\theta$', 'fontsize', 16, 'interpreter', 'latex');
ylabel('$n (\times 10^{3} \mu m^{-2})$', 'fontsize', 16, 'interpreter', 'latex');
grid on
legend('location', 'northeast', 'fontsize', 10, 'Interpreter', 'latex'); % Reduced font size
% Fourth subplot
nexttile; % Equivalent to subplot(2, 2, 4)
for idx = 1:length(data_struct)
theta_values = data_struct(idx).theta_values;
k_roton_values = data_struct(idx).k_roton_values;
plot(theta_values, k_roton_values * 1E-6, '-o', LineWidth=2.0, DisplayName=['$w_z = 2 \pi \times $', num2str(data_struct(idx).wz_value), ' Hz']);
hold on
end
xlabel('$\theta$','fontsize',16,'interpreter','latex');
ylabel('$k_{roton} (\mu m^{-1})$','fontsize',16,'interpreter','latex');
grid on
legend('location', 'northeast','fontsize', 10, 'Interpreter','latex')
% Adjust layout to minimize space
t.TileSpacing = 'compact'; % Minimize space between tiles
t.Padding = 'compact'; % Minimize padding around the layout

195
Estimations/DipolarDispersionAndRotonInstabilityBoundary/ExtractingKRoton.m
View File

@ -36,13 +36,21 @@ DyMagneticMoment = 9.93*BohrMagneton;
%-------DEPLOY-------%
nadd2s = 0.005:0.005:0.5;
as_to_add = 0.35:0.001:1.15;
as_to_add = 0.250:0.001:1.15;
data_struct = struct;
wz_values = [150, 300, 500, 750];
% wz_values = [150, 300, 500, 750];
% kvec = linspace(0, 5e6, 1000); % Vector of magnitudes of k vector
wz_values = [1000, 3000, 5000, 7000];
kvec = linspace(0, 15e6, 1000); % Vector of magnitudes of k vector
% wz_values = [10000, 13000, 15000];
% kvec = linspace(0, 25e6, 1000); % Vector of magnitudes of k vector
theta_values = 0:5:45; % Range of theta values
phi = 0; % Azimuthal angle of momentum vector
kvec = linspace(0, 5e6, 1000); % Vector of magnitudes of k vector
for mainloop_idx = 1:length(wz_values)
format long
@ -135,8 +143,80 @@ for mainloop_idx = 1:length(wz_values)
%}
end
save('ExtractingKRoton_Result.mat', 'data_struct');
save('.\Results\ExtractingKRoton_Result.mat', 'data_struct');
%% Extracting values from the roton instability boundary for tilted dipoles - fixed atom number, trap frequency
%-------DEPLOY-------%
N = 1E5;
add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
area = 100; % in square micrometers
ppmum = N / area;
nadd2s = ppmum*1E12*add^2;
as_to_add = 0.150:0.001:1.15;
data_struct = struct;
wz_values = [500, 750, 1000, 2000];
kvec = linspace(0, 15e6, 1000); % Vector of magnitudes of k vector
theta_values = 0:5:90; % Range of theta values
phi = 90; % Azimuthal angle of momentum vector
for mainloop_idx = 1:length(wz_values)
format long
PlanckConstantReduced = 6.62607015E-34/(2*pi);
AtomicMassUnit = 1.660539066E-27;
Dy164Mass = 163.929174751*AtomicMassUnit;
VacuumPermeability = 1.25663706212E-6;
BohrMagneton = 9.274009994E-24;
DyMagneticMoment = 9.93*BohrMagneton;
wz = 2 * pi * wz_values(mainloop_idx); % Trap frequency in the tight confinement direction
lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length
add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
gdd = VacuumPermeability*DyMagneticMoment^2/3;
var_widths = zeros(length(as_to_add), length(nadd2s));
x0 = 5;
Aineq = [];
Bineq = [];
Aeq = [];
Beq = [];
lb = [1];
ub = [10];
nonlcon = [];
fminconopts = optimoptions(@fmincon,'Display','off', 'StepTolerance', 1.0000e-11, 'MaxIterations',1500);
for idx = 1:length(nadd2s)
for jdx = 1:length(as_to_add)
AtomNumberDensity = nadd2s(idx) / add^2; % Areal density of atoms
as = (as_to_add(jdx) * add); % Scattering length
gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength
TotalEnergyPerParticle = @(x) computeTotalEnergyPerParticle(x, as, AtomNumberDensity, wz, lz, gs, add, gdd, PlanckConstantReduced);
sigma = fmincon(TotalEnergyPerParticle, x0, Aineq, Bineq, Aeq, Beq, lb, ub, nonlcon, fminconopts);
var_widths(jdx, idx) = sigma;
end
end
eps_dd_values = zeros(length(theta_values), 1);
k_roton_values = zeros(length(theta_values), 1);
for idx = 1:length(theta_values)
theta = theta_values(idx);
[eps_dd_values(idx), k_roton_values(idx)] = extractFromBoundaryPoint(theta, phi, nadd2s, as_to_add, var_widths, wz, lz, kvec);
end
data_struct(mainloop_idx).wz_value = wz / (2 * pi);
data_struct(mainloop_idx).theta_values = theta_values;
data_struct(mainloop_idx).eps_dd_values = eps_dd_values;
data_struct(mainloop_idx).k_roton_values = k_roton_values;
end
save('.\Results\ExtractingKRoton_Result_FixedDensity_phi90.mat', 'data_struct');
%%
function [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, theta, phi)
Go = sqrt(pi) * (k * MeanWidth/sqrt(2)) .* exp((k * MeanWidth/sqrt(2)).^2) .* erfc((k * MeanWidth/sqrt(2)));
@ -314,3 +394,110 @@ function [eps_dd, AtomNumberDensity, k_roton] = extractFromBoundaryCurve(theta,
k_roton = NaN;
end
end
function [eps_dd, k_roton] = extractFromBoundaryPoint(theta, phi, nadd2s, as_to_add, var_widths, wz, lz, kvec)
format long
PlanckConstantReduced = 6.62607015E-34/(2*pi);
AtomicMassUnit = 1.660539066E-27;
Dy164Mass = 163.929174751*AtomicMassUnit;
VacuumPermeability = 1.25663706212E-6;
BohrMagneton = 9.274009994E-24;
DyMagneticMoment = 9.93*BohrMagneton;
add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
gdd = VacuumPermeability*DyMagneticMoment^2/3;
phase_diagram = zeros(length(as_to_add), length(nadd2s));
w0 = 2 * pi * 61.6316; % Trap frequency in the tight confinement direction
l0 = sqrt(PlanckConstantReduced/(Dy164Mass * w0)); % Defining a harmonic oscillator length
for idx = 1:length(nadd2s)
for jdx = 1:length(as_to_add)
AtomNumberDensity = nadd2s(idx) / add^2; % Areal density of atoms
as = (as_to_add(jdx) * add); % Scattering length
eps_dd = add/as; % Relative interaction strength
gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength
gdd = VacuumPermeability*DyMagneticMoment^2/3;
MeanWidth = var_widths(jdx, idx) * lz; % Mean width of Gaussian ansatz
[Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(kvec, gs, gdd, MeanWidth, theta, phi); % DDI potential in k-space
% == Quantum Fluctuations term == %
gammaQF = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2));
gamma5 = sqrt(2/5) / (sqrt(pi) * MeanWidth)^(3/2);
gQF = gamma5 * gammaQF;
% == Dispersion relation == %
DeltaK = ((PlanckConstantReduced^2 .* kvec.^2) ./ (2 * Dy164Mass)) + ((2 * AtomNumberDensity) .* Ukk) + (3 * gQF * AtomNumberDensity^(3/2));
EpsilonK = sqrt(((PlanckConstantReduced^2 .* kvec.^2) ./ (2 * Dy164Mass)) .* DeltaK);
phase_diagram(jdx, idx) = ~isreal(EpsilonK);
end
end
matrix = phase_diagram;
% Initialize arrays to store row and column indices of transitions
row_indices = [];
col_indices = [];
% Loop through the matrix to find transitions from 0 to 1
[rows, cols] = size(matrix);
for j = 1:cols
for i = 2:rows
if matrix(i-1, j) == 1 && matrix(i, j) == 0
row_indices = [row_indices; i-1];
col_indices = [col_indices; j];
break; % Stop after the first transition in the column
end
end
end
% Now extract the values from the corresponding vectors
xvals = zeros(length(col_indices), 1);
yvals = zeros(length(row_indices), 1);
for k = 1:length(row_indices)
row = row_indices(k);
col = col_indices(k);
xvals(k) = nadd2s(col);
yvals(k) = as_to_add(row);
end
instability_boundary = [xvals, yvals];
if ~isempty(instability_boundary)
val = instability_boundary(2);
AtomNumberDensity = instability_boundary(1) / add^2; % Areal density of atoms
as = val * add; % Scattering length
eps_dd = 1/val; % Relative interaction strength
gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength
x0 = 5;
Aineq = [];
Bineq = [];
Aeq = [];
Beq = [];
lb = [1];
ub = [10];
nonlcon = [];
fminconopts = optimoptions(@fmincon,'Display','off', 'StepTolerance', 1.0000e-11, 'MaxIterations',1500);
TotalEnergyPerParticle = @(x) computeTotalEnergyPerParticle(x, as, AtomNumberDensity, wz, lz, gs, add, gdd, PlanckConstantReduced);
sigma = fmincon(TotalEnergyPerParticle, x0, Aineq, Bineq, Aeq, Beq, lb, ub, nonlcon, fminconopts);
MeanWidth = sigma * lz; % Mean width of Gaussian ansatz
[Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(kvec, gs, gdd, MeanWidth, theta, phi); % DDI potential in k-space
% == Quantum Fluctuations term == %
gammaQF = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2));
gamma5 = sqrt(2/5) / (sqrt(pi) * MeanWidth)^(3/2);
gQF = gamma5 * gammaQF;
DeltaK = ((PlanckConstantReduced^2 .* kvec.^2) ./ (2 * Dy164Mass)) + ((2 * AtomNumberDensity) .* Ukk) + (3 * gQF * AtomNumberDensity^(3/2));
EpsilonK = sqrt(((PlanckConstantReduced^2 .* kvec.^2) ./ (2 * Dy164Mass)) .* DeltaK);
k_roton_indices = find(imag(EpsilonK) ~= 0);
if ~isempty(k_roton_indices)
k_roton = median(kvec(k_roton_indices));
else
k_roton = NaN;
end
else
eps_dd = NaN;
k_roton = NaN;
end
end

1526
Estimations/DipolarDispersionAndRotonInstabilityBoundary/FeshbachResonances.nb Normal file
View File

File diff suppressed because it is too large Load Diff

38
Estimations/DipolarDispersionAndRotonInstabilityBoundary/bwhpc_matlab_gpe_sim_cpu.slurm Normal file
View File

@ -0,0 +1,38 @@
#!/bin/bash
########### Begin SLURM header ###########
#Partition
#SBATCH --partition=cpu-single
# Request number of nodes and CPU cores per node for job
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=1
#SBATCH --cpus-per-task=10
#SBATCH --mem=2G
# Estimated wallclock time for job
#SBATCH --time=00:30:00
#SBATCH --job-name=simulation
#SBATCH --error=simulation.err
#SBATCH --output=simulation.out
########### End SLURM header ##########
echo "Working Directory: $PWD"
echo "Running on host $HOSTNAME"
echo "Job id: $SLURM_JOB_ID"
echo "Job name: $SLURM_JOB_NAME"
echo "Number of nodes allocated to job: $SLURM_JOB_NUM_NODES"
echo "Number of cores allocated to job: $SLURM_JOB_CPUS_PER_NODE"
# Load module
module load math/matlab/R2023a
echo Directory is `pwd`
echo "Initiating Job..."
# Start a Matlab program
matlab -nodisplay -nosplash -r "ExtractingKRoton"
# notice for tests
echo "Job terminated successfully"
exit