karthik/Calculations
1
0
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity

Clean up, minor corrections to keep notation consistent.

Browse Source
This commit is contained in:
Karthik 2025-01-14 17:30:45 +01:00
parent 6093c700ab
commit 0833bddee8
5 changed files with 453 additions and 62 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

83
Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersion2D.m
View File

@ -29,9 +29,12 @@ Dy164IsotopicAbundance = 0.2826;
DyMagneticMoment = 9.93*BohrMagneton; DyMagneticMoment = 9.93*BohrMagneton;
%% 2-D DDI Potential in k-space, with Gaussian ansatz width determined by constrained minimization %% 2-D DDI Potential in k-space, with Gaussian ansatz width determined by constrained minimization
% With user-defined values of interaction, density and tilt
w0 = 2*pi*61.6316; % Angular frequency unit [s^-1] % w0 = 2*pi*61.6316; % Angular frequency unit [s^-1]
l0 = sqrt(PlanckConstantReduced/(Dy164Mass*w0)); % Defining a harmonic oscillator length % l0 = sqrt(PlanckConstantReduced/(Dy164Mass*w0));
% % Defining a harmonic oscillator length - heredue to the choice of w0, l0
% is 1 micrometer
wz = 2 * pi * 300; % Trap frequency in the tight confinement direction wz = 2 * pi * 300; % Trap frequency in the tight confinement direction
lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length
@ -39,14 +42,14 @@ lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz));
% Number of grid points in each direction % Number of grid points in each direction
Params.Nx = 128; Params.Nx = 128;
Params.Ny = 128; Params.Ny = 128;
Params.Lx = 150*l0; Params.Lx = 150*1e-6;
Params.Ly = 150*l0; Params.Ly = 150*1e-6;
[Transf] = setupSpace(Params); [Transf] = setupSpace(Params);
nadd2s = 0.110; nadd2s = 0.110; % Number density * add^2
as_to_add = 0.782; as_to_add = 0.782; % 1/edd
Params.alpha = 10; Params.theta = 57; % Polar angle of dipole moment
Params.phi = 0; Params.eta = 0; % Azimuthal angle of dipole moment
add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
gdd = VacuumPermeability*DyMagneticMoment^2/3; gdd = VacuumPermeability*DyMagneticMoment^2/3;
@ -61,15 +64,13 @@ ub = [10];
nonlcon = []; nonlcon = [];
fminconopts = optimoptions(@fmincon,'Display','off', 'StepTolerance', 1.0000e-11, 'MaxIterations',1500); fminconopts = optimoptions(@fmincon,'Display','off', 'StepTolerance', 1.0000e-11, 'MaxIterations',1500);
AtomNumberDensity = nadd2s / add^2; % Areal density of atoms AtomNumberDensity = nadd2s / add^2; % Number density of atoms
as = as_to_add * add; % Scattering length as = as_to_add * add; % Scattering length
eps_dd = add/as; % Relative interaction strength eps_dd = add/as; % Relative interaction strength
gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength
TotalEnergyPerParticle = @(x) computeTotalEnergyPerParticle(x, as, AtomNumberDensity, wz, lz, gs, add, gdd, PlanckConstantReduced); TotalEnergyPerParticle = @(x) computeTotalEnergyPerParticle(x, as, AtomNumberDensity, wz, lz, gs, add, gdd, PlanckConstantReduced);
sigma = fmincon(TotalEnergyPerParticle, x0, Aineq, Bineq, Aeq, Beq, lb, ub, nonlcon, fminconopts); sigma = fmincon(TotalEnergyPerParticle, x0, Aineq, Bineq, Aeq, Beq, lb, ub, nonlcon, fminconopts);
%% Chosen values of interaction, density and tilt
MeanWidth = sigma * lz; MeanWidth = sigma * lz;
% == 2-D DDI Potential in k-space == % % == 2-D DDI Potential in k-space == %
@ -79,13 +80,13 @@ VDk_fftshifted = fftshift(VDk);
figure(11) figure(11)
clf clf
set(gcf,'Position',[50 50 950 750]) 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 imagesc(fftshift(Transf.kx)*1e-6, fftshift(Transf.ky)*1e-6, VDk_fftshifted); % Specify x and y data for axes
set(gca, 'YDir', 'normal'); % Correct the y-axis direction set(gca, 'YDir', 'normal'); % Correct the y-axis direction
cbar1 = colorbar; cbar1 = colorbar;
cbar1.Label.Interpreter = 'latex'; cbar1.Label.Interpreter = 'latex';
xlabel('$k_x l_o$','fontsize',16,'interpreter','latex'); xlabel('$k_x l_o$','fontsize',16,'interpreter','latex');
ylabel('$k_y 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') title(['2-D DDI Potential: $\theta = ',num2str(Params.theta), '; \eta = ', num2str(Params.eta),'$'],'fontsize',16,'interpreter','latex')
% == Quantum Fluctuations term == % % == Quantum Fluctuations term == %
gammaQF = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2)); gammaQF = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2));
@ -108,33 +109,33 @@ EpsilonK = double(imag(EpsilonK) ~= 0); % 'is
figure(12) figure(12)
clf clf
set(gcf,'Position',[50 50 950 750]) set(gcf,'Position',[50 50 950 750])
imagesc(fftshift(Transf.kx)*l0, fftshift(Transf.ky)*l0, EpsilonK); % Specify x and y data for axes imagesc(fftshift(Transf.kx)*1e-6, fftshift(Transf.ky)*1e-6, EpsilonK); % Specify x and y data for axes
set(gca, 'YDir', 'normal'); % Correct the y-axis direction set(gca, 'YDir', 'normal'); % Correct the y-axis direction
cbar1 = colorbar; cbar1 = colorbar;
cbar1.Label.Interpreter = 'latex'; cbar1.Label.Interpreter = 'latex';
xlabel('$k_x l_o$','fontsize',16,'interpreter','latex'); xlabel('$k_x l_o$','fontsize',16,'interpreter','latex');
ylabel('$k_y 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') title(['2-D Dispersion: $\theta = ',num2str(Params.theta), '; \eta = ', num2str(Params.eta),'$'],'fontsize',16,'interpreter','latex')
%% %% Cycle through angles
% Define values for alpha and phi % Define values for theta and eta
alpha_values = 0:5:90; % Range of alpha values (you can modify this) theta_values = 0:10:90; % Range of theta values (you can modify this)
phi_values = 0:2:90; % Range of phi values (you can modify this) eta_values = 0:10:90; % Range of eta values (you can modify this)
% Set up VideoWriter object % Set up VideoWriter object to produce a movie
v = VideoWriter('potential_movie', 'MPEG-4'); % Create a video object % v = VideoWriter('potential_movie', 'MPEG-4'); % Create a video object
v.FrameRate = 5; % Frame rate of the video % v.FrameRate = 5; % Frame rate of the video
open(v); % Open the video file % open(v); % Open the video file
% Loop over alpha and phi values % Loop over theta and eta values
for alpha = alpha_values for theta = theta_values
for phi = phi_values for eta = eta_values
% Update Params with current alpha and phi % Update Params with current theta and eta
Params.alpha = alpha; Params.theta = theta;
Params.phi = phi; Params.eta = eta;
% Compute the potential for the current alpha and phi % Compute the potential for the current theta and eta
% == 2-D DDI Potential in k-space == % % == 2-D DDI Potential in k-space == %
VDk = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth); VDk = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth);
VDk_fftshifted = fftshift(VDk); VDk_fftshifted = fftshift(VDk);
@ -158,25 +159,25 @@ for alpha = alpha_values
EpsilonK = double(imag(EpsilonK) ~= 0); % 'isreal' returns 0 for complex numbers and 1 for real numbers, so we negate it EpsilonK = double(imag(EpsilonK) ~= 0); % 'isreal' returns 0 for complex numbers and 1 for real numbers, so we negate it
% Plot the result % Plot the result
figure(12) figure(13)
clf clf
set(gcf,'Position',[50 50 950 750]) set(gcf,'Position',[50 50 950 750])
imagesc(fftshift(Transf.kx)*l0, fftshift(Transf.ky)*l0, EpsilonK); % Specify x and y data for axes imagesc(fftshift(Transf.kx)*1e-6, fftshift(Transf.ky)*1e-6, EpsilonK); % Specify x and y data for axes
set(gca, 'YDir', 'normal'); % Correct the y-axis direction set(gca, 'YDir', 'normal'); % Correct the y-axis direction
cbar1 = colorbar; cbar1 = colorbar;
cbar1.Label.Interpreter = 'latex'; cbar1.Label.Interpreter = 'latex';
xlabel('$k_x$','fontsize',16,'interpreter','latex'); xlabel('$k_x l_o$','fontsize',16,'interpreter','latex');
ylabel('$k_y$','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') title(['$\theta = ',num2str(Params.theta), '; \eta = ', num2str(Params.eta),'$'],'fontsize',16,'interpreter','latex')
% Capture the frame and write to video % Capture the frame and write to video
frame = getframe(gcf); % Capture the current figure % frame = getframe(gcf); % Capture the current figure
writeVideo(v, frame); % Write the frame to the video % writeVideo(v, frame); % Write the frame to the video
end end
end end
% Close the video file % Close the video file
close(v); % close(v);
%% %%
function [Transf] = setupSpace(Params) function [Transf] = setupSpace(Params)
@ -216,14 +217,14 @@ end
function VDk = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth) function VDk = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth)
% == Calculating the DDI potential in Fourier space with appropriate cutoff == % % == Calculating the DDI potential in Fourier space with appropriate cutoff == %
% Angles of the dipole moment are defined in and away from the X-Z plane
% Interaction in K space % Interaction in K space
QX = Transf.KX*MeanWidth/sqrt(2); QX = Transf.KX*MeanWidth/sqrt(2);
QY = Transf.KY*MeanWidth/sqrt(2); QY = Transf.KY*MeanWidth/sqrt(2);
Qsq = QX.^2 + QY.^2; Qsq = QX.^2 + QY.^2;
absQ = sqrt(Qsq); absQ = sqrt(Qsq);
QDsq = QX.^2*cos(Params.phi)^2 + QY.^2*sin(Params.phi)^2; QDsq = QX.^2*cos(Params.eta)^2 + QY.^2*sin(Params.eta)^2; % eta here is the azimuthal angle of the dipole moment (angle from the x-axis)
% Bare interaction % Bare interaction
Fpar = -1 + 3*sqrt(pi)*QDsq.*erfcx(absQ)./absQ; % Scaled complementary error function erfcx(x) = e^(x^2) * erfc(x) Fpar = -1 + 3*sqrt(pi)*QDsq.*erfcx(absQ)./absQ; % Scaled complementary error function erfcx(x) = e^(x^2) * erfc(x)
@ -231,7 +232,7 @@ function VDk = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth)
Fpar(absQ == 0) = -1; Fpar(absQ == 0) = -1;
% Full DDI % Full DDI
VDk = (Fpar*sin(Params.alpha)^2 + Fperp*cos(Params.alpha)^2); VDk = (Fpar*sin(Params.theta)^2 + Fperp*cos(Params.theta)^2); % theta here is the polar angle of the dipole moment (angle from the z-axis)
end end
function ret = computeTotalEnergyPerParticle(x, as, AtomNumberDensity, wz, lz, gs, add, gdd, PlanckConstantReduced) function ret = computeTotalEnergyPerParticle(x, as, AtomNumberDensity, wz, lz, gs, add, gdd, PlanckConstantReduced)

