# Set of three self-consistent coupled 'gap' equations

For a function of the kind of 5 variables, where this function must satisfy the following self-consistent coupled 'gap' equations.

I want to reproduce the graphics in the paper "Phases of QCD: lattice thermodynamics and a field theoretical model" by Claudia Ratti, Michael A. Thaler and Wolfram Weise.

For the coding in Mathematica, I've replaced $$\Phi$$ by "Phi" and $$\Phi^{*}$$ by "phi".

Nc = 3;
Nf = 2;

\[CapitalLambda] = 587.9;    (*MeV*)
m = 5.6;         (*MeV*)
Tc = 215;
G = 2.44/\[CapitalLambda]^2;
a0 = 3.51;
a1 = -2.47;
a2 = 15.22;
b4 = -1.75;
T0 = 270;

Ep[p_, M_] := Sqrt[p^2 + M^2];

B2[T_] := a0 + a1 (T0/T) + a2 (T0/T)^2;
B4[T_] := b4 (T0/T)^3;

U[\[Phi]_, \[CapitalPhi]_, T_] = (-(1/2)*B2[T]*\[Phi]*\[CapitalPhi] +  B4[T]*Log[1 - 6*\[Phi]*\[CapitalPhi] + 4*(\[CapitalPhi]^3 + \[Phi]^3) - 3*(\[Phi]*\[CapitalPhi])^2])*T^4;

\[CapitalOmega][\[Mu]_, T_, M_, \[CapitalPhi]_, \[Phi]_] := \[CapitalOmega]0[M] + U[\[Phi], \[CapitalPhi], T] - 2*Nf*NIntegrate[p^2/(2 \[Pi]^2)T (Log[1 + E^(-((3 Ep[p, M])/T) - 3 \[Mu]/T) + 3 E^(-((2 Ep[p, M])/T) - 2 \[Mu]/T) \[CapitalPhi] + 3 E^(-(Ep[p, M]/T) - \[Mu]/T) \[Phi]] + Log[1 + E^(-((3 Ep[p, M])/T) + 3 \[Mu]/T) + 3 E^(-(Ep[p, M]/T) + \[Mu]/T) \[CapitalPhi] + 3 E^(-((2 Ep[p, M])/T) + 2 \[Mu]/T) \[Phi]]), {p, 0, Infinity}];


and the gap equations are written as:


\[Delta]\[CapitalOmega]\[Delta]M[\[Mu]_, T_, M_, \[CapitalPhi]_, \[Phi]_] := (M - m)/G - 2*Nf*Nc*(4 \[Pi])/(2 \[Pi])^3*M (1/2 \[CapitalLambda] Sqrt[M^2 + \[CapitalLambda]^2] - 1/2 M^2 ArcTanh[\[CapitalLambda]/Sqrt[M^2 + \[CapitalLambda]^2]]) - 2*Nf*NIntegrate[-((3 p^2 (1 + E^((2 (\[Mu] + Ep[p, M]))/T) \[Phi] + 2 E^((\[Mu] + Ep[p, M])/T) \[CapitalPhi])\!$$\*SuperscriptBox[\(Ep$$, TagBox[RowBox[{"(", RowBox[{"0", ",", "1"}], ")"}],Derivative],MultilineFunction->None]\)[p, M])/(2 \[Pi]^2 (1 + E^((3 (\[Mu] + Ep[p, M]))/T) + 3 E^((2 (\[Mu] + Ep[p, M]))/T) \[Phi] + 3 E^((\[Mu] + Ep[p, M])/T) \[CapitalPhi]))) - (3 E^(\[Mu]/T) p^2 (E^((2 \[Mu])/T) + 2 E^((\[Mu] + Ep[p, M])/T) \[Phi] +  E^((2 Ep[p, M])/T) \[CapitalPhi])\!$$\*SuperscriptBox[\(Ep$$, TagBox[RowBox[{"(", RowBox[{"0", ",", "1"}], ")"}],Derivative],MultilineFunction->None]\)[p, M])/( 2 \[Pi]^2 (E^((3 \[Mu])/T) + E^((3 Ep[p, M])/T) +  3 E^((2 \[Mu] + Ep[p, M])/T) \[Phi] +  3 E^((\[Mu] + 2 Ep[p, M])/T) \[CapitalPhi])), {p, 0, Infinity}];

