Calculations/Quasi2DBogoliubovSpectrum.m

%% Physical constants
PlanckConstant               = 6.62607015E-34;
PlanckConstantReduced        = 6.62607015E-34/(2*pi);
FineStructureConstant        = 7.2973525698E-3;
ElectronMass                 = 9.10938291E-31;
GravitationalConstant        = 6.67384E-11;
ProtonMass                   = 1.672621777E-27;
AtomicMassUnit               = 1.660539066E-27; 
BohrRadius                   = 5.2917721067E-11; 
BohrMagneton                 = 9.274009994E-24; 
BoltzmannConstant            = 1.38064852E-23;
StandardGravityAcceleration  = 9.80665;
SpeedOfLight                 = 299792458;
StefanBoltzmannConstant      = 5.670373E-8;
ElectronCharge               = 1.602176634E-19;
VacuumPermeability           = 1.25663706212E-6; 
DielectricConstant           = 8.8541878128E-12;
ElectronGyromagneticFactor   = -2.00231930436153;
AvogadroConstant             = 6.02214076E23;
ZeroKelvin                   = 273.15;
GravitationalAcceleration    = 9.80553;
VacuumPermittivity           = 1 / (SpeedOfLight^2 * VacuumPermeability);
HartreeEnergy                = ElectronCharge^2 / (4 * pi * VacuumPermittivity * BohrRadius);
AtomicUnitOfPolarizability   = (ElectronCharge^2 * BohrRadius^2) / HartreeEnergy; % Or simply 4*pi*VacuumPermittivity*BohrRadius^3

% Dy specific constants
Dy164Mass                    = 163.929174751*1.660539066E-27;
Dy164IsotopicAbundance       = 0.2826;
DyMagneticMoment             = 9.93*9.274009994E-24;

%% 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
as                = 102.515 * BohrRadius;                                                            % Scattering 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
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); 

figure(1)
set(gcf,'Position',[50 50 950 750])
xvals = (k .* add);
yvals = EpsilonK ./ PlanckConstant;
plot(xvals, yvals,LineWidth=2.0)
title(horzcat(['$a_s = ',num2str(round(1/eps_dd,3)),'a_{dd}, '], ['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

%%

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);
   Fka                       = (3 * cos(deg2rad(alpha))^2 - 1) + ((3 * Go) .* ((sin(deg2rad(alpha))^2 .* sin(deg2rad(phi))^2) - cos(deg2rad(alpha))^2));
   Ukk                       = (gs + (gdd * Fka)) * gamma4;
end