Physical Constants

PlanckConstant=6.62606957 (10^(-34));
PlanckConstantReduced=(PlanckConstant/(2  π));
FineStructureConstant=7.2973525698 (10^(-3));
ElectronMass=9.10938291 (10^(-31));
GravitationalConstant=6.67384 (10^(-11));
ProtonMass=1.672621777 (10^(-27));
AtomicMassUnit=1.66053878283 (10^(-27));
BohrRadius=0.52917721092 (10^(-10));
BohrMagneton=927.400968 (10^(-26));
BoltzmannConstant=1.3806488 (10^(-23));
StandardGravityAcceleration=9.80665;
SpeedOfLight=299792458;
StefanBoltzmannConstant=5.670373 (10^(-8));
ElectronCharge=1.602176565 (10^(-19));
VacuumPermeability=4π (10^(-7));
DielectricConstant=(1/((SpeedOfLight^2) VacuumPermeability));
ElectronGyromagneticFactor=-2.00231930436153;
AvogadroConstant=6.02214129 (10^23);

Dysprosium-specific constants

m = 164 AtomicMassUnit;
d=9.93 BohrMagneton;

Definitions

g[as_]:= ((4π (PlanckConstantReduced^2)as)/m) ; (*Contact coupling strength as a function of the s-wave scattering length*)
(Subscript[C, dd])= VacuumPermeability (d^2); (*Dipolar coupling strength*)
(Subscript[ϵ, dd ])[as_]:=((Subscript[C, dd])/(3g[as])); (*Coupling strength ratio*)

(* Analytic solution under the Thomas Fermi approximation - neglecting the kinetic energy term in the time-independent GPE for a dipolar gas in a harmonic trap with cylindrical symmetry*)

(*k below is the aspect ratio, k = Rx/Rz*)
h[κ_]:=((1/(√(1-(κ^2)))))Log[((1+√(1-(κ^2)))/(1-√(1-(κ^2))))];(*Function definition for an prolate spheroid, k < 1*)
j[κ_]:=((2/(√((κ^2)-1))))ArcTan[√((κ^2)-1)];(*Function definition for an oblate spheroid, k > 1*)
Ξ[κ_]:=If[κ<1,h[κ],j[κ]]; 
f[κ_]:= ((2+(κ^2)(4 - 3 Ξ[κ]))/(2(1-(κ^2)))); (*Monotonously decreasing function having values in the range 1 ≥ f(k) ≥ -2 passing through 0 at k=1*)
(*k is first obtained by solving the following transcendental equation*)
AspectRatioEquation[κ_, as_] := 3 (κ^2) (Subscript[ϵ, dd ])[as]((((γ^2)/2)+1)((f[κ])/(1-(κ^2)))-1)+((Subscript[ϵ, dd ])[as]-1)((κ^2)-(γ^2));
(*The numbers then can be plugged in to the following equations to get physical quantities*)
(Subscript[R, x])[κ_,as_]:=(((((15 g[as] n κ)/(4 π m ((Subscript[ω, x])^2)))(1 +(Subscript[ϵ, dd ])[as]((3/2)(((κ^2)f[κ])/(1-(κ^2)))-1))))^(1/5)); (*Radius of cloud in plane*)
(Subscript[R, z])[κ_, as_]:=(((Subscript[R, x])[κ, as])/κ); (*Radius of cloud orthogonal to the plane*)
(Subscript[n, o ])[κ_, as_]:=((15 n)/(8 π (((Subscript[R, x])[κ, as])^2) (Subscript[R, z])[κ, as])); (*Atom Number Density*)
ChemicalPotential[κ_, as_]:=g[as] (Subscript[n, o])[κ, as](1 -(Subscript[ϵ, dd ])[as]f[κ]); (*Chemical Potential*)

(*Dipolar potential*)
(Subscript[ϕ, dd])[κ_, as_]:= (((Subscript[n, o])[κ, as](Subscript[C, dd]))/3)(((ρ^2)/(((Subscript[R, x])[κ, as])^2))-((2(z^2))/((((Subscript[R, z])[κ])^2) ))-f[κ](1-(3/2)(((ρ^2)-2(z^2))/(2 (((Subscript[R, x])[κ, as])^2) - (((Subscript[R, z])[κ])^2)))));

Checks

Plot[Evaluate[(Subscript[ϵ, dd ])[as]] /. as -> sl BohrRadius, {sl,120, 150}, PlotRange->Full, FrameLabel->{HoldForm[(Subscript[a, s])],HoldForm[(Subscript[ϵ, dd ])]},PlotLabel->None,LabelStyle->{16,Black,Bold},GridLines->Automatic,Frame->True,FrameStyle->Directive[Black,Thickness[Large]]]
Plot[Evaluate[Ξ[κ]], {κ, 1, 2}, PlotRange->Full, FrameLabel->{HoldForm[κ], HoldForm[Ξ]},PlotLabel->None,LabelStyle->{16,Black,Bold},GridLines->Automatic,Frame->True,FrameStyle->Directive[Black,Thickness[Large]]]
LogLinearPlot[Evaluate[f[κ]], {κ, 0.01, 100}, PlotRange->Full, FrameLabel->{HoldForm[κ], HoldForm[f]},PlotLabel->None,LabelStyle->{16,Black,Bold},GridLines->Automatic,Frame->True,FrameStyle->Directive[Black,Thickness[Large]]]

Number density and chemical potential from experimental parameters

