diff --git a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersion2D.m b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersion2D.m index a2d26e1..cd7690c 100644 --- a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersion2D.m +++ b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersion2D.m @@ -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 diff --git a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ExtractingKRoton.m b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ExtractingKRoton.m index 437abf1..1f26390 100644 --- a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ExtractingKRoton.m +++ b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ExtractingKRoton.m @@ -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);