%% 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.35:0.001:1.15; data_struct = struct; wz_values = [150, 300, 500, 750]; theta_values = 0:5:45; % Range of theta values phi = 0; % Azimuthal angle of momentum vector kvec = linspace(0, 5e6, 1000); % Vector of magnitudes of k 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('ExtractingKRoton_Result.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