From 7be9f6ddd8cd2cbff479880f8db932f70812d2af Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Fri, 10 Jan 2025 11:17:38 +0100 Subject: [PATCH] Removal of duplicate script --- ...larDispersionAndRotonInstabilityBoundary.m | 335 ------------------ 1 file changed, 335 deletions(-) delete mode 100644 Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersionAndRotonInstabilityBoundary.m diff --git a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersionAndRotonInstabilityBoundary.m b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersionAndRotonInstabilityBoundary.m deleted file mode 100644 index e40489f..0000000 --- a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersionAndRotonInstabilityBoundary.m +++ /dev/null @@ -1,335 +0,0 @@ -%% 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*AtomicMassUnit; -Dy164IsotopicAbundance = 0.2826; -DyMagneticMoment = 9.93*BohrMagneton; - -%% 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, 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 -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 - -%% For different interaction strengths - -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, 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 = [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 = 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; - - [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); - - 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; -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') - -%% For 3 points on the roton instability boundary - -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.7296048721, 5.93178]; - -add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length -EpsilonKs = zeros(length(k), length(nadd2s)); -ScatteringLengths = zeros(length(as_to_add), 1); -AtomNumber = zeros(length(nadd2s), 1); -w0 = 2 * pi * 61.6316; % Trap frequency in the tight confinement direction -l0 = sqrt(PlanckConstantReduced/(Dy164Mass * w0)); % Defining a harmonic oscillator length -tsize = 10 * l0; - -for idx = 1:length(nadd2s) - AtomNumberDensity = nadd2s(idx) / add^2; % Areal density of atoms - AtomNumber(idx) = AtomNumberDensity*tsize^2; - as = (as_to_add(idx) * add); % Scattering length - ScatteringLengths(idx) = as/BohrRadius; - 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') - -%% Mean widths of the variational Gaussian ansatz - extremize the total mean field energy per particle wrt to the variational parameter - -wz = 2 * pi * 72.4; % Trap frequency in the tight confinement direction -lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length -add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length -gdd = VacuumPermeability*DyMagneticMoment^2/3; -AtomNumberDensity = 0.0978 / add^2; -as = 0.784 * add; % Scattering length -gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength -TotalEnergyPerParticle = @(x) computeTotalEnergyPerParticle(x, as, AtomNumberDensity, wz, lz, gs, add, gdd, PlanckConstantReduced); - -x0 = 5; -Aineq = []; -Bineq = []; -Aeq = []; -Beq = []; -lb = [1]; -ub = [7]; -nonlcon = []; -fminconopts = optimoptions(@fmincon,'Display','off', 'StepTolerance', 1.0000e-11, 'MaxIterations',1500); -sigma = fmincon(TotalEnergyPerParticle, x0, Aineq, Bineq, Aeq, Beq, lb, ub, nonlcon, fminconopts); -fprintf(['Variational width of Gaussian ansatz = ' num2str(sigma) ' * lz \n']) - -%% Mean widths of the variational Gaussian ansatz - extremize the total mean field energy per particle wrt to the variational parameter - -wz = 2 * pi * 72.4; % Trap frequency in the tight confinement direction -lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length -add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length -gdd = VacuumPermeability*DyMagneticMoment^2/3; - -nadd2s = 0.05:0.001:0.25; -as_to_add = 0.74:0.001:0.79; -var_widths = zeros(length(as_to_add), length(nadd2s)); - -x0 = 5; -Aineq = []; -Bineq = []; -Aeq = []; -Beq = []; -lb = [1]; -ub = [10]; -nonlcon = []; -fminconopts = optimoptions(@fmincon,'Display','off', 'StepTolerance', 1.0000e-11, 'MaxIterations',1500); - -for idx = 1:length(nadd2s) - for jdx = 1:length(as_to_add) - AtomNumberDensity = nadd2s(idx) / add^2; % Areal density of atoms - as = (as_to_add(jdx) * add); % Scattering length - gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength - TotalEnergyPerParticle = @(x) computeTotalEnergyPerParticle(x, as, AtomNumberDensity, wz, lz, gs, add, gdd, PlanckConstantReduced); - sigma = fmincon(TotalEnergyPerParticle, x0, Aineq, Bineq, Aeq, Beq, lb, ub, nonlcon, fminconopts); - var_widths(jdx, idx) = sigma; - end -end - -figure(4) -clf -set(gcf,'Position',[50 50 950 750]) -imagesc(nadd2s, as_to_add, var_widths); % Specify x and y data for axes -set(gca, 'YDir', 'normal'); % Correct the y-axis direction -colorbar; % Add a colorbar -xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex'); -ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex'); - -% ====================================================================================================================================================== % - -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 -instability_boundary = zeros(length(as_to_add), length(nadd2s)); -ScatteringLengths = zeros(length(as_to_add), 1); -AtomNumber = zeros(length(nadd2s), 1); -w0 = 2 * pi * 61.6316; % Trap frequency in the tight confinement direction -l0 = sqrt(PlanckConstantReduced/(Dy164Mass * w0)); % Defining a harmonic oscillator length -tsize = 10 * l0; - -for idx = 1:length(nadd2s) - for jdx = 1:length(as_to_add) - AtomNumberDensity = nadd2s(idx) / add^2; % Areal density of atoms - AtomNumber(idx) = AtomNumberDensity*tsize^2; - as = (as_to_add(jdx) * add); % Scattering length - ScatteringLengths(jdx) = as/BohrRadius; - 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(jdx, 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); - instability_boundary(jdx, idx) = ~isreal(EpsilonK); - end -end - -nadd2s_from_figure = [0.04974, 0.05383, 0.05655, 0.06609, 0.06916, 0.07291, 0.07836, 0.08517, 0.09063, 0.0978, 0.10459, 0.11345, 0.11822, 0.12231, 0.12674, 0.13117, 0.13560, 0.14003, 0.14548, 0.15127, 0.15775, 0.16660, 0.17546, 0.18364, 0.19557, 0.20579, 0.21839, 0.23850, 0.25144]; -as_to_add_from_figure = [0.76383, 0.76766, 0.76974, 0.77543, 0.77675, 0.77828, 0.78003, 0.78178, 0.78288, 0.7840, 0.78474, 0.78540, 0.78562, 0.78572, 0.78583, 0.78583, 0.78583, 0.78583, 0.78567, 0.78551, 0.78529, 0.78485, 0.78441, 0.78386, 0.78310, 0.78233, 0.78135, 0.77970, 0.77861]; - -figure(5) -clf -set(gcf,'Position',[50 50 950 750]) - - -imagesc(nadd2s, as_to_add, instability_boundary); % Specify x and y data for axes -hold on -plot(nadd2s_from_figure, as_to_add_from_figure, 'r*-', 'LineWidth', 2); % Plot the curve (red line) -set(gca, 'YDir', 'normal'); % Correct the y-axis direction -colorbar; % Add a colorbar -xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex'); -ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex'); - -%{ -imagesc(AtomNumber*1E-5, ScatteringLengths, instability_boundary); % Specify x and y data for axes -set(gca, 'YDir', 'normal'); % Correct the y-axis direction -cbar1 = colorbar; -cbar1.Label.Interpreter = 'latex'; -ylabel(cbar1,'$(\times 10^{-31})$','FontSize',16,'Rotation',270) -xlabel(' Atom number for a trap area of 100$\mu m^2 ~ (\times 10^5)$','fontsize',16,'interpreter','latex'); -ylabel('Scattering length ($\times a_0$)','fontsize',16,'interpreter','latex'); -title('Roton instability boundary','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); - 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 - -function ret = computeTotalEnergyPerParticle(x, as, AtomNumberDensity, wz, lz, gs, add, gdd, PlanckConstantReduced) - eps_dd = add/as; % Relative interaction strength - MeanWidth = x * lz; - gammaQF = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2)); % Quantum Fluctuations term - gamma4 = 1/(sqrt(2*pi) * MeanWidth); - gamma5 = sqrt(2/5) / (sqrt(pi) * MeanWidth)^(3/2); - gQF = gamma5 * gammaQF; - Energy_AxialComponent = (PlanckConstantReduced * wz) * ((lz^2/(4 * MeanWidth^2)) + (MeanWidth^2/(4 * lz^2))); - Energy_TransverseComponent = (0.5 * (gs + (2*gdd)) * gamma4 * AtomNumberDensity) + ((2/5) * gQF * AtomNumberDensity^(3/2)); - ret = (Energy_AxialComponent + Energy_TransverseComponent) / (PlanckConstantReduced * wz); -end - - -