Latest.
This commit is contained in:
parent
0833bddee8
commit
f3ac9e3b12
@ -168,7 +168,7 @@ for theta = theta_values
|
||||
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.theta), '; \eta = ', num2str(Params.eta),'$'],'fontsize',16,'interpreter','latex')
|
||||
title(['2-D Dispersion: $\theta = ',num2str(Params.theta), '; \eta = ', num2str(Params.eta),'$'],'fontsize',16,'interpreter','latex')
|
||||
|
||||
% Capture the frame and write to video
|
||||
% frame = getframe(gcf); % Capture the current figure
|
||||
|
||||
@ -28,19 +28,17 @@ Dy164Mass = 163.929174751*AtomicMassUnit;
|
||||
Dy164IsotopicAbundance = 0.2826;
|
||||
DyMagneticMoment = 9.93*BohrMagneton;
|
||||
|
||||
%% k_roton at the instability boundary for tilted dipoles
|
||||
%% Roton instability boundary for tilted dipoles
|
||||
|
||||
wz = 2 * pi * 72.4; % Trap frequency in the tight confinement direction
|
||||
theta = 0; % Polar angle of dipole moment
|
||||
|
||||
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;
|
||||
|
||||
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;
|
||||
@ -64,13 +62,11 @@ for idx = 1:length(nadd2s)
|
||||
end
|
||||
end
|
||||
|
||||
% ====================================================================================================================================================== %
|
||||
%% ====================================================================================================================================================== %
|
||||
|
||||
theta = 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
|
||||
@ -99,119 +95,58 @@ for idx = 1:length(nadd2s)
|
||||
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
|
||||
|
||||
figure(6)
|
||||
clf
|
||||
set(gcf,'Position',[50 50 950 750])
|
||||
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
|
||||
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')
|
||||
|
||||
|
||||
%%
|
||||
|
||||
matrix = instability_boundary;
|
||||
|
||||
% 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];
|
||||
col_indices = [col_indices; j];
|
||||
break; % Stop after the first transition in the column
|
||||
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
|
||||
% 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
|
||||
|
||||
% 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
|
||||
% Plot the extracted values
|
||||
figure(7);
|
||||
plot(xvals, yvals, '-o');
|
||||
title('Instability Boundary');
|
||||
xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex');
|
||||
ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex')
|
||||
|
||||
%%
|
||||
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);
|
||||
|
||||
Loading…
Reference in New Issue
Block a user