%% 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; % 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 * lz)^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