503 lines
25 KiB
Matlab
503 lines
25 KiB
Matlab
%% 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;
|
|
|
|
%% Extracting values from the roton instability boundary for tilted dipoles
|
|
|
|
%-------TEST-------%
|
|
% nadd2s = 0.05:0.005:0.25;
|
|
% as_to_add = 0.76:0.001:0.81;
|
|
|
|
%-------DEPLOY-------%
|
|
nadd2s = 0.005:0.005:0.5;
|
|
as_to_add = 0.250:0.001:1.15;
|
|
|
|
data_struct = struct;
|
|
|
|
% wz_values = [150, 300, 500, 750];
|
|
% kvec = linspace(0, 5e6, 1000); % Vector of magnitudes of k vector
|
|
|
|
wz_values = [1000, 3000, 5000, 7000];
|
|
kvec = linspace(0, 15e6, 1000); % Vector of magnitudes of k vector
|
|
|
|
% wz_values = [10000, 13000, 15000];
|
|
% kvec = linspace(0, 25e6, 1000); % Vector of magnitudes of k vector
|
|
|
|
theta_values = 0:5:45; % Range of theta values
|
|
phi = 0; % Azimuthal angle of momentum vector
|
|
|
|
for mainloop_idx = 1:length(wz_values)
|
|
format long
|
|
|
|
PlanckConstantReduced = 6.62607015E-34/(2*pi);
|
|
AtomicMassUnit = 1.660539066E-27;
|
|
Dy164Mass = 163.929174751*AtomicMassUnit;
|
|
VacuumPermeability = 1.25663706212E-6;
|
|
BohrMagneton = 9.274009994E-24;
|
|
DyMagneticMoment = 9.93*BohrMagneton;
|
|
|
|
wz = 2 * pi * wz_values(mainloop_idx); % 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;
|
|
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
|
|
|
|
eps_dd_values = zeros(length(theta_values), 1);
|
|
n_values = zeros(length(theta_values), 1);
|
|
k_roton_values = zeros(length(theta_values), 1);
|
|
|
|
for idx = 1:length(theta_values)
|
|
theta = theta_values(idx);
|
|
[eps_dd_values(idx), n_values(idx), k_roton_values(idx)] = extractFromBoundaryCurve(theta, phi, nadd2s, as_to_add, var_widths, wz, lz, kvec);
|
|
end
|
|
|
|
data_struct(mainloop_idx).wz_value = wz / (2 * pi);
|
|
data_struct(mainloop_idx).theta_values = theta_values;
|
|
data_struct(mainloop_idx).eps_dd_values = eps_dd_values;
|
|
data_struct(mainloop_idx).n_values = n_values;
|
|
data_struct(mainloop_idx).k_roton_values = k_roton_values;
|
|
|
|
%{
|
|
figure(13)
|
|
clf
|
|
set(gcf,'Position',[50 50 950 750])
|
|
plot(theta_values, eps_dd_values, '-o', LineWidth=2.0)
|
|
xlabel('$\theta$','fontsize',16,'interpreter','latex');
|
|
ylabel('$\epsilon_{dd}$','fontsize',16,'interpreter','latex');
|
|
% title([''],'fontsize',16,'interpreter','latex')
|
|
grid on
|
|
|
|
figure(14)
|
|
clf
|
|
set(gcf,'Position',[50 50 950 750])
|
|
plot(theta_values, (1./eps_dd_values) * (add/BohrRadius), '-o', LineWidth=2.0)
|
|
xlabel('$\theta$','fontsize',16,'interpreter','latex');
|
|
ylabel('$a_s (\times a_o)$','fontsize',16,'interpreter','latex');
|
|
% title([''],'fontsize',16,'interpreter','latex')
|
|
grid on
|
|
|
|
figure(15)
|
|
clf
|
|
set(gcf,'Position',[50 50 950 750])
|
|
plot(theta_values, n_values * 1E-15, '-o', LineWidth=2.0)
|
|
xlabel('$\theta$','fontsize',16,'interpreter','latex');
|
|
ylabel('$n (\times 10^{3} \mu m^{-2})$','fontsize',16,'interpreter','latex');
|
|
% title([''],'fontsize',16,'interpreter','latex')
|
|
grid on
|
|
|
|
figure(16)
|
|
clf
|
|
set(gcf,'Position',[50 50 950 750])
|
|
plot(theta_values, k_roton_values * 1E-6, '-o', LineWidth=2.0)
|
|
xlabel('$\theta$','fontsize',16,'interpreter','latex');
|
|
ylabel('$k_{roton} (\mu m^{-1})$','fontsize',16,'interpreter','latex');
|
|
% title([''],'fontsize',16,'interpreter','latex')
|
|
grid on
|
|
%}
|
|
end
|
|
|
|
save('.\Results\ExtractingKRoton_Result.mat', 'data_struct');
|
|
|
|
%% Extracting values from the roton instability boundary for tilted dipoles - fixed atom number, trap frequency
|
|
|
|
%-------DEPLOY-------%
|
|
|
|
N = 1E5;
|
|
add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
|
|
area = 100; % in square micrometers
|
|
ppmum = N / area;
|
|
nadd2s = ppmum*1E12*add^2;
|
|
as_to_add = 0.150:0.001:1.15;
|
|
|
|
data_struct = struct;
|
|
|
|
wz_values = [500, 750, 1000, 2000];
|
|
kvec = linspace(0, 15e6, 1000); % Vector of magnitudes of k vector
|
|
|
|
theta_values = 0:5:90; % Range of theta values
|
|
phi = 90; % Azimuthal angle of momentum vector
|
|
|
|
for mainloop_idx = 1:length(wz_values)
|
|
format long
|
|
|
|
PlanckConstantReduced = 6.62607015E-34/(2*pi);
|
|
AtomicMassUnit = 1.660539066E-27;
|
|
Dy164Mass = 163.929174751*AtomicMassUnit;
|
|
VacuumPermeability = 1.25663706212E-6;
|
|
BohrMagneton = 9.274009994E-24;
|
|
DyMagneticMoment = 9.93*BohrMagneton;
|
|
|
|
wz = 2 * pi * wz_values(mainloop_idx); % 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;
|
|
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
|
|
|
|
eps_dd_values = zeros(length(theta_values), 1);
|
|
k_roton_values = zeros(length(theta_values), 1);
|
|
|
|
for idx = 1:length(theta_values)
|
|
theta = theta_values(idx);
|
|
[eps_dd_values(idx), k_roton_values(idx)] = extractFromBoundaryPoint(theta, phi, nadd2s, as_to_add, var_widths, wz, lz, kvec);
|
|
end
|
|
|
|
data_struct(mainloop_idx).wz_value = wz / (2 * pi);
|
|
data_struct(mainloop_idx).theta_values = theta_values;
|
|
data_struct(mainloop_idx).eps_dd_values = eps_dd_values;
|
|
data_struct(mainloop_idx).k_roton_values = k_roton_values;
|
|
|
|
end
|
|
|
|
save('.\Results\ExtractingKRoton_Result_FixedDensity_phi90.mat', 'data_struct');
|
|
%%
|
|
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
|
|
|
|
function [eps_dd, AtomNumberDensity, k_roton] = extractFromBoundaryCurve(theta, phi, nadd2s, as_to_add, var_widths, wz, lz, kvec)
|
|
|
|
format long
|
|
|
|
PlanckConstantReduced = 6.62607015E-34/(2*pi);
|
|
AtomicMassUnit = 1.660539066E-27;
|
|
Dy164Mass = 163.929174751*AtomicMassUnit;
|
|
VacuumPermeability = 1.25663706212E-6;
|
|
BohrMagneton = 9.274009994E-24;
|
|
DyMagneticMoment = 9.93*BohrMagneton;
|
|
add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
|
|
gdd = VacuumPermeability*DyMagneticMoment^2/3;
|
|
phase_diagram = zeros(length(as_to_add), length(nadd2s));
|
|
w0 = 2 * pi * 61.6316; % Trap frequency in the tight confinement direction
|
|
l0 = sqrt(PlanckConstantReduced/(Dy164Mass * w0)); % Defining a harmonic oscillator length
|
|
|
|
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
|
|
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(kvec, 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 .* kvec.^2) ./ (2 * Dy164Mass)) + ((2 * AtomNumberDensity) .* Ukk) + (3 * gQF * AtomNumberDensity^(3/2));
|
|
EpsilonK = sqrt(((PlanckConstantReduced^2 .* kvec.^2) ./ (2 * Dy164Mass)) .* DeltaK);
|
|
phase_diagram(jdx, idx) = ~isreal(EpsilonK);
|
|
end
|
|
end
|
|
%{
|
|
figure(11)
|
|
clf
|
|
set(gcf,'Position',[50 50 950 750])
|
|
imagesc(nadd2s, as_to_add, phase_diagram); % 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');
|
|
title(['$\theta = ',num2str(theta), '; \phi = 0','$', '(Along Y)'],'fontsize',16,'interpreter','latex')
|
|
%}
|
|
%-------------%
|
|
|
|
matrix = phase_diagram;
|
|
|
|
% Initialize arrays to store row and column indices of transitions
|
|
row_indices = [];
|
|
col_indices = [];
|
|
|
|
% Loop through the matrix to find transitions from 0 to 1
|
|
[rows, cols] = size(matrix);
|
|
for j = 1:cols
|
|
for i = 2:rows
|
|
if matrix(i-1, j) == 1 && matrix(i, j) == 0
|
|
row_indices = [row_indices; i-1];
|
|
col_indices = [col_indices; j];
|
|
break; % Stop after the first transition in the column
|
|
end
|
|
end
|
|
end
|
|
|
|
% Now extract the values from the corresponding vectors
|
|
xvals = zeros(length(col_indices), 1);
|
|
yvals = zeros(length(row_indices), 1);
|
|
for k = 1:length(row_indices)
|
|
row = row_indices(k);
|
|
col = col_indices(k);
|
|
xvals(k) = nadd2s(col);
|
|
yvals(k) = as_to_add(row);
|
|
end
|
|
|
|
instability_boundary = [xvals, yvals];
|
|
|
|
%-------------%
|
|
|
|
% Degree of the polynomial to fit
|
|
n = 5; % For a quadratic fit
|
|
|
|
% Fit the polynomial
|
|
p = polyfit(xvals, yvals, n);
|
|
|
|
% Display the polynomial coefficients
|
|
% disp('Polynomial coefficients:');
|
|
% disp(p);
|
|
|
|
% Evaluate the polynomial at points in x
|
|
y_fit = polyval(p, xvals);
|
|
|
|
%{
|
|
% Plot the original data and the fitted polynomial curve
|
|
figure(12);
|
|
clf
|
|
set(gcf,'Position',[50 50 950 750])
|
|
plot(xvals, yvals, 'o', 'LineWidth', 2.0, 'DisplayName', 'Extracted boundary points'); % Original data
|
|
hold on;
|
|
plot(xvals, y_fit, '-r','LineWidth', 2.0, 'DisplayName', ['Polynomial Fit (degree ' num2str(n) ')']); % Fitted curve
|
|
ylim([min(as_to_add) max(as_to_add)])
|
|
xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex');
|
|
ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex')
|
|
title(['$\theta = ',num2str(theta), '; \phi = 0','$', '(Along Y)'],'fontsize',16,'interpreter','latex')
|
|
legend('show');
|
|
grid on;
|
|
%}
|
|
[val, idx] = max(y_fit);
|
|
|
|
% Round down to 4 decimal places
|
|
rounded_val = floor(val * 10^4) / 10^4;
|
|
% Find nearest from original vector of boundary points
|
|
[~, nearest_idx] = min(abs(instability_boundary(:, 2) - rounded_val));
|
|
nearest_val = instability_boundary(nearest_idx, 2);
|
|
% Choose the scalar value between the two
|
|
if ~isscalar(nearest_val)
|
|
val = rounded_val;
|
|
else
|
|
val = nearest_val;
|
|
idx = nearest_idx;
|
|
end
|
|
|
|
AtomNumberDensity = xvals(idx) / add^2; % Areal density of atoms
|
|
as = val * add; % Scattering length
|
|
eps_dd = 1/val; % Relative interaction strength
|
|
gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength
|
|
x0 = 5;
|
|
Aineq = [];
|
|
Bineq = [];
|
|
Aeq = [];
|
|
Beq = [];
|
|
lb = [1];
|
|
ub = [10];
|
|
nonlcon = [];
|
|
fminconopts = optimoptions(@fmincon,'Display','off', 'StepTolerance', 1.0000e-11, 'MaxIterations',1500);
|
|
TotalEnergyPerParticle = @(x) computeTotalEnergyPerParticle(x, as, AtomNumberDensity, wz, lz, gs, add, gdd, PlanckConstantReduced);
|
|
sigma = fmincon(TotalEnergyPerParticle, x0, Aineq, Bineq, Aeq, Beq, lb, ub, nonlcon, fminconopts);
|
|
MeanWidth = sigma * lz; % Mean width of Gaussian ansatz
|
|
[Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(kvec, 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;
|
|
DeltaK = ((PlanckConstantReduced^2 .* kvec.^2) ./ (2 * Dy164Mass)) + ((2 * AtomNumberDensity) .* Ukk) + (3 * gQF * AtomNumberDensity^(3/2));
|
|
EpsilonK = sqrt(((PlanckConstantReduced^2 .* kvec.^2) ./ (2 * Dy164Mass)) .* DeltaK);
|
|
k_roton_indices = find(imag(EpsilonK) ~= 0);
|
|
if ~isempty(k_roton_indices)
|
|
k_roton = median(kvec(k_roton_indices));
|
|
else
|
|
k_roton = NaN;
|
|
end
|
|
end
|
|
|
|
function [eps_dd, k_roton] = extractFromBoundaryPoint(theta, phi, nadd2s, as_to_add, var_widths, wz, lz, kvec)
|
|
|
|
format long
|
|
|
|
PlanckConstantReduced = 6.62607015E-34/(2*pi);
|
|
AtomicMassUnit = 1.660539066E-27;
|
|
Dy164Mass = 163.929174751*AtomicMassUnit;
|
|
VacuumPermeability = 1.25663706212E-6;
|
|
BohrMagneton = 9.274009994E-24;
|
|
DyMagneticMoment = 9.93*BohrMagneton;
|
|
add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
|
|
gdd = VacuumPermeability*DyMagneticMoment^2/3;
|
|
phase_diagram = zeros(length(as_to_add), length(nadd2s));
|
|
w0 = 2 * pi * 61.6316; % Trap frequency in the tight confinement direction
|
|
l0 = sqrt(PlanckConstantReduced/(Dy164Mass * w0)); % Defining a harmonic oscillator length
|
|
|
|
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
|
|
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(kvec, 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 .* kvec.^2) ./ (2 * Dy164Mass)) + ((2 * AtomNumberDensity) .* Ukk) + (3 * gQF * AtomNumberDensity^(3/2));
|
|
EpsilonK = sqrt(((PlanckConstantReduced^2 .* kvec.^2) ./ (2 * Dy164Mass)) .* DeltaK);
|
|
phase_diagram(jdx, idx) = ~isreal(EpsilonK);
|
|
end
|
|
end
|
|
|
|
matrix = phase_diagram;
|
|
|
|
% Initialize arrays to store row and column indices of transitions
|
|
row_indices = [];
|
|
col_indices = [];
|
|
|
|
% Loop through the matrix to find transitions from 0 to 1
|
|
[rows, cols] = size(matrix);
|
|
for j = 1:cols
|
|
for i = 2:rows
|
|
if matrix(i-1, j) == 1 && matrix(i, j) == 0
|
|
row_indices = [row_indices; i-1];
|
|
col_indices = [col_indices; j];
|
|
break; % Stop after the first transition in the column
|
|
end
|
|
end
|
|
end
|
|
|
|
% Now extract the values from the corresponding vectors
|
|
xvals = zeros(length(col_indices), 1);
|
|
yvals = zeros(length(row_indices), 1);
|
|
for k = 1:length(row_indices)
|
|
row = row_indices(k);
|
|
col = col_indices(k);
|
|
xvals(k) = nadd2s(col);
|
|
yvals(k) = as_to_add(row);
|
|
end
|
|
|
|
instability_boundary = [xvals, yvals];
|
|
|
|
if ~isempty(instability_boundary)
|
|
val = instability_boundary(2);
|
|
AtomNumberDensity = instability_boundary(1) / add^2; % Areal density of atoms
|
|
as = val * add; % Scattering length
|
|
eps_dd = 1/val; % Relative interaction strength
|
|
gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength
|
|
x0 = 5;
|
|
Aineq = [];
|
|
Bineq = [];
|
|
Aeq = [];
|
|
Beq = [];
|
|
lb = [1];
|
|
ub = [10];
|
|
nonlcon = [];
|
|
fminconopts = optimoptions(@fmincon,'Display','off', 'StepTolerance', 1.0000e-11, 'MaxIterations',1500);
|
|
TotalEnergyPerParticle = @(x) computeTotalEnergyPerParticle(x, as, AtomNumberDensity, wz, lz, gs, add, gdd, PlanckConstantReduced);
|
|
sigma = fmincon(TotalEnergyPerParticle, x0, Aineq, Bineq, Aeq, Beq, lb, ub, nonlcon, fminconopts);
|
|
MeanWidth = sigma * lz; % Mean width of Gaussian ansatz
|
|
[Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(kvec, 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;
|
|
DeltaK = ((PlanckConstantReduced^2 .* kvec.^2) ./ (2 * Dy164Mass)) + ((2 * AtomNumberDensity) .* Ukk) + (3 * gQF * AtomNumberDensity^(3/2));
|
|
EpsilonK = sqrt(((PlanckConstantReduced^2 .* kvec.^2) ./ (2 * Dy164Mass)) .* DeltaK);
|
|
k_roton_indices = find(imag(EpsilonK) ~= 0);
|
|
if ~isempty(k_roton_indices)
|
|
k_roton = median(kvec(k_roton_indices));
|
|
else
|
|
k_roton = NaN;
|
|
end
|
|
else
|
|
eps_dd = NaN;
|
|
k_roton = NaN;
|
|
end
|
|
end |