From 2cd562ab920058e1b71179b5fddc56300cb023c0 Mon Sep 17 00:00:00 2001 From: Karthik Date: Sun, 29 Jun 2025 19:09:49 +0200 Subject: [PATCH] New script to generate an animation of the roton softening. --- .../RotonInstability/AnimateRotonSoftening.m | 215 ++++++++++++++++++ 1 file changed, 215 insertions(+) create mode 100644 Estimations/RotonInstability/AnimateRotonSoftening.m diff --git a/Estimations/RotonInstability/AnimateRotonSoftening.m b/Estimations/RotonInstability/AnimateRotonSoftening.m new file mode 100644 index 0000000..9a822c0 --- /dev/null +++ b/Estimations/RotonInstability/AnimateRotonSoftening.m @@ -0,0 +1,215 @@ +% === Preliminaries === +clear; clc; + +%% 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; + +% Dy specific constants +Dy164Mass = 163.929174751*AtomicMassUnit; +DyMagneticMoment = 9.93 * BohrMagneton; + +% === System Parameters === +AtomNumber = 1E5; +wz = 2 * pi * 72.4; +lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); +Trapsize = 7.5815 * lz; +theta = 0; +phi = 0; +MeanWidth = 5.7304888515 * lz; +k = linspace(0, 3e6, 1000); +AtomNumberDensity = AtomNumber / Trapsize^2; +add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); +gdd = VacuumPermeability*DyMagneticMoment^2/3; + +% === eps_dd sweep parameters === +eps_dd_values = linspace(1.15, 1.2755, 40); +max_kval = max(k .* add); +max_yval = 200; % Hz - manually set this based on known safe upper bound +frames(length(eps_dd_values)) = struct('cdata', [], 'colormap', []); +pauseFrameIdx = 0; + +% History for trail +xHist = {}; +yHist = {}; + +% === Plotting === +figure(1) +set(gcf,'Position',[100 100 950 750]) +set(gcf, 'Color', 'w'); % Set figure background to white +set(gca, 'Color', 'w'); % Set axes background to white + + +for idx = 1:length(eps_dd_values) + eps_dd = eps_dd_values(idx); + as = add / eps_dd; + gs = 4 * pi * PlanckConstantReduced^2 / Dy164Mass * as; + + % Compute DDI potential + [~, ~, ~, Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, theta, phi); + + % Quantum fluctuations + 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)); + + % Check for instability + if any(~isreal(DeltaK)) + disp(['Instability reached at eps_dd = ', num2str(eps_dd)]); + break + end + + EpsilonK = sqrt(((PlanckConstantReduced^2 .* k.^2) ./ (2 * Dy164Mass)) .* DeltaK); + + % Store current curve + xvals = k .* add; + yvals = EpsilonK ./ PlanckConstant; + xHist{end+1} = xvals; + yHist{end+1} = yvals; + + % Plot all past curves + clf; + hold on; + Nfade = length(xHist); + cmap = coolwarm(Nfade); % Use custom coolwarm colormap + for n = 1:5:Nfade-1 + plot(xHist{n}, yHist{n}, 'Color', [cmap(n,:), 0.2 + 0.6*(n/Nfade)], 'LineWidth', 1.5); + end + + % Highlight current curve + minDelta = min(DeltaK); + [~, minIdx] = min(DeltaK); + + if (minDelta < 0.0001 * max(DeltaK)) && (minIdx > 10) + plot(xvals, yvals, 'Color', [180, 4, 38]/255, 'LineWidth', 5); % Orange + + % Arrow + x_arrow = xvals(minIdx); + y_arrow = yvals(minIdx); + dx = 0; + dy = -0.1 * max_yval; + + quiver(x_arrow, y_arrow - dy/2, dx, dy, ... + 'AutoScale', 'off', 'MaxHeadSize', 1.5, ... + 'Color', 'k', 'LineWidth', 2); + + text(x_arrow, y_arrow - dy/2 + 0.05 * max_yval, ... + '$2\pi/\lambda_{\mathrm{rot}}$', ... + 'Interpreter', 'latex', 'FontSize', 16, 'Color', 'k', ... + 'HorizontalAlignment', 'center'); + + pauseFrameIdx = idx; + else + plot(xvals, yvals, 'Color', [59, 76, 192]/255, 'LineWidth', 5); % Blue + end + + % Axes + % title(['$a_s = ', num2str(round(1/eps_dd,3)), 'a_{dd}, \quad na_{dd}^2 = ', num2str(round(AtomNumberDensity * add^2,4)), '$'], ... + % 'fontsize', 16, 'interpreter', 'latex') + % xlabel('$k_{\rho}a_{dd}$','fontsize',26,'interpreter','latex') + % ylabel('$\epsilon(k_{\rho})/h$ (Hz)','fontsize',26,'interpreter','latex') + + xlim([0, max_kval]) + ylim([0, max_yval]) + set(gca, 'XTick', [], 'YTick', [], 'Color', 'w', 'XColor', 'none', 'YColor', 'none'); + % Arrow positions (in normalized figure coordinates) + ax = gca; + ax.Units = 'normalized'; + pos = ax.Position; + + % X-axis arrow (horizontal) + annotation('arrow', ... + [pos(1), pos(1) + pos(3)*0.95], ... + [pos(2), pos(2)], ... + 'Color', 'k', 'LineWidth', 2); + + % Y-axis arrow (vertical) + annotation('arrow', ... + [pos(1), pos(1)], ... + [pos(2), pos(2) + pos(4)*0.95], ... + 'Color', 'k', 'LineWidth', 2); + + % x-label + text(pos(1) + pos(3)*0.5, pos(2) - 0.15, 'Momentum', ... + 'Units', 'normalized', ... + 'HorizontalAlignment', 'center', ... + 'VerticalAlignment', 'top', ... + 'FontName', 'Bahnschrift' , ... + 'FontSize', 26); + + % y-label + text(pos(1) - 0.2, pos(2) + pos(4)*0.5, 'Energy', ... + 'Units', 'normalized', ... + 'HorizontalAlignment', 'center', ... + 'VerticalAlignment', 'middle', ... + 'FontName', 'Bahnschrift' , ... + 'FontSize', 26, ... + 'Rotation', 90); + + + drawnow; + + frames(idx) = getframe(gcf); +end + +%% === Save as GIF === +filename = 'roton_dispersion.gif'; +delayTime = 0.15; +pauseRepeats = 10; + +for i = 1:length(frames) + [A, map] = rgb2ind(frame2im(frames(i)), 256); + if i == 1 + imwrite(A, map, filename, 'gif', 'LoopCount', Inf, 'DelayTime', delayTime); + elseif i == pauseFrameIdx + for j = 1:pauseRepeats + imwrite(A, map, filename, 'gif', 'WriteMode', 'append', 'DelayTime', delayTime); + end + else + imwrite(A, map, filename, 'gif', 'WriteMode', 'append', 'DelayTime', delayTime); + end +end + +disp(['GIF saved as ', filename]); + +%% === Helper Function === +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 cmap = coolwarm(n) + % Generates a coolwarm colormap with n colors + c1 = [59, 76, 192]/255; + c2 = [180, 4, 38]/255; + r = linspace(c1(1), c2(1), n)'; + g = linspace(c1(2), c2(2), n)'; + b = linspace(c1(3), c2(3), n)'; + cmap = [r g b]; +end