3
$\begingroup$

This question is also asked here and here.

I would like to reproduce the two solid curves in Fig. 1a of this paper. The total energy is given in Appendix A, and the curve is obtained by minimizing the energy with respect to the variational parameters $\sigma_k$ and $d$. I'm not sure what dimensionless variables they used on their code, but I'll expresses the energy using $(x,y,z)=(\sigma_x,\sigma_y,\sigma_z)/\bar{a}$; $(wx,wy,wz)=(\omega_x,\omega_y,\omega_z)/\bar{\omega}$; $as=a_s/\bar{a}$ and so on (Notice that $\bar{a}$ has length dimension and $\bar{\omega}$ has the dimension inverse time).

ClearAll["Global`*"] Idip[X_?NumericQ, Y_?NumericQ, u_?NumericQ] := Re[Exp[-u^2/2] - 3 X*Y/(1 - X^2)^(3/2) NIntegrate[ v^2 Exp[-u^2 X^2 v^2/(2 (1 - X^2) (1 - v^2))]/(Sqrt[1 - v^2] Sqrt[1 - v^2 (1 - Y^2)/(1 - X^2)]), {v, 0, Sqrt[1 - X^2]}]]; f[X_?NumericQ, Y_?NumericQ] := Re[Idip[X, Y, 0]]; Iqf[u_?NumericQ] := Re[2 Exp[-5 u^2/8]/Sqrt[Pi] NIntegrate[ Exp[-l^2] Cosh[Sqrt[2/5] u*l]^(5/2), {l, 0, Infinity}]]; ettot[x_, y_, z_, u_, wy_] := h*w[wy]/Kb (1/4 (1/x^2 + 1/y^2 + 1/z^2) + 1/4 (wx[wy]^2 x^2 (1 + u^2/2) + wy^2 y^2 + wz[wy]^2 z^2) + as[wy]*n/(2 Sqrt[2 Pi] (x*y*z)) (1 + Exp[-u^2/2]) + add[wy]*n/(2 Sqrt[2 Pi] (x*y*z)) (-f[x/z, y/z] - Idip[x/z, y/z, u]) + 512*Sqrt[2] as[wy]^(5/2) n^(3/2)/(75 Sqrt[5] Pi^(7/4) (x*y*z)^(3/2)) (1 + 3 add[wy]^2/as[wy]^2) Iqf[u]);

Obs: In the code, $wx$ and $wz$ are dimensionless, but $wy$ has the dimension inverse time. Morover, since they have used Kelvin units, I divided the total energy by $Kb$.

In order to find the local minima I have used the following code based on the one found here:

