157 lines
7.9 KiB
Mathematica
157 lines
7.9 KiB
Mathematica
|
%% 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
|