Latest - complete script for extraction of edd, n, k_roton for different tilts of the dipole moment and different trapping frequencies and minor modifications to other scripts.

This commit is contained in:
Karthik 2025-01-15 18:52:24 +01:00
parent f3ac9e3b12
commit bc6234e046
4 changed files with 261 additions and 113 deletions

View File

@ -77,7 +77,7 @@ MeanWidth = sigma * lz;
VDk = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth);
VDk_fftshifted = fftshift(VDk);
figure(11)
figure(8)
clf
set(gcf,'Position',[50 50 950 750])
imagesc(fftshift(Transf.kx)*1e-6, fftshift(Transf.ky)*1e-6, VDk_fftshifted); % Specify x and y data for axes
@ -106,7 +106,7 @@ end
EpsilonK = double(imag(EpsilonK) ~= 0); % 'isreal' returns 0 for complex numbers and 1 for real numbers, so we negate it
figure(12)
figure(9)
clf
set(gcf,'Position',[50 50 950 750])
imagesc(fftshift(Transf.kx)*1e-6, fftshift(Transf.ky)*1e-6, EpsilonK); % Specify x and y data for axes
@ -159,7 +159,7 @@ for theta = theta_values
EpsilonK = double(imag(EpsilonK) ~= 0); % 'isreal' returns 0 for complex numbers and 1 for real numbers, so we negate it
% Plot the result
figure(13)
figure(10)
clf
set(gcf,'Position',[50 50 950 750])
imagesc(fftshift(Transf.kx)*1e-6, fftshift(Transf.ky)*1e-6, EpsilonK); % Specify x and y data for axes

View File