81
Estimations/DipolarDispersionAndRotonInstabilityBoundary/ExtractingKRoton.m
View File

@ -66,7 +66,7 @@ end
% ====================================================================================================================================================== % % ====================================================================================================================================================== %
alpha = 0; % Polar angle of dipole moment theta = 0; % Polar angle of dipole moment
phi = 0; % Azimuthal angle of momentum vector phi = 0; % Azimuthal angle of momentum vector
k = linspace(0, 2.25e6, 1000); % Vector of magnitudes of k vector k = linspace(0, 2.25e6, 1000); % Vector of magnitudes of k vector
instability_boundary = zeros(length(as_to_add), length(nadd2s)); instability_boundary = zeros(length(as_to_add), length(nadd2s));
@ -88,7 +88,7 @@ for idx = 1:length(nadd2s)
gdd = VacuumPermeability*DyMagneticMoment^2/3; gdd = VacuumPermeability*DyMagneticMoment^2/3;
MeanWidth = var_widths(jdx, idx) * lz; % Mean width of Gaussian ansatz 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 [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, theta, phi); % DDI potential in k-space
% == Quantum Fluctuations term == % % == Quantum Fluctuations term == %
gammaQF = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2)); gammaQF = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2));
@ -154,4 +154,79 @@ 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'); 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') ylabel('$k_{\rho}a_{dd}$','fontsize',16,'interpreter','latex')
title('$k_{roton}$ at the instability boundary','fontsize',16,'interpreter','latex') title('$k_{roton}$ at the instability boundary','fontsize',16,'interpreter','latex')
grid on grid on
%%
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 [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, theta, 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(theta))^2 - 1) + ((3 * Go) .* ((sin(deg2rad(theta))^2 .* sin(deg2rad(phi))^2) - cos(deg2rad(theta))^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

13
Estimations/DipolarDispersionAndRotonInstabilityBoundary/RIBForTiltedDipoles.m
View File

@ -62,7 +62,7 @@ end
% ====================================================================================================================================================== % % ====================================================================================================================================================== %
alpha = 0; % Polar angle of dipole moment theta = 0; % Polar angle of dipole moment
phi = 0; % Azimuthal angle of momentum vector phi = 0; % Azimuthal angle of momentum vector
k = linspace(0, 2.25e6, 1000); % Vector of magnitudes of k vector k = linspace(0, 2.25e6, 1000); % Vector of magnitudes of k vector
instability_boundary = zeros(length(as_to_add), length(nadd2s)); instability_boundary = zeros(length(as_to_add), length(nadd2s));
@ -83,7 +83,7 @@ for idx = 1:length(nadd2s)
gdd = VacuumPermeability*DyMagneticMoment^2/3; gdd = VacuumPermeability*DyMagneticMoment^2/3;
MeanWidth = var_widths(jdx, idx) * lz; % Mean width of Gaussian ansatz 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 [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, theta, phi); % DDI potential in k-space
% == Quantum Fluctuations term == % % == Quantum Fluctuations term == %
gammaQF = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2)); gammaQF = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2));
@ -126,10 +126,10 @@ ylabel('Scattering length ($\times a_0$)','fontsize',16,'interpreter','latex');
title('Roton instability boundary','fontsize',16,'interpreter','latex') title('Roton instability boundary','fontsize',16,'interpreter','latex')
%% %%
function [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, alpha, phi) function [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, theta, phi)
Go = sqrt(pi) * (k * MeanWidth/sqrt(2)) .* exp((k * MeanWidth/sqrt(2)).^2) .* erfc((k * MeanWidth/sqrt(2))); Go = sqrt(pi) * (k * MeanWidth/sqrt(2)) .* exp((k * MeanWidth/sqrt(2)).^2) .* erfc((k * MeanWidth/sqrt(2)));
gamma4 = 1/(sqrt(2*pi) * MeanWidth); 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)); Fka = (3 * cos(deg2rad(theta))^2 - 1) + ((3 * Go) .* ((sin(deg2rad(theta))^2 .* sin(deg2rad(phi))^2) - cos(deg2rad(theta))^2));
Ukk = (gs + (gdd * Fka)) * gamma4; Ukk = (gs + (gdd * Fka)) * gamma4;
end end
@ -143,7 +143,4 @@ function ret = computeTotalEnergyPerParticle(x, as, AtomNumberDensity, wz, lz, g
Energy_AxialComponent = (PlanckConstantReduced * wz) * ((lz^2/(4 * MeanWidth^2)) + (MeanWidth^2/(4 * lz^2))); 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)); Energy_TransverseComponent = (0.5 * (gs + (2*gdd)) * gamma4 * AtomNumberDensity) + ((2/5) * gQF * AtomNumberDensity^(3/2));
ret = (Energy_AxialComponent + Energy_TransverseComponent) / (PlanckConstantReduced * wz); ret = (Energy_AxialComponent + Energy_TransverseComponent) / (PlanckConstantReduced * wz);
end end