\[Delta]\[CapitalOmega]\[Delta]\[CapitalPhi][\[Mu]_, T_, M_, \[CapitalPhi]_, \[Phi]_] :=  T^4 (-(1/2) \[Phi] B2[T] - (6 (\[Phi] + \[Phi]^2 \[CapitalPhi] - 2 \[CapitalPhi]^2) B4[T])/(1 + 4 \[Phi]^3 - 6 \[Phi] \[CapitalPhi] -  3 \[Phi]^2 \[CapitalPhi]^2 + 4 \[CapitalPhi]^3)) -  2*Nf*NIntegrate[(3 E^((\[Mu] + Ep[p, M])/T) p^2 T (1/(E^((3 \[Mu] - Ep[p, M])/T) + E^((2 Ep[p, M])/T) + 3 E^((2 \[Mu])/T) \[Phi] +  3 E^((\[Mu] + Ep[p, M])/T) \[CapitalPhi]) + 1/( 1 + E^((3 (\[Mu] + Ep[p, M]))/T) +  3 E^((2 (\[Mu] + Ep[p, M]))/T) \[Phi] +  3 E^((\[Mu] + Ep[p, M])/T) \[CapitalPhi])))/( 2 \[Pi]^2), {p, 0, Infinity}];

\[Delta]\[CapitalOmega]\[Delta]\[Phi][\[Mu]_, T_,  M_, \[CapitalPhi]_, \[Phi]_] := T^4 (-(1/2) \[CapitalPhi] B2[T] - (6 (\[CapitalPhi] + \[Phi] (-2 \[Phi] + \[CapitalPhi]^2)) B4[T])/(1 + 4 \[Phi]^3 - 6 \[Phi] \[CapitalPhi] - 3 \[Phi]^2 \[CapitalPhi]^2 + 4 \[CapitalPhi]^3)) - 2*Nf*NIntegrate[(3 E^((2 \[Mu] + Ep[p, M])/T) p^2 T (E^(Ep[p, M]/T)/( 1 + E^((3 (\[Mu] + Ep[p, M]))/T) +  3 E^((2 (\[Mu] + Ep[p, M]))/T) \[Phi] + 3 E^((\[Mu] + Ep[p, M])/T) \[CapitalPhi]) + 1/(E^((3 \[Mu])/T) + E^((3 Ep[p, M])/T) +  3 E^((2 \[Mu] + Ep[p, M])/T) \[Phi] + 3 E^((\[Mu] + 2 Ep[p, M])/T) \[CapitalPhi])))/(2 \[Pi]^2), {p, 0, Infinity}];



I tried to plot these graphs by using FindRoot of all these 3 gap equations and defining the values of Phi, phi and M, and the result will be a function of mu and T. Later I take a desired value for mu and do a Table varying T. For this, I wrote a solution of the type:

sol\[Phi][\[Mu]_, T_] :=  x /. FindRoot[{Re[\[Delta]\[CapitalOmega]\[Delta]\[CapitalPhi][\[Mu]\, T, w, y, x]] == 0, Re[\[Delta]\[CapitalOmega]\[Delta]\[Phi][\[Mu], T, w, y, x]] == 0, Re[\[Delta]\[CapitalOmega]\[Delta]M[\[Mu],T, w, y, x]] == 0}, {{x, 0, 0.8}, {y, 0, 0.8}, {w, 10, 350}}];


Where "dfdPhi", "dfdphi" and "dfdM" are the gap equations above, but they are simply derivatives of $$f$$ with respect to $$\Phi$$, $$\Phi^{*}$$ and $$M$$ respectively. In the FindRoot, for solutions of "Phi" I defined a variable "x" from $$0.35$$ to $$0.75$$, for "phi" I defined a variable "y" from $$0.35$$ to $$0.75$$ as well, and for "M" I defined a variable "w" from $$0$$ to $$350$$, which are approximately the range of values the figure.

I tried, as said before, to use Table for plotting Phi vs T, for a specific mu, e.g., mu=340 (the problematic curve), as follows:

