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

Plotting for different interaction strengths added.

Browse Source
This commit is contained in:
Karthik 2024-11-05 12:07:15 +01:00
parent cc446a99bb
commit cbf817e7e9
1 changed files with 55 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

56
Quasi2DBogoliubovSpectrum.m
View File

@ -66,8 +66,62 @@ xlabel('$k_{\rho}a_{dd}$','fontsize',16,'interpreter','latex')
ylabel('$\epsilon(k_{\rho})/h$ (Hz)','fontsize',16,'interpreter','latex')
grid on
%%
%% Bogoliubov excitation spectrum for quasi-2D dipolar gas with QF correction
AtomNumber = 1E5; % Total atom number in the system
wz = 2 * pi * 72.4; % Trap frequency in the tight confinement direction
lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length
Trapsize = 7.5815 * lz; % Trap is assumed to be a box of finite extent , given here in units of the harmonic oscillator length
alpha = 0; % Polar angle of dipole moment
phi = 0; % Azimuthal angle of momentum vector
MeanWidth = 5.7304888515 * lz; % Mean width of Gaussian ansatz
k = linspace(0, 3e6, 1000); % Vector of magnitudes of k vector
AtomNumberDensity = AtomNumber / Trapsize^2; % Areal density of atoms
add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
ScatteringLengths = [];
eps_dds = [];
EpsilonKs = [];
for a = linspace(131,102.515,5)
as = a * BohrRadius; % Scattering length
eps_dd = add/as; % Relative interaction strength
gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength
gdd = VacuumPermeability*DyMagneticMoment^2/3;
[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);
ScatteringLengths(end+1) = as;
eps_dds(end+1) = eps_dd;
EpsilonKs(end+1,:) = EpsilonK;
end
figure(2)
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(1/eps_dds(1),3)),'a_{dd}$'])
hold on
for idx = 2:length(ScatteringLengths)
yvals = EpsilonKs(idx, :) ./ PlanckConstant;
plot(xvals, yvals,LineWidth=2.0, DisplayName=['$a_s = ',num2str(round(1/eps_dds(idx),3)),'a_{dd}$'])
end
title(['$na_{dd}^2 = ',num2str(round(AtomNumberDensity * add^2,4)),'$'],'fontsize',16,'interpreter','latex')
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)
Go = sqrt(pi) * (k * MeanWidth/sqrt(2)) .* exp((k * MeanWidth/sqrt(2)).^2) .* erfc((k * MeanWidth/sqrt(2)));
gamma4 = 1/(sqrt(2*pi) * MeanWidth);