Added code to determine roton instability boundary.
This commit is contained in:
parent
a8a0901eaf
commit
4ad098ae9a
@ -134,7 +134,7 @@ k = linspace(0, 2.25e6, 1000);
|
|||||||
|
|
||||||
nadd2s = [0.0844, 0.0978, 0.123];
|
nadd2s = [0.0844, 0.0978, 0.123];
|
||||||
as_to_add = [0.7730, 0.7840, 0.7819];
|
as_to_add = [0.7730, 0.7840, 0.7819];
|
||||||
var_widths = [4.97165, 5.72960, 5.93178];
|
var_widths = [4.97165, 5.7296048721, 5.93178];
|
||||||
|
|
||||||
add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
|
add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
|
||||||
EpsilonKs = zeros(length(k), length(nadd2s));
|
EpsilonKs = zeros(length(k), length(nadd2s));
|
||||||
@ -179,6 +179,112 @@ 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
|
%% 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
|
||||||
|
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;
|
||||||
|
AtomNumberDensity = 0.0978 / add^2;
|
||||||
|
as = 0.784 * add; % Scattering length
|
||||||
|
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
|
||||||
|
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.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
|
||||||
|
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));
|
||||||
|
|
||||||
|
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(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');
|
||||||
|
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, alpha, phi)
|
||||||
@ -187,3 +293,15 @@ function [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanW
|
|||||||
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(alpha))^2 - 1) + ((3 * Go) .* ((sin(deg2rad(alpha))^2 .* sin(deg2rad(phi))^2) - cos(deg2rad(alpha))^2));
|
||||||
Ukk = (gs + (gdd * Fka)) * gamma4;
|
Ukk = (gs + (gdd * Fka)) * gamma4;
|
||||||
end
|
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
|
||||||
1
roton_instability_project.json
Normal file
1
roton_instability_project.json
Normal file
@ -0,0 +1 @@
|
|||||||
|
{"version":[4,2],"axesColl":[{"name":"XY","type":"XYAxes","isLogX":false,"isLogY":false,"noRotation":true,"calibrationPoints":[{"px":304.86283783783784,"py":633.6364864864864,"dx":"0.1","dy":"0.745","dz":null},{"px":694.7418918918919,"py":633.6364864864864,"dx":"0.2","dy":"0.745","dz":null},{"px":107.59864864864865,"py":573.8594594594595,"dx":"0.1","dy":"0.745","dz":null},{"px":107.59864864864865,"py":27.895945945945947,"dx":"0.2","dy":"0.79","dz":null}]}],"datasetColl":[{"name":"Default Dataset","axesName":"XY","colorRGB":[200,0,0,255],"metadataKeys":[],"data":[{"x":108.92702702702702,"y":345.3783783783784,"value":[0.04974446337308348,0.7638321167883212]},{"x":124.86756756756756,"y":298.88513513513516,"value":[0.053833049403747876,0.7676642335766424]},{"x":135.4945945945946,"y":273.6459459459459,"value":[0.05655877342419081,0.7697445255474452]},{"x":172.6891891891892,"y":204.57027027027027,"value":[0.06609880749574107,0.7754379562043796]},{"x":564.5608108108108,"y":90.32972972972973,"value":[0.1666098807495741,0.7848540145985402]},{"x":630.9797297297297,"y":102.28513513513514,"value":[0.1836456558773424,0.7838686131386862]},{"x":677.472972972973,"y":111.58378378378379,"value":[0.19557069846678027,0.7831021897810219]},{"x":717.3243243243244,"y":120.88243243243244,"value":[0.20579216354344126,0.7823357664233577]},{"x":766.4743243243244,"y":132.83783783783784,"value":[0.2183986371379898,0.7813503649635036]},{"x":844.8486486486487,"y":152.76351351351352,"value":[0.2385008517887564,0.7797080291970803]},{"x":895.327027027027,"y":166.0472972972973,"value":[0.25144804088586037,0.7786131386861314]},{"x":184.6445945945946,"y":188.62972972972972,"value":[0.06916524701873936,0.7767518248175183]},{"x":199.25675675675674,"y":170.03243243243244,"value":[0.07291311754684839,0.7782846715328468]},{"x":220.51081081081082,"y":148.7783783783784,"value":[0.07836456558773425,0.780036496350365]},{"x":247.07837837837837,"y":127.52432432432433,"value":[0.08517887563884158,0.7817883211678832]},{"x":268.3324324324324,"y":114.24054054054054,"value":[0.09063032367972744,0.7828832116788321]},{"x":322.79594594594596,"y":91.65810810810811,"value":[0.10459965928449745,0.7847445255474452]},{"x":357.33378378378376,"y":83.68783783783783,"value":[0.11345826235093698,0.7854014598540147]},{"x":375.93108108108106,"y":81.03108108108108,"value":[0.1182282793867121,0.7856204379562044]},{"x":391.8716216216216,"y":79.70270270270271,"value":[0.1223168654173765,0.7857299270072993]},{"x":409.14054054054054,"y":78.37432432432432,"value":[0.12674616695059626,0.7858394160583941]},{"x":426.40945945945947,"y":78.37432432432432,"value":[0.13117546848381603,0.7858394160583941]},{"x":443.6783783783784,"y":78.37432432432432,"value":[0.1356047700170358,0.7858394160583941]},{"x":460.9472972972973,"y":78.37432432432433,"value":[0.14003407155025555,0.7858394160583941]},{"x":482.2013513513514,"y":80.3668918918919,"value":[0.14548551959114142,0.7856751824817518]},{"x":504.7837837837838,"y":82.35945945945946,"value":[0.15127768313458265,0.7855109489051095]},{"x":530.022972972973,"y":85.01621621621622,"value":[0.1577512776831346,0.7852919708029197]},{"x":599.0986486486487,"y":95.64324324324325,"value":[0.17546848381601365,0.7844160583941606]}],"autoDetectionData":null}],"measurementColl":[]}
|
||||||
BIN
roton_instability_project.tar
Normal file
BIN
roton_instability_project.tar
Normal file
Binary file not shown.
Loading…
Reference in New Issue
Block a user