@ -28,123 +28,114 @@ Dy164Mass = 163.929174751*AtomicMassUnit;
Dy164IsotopicAbundance = 0.2826;
DyMagneticMoment = 9.93*BohrMagneton;
%% Roton instability boundary for tilted dipoles
%% Extracting values from the roton instability boundary for tilted dipoles
wz = 2 * pi * 72.4; % Trap frequency in the tight confinement direction
theta = 0; % Polar angle of dipole moment
%-------TEST-------%
% nadd2s = 0.05:0.005:0.25;
% as_to_add = 0.76:0.0001:0.81;
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;
%-------DEPLOY-------%
nadd2s = 0.005:0.005:0.5;
as_to_add = 0.35:0.0001:1.15;
nadd2s = 0.05:0.001:0.25;
as_to_add = 0.76:0.001:0.81;
var_widths = zeros(length(as_to_add), length(nadd2s));
data_struct = struct;
wz_values = [150, 300, 500, 750];
theta_values = 0:5:45; % Range of theta values (you can modify this)
phi = 0; % Azimuthal angle of momentum vector
kvec = linspace(0, 2.25e6, 1000); % Vector of magnitudes of k vector
x0 = 5;
Aineq = [];
Bineq = [];
Aeq = [];
Beq = [];
lb = [1];
ub = [10];
nonlcon = [];
fminconopts = optimoptions(@fmincon,'Display','off', 'StepTolerance', 1.0000e-11, 'MaxIterations',1500);
for mainloop_idx = 1:length(wz_values)
format long
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
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;
%% ====================================================================================================================================================== %
phi = 0; % Azimuthal angle of momentum vector
k = linspace(0, 2.25e6, 1000); % Vector of magnitudes of k vector
instability_boundary = zeros(length(as_to_add), length(nadd2s));
ScatteringLengths = zeros(length(as_to_add), 1);
AtomNumber = zeros(length(nadd2s), 1);
w0 = 2 * pi * 61.6316; % Trap frequency in the tight confinement direction
l0 = sqrt(PlanckConstantReduced/(Dy164Mass * w0)); % Defining a harmonic oscillator length
tsize = 10 * l0;
for idx = 1:length(nadd2s)
for jdx = 1:length(as_to_add)
AtomNumberDensity = nadd2s(idx) / add^2; % Areal density of atoms
AtomNumber(idx) = AtomNumberDensity*tsize^2;
as = (as_to_add(jdx) * add); % Scattering length
ScatteringLengths(jdx) = as/BohrRadius;
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
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));
[Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, 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 .* k.^2) ./ (2 * Dy164Mass)) + ((2 * AtomNumberDensity) .* Ukk) + (3 * gQF * AtomNumberDensity^(3/2));
EpsilonK = sqrt(((PlanckConstantReduced^2 .* k.^2) ./ (2 * Dy164Mass)) .* DeltaK);
instability_boundary(jdx, idx) = ~isreal(EpsilonK);
end
end
figure(6)
clf
set(gcf,'Position',[50 50 950 750])
imagesc(nadd2s, as_to_add, instability_boundary); % Specify x and y data for axes
set(gca, 'YDir', 'normal'); % Correct the y-axis direction
colorbar; % Add a colorbar
xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex');
ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex');
title('Roton instability boundary','fontsize',16,'interpreter','latex')
%%
matrix = instability_boundary;
% 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];
col_indices = [col_indices; j];
break; % Stop after the first transition in the column
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);
n_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), n_values(idx), k_roton_values(idx)] = extractFromBoundaryCurve(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).n_values = n_values;
data_struct(mainloop_idx).k_roton_values = k_roton_values;
%{
figure(13)
clf
set(gcf,'Position',[50 50 950 750])
plot(theta_values, eps_dd_values, '-o', LineWidth=2.0)
xlabel('$\theta$','fontsize',16,'interpreter','latex');
ylabel('$\epsilon_{dd}$','fontsize',16,'interpreter','latex');
% title([''],'fontsize',16,'interpreter','latex')
grid on
figure(14)
clf
set(gcf,'Position',[50 50 950 750])
plot(theta_values, (1./eps_dd_values) * (add/BohrRadius), '-o', LineWidth=2.0)
xlabel('$\theta$','fontsize',16,'interpreter','latex');
ylabel('$a_s (\times a_o)$','fontsize',16,'interpreter','latex');
% title([''],'fontsize',16,'interpreter','latex')
grid on
figure(15)
clf
set(gcf,'Position',[50 50 950 750])
plot(theta_values, n_values * 1E-15, '-o', LineWidth=2.0)
xlabel('$\theta$','fontsize',16,'interpreter','latex');
ylabel('$n (\times 10^{3} \mu m^{-2})$','fontsize',16,'interpreter','latex');
% title([''],'fontsize',16,'interpreter','latex')
grid on
figure(16)
clf
set(gcf,'Position',[50 50 950 750])
plot(theta_values, k_roton_values * 1E-6, '-o', LineWidth=2.0)
xlabel('$\theta$','fontsize',16,'interpreter','latex');
ylabel('$k_{roton} (\mu m^{-1})$','fontsize',16,'interpreter','latex');
% title([''],'fontsize',16,'interpreter','latex')
grid on
%}
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
% Plot the extracted values
figure(7);
plot(xvals, yvals, '-o');
title('Instability Boundary');
xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex');
ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex')
save('ExtractingKRoton_Result.mat', 'data_struct');
%%
function [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, theta, phi)
@ -164,4 +155,161 @@ function ret = computeTotalEnergyPerParticle(x, as, AtomNumberDensity, wz, lz, g
Energy_AxialComponent = (PlanckConstantReduced * wz) * ((lz^2/(4 * MeanWidth^2)) + (MeanWidth^2/(4 * lz^2)));
Energy_TransverseComponent = (0.5 * (gs + (2*gdd)) * gamma4 * AtomNumberDensity) + ((2/5) * gQF * AtomNumberDensity^(3/2));
ret = (Energy_AxialComponent + Energy_TransverseComponent) / (PlanckConstantReduced * wz);
end
end
function [eps_dd, AtomNumberDensity, k_roton] = extractFromBoundaryCurve(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
%{
figure(11)
clf
set(gcf,'Position',[50 50 950 750])
imagesc(nadd2s, as_to_add, phase_diagram); % Specify x and y data for axes
set(gca, 'YDir', 'normal'); % Correct the y-axis direction
colorbar; % Add a colorbar
xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex');
ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex');
title(['$\theta = ',num2str(theta), '; \phi = 0','$', '(Along Y)'],'fontsize',16,'interpreter','latex')
%}
%-------------%
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];
%-------------%
% Degree of the polynomial to fit
n = 5; % For a quadratic fit
% Fit the polynomial
p = polyfit(xvals, yvals, n);
% Display the polynomial coefficients
% disp('Polynomial coefficients:');
% disp(p);
% Evaluate the polynomial at points in x
y_fit = polyval(p, xvals);
%{
% Plot the original data and the fitted polynomial curve
figure(12);
clf
set(gcf,'Position',[50 50 950 750])
plot(xvals, yvals, 'o', 'LineWidth', 2.0, 'DisplayName', 'Extracted boundary points'); % Original data
hold on;
plot(xvals, y_fit, '-r','LineWidth', 2.0, 'DisplayName', ['Polynomial Fit (degree ' num2str(n) ')']); % Fitted curve
ylim([min(as_to_add) max(as_to_add)])
xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex');
ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex')
title(['$\theta = ',num2str(theta), '; \phi = 0','$', '(Along Y)'],'fontsize',16,'interpreter','latex')
legend('show');
grid on;
%}
[val, idx] = max(y_fit);
% Round down to 4 decimal places
rounded_val = floor(val * 10^4) / 10^4;
% Find nearest from original vector of boundary points
[~, nearest_idx] = min(abs(instability_boundary(:, 2) - rounded_val));
nearest_val = instability_boundary(nearest_idx, 2);
% Choose the scalar value between the two
if ~isscalar(nearest_val)
val = rounded_val;
else
val = nearest_val;
end
AtomNumberDensity = xvals(idx) / 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
end

View File

@ -62,7 +62,7 @@ end
%% ====================================================================================================================================================== %
figure(6)
figure(7)
clf
set(gcf,'Position',[50 50 1850 750])
@ -181,7 +181,7 @@ v.FrameRate = 5; % Frame rate of the video
open(v); % Open the video file
for theta = theta_values
figure(6)
figure(7)
clf
set(gcf,'Position',[50 50 1850 750])
phi = 0; % Azimuthal angle of momentum vector

View File

@ -76,7 +76,7 @@ for idx = 1:length(nadd2s)
end
end
figure(7)
figure
clf
set(gcf,'Position',[50 50 950 750])
imagesc(AtomNumber*1E-5, ScatteringLengths, QF * 1E31); % Specify x and y data for axes