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

Added plotting of dispersion on the roton boundary.

Browse Source
This commit is contained in:
Karthik 2024-11-05 20:52:59 +01:00
parent cbf817e7e9
commit 5e130e2d93
1 changed files with 68 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

82
Quasi2DBogoliubovSpectrum.m
View File

@ -24,9 +24,9 @@ HartreeEnergy = ElectronCharge^2 / (4 * pi * VacuumPermittivity *
AtomicUnitOfPolarizability = (ElectronCharge^2 * BohrRadius^2) / HartreeEnergy; % Or simply 4*pi*VacuumPermittivity*BohrRadius^3 AtomicUnitOfPolarizability = (ElectronCharge^2 * BohrRadius^2) / HartreeEnergy; % Or simply 4*pi*VacuumPermittivity*BohrRadius^3
% Dy specific constants % Dy specific constants
Dy164Mass = 163.929174751*1.660539066E-27; Dy164Mass = 163.929174751*AtomicMassUnit;
Dy164IsotopicAbundance = 0.2826; Dy164IsotopicAbundance = 0.2826;
DyMagneticMoment = 9.93*9.274009994E-24; DyMagneticMoment = 9.93*BohrMagneton;
%% Bogoliubov excitation spectrum for quasi-2D dipolar gas with QF correction %% Bogoliubov excitation spectrum for quasi-2D dipolar gas with QF correction
AtomNumber = 1E5; % Total atom number in the system AtomNumber = 1E5; % Total atom number in the system
@ -37,7 +37,9 @@ Trapsize = 7.5815 * lz;
alpha = 0; % Polar angle of dipole moment alpha = 0; % Polar angle of dipole moment
phi = 0; % Azimuthal angle of momentum vector phi = 0; % Azimuthal angle of momentum vector
MeanWidth = 5.7304888515 * lz; % Mean width of Gaussian ansatz MeanWidth = 5.7304888515 * lz; % Mean width of Gaussian ansatz
k = linspace(0, 3e6, 1000); % Vector of magnitudes of k vector k = linspace(0, 2e6, 1000); % Vector of magnitudes of k vector
% no = 2.0429e+15, eps_dd = 1.2755, as = 5.4249e-09
AtomNumberDensity = AtomNumber / Trapsize^2; % Areal density of atoms AtomNumberDensity = AtomNumber / Trapsize^2; % Areal density of atoms
add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
@ -74,17 +76,17 @@ Trapsize = 7.5815 * lz;
alpha = 0; % Polar angle of dipole moment alpha = 0; % Polar angle of dipole moment
phi = 0; % Azimuthal angle of momentum vector phi = 0; % Azimuthal angle of momentum vector
MeanWidth = 5.7304888515 * lz; % Mean width of Gaussian ansatz MeanWidth = 5.7304888515 * lz; % Mean width of Gaussian ansatz
k = linspace(0, 3e6, 1000); % Vector of magnitudes of k vector k = linspace(0, 2e6, 1000); % Vector of magnitudes of k vector
AtomNumberDensity = AtomNumber / Trapsize^2; % Areal density of atoms AtomNumberDensity = AtomNumber / Trapsize^2; % Areal density of atoms
add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
ScatteringLengths = []; ScatteringLengths = [108.5, 105.9, 103.3, 102.5150];
eps_dds = []; eps_dds = zeros(1, length(ScatteringLengths));
EpsilonKs = []; EpsilonKs = zeros(length(k), length(ScatteringLengths));
for a = linspace(131,102.515,5) for idx = 1:length(ScatteringLengths)
as = a * BohrRadius; % Scattering length as = ScatteringLengths(idx) * BohrRadius; % Scattering length
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;
@ -100,19 +102,19 @@ for a = linspace(131,102.515,5)
DeltaK = ((PlanckConstantReduced^2 .* k.^2) ./ (2 * Dy164Mass)) + ((2 * AtomNumberDensity) .* Ukk) + (3 * gQF * AtomNumberDensity^(3/2)); DeltaK = ((PlanckConstantReduced^2 .* k.^2) ./ (2 * Dy164Mass)) + ((2 * AtomNumberDensity) .* Ukk) + (3 * gQF * AtomNumberDensity^(3/2));
EpsilonK = sqrt(((PlanckConstantReduced^2 .* k.^2) ./ (2 * Dy164Mass)) .* DeltaK); EpsilonK = sqrt(((PlanckConstantReduced^2 .* k.^2) ./ (2 * Dy164Mass)) .* DeltaK);
ScatteringLengths(end+1) = as; eps_dds(idx) = eps_dd;
eps_dds(end+1) = eps_dd; EpsilonKs(:,idx) = EpsilonK;
EpsilonKs(end+1,:) = EpsilonK;
end end
figure(2) figure(2)
clf
set(gcf,'Position',[50 50 950 750]) set(gcf,'Position',[50 50 950 750])
xvals = (k .* add); xvals = (k .* add);
yvals = EpsilonKs(1, :) ./ PlanckConstant; yvals = EpsilonKs(:, 1) ./ PlanckConstant;
plot(xvals, yvals,LineWidth=2.0, DisplayName=['$a_s = ',num2str(round(1/eps_dds(1),3)),'a_{dd}$']) plot(xvals, yvals,LineWidth=2.0, DisplayName=['$a_s = ',num2str(round(1/eps_dds(1),3)),'a_{dd}$'])
hold on hold on
for idx = 2:length(ScatteringLengths) for idx = 2:length(ScatteringLengths)
yvals = EpsilonKs(idx, :) ./ PlanckConstant; yvals = EpsilonKs(:, idx) ./ PlanckConstant;
plot(xvals, yvals,LineWidth=2.0, DisplayName=['$a_s = ',num2str(round(1/eps_dds(idx),3)),'a_{dd}$']) plot(xvals, yvals,LineWidth=2.0, DisplayName=['$a_s = ',num2str(round(1/eps_dds(idx),3)),'a_{dd}$'])
end end
title(['$na_{dd}^2 = ',num2str(round(AtomNumberDensity * add^2,4)),'$'],'fontsize',16,'interpreter','latex') title(['$na_{dd}^2 = ',num2str(round(AtomNumberDensity * add^2,4)),'$'],'fontsize',16,'interpreter','latex')
@ -121,6 +123,58 @@ ylabel('$\epsilon(k_{\rho})/h$ (Hz)','fontsize',16,'interpreter','latex')
grid on grid on
legend('location', 'northwest','fontsize',16, 'Interpreter','latex') legend('location', 'northwest','fontsize',16, 'Interpreter','latex')
%% Bogoliubov excitation spectrum for quasi-2D dipolar gas with QF correction
wz = 2 * pi * 72.4; % Trap frequency in the tight confinement direction
lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length
alpha = 0; % Polar angle of dipole moment
phi = 0; % Azimuthal angle of momentum vector
k = linspace(0, 2.25e6, 1000); % Vector of magnitudes of k vector
nadd2s = [0.0844, 0.0978, 0.123];
as_to_add = [0.7730, 0.7840, 0.7819];
var_widths = [4.97165, 5.72960, 5.93178];
add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
EpsilonKs = zeros(length(k), length(nadd2s));
for idx = 1:length(nadd2s)
AtomNumberDensity = nadd2s(idx) / add^2; % Areal density of atoms
as = (as_to_add(idx) * 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(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);
EpsilonKs(:,idx) = EpsilonK;
end
figure(3)
clf
set(gcf,'Position',[50 50 950 750])
xvals = (k .* add);
yvals = EpsilonKs(:, 1) ./ PlanckConstant;
plot(xvals, yvals,LineWidth=2.0, DisplayName=['$a_s = ',num2str(round(as_to_add(1),4)),'a_{dd}, na_{dd}^2 = ',num2str(round(nadd2s(1),4)),'$'])
hold on
for idx = 2:length(nadd2s)
yvals = EpsilonKs(:, idx) ./ PlanckConstant;
plot(xvals, yvals,LineWidth=2.0, DisplayName=['$a_s = ',num2str(round(as_to_add(idx),4)),'a_{dd}, na_{dd}^2 = ',num2str(round(nadd2s(idx),4)),'$'])
end
xlabel('$k_{\rho}a_{dd}$','fontsize',16,'interpreter','latex')
ylabel('$\epsilon(k_{\rho})/h$ (Hz)','fontsize',16,'interpreter','latex')
grid on
legend('location', 'northwest','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)
Go = sqrt(pi) * (k * MeanWidth/sqrt(2)) .* exp((k * MeanWidth/sqrt(2)).^2) .* erfc((k * MeanWidth/sqrt(2))); Go = sqrt(pi) * (k * MeanWidth/sqrt(2)) .* exp((k * MeanWidth/sqrt(2)).^2) .* erfc((k * MeanWidth/sqrt(2)));