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

Name change

Browse Source
This commit is contained in:
Karthik 2025-01-10 11:14:41 +01:00
parent ac47adcce3
commit de6940f2e1
1 changed files with 42 additions and 91 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

107
Estimations/UniformDipolarBECRotonInstabilityBoundary.m → Estimations/ReproduceBlairBlakieResults.m
View File

@ -138,11 +138,17 @@ var_widths = [4.97165, 5.7296048721, 5.93178];
add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
EpsilonKs = zeros(length(k), length(nadd2s)); EpsilonKs = zeros(length(k), 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 idx = 1:length(nadd2s)
AtomNumberDensity = nadd2s(idx) / add^2; % Areal density of atoms AtomNumberDensity = nadd2s(idx) / add^2; % Areal density of atoms
AtomNumber(idx) = AtomNumberDensity*tsize^2;
as = (as_to_add(idx) * add); % Scattering length as = (as_to_add(idx) * add); % Scattering length
ScatteringLengths(idx) = as/BohrRadius;
eps_dd = add/as; % Relative interaction strength eps_dd = add/as; % Relative interaction strength
gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength
gdd = VacuumPermeability*DyMagneticMoment^2/3; gdd = VacuumPermeability*DyMagneticMoment^2/3;
@ -181,11 +187,11 @@ legend('location', 'northwest','fontsize',16, 'Interpreter','latex')
wz = 2 * pi * 72.4; % Trap frequency in the tight confinement direction wz = 2 * pi * 72.4; % Trap frequency in the tight confinement direction
lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length
gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength
add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
gdd = VacuumPermeability*DyMagneticMoment^2/3; gdd = VacuumPermeability*DyMagneticMoment^2/3;
AtomNumberDensity = 0.0978 / add^2; AtomNumberDensity = 0.0978 / add^2;
as = 0.784 * add; % Scattering length as = 0.784 * 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); TotalEnergyPerParticle = @(x) computeTotalEnergyPerParticle(x, as, AtomNumberDensity, wz, lz, gs, add, gdd, PlanckConstantReduced);
x0 = 5; x0 = 5;
@ -204,7 +210,6 @@ fprintf(['Variational width of Gaussian ansatz = ' num2str(sigma) ' * lz \n'])
wz = 2 * pi * 72.4; % Trap frequency in the tight confinement direction wz = 2 * pi * 72.4; % Trap frequency in the tight confinement direction
lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length
gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength
add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
gdd = VacuumPermeability*DyMagneticMoment^2/3; gdd = VacuumPermeability*DyMagneticMoment^2/3;
@ -226,6 +231,7 @@ for idx = 1:length(nadd2s)
for jdx = 1:length(as_to_add) for jdx = 1:length(as_to_add)
AtomNumberDensity = nadd2s(idx) / add^2; % Areal density of atoms AtomNumberDensity = nadd2s(idx) / add^2; % Areal density of atoms
as = (as_to_add(jdx) * add); % Scattering length 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); TotalEnergyPerParticle = @(x) computeTotalEnergyPerParticle(x, as, AtomNumberDensity, wz, lz, gs, add, gdd, PlanckConstantReduced);
sigma = fmincon(TotalEnergyPerParticle, x0, Aineq, Bineq, Aeq, Beq, lb, ub, nonlcon, fminconopts); sigma = fmincon(TotalEnergyPerParticle, x0, Aineq, Bineq, Aeq, Beq, lb, ub, nonlcon, fminconopts);
var_widths(jdx, idx) = sigma; var_widths(jdx, idx) = sigma;
@ -247,11 +253,18 @@ alpha = 0;
phi = 0; % Azimuthal angle of momentum vector phi = 0; % Azimuthal angle of momentum vector
k = linspace(0, 2.25e6, 1000); % Vector of magnitudes of k vector k = linspace(0, 2.25e6, 1000); % Vector of magnitudes of k vector
instability_boundary = zeros(length(as_to_add), length(nadd2s)); 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 idx = 1:length(nadd2s)
for jdx = 1:length(as_to_add) for jdx = 1:length(as_to_add)
AtomNumberDensity = nadd2s(idx) / add^2; % Areal density of atoms AtomNumberDensity = nadd2s(idx) / add^2; % Areal density of atoms
AtomNumber(idx) = AtomNumberDensity*tsize^2;
as = (as_to_add(jdx) * add); % Scattering length as = (as_to_add(jdx) * add); % Scattering length
ScatteringLengths(jdx) = as/BohrRadius;
eps_dd = add/as; % Relative interaction strength eps_dd = add/as; % Relative interaction strength
gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength
gdd = VacuumPermeability*DyMagneticMoment^2/3; gdd = VacuumPermeability*DyMagneticMoment^2/3;
@ -277,6 +290,8 @@ as_to_add_from_figure = [0.76383, 0.76766, 0.76974, 0.77543, 0.77675, 0.77828,
figure(5) figure(5)
clf clf
set(gcf,'Position',[50 50 950 750]) set(gcf,'Position',[50 50 950 750])
imagesc(nadd2s, as_to_add, instability_boundary); % Specify x and y data for axes imagesc(nadd2s, as_to_add, instability_boundary); % Specify x and y data for axes
hold on hold on
plot(nadd2s_from_figure, as_to_add_from_figure, 'r*-', 'LineWidth', 2); % Plot the curve (red line) plot(nadd2s_from_figure, as_to_add_from_figure, 'r*-', 'LineWidth', 2); % Plot the curve (red line)
@ -284,84 +299,17 @@ set(gca, 'YDir', 'normal'); % Correct the y-axis direction
colorbar; % Add a colorbar colorbar; % Add a colorbar
xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex'); xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex');
ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex'); ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex');
title('Roton instability boundary','fontsize',16,'interpreter','latex')
%% Roton instability boundary for tilted dipoles %{
imagesc(AtomNumber*1E-5, ScatteringLengths, instability_boundary); % Specify x and y data for axes
wz = 2 * pi * 72.4; % Trap frequency in the tight confinement direction
lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length
gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength
add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
gdd = VacuumPermeability*DyMagneticMoment^2/3;
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));
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
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
% ====================================================================================================================================================== %
alpha = 45; % Polar angle of dipole moment
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));
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(k, gs, gdd, MeanWidth, alpha, 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
nadd2s_from_figure = [0.04974, 0.05383, 0.05655, 0.06609, 0.06916, 0.07291, 0.07836, 0.08517, 0.09063, 0.0978, 0.10459, 0.11345, 0.11822, 0.12231, 0.12674, 0.13117, 0.13560, 0.14003, 0.14548, 0.15127, 0.15775, 0.16660, 0.17546, 0.18364, 0.19557, 0.20579, 0.21839, 0.23850, 0.25144];
as_to_add_from_figure = [0.76383, 0.76766, 0.76974, 0.77543, 0.77675, 0.77828, 0.78003, 0.78178, 0.78288, 0.7840, 0.78474, 0.78540, 0.78562, 0.78572, 0.78583, 0.78583, 0.78583, 0.78583, 0.78567, 0.78551, 0.78529, 0.78485, 0.78441, 0.78386, 0.78310, 0.78233, 0.78135, 0.77970, 0.77861];
figure(6)
clf
set(gcf,'Position',[50 50 950 750])
imagesc(nadd2s, as_to_add, instability_boundary); % Specify x and y data for axes
hold on
plot(nadd2s_from_figure, as_to_add_from_figure, 'r*-', 'LineWidth', 2); % Plot the curve (red line)
set(gca, 'YDir', 'normal'); % Correct the y-axis direction set(gca, 'YDir', 'normal'); % Correct the y-axis direction
colorbar; % Add a colorbar cbar1 = colorbar;
xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex'); cbar1.Label.Interpreter = 'latex';
ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex'); ylabel(cbar1,'$(\times 10^{-31})$','FontSize',16,'Rotation',270)
xlabel(' Atom number for a trap area of 100$\mu m^2 ~ (\times 10^5)$','fontsize',16,'interpreter','latex');
ylabel('Scattering length ($\times a_0$)','fontsize',16,'interpreter','latex');
title('Roton instability boundary','fontsize',16,'interpreter','latex') title('Roton instability boundary','fontsize',16,'interpreter','latex')
%}
%% %%
function [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, alpha, phi) function [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, alpha, phi)
@ -382,3 +330,6 @@ function ret = computeTotalEnergyPerParticle(x, as, AtomNumberDensity, wz, lz, g
Energy_TransverseComponent = (0.5 * (gs + (2*gdd)) * gamma4 * AtomNumberDensity) + ((2/5) * gQF * AtomNumberDensity^(3/2)); Energy_TransverseComponent = (0.5 * (gs + (2*gdd)) * gamma4 * AtomNumberDensity) + ((2/5) * gQF * AtomNumberDensity^(3/2));
ret = (Energy_AxialComponent + Energy_TransverseComponent) / (PlanckConstantReduced * wz); ret = (Energy_AxialComponent + Energy_TransverseComponent) / (PlanckConstantReduced * wz);
end end