Table[{T, sol\[CapitalPhi][340, T]}, {T, 1.1, 302, 50}]


The result is for a few points is:

Out=  {{1.1, -2.91187*10^-37}, {51.1, 0.0217319}, {101.1,
0.0569066}, {151.1, 0.159099}, {201.1, 0.421985}, {251.1,
0.62148}, {301.1, 0.7305}}


The plotted points looks like this:

And for more points::

Table[{T, solPhi[340, T]}, {T, 1, 301, 10}]


The result:

Out=  {{1.1, -2.91187*10^-37}, {11.1, 0.00453317}, {21.1,
0.00861518}, {31.1, 0.0127729}, {41.1, 0.0171036}, {51.1,
0.0217319}, {61.1, 0.026818}, {71.1, 0.0325687}, {81.1,
0.0392518}, {91.1, 0.0472141}, {101.1, 0.0569066}, {111.1,
0.0689161}, {121.1, 0.0840072}, {131.1, 0.103173}, {141.1,
0.127685}, {151.1, 0.159099}, {161.1, 0.19903}, {171.1,
0.248304}, {181.1, 0.305262}, {191.1, 0.364994}, {201.1,
0.421985}, {211.1, 0.473195}, {221.1, 0.694858}, {231.1,
0.557162}, {241.1, -0.0897173}, {251.1, 0.62148}, {261.1,
0.648192}, {271.1, 0.704975}, {281.1, -0.3734}, {291.1,
0.712885}, {301.1, 0.7305}}


with plotting:

As you can see, nothing like the original figure for $$\mu=340$$ from the paper. It looks like more for $$\mu=200$$.

Am I doing something wrong or any suggestions? Anything will be welcome!! Thanks.

• Short, focused questions that present a minimal working example get more attention and replies in general. Your question is physics-specific so you will probably have to wait a bit longer for an interested and knowledgeable respondent to take it up.
– Syed
Commented Mar 18, 2022 at 7:47
• Your code does not evaluate successfully with v13.0.1 (NIntegrate::izero and FindRoot::nlnum errors which then produce ReplaceAll::reps errors) Commented Mar 18, 2022 at 15:54
• Parameter T0 is not defined. Also in your code there is defined function dfPhi while in FindRoot it used with name dfdPhi. Function dfdM diverges, probably it defined with typo. Commented Mar 18, 2022 at 16:45
• @AlexTrounev I gonna check it and edit the question...I simplified and changed the notation to make it shorter here. Commented Mar 18, 2022 at 22:54
• @BobHanlon yes, probably. My version is 12.0, I noticed also it doesn't work in version 11. Commented Mar 18, 2022 at 22:57

With this code we can reproduce picture for mass and $$\Phi$$.

Nc = 3;
Nf = 2;

Lambda = \[CapitalLambda] = 651;(*MeV*)m = 5.6;(*MeV*)Tc = 215;
G = 10.08 10^-6;
a0 = 3.51;
a1 = -2.47;
a2 = 15.22;
b4 = -1.75;
T0 = 270;

Ep[p_, M_] := Sqrt[p^2 + M^2];
\[CapitalOmega]0[M_] := (M - m)^2/(2 G);

B2[T_] := a0 + a1 (T0/T) + a2 (T0/T)^2;
B4[T_] := b4 (T0/T)^3;

U[\[Phi]_, \[CapitalPhi]_,
T_] = (-(1/2)*B2[T]*\[Phi]*\[CapitalPhi] +
B4[T]*Log[
1 - 6*\[Phi]*\[CapitalPhi] + 4*(\[CapitalPhi]^3 + \[Phi]^3) -
3*(\[Phi]*\[CapitalPhi])^2])*T^4;

\[CapitalOmega][\[Mu]_, T_,
M_, \[CapitalPhi]_, \[Phi]_] := \[CapitalOmega]0[M] +
U[\[Phi], \[CapitalPhi], T] -
2*Nf/(2 \[Pi]^2) T *
NIntegrate[
p^2 (Log[
1 + E^(-3 ( Ep[p, M] - \[Mu])/T) +
3 E^(-(Ep[p, M] - \[Mu])/T) \[CapitalPhi] +
3 E^(-2 (Ep[p, M] - \[Mu])/T) \[Phi]] +
Log[1 + E^(-3 ( Ep[p, M] + \[Mu])/T) +
3 \[Phi] E^(-(Ep[p, M] + \[Mu])/T) +
3 E^(-2 (Ep[p, M] + \[Mu])/T) \[CapitalPhi] ]), {p, 0,
Infinity}];