318
Estimations/DipolarDispersionAndRotonInstabilityBoundary/RIBForTiltedDipoles_TwoDirections.m Normal file
View File

@ -0,0 +1,318 @@
%% 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.005:0.0025:0.5;
as_to_add = 0.3:0.025:0.95;
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(6)
clf
set(gcf,'Position',[50 50 1850 750])
theta = 66; % 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, theta, 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
subplot(1, 2, 1); % 1 row, 2 columns, first subplot
imagesc(nadd2s, as_to_add, instability_boundary); % Specify x and y data for axes
set(gca, 'YDir', 'normal'); % Correct the y-axis direction
colorbar; % Add a colorbar
caxis([0 1])
xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex');
ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex');
title(['Along Y: $\theta = ',num2str(theta), '; \phi = ', num2str(phi),'$'],'fontsize',16,'interpreter','latex')
theta = 66; % Polar angle of dipole moment
phi = 90; % 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, theta, 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
% set(gcf,'Position',[50 50 950 750])
subplot(1, 2, 2); % 1 row, 2 columns, first subplot
imagesc(nadd2s, as_to_add, instability_boundary); % Specify x and y data for axes
set(gca, 'YDir', 'normal'); % Correct the y-axis direction
colorbar; % Add a colorbar
caxis([0 1])
xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex');
ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex');
title(['Along X: $\theta = ',num2str(theta), '; \phi = ', num2str(phi),'$'],'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';
caxis([0 1])
% 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');
%}
sgtitle('Mean-field instability boundary','fontsize',16,'interpreter','latex')
%% Cycle through angles
% Define values for theta and phi
theta_values = 0:2:90; % Range of theta values (you can modify this)
% Set up VideoWriter object to produce a movie
v = VideoWriter('rib_movie', 'MPEG-4'); % Create a video object
v.FrameRate = 5; % Frame rate of the video
open(v); % Open the video file
for theta = theta_values
figure(6)
clf
set(gcf,'Position',[50 50 1850 750])
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, theta, 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
subplot(1, 2, 1); % 1 row, 2 columns, first subplot
imagesc(nadd2s, as_to_add, instability_boundary); % Specify x and y data for axes
set(gca, 'YDir', 'normal'); % Correct the y-axis direction
colorbar; % Add a colorbar
caxis([0 1])
xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex');
ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex');
title(['Along Y: $\theta = ',num2str(theta), '; \phi = ', num2str(phi),'$'],'fontsize',16,'interpreter','latex')
phi = 90; % 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, theta, 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
% set(gcf,'Position',[50 50 950 750])
subplot(1, 2, 2); % 1 row, 2 columns, first subplot
imagesc(nadd2s, as_to_add, instability_boundary); % Specify x and y data for axes
set(gca, 'YDir', 'normal'); % Correct the y-axis direction
colorbar; % Add a colorbar
caxis([0 1])
xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex');
ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex');
title(['Along X: $\theta = ',num2str(theta), '; \phi = ', num2str(phi),'$'],'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';
caxis([0 1])
% 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');
%}
% Capture the frame and write to video
frame = getframe(gcf); % Capture the current figure
writeVideo(v, frame); % Write the frame to the video
% sgtitle('Mean-field instability boundary','fontsize',16,'interpreter','latex')
end
% Close the video file
close(v);
%%
function [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, theta, 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(theta))^2 - 1) + ((3 * Go) .* ((sin(deg2rad(theta))^2 .* sin(deg2rad(phi))^2) - cos(deg2rad(theta))^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

20
Estimations/DipolarDispersionAndRotonInstabilityBoundary/ReproduceBlairBlakieResults.m
View File

@ -34,7 +34,7 @@ wz = 2 * pi * 72.4;
lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length
as = 102.515 * BohrRadius; % Scattering 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 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 theta = 0; % Polar angle of dipole moment
phi = 0; % Azimuthal angle of momentum vector phi = 0; % Azimuthal angle of momentum vector
MeanWidth = 5.7304888515 * lz; % Mean width of Gaussian ansatz MeanWidth = 5.7304888515 * lz; % Mean width of Gaussian ansatz
k = linspace(0, 2e6, 1000); % Vector of magnitudes of k vector k = linspace(0, 2e6, 1000); % Vector of magnitudes of k vector
@ -47,7 +47,7 @@ eps_dd = add/as;
gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength
gdd = VacuumPermeability*DyMagneticMoment^2/3; gdd = VacuumPermeability*DyMagneticMoment^2/3;
[Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, alpha, phi); % DDI potential in k-space [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, theta, phi); % DDI potential in k-space
% == Quantum Fluctuations term == % % == Quantum Fluctuations term == %
gammaQF = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2)); gammaQF = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2));
@ -74,7 +74,7 @@ AtomNumber = 1E5;
wz = 2 * pi * 72.4; % Trap frequency in the tight confinement direction wz = 2 * pi * 72.4; % Trap frequency in the tight confinement direction
lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length 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 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 theta = 0; % Polar angle of dipole moment
phi = 0; % Azimuthal angle of momentum vector phi = 0; % Azimuthal angle of momentum vector
MeanWidth = 5.7304888515 * lz; % Mean width of Gaussian ansatz MeanWidth = 5.7304888515 * lz; % Mean width of Gaussian ansatz
k = linspace(0, 2e6, 1000); % Vector of magnitudes of k vector k = linspace(0, 2e6, 1000); % Vector of magnitudes of k vector
@ -92,7 +92,7 @@ for idx = 1:length(ScatteringLengths)
gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength
gdd = VacuumPermeability*DyMagneticMoment^2/3; gdd = VacuumPermeability*DyMagneticMoment^2/3;
[Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, alpha, phi); % DDI potential in k-space [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, theta, phi); % DDI potential in k-space
% == Quantum Fluctuations term == % % == Quantum Fluctuations term == %
gammaQF = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2)); gammaQF = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2));
@ -128,7 +128,7 @@ legend('location', 'northwest','fontsize',16, 'Interpreter','latex')
wz = 2 * pi * 72.4; % Trap frequency in the tight confinement direction wz = 2 * pi * 72.4; % Trap frequency in the tight confinement direction
lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length
alpha = 0; % Polar angle of dipole moment theta = 0; % Polar angle of dipole moment
phi = 0; % Azimuthal angle of momentum vector phi = 0; % Azimuthal angle of momentum vector
k = linspace(0, 2.25e6, 1000); % Vector of magnitudes of k vector k = linspace(0, 2.25e6, 1000); % Vector of magnitudes of k vector
@ -154,7 +154,7 @@ for idx = 1:length(nadd2s)
gdd = VacuumPermeability*DyMagneticMoment^2/3; gdd = VacuumPermeability*DyMagneticMoment^2/3;
MeanWidth = var_widths(idx) * lz; % Mean width of Gaussian ansatz 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 [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, theta, phi); % DDI potential in k-space
% == Quantum Fluctuations term == % % == Quantum Fluctuations term == %
gammaQF = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2)); gammaQF = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2));
@ -249,7 +249,7 @@ ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex');
% ====================================================================================================================================================== % % ====================================================================================================================================================== %
alpha = 0; % Polar angle of dipole moment theta = 0; % Polar angle of dipole moment
phi = 0; % Azimuthal angle of momentum vector phi = 0; % Azimuthal angle of momentum vector
k = linspace(0, 2.25e6, 1000); % Vector of magnitudes of k vector k = linspace(0, 2.25e6, 1000); % Vector of magnitudes of k vector
instability_boundary = zeros(length(as_to_add), length(nadd2s)); instability_boundary = zeros(length(as_to_add), length(nadd2s));
@ -270,7 +270,7 @@ for idx = 1:length(nadd2s)
gdd = VacuumPermeability*DyMagneticMoment^2/3; gdd = VacuumPermeability*DyMagneticMoment^2/3;
MeanWidth = var_widths(jdx, idx) * lz; % Mean width of Gaussian ansatz 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 [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, theta, phi); % DDI potential in k-space
% == Quantum Fluctuations term == % % == Quantum Fluctuations term == %
gammaQF = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2)); gammaQF = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2));
@ -312,10 +312,10 @@ title('Roton instability boundary','fontsize',16,'interpreter','latex')
%} %}
%% %%
function [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, alpha, phi) function [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, theta, phi)
Go = sqrt(pi) * (k * MeanWidth/sqrt(2)) .* exp((k * MeanWidth/sqrt(2)).^2) .* erfc((k * MeanWidth/sqrt(2))); Go = sqrt(pi) * (k * MeanWidth/sqrt(2)) .* exp((k * MeanWidth/sqrt(2)).^2) .* erfc((k * MeanWidth/sqrt(2)));
gamma4 = 1/(sqrt(2*pi) * MeanWidth); 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)); Fka = (3 * cos(deg2rad(theta))^2 - 1) + ((3 * Go) .* ((sin(deg2rad(theta))^2 .* sin(deg2rad(phi))^2) - cos(deg2rad(theta))^2));
Ukk = (gs + (gdd * Fka)) * gamma4; Ukk = (gs + (gdd * Fka)) * gamma4;
end end