diff --git a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersion2D.m b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersion2D.m new file mode 100644 index 0000000..7113df3 --- /dev/null +++ b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersion2D.m @@ -0,0 +1,290 @@ +%% 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; + +%% 2-D DDI Potential in k-space, with Gaussian ansatz width determined by constrained minimization + +wz = 2 * pi * 300; % Trap frequency in the tight confinement direction +lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length + +% Number of grid points in each direction +Params.Nx = 128; +Params.Ny = 128; + +% Dimensions (in units of l0) +w0 = 2*pi*61.6316; % Angular frequency unit [s^-1] +l0 = sqrt(PlanckConstantReduced/(Dy164Mass*w0)); % Defining a harmonic oscillator length +Params.Lx = 150*l0; +Params.Ly = 150*l0; +[Transf] = setupSpace(Params); + +nadd2s = 0.05:0.001:0.25; +as_to_add = 0.74:0.001:0.79; + +add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length +gdd = VacuumPermeability*DyMagneticMoment^2/3; + +%% +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(10) +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'); + +%% Chosen values of interaction, density and tilt + +nadd = 0.110; +asadd = 0.782; +Params.alpha = 0; +Params.phi = 0; + +[~, nadd2sidx] = min(abs(nadd2s - nadd)); +[~, asaddidx] = min(abs(as_to_add - asadd)); + +AtomNumberDensity = nadd2s(nadd2sidx) / add^2; % Areal density of atoms +as = (as_to_add(asaddidx) * 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(asaddidx, nadd2sidx)*lz; + +% == 2-D DDI Potential in k-space == % +VDk = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth); +VDk_fftshifted = fftshift(VDk); + +figure(11) +clf +set(gcf,'Position',[50 50 950 750]) +imagesc(fftshift(Transf.kx)*l0, fftshift(Transf.ky)*l0, VDk_fftshifted); % Specify x and y data for axes +set(gca, 'YDir', 'normal'); % Correct the y-axis direction +cbar1 = colorbar; +cbar1.Label.Interpreter = 'latex'; +xlabel('$k_x l_o$','fontsize',16,'interpreter','latex'); +ylabel('$k_y l_o$','fontsize',16,'interpreter','latex'); +title(['$\theta = ',num2str(Params.alpha), '; \phi = ', num2str(Params.phi),'$'],'fontsize',16,'interpreter','latex') + +% == 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; + +EpsilonK = zeros(length(Transf.ky), length(Transf.kx)); +gs_tilde = gs / (sqrt(2*pi) * MeanWidth); + +% == Dispersion relation == % +for idx = 1:length(Transf.kx) + for jdx = 1:length(Transf.ky) + DeltaK = ((PlanckConstantReduced^2 .* (Transf.kx(idx).^2 + Transf.ky(jdx).^2)) ./ (2 * Dy164Mass)) + (2 * AtomNumberDensity * gs_tilde) + ((2 * AtomNumberDensity) .* VDk_fftshifted(jdx, idx)) + (3 * gQF * AtomNumberDensity^(3/2)); + EpsilonK(jdx, idx) = sqrt(((PlanckConstantReduced^2 .* (Transf.kx(idx).^2 + Transf.ky(jdx).^2)) ./ (2 * Dy164Mass)) .* DeltaK); + end +end + +EpsilonK = double(imag(EpsilonK) ~= 0); % 'isreal' returns 0 for complex numbers and 1 for real numbers, so we negate it + +figure(12) +clf +set(gcf,'Position',[50 50 950 750]) +imagesc(fftshift(Transf.kx)*l0, fftshift(Transf.ky)*l0, EpsilonK); % Specify x and y data for axes +set(gca, 'YDir', 'normal'); % Correct the y-axis direction +cbar1 = colorbar; +cbar1.Label.Interpreter = 'latex'; +xlabel('$k_x l_o$','fontsize',16,'interpreter','latex'); +ylabel('$k_y l_o$','fontsize',16,'interpreter','latex'); +title(['$\theta = ',num2str(Params.alpha), '; \phi = ', num2str(Params.phi),'$'],'fontsize',16,'interpreter','latex') + +%% + +nadd = 0.110; +asadd = 0.782; +% Define values for alpha and phi +alpha_values = 0:5:90; % Range of alpha values (you can modify this) +phi_values = 0:2:90; % Range of phi values (you can modify this) + +[~, nadd2sidx] = min(abs(nadd2s - nadd)); +[~, asaddidx] = min(abs(as_to_add - asadd)); + +AtomNumberDensity = nadd2s(nadd2sidx) / add^2; % Areal density of atoms +as = (as_to_add(asaddidx) * 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(asaddidx, nadd2sidx)*lz; + +% Set up VideoWriter object +v = VideoWriter('potential_movie', 'MPEG-4'); % Create a video object +v.FrameRate = 5; % Frame rate of the video +open(v); % Open the video file + +% Loop over alpha and phi values +for alpha = alpha_values + for phi = phi_values + % Update Params with current alpha and phi + Params.alpha = alpha; + Params.phi = phi; + + % Compute the potential for the current alpha and phi + % == 2-D DDI Potential in k-space == % + VDk = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth); + VDk_fftshifted = fftshift(VDk); + + % == 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; + + EpsilonK = zeros(length(Transf.ky), length(Transf.kx)); + gs_tilde = gs / (sqrt(2*pi) * MeanWidth); + + % == Dispersion relation == % + for idx = 1:length(Transf.kx) + for jdx = 1:length(Transf.ky) + DeltaK = ((PlanckConstantReduced^2 .* (Transf.kx(idx).^2 + Transf.ky(jdx).^2)) ./ (2 * Dy164Mass)) + (2 * AtomNumberDensity * gs_tilde) + ((2 * AtomNumberDensity) .* VDk_fftshifted(jdx, idx)) + (3 * gQF * AtomNumberDensity^(3/2)); + EpsilonK(jdx, idx) = sqrt(((PlanckConstantReduced^2 .* (Transf.kx(idx).^2 + Transf.ky(jdx).^2)) ./ (2 * Dy164Mass)) .* DeltaK); + end + end + + EpsilonK = double(imag(EpsilonK) ~= 0); % 'isreal' returns 0 for complex numbers and 1 for real numbers, so we negate it + + % Plot the result + figure(12) + clf + set(gcf,'Position',[50 50 950 750]) + imagesc(fftshift(Transf.kx)*l0, fftshift(Transf.ky)*l0, EpsilonK); % Specify x and y data for axes + set(gca, 'YDir', 'normal'); % Correct the y-axis direction + cbar1 = colorbar; + cbar1.Label.Interpreter = 'latex'; + xlabel('$k_x$','fontsize',16,'interpreter','latex'); + ylabel('$k_y$','fontsize',16,'interpreter','latex'); + title(['$\theta = ',num2str(Params.alpha), '; \phi = ', num2str(Params.phi),'$'],'fontsize',16,'interpreter','latex') + + % Capture the frame and write to video + frame = getframe(gcf); % Capture the current figure + writeVideo(v, frame); % Write the frame to the video + end +end + +% Close the video file +close(v); + +%% +function [Transf] = setupSpace(Params) + + Transf.Xmax = 0.5*Params.Lx; + Transf.Ymax = 0.5*Params.Ly; + + Nx = Params.Nx; + Ny = Params.Ny; + + % Fourier grids + x = linspace(-0.5*Params.Lx,0.5*Params.Lx-Params.Lx/Params.Nx,Params.Nx); + Kmax = pi*Params.Nx/Params.Lx; + kx = linspace(-Kmax,Kmax,Nx+1); + kx = kx(1:end-1); + dkx = kx(2)-kx(1); + kx = fftshift(kx); + + y = linspace(-0.5*Params.Ly,0.5*Params.Ly-Params.Ly/Params.Ny,Params.Ny); + Kmax = pi*Params.Ny/Params.Ly; + ky = linspace(-Kmax,Kmax,Ny+1); + ky = ky(1:end-1); + dky = ky(2)-ky(1); + ky = fftshift(ky); + + [Transf.X,Transf.Y] = ndgrid(x,y); + [Transf.KX,Transf.KY] = ndgrid(kx,ky); + Transf.x = x; + Transf.y = y; + Transf.kx = kx; + Transf.ky = ky; + Transf.dx = x(2)-x(1); + Transf.dy = y(2)-y(1); + Transf.dkx = dkx; + Transf.dky = dky; +end + +function VDk = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth) +% == Calculating the DDI potential in Fourier space with appropriate cutoff == % + + % Interaction in K space + QX = Transf.KX*MeanWidth/sqrt(2); + QY = Transf.KY*MeanWidth/sqrt(2); + + Qsq = QX.^2 + QY.^2; + absQ = sqrt(Qsq); + QDsq = QX.^2*cos(Params.phi)^2 + QY.^2*sin(Params.phi)^2; + + % Bare interaction + Fpar = -1 + 3*sqrt(pi)*QDsq.*erfcx(absQ)./absQ; % Scaled complementary error function erfcx(x) = e^(x^2) * erfc(x) + Fperp = 2 - 3*sqrt(pi).*absQ.*erfcx(absQ); + Fpar(absQ == 0) = -1; + + % Full DDI + VDk = (Fpar*sin(Params.alpha)^2 + Fperp*cos(Params.alpha)^2); +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 + + + diff --git a/Estimations/ReproduceBlairBlakieResults.m b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersionAndRotonInstabilityBoundary.m similarity index 100% rename from Estimations/ReproduceBlairBlakieResults.m rename to Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersionAndRotonInstabilityBoundary.m diff --git a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ExtractingKRoton.m b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ExtractingKRoton.m new file mode 100644 index 0000000..073acf3 --- /dev/null +++ b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ExtractingKRoton.m @@ -0,0 +1,157 @@ +%% 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; + +%% k_roton at the instability boundary for tilted dipoles + +wz = 2 * pi * 500; % 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.2:0.005:0.75; +% as_to_add = 0.4:0.002:0.5; + +nadd2s = 0.05:0.005:0.25; +as_to_add = 0.50:0.001:0.80; + +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 + +% ====================================================================================================================================================== % + +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)); +k_roton = 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); + k_roton_indices = find(imag(EpsilonK) ~= 0); + if ~isempty(k_roton_indices) + k_roton(jdx, idx) = k(k_roton_indices(1)); + else + k_roton(jdx, idx) = NaN; + end + end +end +% +k_roton_vals = (k_roton .* add); +% +figure(8) +clf +set(gcf,'Position',[50 50 950 750]) +imagesc(AtomNumber*1E-5, ScatteringLengths, k_roton_vals); % Specify x and y data for axes +set(gca, 'YDir', 'normal'); % Correct the y-axis direction +cbar1 = colorbar; +cbar1.Label.Interpreter = 'latex'; +% ylabel(cbar1,'$$','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') +% +% Get the size of the matrix +k_roton_vals = flipud(k_roton_vals); +[rows, cols] = size(k_roton_vals); + +first_nonnan_row = zeros(1, cols); + +% Loop through each column +for col = 1:cols + nonnan_rows = find(~isnan(k_roton_vals(:, col))); + + if ~isempty(nonnan_rows) + first_nonnan_row(col) = nonnan_rows(1); + else + first_nonnan_row(col) = NaN; % Use NaN to represent no non-zero elements in this column + end +end + +% Create column indices (1 to number of columns) +column_indices = 1:cols; +% +% Use row and column indices to extract the first non-zero elements +k_roton_instability_boundary = arrayfun(@(r, c) k_roton_vals(r, c), first_nonnan_row(~isnan(first_nonnan_row)), column_indices(~isnan(first_nonnan_row))); + +figure(9) +clf +set(gcf,'Position',[50 50 950 750]) +xvals = AtomNumber*1E-5; +yvals = k_roton_instability_boundary; +plot(xvals', yvals,LineWidth=2.0) +xlabel(' Atom number for a trap area of 100$\mu m^2 ~ (\times 10^5)$','fontsize',16,'interpreter','latex'); +ylabel('$k_{\rho}a_{dd}$','fontsize',16,'interpreter','latex') +title('$k_{roton}$ at the instability boundary','fontsize',16,'interpreter','latex') +grid on \ No newline at end of file diff --git a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/RIBForTiltedDipoles.m b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/RIBForTiltedDipoles.m new file mode 100644 index 0000000..c967f59 --- /dev/null +++ b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/RIBForTiltedDipoles.m @@ -0,0 +1,149 @@ +%% 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; + +%% Roton instability boundary for tilted dipoles + +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.76:0.001:0.81; +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 + +% ====================================================================================================================================================== % + +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(6) +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 + + + diff --git a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ReproduceBlairBlakieResults.m b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ReproduceBlairBlakieResults.m new file mode 100644 index 0000000..e40489f --- /dev/null +++ b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ReproduceBlairBlakieResults.m @@ -0,0 +1,335 @@ +%% 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 + + + diff --git a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ScalingOfTheQFTerm.m b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ScalingOfTheQFTerm.m new file mode 100644 index 0000000..36bca94 --- /dev/null +++ b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ScalingOfTheQFTerm.m @@ -0,0 +1,89 @@ +%% 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; + +%% Scaling of the QF term + +wz = 2 * pi * 72.4; % Trap frequency in the tight confinement direction +lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length +gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength +add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length +gdd = VacuumPermeability*DyMagneticMoment^2/3; + +nadd2s = 0.05:0.01:0.25; +as_to_add = 0.76:0.01:0.81; + +QF = 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; + +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 + AtomNumber(idx) = AtomNumberDensity*tsize^2; + as = (as_to_add(jdx) * add); % Scattering length + gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength + ScatteringLengths(jdx) = as/BohrRadius; + TotalEnergyPerParticle = @(x) computeTotalEnergyPerParticle(x, as, AtomNumberDensity, wz, lz, gs, add, gdd, PlanckConstantReduced); + sigma = fmincon(TotalEnergyPerParticle, x0, Aineq, Bineq, Aeq, Beq, lb, ub, nonlcon, fminconopts); + eps_dd = add/as; % Relative interaction strength + + % == Quantum Fluctuations term == % + MeanWidth = sigma * lz; + 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; + QF(jdx, idx) = 3 * gQF * AtomNumberDensity^(3/2); + end +end + +figure(7) +clf +set(gcf,'Position',[50 50 950 750]) +imagesc(AtomNumber*1E-5, ScatteringLengths, QF * 1E31); % 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('Scaling of the quantum fluctuations term','fontsize',16,'interpreter','latex') \ No newline at end of file diff --git a/Estimations/roton_instability_project.json b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/roton_instability_project.json similarity index 100% rename from Estimations/roton_instability_project.json rename to Estimations/DipolarDispersionAndRotonInstabilityBoundary/roton_instability_project.json diff --git a/Estimations/roton_instability_project.tar b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/roton_instability_project.tar similarity index 100% rename from Estimations/roton_instability_project.tar rename to Estimations/DipolarDispersionAndRotonInstabilityBoundary/roton_instability_project.tar