# Finding local minima of an energy defined by integration

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 as[wy]^(5/2) n^(3/2)/(75 Sqrt Pi^(7/4) (x*y*z)^(3/2)) (1 +


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[]} /. 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[]} /. 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? 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.

• Recommend that you use NumericQ rather than NumberQ. Try (# /@ {2, 2., Pi, E, x} )& /@ {NumberQ, NumericQ} – Bob Hanlon Oct 13 at 15:47
• 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 :) – halirutan Oct 16 at 0:53
• 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. – halirutan Oct 16 at 9:14
• @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$ – Alex Trounev Oct 16 at 21:02
• @DineshShankar If there is one typo, then it is possible there is another. – Alex Trounev Oct 17 at 18:29

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 Pi^(7/4) (x*y*z)^(3/2)) (1 +
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*)
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[]} /. 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[]} /. 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 == b2 + k2 et2,
b1 + k1 et1 == 0, b2 + k2 et2 == 0,
b2 + k2 et2 == 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[]} /. 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[]} /. 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 • 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? – Dinesh Shankar Oct 22 at 8:29
• @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. – Alex Trounev Oct 22 at 9:47
• I got it! Thank you for your help. – Dinesh Shankar Oct 22 at 10:39
• @DineshShankar You're welcome! – Alex Trounev Oct 22 at 11:10