diff --git a/Quasi2DBogoliubovSpectrum.m b/Quasi2DBogoliubovSpectrum.m index 64810b1..0b27b1b 100644 --- a/Quasi2DBogoliubovSpectrum.m +++ b/Quasi2DBogoliubovSpectrum.m @@ -24,9 +24,9 @@ HartreeEnergy = ElectronCharge^2 / (4 * pi * VacuumPermittivity * AtomicUnitOfPolarizability = (ElectronCharge^2 * BohrRadius^2) / HartreeEnergy; % Or simply 4*pi*VacuumPermittivity*BohrRadius^3 % Dy specific constants -Dy164Mass = 163.929174751*1.660539066E-27; +Dy164Mass = 163.929174751*AtomicMassUnit; Dy164IsotopicAbundance = 0.2826; -DyMagneticMoment = 9.93*9.274009994E-24; +DyMagneticMoment = 9.93*BohrMagneton; %% Bogoliubov excitation spectrum for quasi-2D dipolar gas with QF correction AtomNumber = 1E5; % Total atom number in the system @@ -37,7 +37,9 @@ Trapsize = 7.5815 * lz; 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 +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 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 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 +k = linspace(0, 2e6, 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) +ScatteringLengths = [108.5, 105.9, 103.3, 102.5150]; +eps_dds = zeros(1, length(ScatteringLengths)); +EpsilonKs = zeros(length(k), length(ScatteringLengths)); +for idx = 1:length(ScatteringLengths) - as = a * BohrRadius; % Scattering length + as = ScatteringLengths(idx) * 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; @@ -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)); EpsilonK = sqrt(((PlanckConstantReduced^2 .* k.^2) ./ (2 * Dy164Mass)) .* DeltaK); - ScatteringLengths(end+1) = as; - eps_dds(end+1) = eps_dd; - EpsilonKs(end+1,:) = EpsilonK; + eps_dds(idx) = eps_dd; + EpsilonKs(:,idx) = EpsilonK; end figure(2) +clf set(gcf,'Position',[50 50 950 750]) 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}$']) hold on 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}$']) end 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 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) Go = sqrt(pi) * (k * MeanWidth/sqrt(2)) .* exp((k * MeanWidth/sqrt(2)).^2) .* erfc((k * MeanWidth/sqrt(2)));