a0 = 5.29*10^(-11); (* Bohr radius *) h = 1.054*10^(-34); (* Reduced Planck's constant *) M = 163.9*1.66*10^(-27); (* Dy-164 mass in kg *) wxx = 2 Pi*70; (* Experimental value of x-frequency *) wzz = 2 Pi*1000; (* Experimental value of z-frequency *) w[wy_] := (wxx*wy*wzz)^(1/3); a[wy_] := Sqrt[h/(M*w[wy])]; wx[wy_] := wxx/w[wy]; (* Dimensionless x-frenquency*) wz[wy_] := wzz/w[wy]; (* Dimensionless w-frenquency*) add[wy_] := 131*a0/a[wy]; (* Dimensionless Dipolar length*) as[wy_] := 70*a0/a[wy]; (* Dimensionless contact length*) n = 10^4; Kb = 1.38*10^(-23); (* Boltzmann constant *) ddata1 = Table[minsol1 =FindMinimum[{ettot[x, y, z, 0, wy/w[wy]], x > 0 && y > 0 && z > 0}, {{x, 1.01}, {y, 1.012}, {z, 0.14}}, Method -> Automatic, PrecisionGoal -> Automatic, AccuracyGoal -> Automatic]; Re[{x, y, z, wy/(2 Pi) , minsol1[[1]]} /. Last[minsol1]], {wy, 2 Pi*100, 2 Pi*250, 2 Pi*20}]; ddata2 = Table[minsol2 =FindMinimum[{ettot[x, y, z, d/x, wy/w[wy]], x > 0 && y > 0 && z > 0 && d > x}, {{x, 1.01}, {y, 1.012}, {z, 0.14}, {d, 1.02}}, Method -> Automatic, PrecisionGoal -> Automatic, AccuracyGoal -> Automatic]; Re[{x, y, z, d, wy/(2 Pi) , minsol2[[1]]} /. Last[minsol2]], {wy, 2 Pi*100, 2 Pi*250, 2 Pi*20}]; ListPlot[{ddata1[[All, {4, 5}]], ddata2[[All, {5, 6}]]}, Joined -> True, PlotRange -> {Automatic, Automatic}]

I get two pretty nice curves using this code. However, I expected the curves to match at $f_y = 200$ Hz Moreover, at $f_y < 200$ Hz it had to be negative. Could helpe me?

enter image description here

Additional comment: There is probably an error in the paper, since neither Figure 1a nor Figure 6a can be obtained from the total energy. This is very strange, however, because the group is well respected in the area. Anyway, should I send an email to authors asking about it? They used probably C++ instead of Mathematica.

Improve this question
$\endgroup$
16
  • 1
    $\begingroup$ Recommend that you use NumericQ rather than NumberQ. Try (# /@ {2, 2., Pi, E, x} )& /@ {NumberQ, NumericQ} $\endgroup$ Commented Oct 13, 2019 at 15:47
  • 1
    $\begingroup$ There are many more errors. For instance, you define a function add[yw_] but inside ettot you simply use it as symbol add without arguments. Same for as and wz. I'm not even sure how you could calculate any results at all :) $\endgroup$ Commented Oct 16, 2019 at 0:53
  • 1
    $\begingroup$ Please copy your code into Mathematica and verify that it is syntactically correct. For instance, there is (2 Pi]) in your code (a line starting with Re) which cannot be correct. $\endgroup$ Commented Oct 16, 2019 at 9:14
  • 1
    $\begingroup$ @DineshShankar Why did you expect that ddata1[[All, {4, 5}]], ddata2[[All, {5, 6}]] should match when $f_y=250 Hz$? In Fig. 1, the curves intersect at $f_y=200$ $\endgroup$ Commented Oct 16, 2019 at 21:02
  • 1
    $\begingroup$ @DineshShankar If there is one typo, then it is possible there is another. $\endgroup$ Commented Oct 17, 2019 at 18:29

1 Answer 1

Reset to default
5
+100
$\begingroup$

After studying the article and theses, as well as article Ground-state properties and elementary excitations of quantum droplets in dipolar Bose-Einstein condensates from which the ansatz was taken, I came to the conclusion that the ansatz may not have been calibrated to find the magnitude of the energy. Therefore, we calibrate the ansatz using the data in Fig. 1:

ClearAll["Global`*"] Idip[X_?NumericQ, Y_?NumericQ, u_?NumericQ] := Re[Exp[-u^2/2] - 3 X*Y/(1 - X^2)^(3/2) NIntegrate[ v^2 Exp[-u^2 X^2 v^2/(2 (1 - X^2) (1 - v^2))]/(Sqrt[ 1 - v^2] Sqrt[1 - v^2 (1 - Y^2)/(1 - X^2)]), {v, 0, Sqrt[1 - X^2]}]]; f[X_?NumericQ, Y_?NumericQ] := Re[Idip[X, Y, 0]]; Iqf[u_?NumericQ] := Re[2 Exp[-5 u^2/8]/Sqrt[Pi] NIntegrate[ Exp[-l^2] Cosh[Sqrt[2/5] u*l]^(5/2), {l, 0, Infinity}]]; ettot[x_, y_, z_, u_, wy_] := h*w[wy]/Kb (1/4 (1/x^2 + 1/y^2 + 1/z^2) + 1/4 (wx[wy]^2 x^2 (1 + u^2/2) + wy^2 y^2 + wz[wy]^2 z^2) + as[wy]*n/(2 Sqrt[2 Pi] (x*y*z)) (1 + Exp[-u^2/2]) + add[wy]*n/(2 Sqrt[2 Pi] (x*y*z)) (-f[x/z, y/z] - Idip[x/z, y/z, u]) + 512*Sqrt[ 2] as[wy]^(5/ 2) n^(3/2)/(75 Sqrt[5] Pi^(7/4) (x*y*z)^(3/2)) (1 + 3/2 add[wy]^2/as[wy]^2) Iqf[u]); a0 = 5.29*10^(-11);(*Bohr radius*)h = 1.054*10^(-34);(*Reduced Planck's constant*)M = 163.9*1.66*10^(-27);(*Dy-164 mass in kg*)wxx = 2 Pi*70;(*Experimental value of x-frequency*)wzz = 2 Pi*1000;(*Experimental value of z-frequency*) w[wy_] := (wxx*wy*wzz)^(1/3); a[wy_] := Sqrt[h/(M*w[wy])]; wx[wy_] := wxx/w[wy];(*Dimensionless x-frenquency*) wz[wy_] := wzz/w[wy];(*Dimensionless w-frenquency*) add[wy_] := 131*a0/a[wy];(*Dimensionless Dipolar length*) as[wy_] := 70*a0/a[wy];(*Dimensionless contact length*)n = 10^4; Kb = 1.38*10^(-23)*10^-9(*nK*);(*Boltzmann constant*)ddata1 = Table[minsol1 = FindMinimum[{ettot[x, y, z, 0, wy/w[wy]], x > 0 && y > 0 && z > 0}, {{x, 1.01}, {y, 1.012}, {z, 0.14}}, Method -> Automatic, PrecisionGoal -> Automatic, AccuracyGoal -> Automatic]; Re[{x, y, z, wy/(2 Pi), minsol1[[1]]} /. Last[minsol1]], {wy, 2 Pi*100, 2 Pi*250, 2 Pi*20}]; ddata2 = Table[ minsol2 = FindMinimum[{ettot[x, y, z, d/x, wy/w[wy]], x > 0 && y > 0 && z > 0 && d > x}, {{x, 1.01}, {y, 1.012}, {z, 0.14}, {d, 1.02}}, Method -> Automatic, PrecisionGoal -> Automatic, AccuracyGoal -> Automatic]; Re[{x, y, z, d, wy/(2 Pi), minsol2[[1]]} /. Last[minsol2]], {wy, 2 Pi*100, 2 Pi*250, 2 Pi*20}];

Compose a new ansatz $E_{total}=k E_{tot}+b$ and find the constants $b,k$

et1 = Interpolation[ddata1[[All, {4, 5}]]]; et2 = Interpolation[ddata2[[All, {5, 6}]]]; sol = NSolve[{b1 + k1 et1[200] == b2 + k2 et2[200], b1 + k1 et1[160] == 0, b2 + k2 et2[144] == 0, b2 + k2 et2[220] == 20}, {b1, b2, k1, k2}]

Now we optimize a new ansatz

{b2, k2} = First[{b2, k2} /. sol] ddata2 = Table[ minsol2 = FindMinimum[{b2 + k2 ettot[x, y, z, d/x, wy/w[wy]], x > 0 && y > 0 && z > 0 && d > x}, {{x, 1.01}, {y, 1.012}, {z, 0.14}, {d, 1.02}}, Method -> Automatic, PrecisionGoal -> Automatic, AccuracyGoal -> Automatic]; Re[{x, y, z, d, wy/(2 Pi), minsol2[[1]]} /. Last[minsol2]], {wy, 2 Pi*100, 2 Pi*250, 2 Pi*20}]; {b1, k1} = First[{b1, k1} /. sol] ddata1 = Table[ minsol1 = FindMinimum[{b1 + k1 ettot[x, y, z, 0, wy/w[wy]], x > 0 && y > 0 && z > 0}, {{x, 1.01}, {y, 1.012}, {z, 0.14}}, Method -> Automatic, PrecisionGoal -> Automatic, AccuracyGoal -> Automatic]; Re[{x, y, z, wy/(2 Pi), minsol1[[1]]} /. Last[minsol1]], {wy, 2 Pi*100, 2 Pi*250, 2 Pi*20}];

Finally, build the curves and plot the data. It looks similar to Fig. 1, but there are differences in curvature. Note that you can not do re-optimization, because the data of re-optimization accurately lie on the curves obtained in the first optimization Figure 1

Improve this answer
$\endgroup$
4
  • $\begingroup$ This is a pretty nice answer. But let me ask you one quick question regarding your answer. It is possible only because you habe already knowed some information about the original Figure, is it? $\endgroup$ Commented Oct 22, 2019 at 8:29
  • $\begingroup$ @DineshShankar We can use numerical data obtained by solving the equation (1). It is only necessary to determine the shift and the intersection point $f_y=200$ Hz.This is determined from an experiment I think. $\endgroup$ Commented Oct 22, 2019 at 9:47
  • $\begingroup$ I got it! Thank you for your help. $\endgroup$ Commented Oct 22, 2019 at 10:39
  • $\begingroup$ @DineshShankar You're welcome! $\endgroup$ Commented Oct 22, 2019 at 11:10

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.