Energia libera

ϕ[m_, T_, 
        q_] := -T/
          2 (m/2/T^2 (1 + (m - 1) q^3) + (m - 1) Log[
                1 - q] +  Log[1 + (m - 1) q])

ϕ[m, T, q]

FullSimplify[D[ϕ[m, T, q], q]]

Plot[ϕ[1.1, 0.62, q], {q, 0, 1}]

Equazione per T*(m) o m*(T)

f[q_, m_] := 
      Sqrt[3/2  q (1 - q) (1 + (m - 1) q)]

Plot[f[q, 0.9], {q, 0, 1}]

Solve[(-3 (m^2) - 
            4 (-2 + m) (1 + m) (-1 + 2 m) T2 + 
            36 (((-1 + m))^2) T2^2) == 0, T2]

Tstar[m_] := 
      Sqrt[((4 (-2 + m) (1 + m) (-1 + 2 m) + 8 (Sqrt[1 - 3 m + 6 (m^2) - 7 (m^3) + 6 (m^4) - 3 (m^5) + (m^6)]))/(72 (((-1 + m))^2)))]

N[Limit[Tstar[m], m -> 1]]

Limit[Tstar[m]/Sqrt[m], m -> Infinity]

Tstar[1.01]

N[Tstar[0]]

Plot[Tstar[m], {m, 0, 10}]

tt := Table[{Tstar[m], m}, {m, 0.0011, 2, 0.001}]

mstar = Interpolation[tt]

mstar[0.472]

Plot[Max[0, mstar[T]], {T, 0, 0.7}]

Equazione per q

FullSimplify[Solve[
          q (1 - q) (1 + (m - 1) q) == 2 T^2/3, 
          q][[1]]]

q1[m_, 
        T_] := (1/(6 (-1 + m))) (2 (-2 + 
                  m) + ((2 (2^(1/3)) (1 + (-1 + m) m))/(((2 - 3 m - 3 (m^2) + 2 (m^3) - 18 (((-1 + m))^2) (T^2) + 3 (Sqrt[(((-1 + m))^2) (-3 (m^2) - 4 (-2 + m) (1 + m) (-1 + 2 m) (T^2) + 36 (((-1 + m))^2) (T^4))])))^(1/3))) + 
            (2^(2/3)) (((2 - 3 m - 3 (m^2) + 2 (m^3) - 18 (((-1 + m))^2) (T^2) + 3 (Sqrt[(((-1 + m))^2) (-3 (m^2) - 4 (-2 + m) (1 + m) (-1 + 2 m) (T^2) + 36 (((-1 + m))^2) (T^4))])))^(1/3)))

Solve[q (1 - q) (1 + (m - 1) q) == 2 T^2/3, 
        q][[3]]

q3[m_, 
        T_] := ((-2 + m)/(3 (-1 + m))) + ((1 - I (Sqrt[3])) (-9 + 9 m - 
                  9 (m^2)))/(9 (2^(2/3)) (-1 + 
                  m) (((54 - 81 m - 81 (m^2) + 54 (m^3) - 486 (T^2) + 972 m (T^2) - 486 (m^2) (T^2) + (Sqrt[4 (((-9 + 9 m - 9 (m^2)))^3) + (((54 - 81 m - 81 (m^2) + 54 (m^3) - 486 (T^2) + 972 m (T^2) - 486 (m^2) (T^2)))^2)])))^(1/3))) - (1/(18 (2^(1/3)) (-1 + 
                    m))) ((1 + I (Sqrt[3])) (((54 - 81 m - 81 (m^2) + 54 (m^3) - 486 (T^2) + 972 m (T^2) - 486 (m^2) (T^2) + (Sqrt[4 (((-9 + 9 m - 9 (m^2)))^3) + (((54 - 81 m - 81 (m^2) + 54 (m^3) - 486 (T^2) + 972 m (T^2) - 486 (m^2) (T^2)))^2)])))^(1/3)))

q1[0.5, 0.1]

q3[1.1, 0.1]

Plot[{q1[1.5, T], q3[1.5, T]}, {T, 0, Tstar[1.5]}, 
      PlotStyle -> {Hue[0], Hue[0.6]}]

q[m_, T_] := If[m > 1, q1[m, T], q3[m, T]]

Plot[q[0.9, T], {T, 0, Tstar[0.9]}]

Termodinamica

F[m_, 
        T_] := (-T/
              2 (1/2/T^2 (1 + (2  m - 1) q^3) + Log[1 - q] + 
                q/(1 + (m - 1) q))) /. 
        q -> q[m, T]

FF[m_, T_] := ϕ[m, T, q[m, T]]/m/T

Σ[m_, T_] := 
      m^2 D[FF[x, T], x] /. x -> m

Plot[Σ[m, 0.4], {m, 0.001, 0.6}]

Plot[F[m, 0.4], {m, 0.001, 0.6}]

Static line

Ts[m_] := 
      T /. FindRoot[Σ[m, T] == 0, {T, 
            0.59 m - 0.4 m (m - 1), 
            0.61 m - 0.4 m (m - 1)}]

Ts[0.01]

ts := Table[{Re[Ts[m]], m}, {m, 0.0075, 1.4, 0.01}]

Export["sta.dat", ts]

ms = Interpolation[ts]

Plot[ms[T], {T, 0.001, 0.67}]

Dynamic line

md[T_] := 
      m /. FindMaximum[
            F[m, T], {m, Max[0.0001, mstar[T] + 0.0001], 
              ms[T]}][[2]]

md[0.6]

Plot[md[T], {T, 0.01, 0.7}]

Final phase diagram

Plot[{ms[T], md[T], Max[0, mstar[T]]}, {T, 0.01, 0.67}, 
      PlotStyle -> {Hue[0], Hue[0.6], Hue[0.4]}]

Export["dias.dat", 
      Table[{T, ms[T]}, {T, 0.01, 0.67, 0.001}]]

Export["diad.dat", 
      Table[{T, md[T]}, {T, 0.01, 0.605, 0.001}]]

Export["diastar.dat", 
      Table[{T, Max[0, mstar[T]]}, {T, 0.01, 0.67, 0.001}]]