\[Delta]\[CapitalOmega]\[Delta]M[\[Mu]_, T_,
M_, \[CapitalPhi]_, \[Phi]_] := (M - m)/G -
6*Nf*(4 \[Pi])/(2 \[Pi])^3*(1/4 Lambda M Sqrt[Lambda^2 + M^2] + (
Lambda M (2 Lambda^2 + M^2))/(8 Sqrt[Lambda^2 + M^2]) + (
Lambda M^5)/(
8 (Lambda^2 + M^2)^(3/2) (1 - Lambda^2/(Lambda^2 + M^2))) -
1/2 M^3 ArcTanh[Lambda/Sqrt[Lambda^2 + M^2]]) -
2*Nf/(2 \[Pi]^2) T*
NIntegrate[((-((3 E^(-((3 (Sqrt[M^2 + p^2] - \[Mu]))/T)) M)/(
Sqrt[M^2 + p^2] T)) - (
6 E^(-((2 (Sqrt[M^2 + p^2] - \[Mu]))/T)) M \[Phi])/(
Sqrt[M^2 + p^2] T) - (
3 E^((-Sqrt[M^2 + p^2] + \[Mu])/T) M \[CapitalPhi])/(
Sqrt[M^2 + p^2] T))/(
1 + E^(-((3 (Sqrt[M^2 + p^2] - \[Mu]))/T)) +
3 E^(-((2 (Sqrt[M^2 + p^2] - \[Mu]))/T)) \[Phi] +
3 E^((-Sqrt[M^2 + p^2] + \[Mu])/T) \[CapitalPhi]) + (-((
3 E^(-((3 (Sqrt[M^2 + p^2] + \[Mu]))/T)) M)/(
Sqrt[M^2 + p^2] T)) - (
3 E^((-Sqrt[M^2 + p^2] - \[Mu])/T) M \[Phi])/(
Sqrt[M^2 + p^2] T) - (
6 E^(-((2 (Sqrt[M^2 + p^2] + \[Mu]))/T)) M \[CapitalPhi])/(
Sqrt[M^2 + p^2] T))/(
1 + E^(-((3 (Sqrt[M^2 + p^2] + \[Mu]))/T)) +
3 E^((-Sqrt[M^2 + p^2] - \[Mu])/T) \[Phi] +
3 E^(-((2 (Sqrt[M^2 + p^2] + \[Mu]))/
T)) \[CapitalPhi])) p^2, {p, 0, Infinity},
AccuracyGoal -> 5, PrecisionGoal -> 2];

\[Delta]\[CapitalOmega]\[Delta]\[CapitalPhi][\[Mu]_, T_,
M_, \[CapitalPhi]_, \[Phi]_] :=
T^4 (-(1/2) \[Phi] B2[
T] + ((-6 \[Phi] - 6 \[Phi]^2 \[CapitalPhi] +
12 \[CapitalPhi]^2) B4[T])/(
1 - 6 \[Phi] \[CapitalPhi] - 3 \[Phi]^2 \[CapitalPhi]^2 +
4 (\[Phi]^3 + \[CapitalPhi]^3))) -
2*Nf/(2 \[Pi]^2) T*
NIntegrate[((3 E^((-Sqrt[M^2 + p^2] + \[Mu])/T))/(
1 + E^(-((3 (Sqrt[M^2 + p^2] - \[Mu]))/T)) +
3 E^(-((2 (Sqrt[M^2 + p^2] - \[Mu]))/T)) \[Phi] +
3 E^((-Sqrt[M^2 + p^2] + \[Mu])/T) \[CapitalPhi]) + (
3 E^(-((2 (Sqrt[M^2 + p^2] + \[Mu]))/T)))/(
1 + E^(-((3 (Sqrt[M^2 + p^2] + \[Mu]))/T)) +
3 E^((-Sqrt[M^2 + p^2] - \[Mu])/T) \[Phi] +
3 E^(-((2 (Sqrt[M^2 + p^2] + \[Mu]))/
T)) \[CapitalPhi])) p^2, {p, 0, Infinity},
AccuracyGoal -> 5, PrecisionGoal -> 2];

