From 0833bddee8b48cc8df2385f20ab26f6191383962 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Tue, 14 Jan 2025 17:30:45 +0100 Subject: [PATCH] Clean up, minor corrections to keep notation consistent. --- .../DipolarDispersion2D.m | 83 ++--- .../ExtractingKRoton.m | 81 ++++- .../RIBForTiltedDipoles.m | 13 +- .../RIBForTiltedDipoles_TwoDirections.m | 318 ++++++++++++++++++ .../ReproduceBlairBlakieResults.m | 20 +- 5 files changed, 453 insertions(+), 62 deletions(-) create mode 100644 Estimations/DipolarDispersionAndRotonInstabilityBoundary/RIBForTiltedDipoles_TwoDirections.m diff --git a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersion2D.m b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersion2D.m index 95a0a79..a2d26e1 100644 --- a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersion2D.m +++ b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersion2D.m @@ -29,9 +29,12 @@ Dy164IsotopicAbundance = 0.2826; DyMagneticMoment = 9.93*BohrMagneton; %% 2-D DDI Potential in k-space, with Gaussian ansatz width determined by constrained minimization +% With user-defined values of interaction, density and tilt -w0 = 2*pi*61.6316; % Angular frequency unit [s^-1] -l0 = sqrt(PlanckConstantReduced/(Dy164Mass*w0)); % Defining a harmonic oscillator length +% w0 = 2*pi*61.6316; % Angular frequency unit [s^-1] +% l0 = sqrt(PlanckConstantReduced/(Dy164Mass*w0)); +% % Defining a harmonic oscillator length - heredue to the choice of w0, l0 +% is 1 micrometer wz = 2 * pi * 300; % Trap frequency in the tight confinement direction lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length @@ -39,14 +42,14 @@ lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Number of grid points in each direction Params.Nx = 128; Params.Ny = 128; -Params.Lx = 150*l0; -Params.Ly = 150*l0; +Params.Lx = 150*1e-6; +Params.Ly = 150*1e-6; [Transf] = setupSpace(Params); -nadd2s = 0.110; -as_to_add = 0.782; -Params.alpha = 10; -Params.phi = 0; +nadd2s = 0.110; % Number density * add^2 +as_to_add = 0.782; % 1/edd +Params.theta = 57; % Polar angle of dipole moment +Params.eta = 0; % Azimuthal angle of dipole moment add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length gdd = VacuumPermeability*DyMagneticMoment^2/3; @@ -61,15 +64,13 @@ ub = [10]; nonlcon = []; fminconopts = optimoptions(@fmincon,'Display','off', 'StepTolerance', 1.0000e-11, 'MaxIterations',1500); -AtomNumberDensity = nadd2s / add^2; % Areal density of atoms +AtomNumberDensity = nadd2s / add^2; % Number density of atoms as = as_to_add * add; % Scattering length eps_dd = add/as; % Relative interaction strength 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); -%% Chosen values of interaction, density and tilt - MeanWidth = sigma * lz; % == 2-D DDI Potential in k-space == % @@ -79,13 +80,13 @@ VDk_fftshifted = fftshift(VDk); figure(11) clf set(gcf,'Position',[50 50 950 750]) -imagesc(fftshift(Transf.kx)*l0, fftshift(Transf.ky)*l0, VDk_fftshifted); % Specify x and y data for axes +imagesc(fftshift(Transf.kx)*1e-6, fftshift(Transf.ky)*1e-6, VDk_fftshifted); % Specify x and y data for axes set(gca, 'YDir', 'normal'); % Correct the y-axis direction cbar1 = colorbar; 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.alpha), '; \phi = ', num2str(Params.phi),'$'],'fontsize',16,'interpreter','latex') +title(['2-D DDI Potential: $\theta = ',num2str(Params.theta), '; \eta = ', num2str(Params.eta),'$'],'fontsize',16,'interpreter','latex') % == Quantum Fluctuations term == % gammaQF = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2)); @@ -108,33 +109,33 @@ EpsilonK = double(imag(EpsilonK) ~= 0); % 'is figure(12) clf set(gcf,'Position',[50 50 950 750]) -imagesc(fftshift(Transf.kx)*l0, fftshift(Transf.ky)*l0, EpsilonK); % Specify x and y data for axes +imagesc(fftshift(Transf.kx)*1e-6, fftshift(Transf.ky)*1e-6, EpsilonK); % Specify x and y data for axes set(gca, 'YDir', 'normal'); % Correct the y-axis direction cbar1 = colorbar; 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.alpha), '; \phi = ', num2str(Params.phi),'$'],'fontsize',16,'interpreter','latex') +title(['2-D Dispersion: $\theta = ',num2str(Params.theta), '; \eta = ', num2str(Params.eta),'$'],'fontsize',16,'interpreter','latex') -%% +%% Cycle through angles -% Define values for alpha and phi -alpha_values = 0:5:90; % Range of alpha values (you can modify this) -phi_values = 0:2:90; % Range of phi values (you can modify this) +% Define values for theta and eta +theta_values = 0:10:90; % Range of theta values (you can modify this) +eta_values = 0:10:90; % Range of eta values (you can modify this) -% Set up VideoWriter object -v = VideoWriter('potential_movie', 'MPEG-4'); % Create a video object -v.FrameRate = 5; % Frame rate of the video -open(v); % Open the video file +% Set up VideoWriter object to produce a movie +% v = VideoWriter('potential_movie', 'MPEG-4'); % Create a video object +% v.FrameRate = 5; % Frame rate of the video +% open(v); % Open the video file -% Loop over alpha and phi values -for alpha = alpha_values - for phi = phi_values - % Update Params with current alpha and phi - Params.alpha = alpha; - Params.phi = phi; +% Loop over theta and eta values +for theta = theta_values + for eta = eta_values + % Update Params with current theta and eta + Params.theta = theta; + Params.eta = eta; - % Compute the potential for the current alpha and phi + % Compute the potential for the current theta and eta % == 2-D DDI Potential in k-space == % VDk = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth); VDk_fftshifted = fftshift(VDk); @@ -158,25 +159,25 @@ for alpha = alpha_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(12) + figure(13) clf set(gcf,'Position',[50 50 950 750]) - imagesc(fftshift(Transf.kx)*l0, fftshift(Transf.ky)*l0, EpsilonK); % Specify x and y data for axes + imagesc(fftshift(Transf.kx)*1e-6, fftshift(Transf.ky)*1e-6, EpsilonK); % Specify x and y data for axes set(gca, 'YDir', 'normal'); % Correct the y-axis direction cbar1 = colorbar; cbar1.Label.Interpreter = 'latex'; - xlabel('$k_x$','fontsize',16,'interpreter','latex'); - ylabel('$k_y$','fontsize',16,'interpreter','latex'); - title(['$\theta = ',num2str(Params.alpha), '; \phi = ', num2str(Params.phi),'$'],'fontsize',16,'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') % Capture the frame and write to video - frame = getframe(gcf); % Capture the current figure - writeVideo(v, frame); % Write the frame to the video + % frame = getframe(gcf); % Capture the current figure + % writeVideo(v, frame); % Write the frame to the video end end % Close the video file -close(v); +% close(v); %% function [Transf] = setupSpace(Params) @@ -216,14 +217,14 @@ end function VDk = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth) % == Calculating the DDI potential in Fourier space with appropriate cutoff == % - + % Angles of the dipole moment are defined in and away from the X-Z plane % 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; + QDsq = QX.^2*cos(Params.eta)^2 + QY.^2*sin(Params.eta)^2; % eta here is the azimuthal angle of the dipole moment (angle from the x-axis) % Bare interaction Fpar = -1 + 3*sqrt(pi)*QDsq.*erfcx(absQ)./absQ; % Scaled complementary error function erfcx(x) = e^(x^2) * erfc(x) @@ -231,7 +232,7 @@ function VDk = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth) Fpar(absQ == 0) = -1; % Full DDI - VDk = (Fpar*sin(Params.alpha)^2 + Fperp*cos(Params.alpha)^2); + VDk = (Fpar*sin(Params.theta)^2 + Fperp*cos(Params.theta)^2); % theta here is the polar angle of the dipole moment (angle from the z-axis) end function ret = computeTotalEnergyPerParticle(x, as, AtomNumberDensity, wz, lz, gs, add, gdd, PlanckConstantReduced) diff --git a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ExtractingKRoton.m b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ExtractingKRoton.m index 073acf3..437abf1 100644 --- a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ExtractingKRoton.m +++ b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ExtractingKRoton.m @@ -66,7 +66,7 @@ end % ====================================================================================================================================================== % -alpha = 0; % Polar angle of dipole moment +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)); @@ -88,7 +88,7 @@ for idx = 1:length(nadd2s) gdd = VacuumPermeability*DyMagneticMoment^2/3; MeanWidth = var_widths(jdx, idx) * lz; % Mean width of Gaussian ansatz - [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, alpha, phi); % DDI potential in k-space + [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)); @@ -154,4 +154,79 @@ 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 \ No newline at end of file +grid on + +%% +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); + 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 \ No newline at end of file diff --git a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/RIBForTiltedDipoles.m b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/RIBForTiltedDipoles.m index c967f59..d5e398a 100644 --- a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/RIBForTiltedDipoles.m +++ b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/RIBForTiltedDipoles.m @@ -62,7 +62,7 @@ end % ====================================================================================================================================================== % -alpha = 0; % Polar angle of dipole moment +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)); @@ -83,7 +83,7 @@ for idx = 1:length(nadd2s) gdd = VacuumPermeability*DyMagneticMoment^2/3; MeanWidth = var_widths(jdx, idx) * lz; % Mean width of Gaussian ansatz - [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, alpha, phi); % DDI potential in k-space + [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)); @@ -126,10 +126,10 @@ ylabel('Scattering length ($\times a_0$)','fontsize',16,'interpreter','latex'); title('Roton instability boundary','fontsize',16,'interpreter','latex') %% -function [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, alpha, phi) +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(alpha))^2 - 1) + ((3 * Go) .* ((sin(deg2rad(alpha))^2 .* sin(deg2rad(phi))^2) - cos(deg2rad(alpha))^2)); + 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 @@ -143,7 +143,4 @@ 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 - - - +end \ No newline at end of file diff --git a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/RIBForTiltedDipoles_TwoDirections.m b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/RIBForTiltedDipoles_TwoDirections.m new file mode 100644 index 0000000..7cc32ff --- /dev/null +++ b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/RIBForTiltedDipoles_TwoDirections.m @@ -0,0 +1,318 @@ +%% 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; + +%% Roton instability boundary for tilted dipoles + +wz = 2 * pi * 72.4; % 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.005:0.0025:0.5; +as_to_add = 0.3:0.025:0.95; +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 + +%% ====================================================================================================================================================== % + +figure(6) +clf +set(gcf,'Position',[50 50 1850 750]) + +theta = 66; % 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)); +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 + + [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 + +subplot(1, 2, 1); % 1 row, 2 columns, first subplot +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 +caxis([0 1]) +xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex'); +ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex'); + +title(['Along Y: $\theta = ',num2str(theta), '; \phi = ', num2str(phi),'$'],'fontsize',16,'interpreter','latex') + +theta = 66; % Polar angle of dipole moment +phi = 90; % 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 + + [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 + +% set(gcf,'Position',[50 50 950 750]) +subplot(1, 2, 2); % 1 row, 2 columns, first subplot +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 +caxis([0 1]) +xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex'); +ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex'); + +title(['Along X: $\theta = ',num2str(theta), '; \phi = ', num2str(phi),'$'],'fontsize',16,'interpreter','latex') + +%{ +imagesc(AtomNumber*1E-5, ScatteringLengths, instability_boundary); % Specify x and y data for axes +set(gca, 'YDir', 'normal'); % Correct the y-axis direction +cbar1 = colorbar; +cbar1.Label.Interpreter = 'latex'; +caxis([0 1]) +% ylabel(cbar1,'$(\times 10^{-31})$','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'); +%} + +sgtitle('Mean-field instability boundary','fontsize',16,'interpreter','latex') + +%% Cycle through angles + +% Define values for theta and phi +theta_values = 0:2:90; % Range of theta values (you can modify this) + +% Set up VideoWriter object to produce a movie +v = VideoWriter('rib_movie', 'MPEG-4'); % Create a video object +v.FrameRate = 5; % Frame rate of the video +open(v); % Open the video file + +for theta = theta_values + figure(6) + clf + set(gcf,'Position',[50 50 1850 750]) + 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 + + [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 + + subplot(1, 2, 1); % 1 row, 2 columns, first subplot + 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 + caxis([0 1]) + xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex'); + ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex'); + + title(['Along Y: $\theta = ',num2str(theta), '; \phi = ', num2str(phi),'$'],'fontsize',16,'interpreter','latex') + + phi = 90; % 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 + + [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 + + % set(gcf,'Position',[50 50 950 750]) + subplot(1, 2, 2); % 1 row, 2 columns, first subplot + 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 + caxis([0 1]) + xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex'); + ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex'); + + title(['Along X: $\theta = ',num2str(theta), '; \phi = ', num2str(phi),'$'],'fontsize',16,'interpreter','latex') + + %{ + imagesc(AtomNumber*1E-5, ScatteringLengths, instability_boundary); % Specify x and y data for axes + set(gca, 'YDir', 'normal'); % Correct the y-axis direction + cbar1 = colorbar; + cbar1.Label.Interpreter = 'latex'; + caxis([0 1]) + % ylabel(cbar1,'$(\times 10^{-31})$','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'); + %} + + % Capture the frame and write to video + frame = getframe(gcf); % Capture the current figure + writeVideo(v, frame); % Write the frame to the video + + % sgtitle('Mean-field instability boundary','fontsize',16,'interpreter','latex') +end + +% Close the video file +close(v); + +%% +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 + + + diff --git a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ReproduceBlairBlakieResults.m b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ReproduceBlairBlakieResults.m index e40489f..f718de5 100644 --- a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ReproduceBlairBlakieResults.m +++ b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/ReproduceBlairBlakieResults.m @@ -34,7 +34,7 @@ wz = 2 * pi * 72.4; lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length as = 102.515 * BohrRadius; % Scattering length Trapsize = 7.5815 * lz; % Trap is assumed to be a box of finite extent , given here in units of the harmonic oscillator length -alpha = 0; % Polar angle of dipole moment +theta = 0; % Polar angle of dipole moment phi = 0; % Azimuthal angle of momentum vector MeanWidth = 5.7304888515 * lz; % Mean width of Gaussian ansatz k = linspace(0, 2e6, 1000); % Vector of magnitudes of k vector @@ -47,7 +47,7 @@ eps_dd = add/as; gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength gdd = VacuumPermeability*DyMagneticMoment^2/3; -[Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, alpha, phi); % DDI potential in k-space +[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)); @@ -74,7 +74,7 @@ AtomNumber = 1E5; wz = 2 * pi * 72.4; % Trap frequency in the tight confinement direction lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length Trapsize = 7.5815 * lz; % Trap is assumed to be a box of finite extent , given here in units of the harmonic oscillator length -alpha = 0; % Polar angle of dipole moment +theta = 0; % Polar angle of dipole moment phi = 0; % Azimuthal angle of momentum vector MeanWidth = 5.7304888515 * lz; % Mean width of Gaussian ansatz k = linspace(0, 2e6, 1000); % Vector of magnitudes of k vector @@ -92,7 +92,7 @@ for idx = 1:length(ScatteringLengths) gs = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as; % Contact interaction strength gdd = VacuumPermeability*DyMagneticMoment^2/3; - [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, alpha, phi); % DDI potential in k-space + [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)); @@ -128,7 +128,7 @@ legend('location', 'northwest','fontsize',16, 'Interpreter','latex') wz = 2 * pi * 72.4; % Trap frequency in the tight confinement direction lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length -alpha = 0; % Polar angle of dipole moment +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 @@ -154,7 +154,7 @@ for idx = 1:length(nadd2s) gdd = VacuumPermeability*DyMagneticMoment^2/3; MeanWidth = var_widths(idx) * lz; % Mean width of Gaussian ansatz - [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, alpha, phi); % DDI potential in k-space + [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)); @@ -249,7 +249,7 @@ ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex'); % ====================================================================================================================================================== % -alpha = 0; % Polar angle of dipole moment +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)); @@ -270,7 +270,7 @@ for idx = 1:length(nadd2s) gdd = VacuumPermeability*DyMagneticMoment^2/3; MeanWidth = var_widths(jdx, idx) * lz; % Mean width of Gaussian ansatz - [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, alpha, phi); % DDI potential in k-space + [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)); @@ -312,10 +312,10 @@ title('Roton instability boundary','fontsize',16,'interpreter','latex') %} %% -function [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, alpha, phi) +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(alpha))^2 - 1) + ((3 * Go) .* ((sin(deg2rad(alpha))^2 .* sin(deg2rad(phi))^2) - cos(deg2rad(alpha))^2)); + 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