AtomNumber= (10^4);
(Subscript[v, x]) = 50;
(Subscript[v, z]) = 500;

SL = Table[i BohrRadius,{i,130, 200}]; (*Range of scattering lengths*)
Plot[AspectRatioEquation[κ, as]/. {γ-> ((Subscript[v, z])/(Subscript[v, x])), as-> 130 BohrRadius}, {κ, 0, 250.0}, PlotRange->Full,FrameLabel->{HoldForm[κ],HoldForm[Values]},PlotLabel->None,LabelStyle->{16,Black,Bold},GridLines->Automatic,Frame->True,FrameStyle->Directive[Black,Thickness[Large]]]
FindRoot[AspectRatioEquation[κ,as]/.{γ-> ((Subscript[v, z])/(Subscript[v, x])),as->130 BohrRadius },{κ,100}]

AR=Table[κ/.FindRoot[AspectRatioEquation[κ,as]/.{γ-> ((Subscript[v, z])/(Subscript[v, x])),as->sl },{κ,100}],{sl,SL}]; (*Corresponding aspect ratios, obtained by solving the transcendental equation*)

ND = MapThread[(Subscript[n, o ]),{AR,SL}]/. {(Subscript[ω, x]) -> 2 π (Subscript[v, x]), n-> AtomNumber};  (*Corresponding number density values*)

CPInKelvin = ((MapThread[ChemicalPotential,{AR,SL}](10^9))/BoltzmannConstant)/. {(Subscript[ω, x]) -> 2 π (Subscript[v, x]), n-> AtomNumber} ; (*Corresponding chemical potential values in Kelvin*)
CPInHz = ((MapThread[ChemicalPotential,{AR,SL}](10^(-3)))/PlanckConstant)/. {(Subscript[ω, x]) -> 2 π (Subscript[v, x]), n-> AtomNumber} ; (*Corresponding chemical potential values in Hz*)

Epsilondd = Evaluate[(Subscript[ϵ, dd ])[as]] /. as -> SL;

data=Transpose@{Epsilondd,AR};
ListLinePlot[data,PlotRange->Full,FrameLabel->{HoldForm[(Subscript[ϵ, dd])],HoldForm[κ]},PlotLabel->None,LabelStyle->{16,Black,Bold},GridLines->Automatic,Frame->True,FrameStyle->Directive[Black,Thickness[Large]]]

(*data=Transpose@{(SL/BohrRadius),ND};*)
data=Transpose@{Epsilondd,ND (10^(-20))};
ListLinePlot[data,PlotRange->Full,FrameLabel->{HoldForm[(Subscript[ϵ, dd])],HoldForm[(Subscript[n, o ])(x (10^20)(m^(-3)))]},PlotLabel->None,LabelStyle->{16,Black,Bold},GridLines->Automatic,Frame->True,FrameStyle->Directive[Black,Thickness[Large]]]

(*data=Transpose@{(SL/BohrRadius),CP};*)
data=Transpose@{Epsilondd,CPInKelvin};
ListLinePlot[data,PlotRange->Full,FrameLabel->{HoldForm[(Subscript[ϵ, dd])],HoldForm[μ[nK]]},PlotLabel->None,LabelStyle->{16,Black,Bold},GridLines->Automatic,Frame->True,FrameStyle->Directive[Black,Thickness[Large]]]

data=Transpose@{Epsilondd,CPInHz};
ListLinePlot[data,PlotRange->Full,FrameLabel->{HoldForm[(Subscript[ϵ, dd])],HoldForm[μ[kHz]]},PlotLabel->None,LabelStyle->{16,Black,Bold},GridLines->Automatic,Frame->True,FrameStyle->Directive[Black,Thickness[Large]]]

Cloud size

AtomNumber= (10^4);
(Subscript[v, x]) = 50;
(Subscript[v, z]) = 500;

SL = Table[i BohrRadius,{i,130, 200}]; (*Range of scattering lengths*)
AR=Table[κ/.FindRoot[AspectRatioEquation[κ,as]/.{γ-> ((Subscript[v, z])/(Subscript[v, x])),as->sl },{κ,100}],{sl,SL}]; (*Corresponding aspect ratios, obtained by solving the transcendental equation*)
RADX = MapThread[(Subscript[R, x]),{AR,SL}]/. {(Subscript[ω, x]) -> 2 π (Subscript[v, x]), n-> AtomNumber};  (*Corresponding radii in x*)
RADZ = MapThread[(Subscript[R, z]),{AR,SL}]/. {(Subscript[ω, x]) -> 2 π (Subscript[v, x]), n-> AtomNumber}; (*Corresponding radii in z*)
Epsilondd = Evaluate[(Subscript[ϵ, dd ])[as]] /. as -> SL;

datax=Transpose@{Epsilondd,RADX (10^6)};
dataz=Transpose@{Epsilondd,RADZ (10^9)};
ListLinePlot[datax,PlotRange->Full,FrameLabel->{HoldForm[(Subscript[ϵ, dd])],HoldForm[Radius[μm]]},PlotLabel->None,LabelStyle->{16,Black,Bold},GridLines->Automatic,Frame->True,FrameStyle->Directive[Black,Thickness[Large]]]
ListLinePlot[dataz,PlotRange->Full,FrameLabel->{HoldForm[(Subscript[ϵ, dd])],HoldForm[Radius[nm]]},PlotLabel->None,LabelStyle->{16,Black,Bold},GridLines->Automatic,Frame->True,FrameStyle->Directive[Black,Thickness[Large]]]