\[Delta]\[CapitalOmega]\[Delta]\[Phi][\[Mu]_, T_,
M_, \[CapitalPhi]_, \[Phi]_] :=
T^4 (-(1/2) \[CapitalPhi] B2[
T] + ((12 \[Phi]^2 - 6 \[CapitalPhi] -
6 \[Phi] \[CapitalPhi]^2) B4[T])/(
1 - 6 \[Phi] \[CapitalPhi] - 3 \[Phi]^2 \[CapitalPhi]^2 +
4 (\[Phi]^3 + \[CapitalPhi]^3)))  -
2*Nf/(2 \[Pi]^2) T*
NIntegrate[((3 E^(-((2 (Sqrt[M^2 + p^2] - \[Mu]))/T)))/(
1 + E^(-((3 (Sqrt[M^2 + p^2] - \[Mu]))/T)) +
3 E^(-((2 (Sqrt[M^2 + p^2] - \[Mu]))/T)) \[Phi] +
3 E^((-Sqrt[M^2 + p^2] + \[Mu])/T) \[CapitalPhi]) + (
3 E^((-Sqrt[M^2 + p^2] - \[Mu])/T))/(
1 + E^(-((3 (Sqrt[M^2 + p^2] + \[Mu]))/T)) +
3 E^((-Sqrt[M^2 + p^2] - \[Mu])/T) \[Phi] +
3 E^(-((2 (Sqrt[M^2 + p^2] + \[Mu]))/
T)) \[CapitalPhi])) p^2, {p, 0, Infinity},
AccuracyGoal -> 5, PrecisionGoal -> 2];


To compute M and $$\Phi$$ we organize Module as follows

s[mu_?NumericQ, T_?NumericQ, x0_, y0_, w0_] :=
Module[{sol, rut, x, y, w},
sol = FindRoot[{\[Delta]\[CapitalOmega]\[Delta]M[mu, T, w, x, y] ==
0, \[Delta]\[CapitalOmega]\[Delta]\[CapitalPhi][mu, T, w, x,
y] == 0, \[Delta]\[CapitalOmega]\[Delta]\[Phi][mu, T, w, x,
y] == 0}, {{x, x0}, {y, y0}, {w, w0}},
MaxIterations -> 1000] // Quiet; root = {x, y, w} /. sol; root]


We use Do lope to compute m and $$\Phi$$ in the range $$\mu$$ and T

Do[sols[mu][0] = s[mu, 1, 0.1, 0.1, 300];
Do[sols[mu][i] =
s[mu, 1 + i 5, sols[mu][i - 1][[1]], sols[mu][i - 1][[2]],
sols[mu][i - 1][[3]]];, {i, 1, 62}];, {mu, {0, 200, 270, 340}}]


Visualization

ListPlot[Table[
Table[{1 + i 5, sols[mu][i][[3]]}, {i, 0, 62}], {mu, {0, 200, 270,
340}}], PlotLegends ->
Table[Row[{"\[Mu] = ", mu}], {mu, {0, 200, 270, 340}}],
Frame -> True, FrameLabel -> {"T", "m"}]

ListPlot[Table[
Table[{1 + i 5, sols[mu][i][[1]]}, {i, 0, 62}], {mu, {0, 200, 270,
340}}], PlotLegends ->
Table[Row[{"\[Mu] = ", mu}], {mu, {0, 200, 270, 340}}],
Frame -> True, FrameLabel -> {"T", "\[CapitalPhi]"}]


• Yes, it worked! They are very similar to the original ones!! You solved my problem. Now I need to study your code. I hope this can help many others too. I can't thank you enough. Anyhow, many thanks, Alex!! God bless you! Commented Mar 20, 2022 at 8:59
• @EverlinMartins You are welcome! Commented Mar 20, 2022 at 10:43