diff --git a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersion2D.m b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersion2D.m index cd7690c..b487963 100644 --- a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersion2D.m +++ b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersion2D.m @@ -77,7 +77,7 @@ MeanWidth = sigma * lz; VDk = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth); VDk_fftshifted = fftshift(VDk); -figure(11) +figure(8) clf set(gcf,'Position',[50 50 950 750]) imagesc(fftshift(Transf.kx)*1e-6, fftshift(Transf.ky)*1e-6, VDk_fftshifted); % Specify x and y data for axes @@ -106,7 +106,7 @@ end EpsilonK = double(imag(EpsilonK) ~= 0); % 'isreal' returns 0 for complex numbers and 1 for real numbers, so we negate it -figure(12) +figure(9) clf set(gcf,'Position',[50 50 950 750]) imagesc(fftshift(Transf.kx)*1e-6, fftshift(Transf.ky)*1e-6, EpsilonK); % Specify x and y data for axes @@ -159,7 +159,7 @@ for theta = theta_values EpsilonK = double(imag(EpsilonK) ~= 0); % 'isreal' returns 0 for complex numbers and 1 for real numbers, so we negate it % Plot the result - figure(13) + figure(10) clf set(gcf,'Position',[50 50 950 750]) imagesc(fftshift(Transf.kx)*1e-6, fftshift(Transf.ky)*1e-6, EpsilonK); % Specify x and y data for axes diff --git a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ExtractingKRoton.m b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ExtractingKRoton.m index 1f26390..ce83121 100644 --- a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ExtractingKRoton.m +++ b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ExtractingKRoton.m @@ -28,123 +28,114 @@ Dy164Mass = 163.929174751*AtomicMassUnit; Dy164IsotopicAbundance = 0.2826; DyMagneticMoment = 9.93*BohrMagneton; -%% Roton instability boundary for tilted dipoles +%% Extracting values from the 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 +%-------TEST-------% +% nadd2s = 0.05:0.005:0.25; +% as_to_add = 0.76:0.0001:0.81; -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; +%-------DEPLOY-------% +nadd2s = 0.005:0.005:0.5; +as_to_add = 0.35:0.0001:1.15; -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)); +data_struct = struct; +wz_values = [150, 300, 500, 750]; +theta_values = 0:5:45; % Range of theta values (you can modify this) +phi = 0; % Azimuthal angle of momentum vector +kvec = linspace(0, 2.25e6, 1000); % Vector of magnitudes of k vector -x0 = 5; -Aineq = []; -Bineq = []; -Aeq = []; -Beq = []; -lb = [1]; -ub = [10]; -nonlcon = []; -fminconopts = optimoptions(@fmincon,'Display','off', 'StepTolerance', 1.0000e-11, 'MaxIterations',1500); +for mainloop_idx = 1:length(wz_values) + format long -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 + 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; -%% ====================================================================================================================================================== % - -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)); -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 + 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)); - [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, 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 .* 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 - -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 + 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 -% 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 - -% 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') +save('ExtractingKRoton_Result.mat', 'data_struct'); %% function [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, theta, phi) @@ -164,4 +155,161 @@ function ret = computeTotalEnergyPerParticle(x, as, AtomNumberDensity, wz, lz, g 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 \ No newline at end of file +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; + 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 diff --git a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/RIBForTiltedDipoles_TwoDirections.m b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/RIBForTiltedDipoles_TwoDirections.m index 7cc32ff..03dfe28 100644 --- a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/RIBForTiltedDipoles_TwoDirections.m +++ b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/RIBForTiltedDipoles_TwoDirections.m @@ -62,7 +62,7 @@ end %% ====================================================================================================================================================== % -figure(6) +figure(7) clf set(gcf,'Position',[50 50 1850 750]) @@ -181,7 +181,7 @@ v.FrameRate = 5; % Frame rate of the video open(v); % Open the video file for theta = theta_values - figure(6) + figure(7) clf set(gcf,'Position',[50 50 1850 750]) phi = 0; % Azimuthal angle of momentum vector diff --git a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ScalingOfTheQFTerm.m b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ScalingOfTheQFTerm.m index 36bca94..4c5ee06 100644 --- a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ScalingOfTheQFTerm.m +++ b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ScalingOfTheQFTerm.m @@ -76,7 +76,7 @@ for idx = 1:length(nadd2s) end end -figure(7) +figure clf set(gcf,'Position',[50 50 950 750]) imagesc(AtomNumber*1E-5, ScatteringLengths, QF * 1E31); % Specify x and y data for axes