xAct meets Asymptotic Safety: the Einstein-Hilbert approximation

Work flow

Preamble

Loading xAct`xTras + Geometry

Block[{Print},<<xAct`xTras`]
$DefInfoQ=False;
$UndefInfoQ=False;

DefManifold[SpaceTime,dim,{μ,ν,α,β,γ,σ,ρ,τ,κ}]
DefMetric[1, metric[-μ,-ν],CD,{";","∇"},PrintAs->"g"]
(** Since we will apply FRG-methods we work will with Euclidean signature **)

(** Defining a second covariant derivative to be used as a 'background derivative'**)
DefCovD[CDbar[-μ],{FromMetric->metric},{":","(OverBar[∇])"}];

Metric Perturbations

DefMetricPerturbation[metric,metricpert,δϵ,PrintAs->"δg"]
ConstantQ[δϵ]^=True;

DefTensor[h[-μ,-ν],SpaceTime,Symmetric[{-μ,-ν}],PrintAs->"h"];
metricpert[LI[n_],μ_,ν_]^:=Which[n===1,h[μ,ν],n>1,0]

From EH-approximation to 2-point functions

Einstein-Hilbert approximation

(** Defining constant symbols to be used as Newton's constant and comological constant **)
DefConstantSymbol[{G,Λ}]

(** Defining the EH action as our truncation for (Subscript[Γ, k]) **)
ΓEH=(1/(16π G))(Sqrt[Detmetric[]])(2Λ - RicciScalarCD[]);

(** Defining constant symbols to be used as gauge-fixing parameters **)
DefConstantSymbol[αgauge,PrintAs->"(Subscript[α, gf])"]
DefConstantSymbol[βgauge,PrintAs->"(Subscript[β, gf])"]

(** Gauge-fixing part of our truncation **)
(**The gauge-fixing part is already in terms of metric fluctuation, so there is not need to expand this part in terms of metric perturbations**)
F[μ_,α_,β_]:=(1/(Sqrt[16π G αgauge]))((1/2)(metric[μ,α]CDbar[β]@#+metric[μ,β]CDbar[α]@#)-((1+βgauge)/dim)metric[α,β]CDbar[μ]@#)&;
ℱ[μ_]:=Module[{α=DummyIn[TangentSpaceTime],β=DummyIn[TangentSpaceTime]},F[μ,α,β]@h[-α,-β]];
Γgf=(1/2)ℱ[-μ]ℱ[μ]//CollectTensors;

(** Fadeev-Popov ghost part of our truncation **)
DefTensor[Gh[-μ],SpaceTime,PrintAs->"c"] (** Defining ghost field **)
DefTensor[GhBar[-μ],SpaceTime,PrintAs->"(OverBar[c])"] (** Defining anti-ghost field **)

Γghost=GhBar[μ]F[-μ,α,β][CDbar[-α]@Gh[-β]+CDbar[-β]@Gh[-α]]//CollectTensors;
 (** Since we are interested in 2-point functions at vanishing fluctuation fields, we are already setting 𝒪(h) to zero **)
(** Warning: FP-ghosts are anti-commuting fields! For the calculation discussed here, the anti-commuting nature of the ghosts fields does not play a major role. But we need to be careful if we extend this beyond bi-linear contributions of anti-commuting fields **)

Quadratic expansion of (Subscript[Γ, k]) and York decomposition

Collecting quadratic terms in the metric fluctuation

(** Computing (((Subscript[h, μν])[(Γ^((2)))])^(μν,αβ))(Subscript[h, αβ]) -- (note: Perturbation[_,2] already takes care of the factor 1/2)**)
hΓ2h=CollectTensors[((ExpandPerturbation@Perturbation[ΓEH,2])/(Sqrt[Detmetric[]]))+2 × Γgf]/.{CD->CDbar,RicciScalarCD->RicciScalarCDbar,RicciCD->RicciCDbar,RiemannCD->RiemannCDbar}
(** We add the gauge-fixing term separately since it is already quadratic in (Subscript[h, μν]) **)

York decomposition

(** Defining the different modes of York decomposition: (Subscript[h, μν])=(Subscript[h, TT,μν]) + (Subscript[∇, μ])(Subscript[ξ, ν]) + (Subscript[∇, μ])(Subscript[ξ, ν]) + (Subscript[∇, μ])(Subscript[∇, ν])σ - (1/d)(Subscript[g, μν])(∇^2)σ + (1/d)(Subscript[g, μν])h **)
DefTensor[hTT[-μ,-ν],SpaceTime,Symmetric[{-μ,-ν}],PrintAs->"(Subscript[h, TT])"]
DefTensor[ξfield[-μ],SpaceTime,PrintAs->"ξ"]
DefTensor[σfield[],SpaceTime,PrintAs->"σ"]
DefTensor[htr[],SpaceTime,PrintAs->"(Subscript[h, tr])"]

(** Defining properties of (Subscript[h, TT]) **)
hTT/:hTT[μ_,ν_]/;(μ+ν===0):=0 (** traceless condition **)
hTT/:CDbar[α_]@hTT[μ_,ν_]/;((α+μ===0)||(α+ν===0)):=0 (** transversality condition **)

(** Defining properties of ξ **)
ξfield/:CDbar[α_]@ξfield[μ_]/;(α+μ===0):=0 (** transversality condition **)

(** Defining a rule to implement the York decomposition **)
RuleYorkDecomp={h[μ_,ν_]:>hTT[μ,ν]+CDbar[μ]@ξfield[ν]+CDbar[ν]@ξfield[μ]+CDbar[μ]@CDbar[ν]@σfield[]-(1/dim)metric[μ,ν]Module[{α=DummyIn[TangentSpaceTime]},CDbar[α]@CDbar[-α]@σfield[]]+(1/dim)metric[μ,ν]htr[]};

(** Computing (((Subscript[h, μν])[(Γ^((2)))])^(μν,αβ))(Subscript[h, αβ]) with York decomposition **)
hΓ2hYorkStep1=hΓ2h/.RuleYorkDecomp//CollectTensors

(** It is useful to introduce constant parameters to keep track of the different 'York sectors' **)
DefConstantSymbol[{ϵTT,ϵξ,ϵσ,ϵtr}];
(** Redefining each York mode according to (Subscript[h, i])->(Subscript[ϵ, i])(Subscript[h, i]) **)
hΓ2hYorkStep2=hΓ2hYorkStep1/.{hTT[μ_,ν_]:>ϵTT hTT[μ,ν],ξfield[μ_]:>ϵξ ξfield[μ],σfield[]:>ϵσ σfield[],htr[]:>ϵtr htr[]}

(** Constructing a Matrix with the different York sectors **)
hΓ2hMatrix=Table[If[x===y,1,1/2]Coefficient[hΓ2hYorkStep2,x y],{x,{ϵTT,ϵξ,ϵσ,ϵtr}},{y,{ϵTT,ϵξ,ϵσ,ϵtr}}]
(** We introduce 'If[x===y,1,1/2]' to avoid the double counting of mixing-terms**)

Bringing terms to a “standard form”

We need to bring all terms of “hΓ2hMatrix” to a standard form given by “field1-Operator-field2”, while eliminating curvature-independent off-diagonal contributions (i.e., we want to convert (Subscript[∇, μ])(Subscript[∇, ν]) into (Subscript[g, μν])(∇^2)).
To accomplish this goal, the general strategy is to perform a sequence of integration by parts combined with sorting of covariant derivatives.

(** First step: define the simplification that we can use at second order in the derivative expansion **) 
DefConstantSymbol[Rbarconst,PrintAs->"(Subscript[(OverBar[R]), 0])"];
EHTruncate[expr_]:=(expr//RiemannToWeyl//RicciToTFRicci)/.WeylCDbar->Zero/.TFRicciCDbar->Zero/.RicciScalarCDbar[]->Rbarconst/.(Rbarconst^(n_/;n>1)):>0//CollectTensors;

(** Let's focus on TT-TT contributions **) 
hΓ2hMatrix[[1,1]]//EHTruncate
(** In this case, field1 = hTT and field2 = hTT **)

(** We can integrate by parts by taking a functional derivative w.r.t 'field2' and multiplying the result by the 'field2' -- the result has the form field2-operator-field1 **)
1/2×hTT[μ,ν]VarD[hTT[μ,ν],CDbar][hΓ2hMatrix[[1,1]]]//CollectTensors (** 1/2 is a kind of 'symmetry factor' **)
%//SortCovDsToDiv[hTT]//SortCovDsToBox[hTT]//CollectTensors (** We sort the CDs to produce terms that will cancell due to properties of hTT **)
ToStdFormStep1=%//EHTruncate

(** The previous step account for simplifications where derivatives acting on 'field1'**)
(** We need to repeate this procedure to bring the result to the form field1-operator-field2, and then implement simplifications to to derivatives acting on field2 **)
1/2×hTT[α,β]VarD[hTT[α,β],CDbar][ToStdFormStep1]//CollectTensors
%//SortCovDsToDiv[hTT]//SortCovDsToBox[hTT]//CollectTensors (** We sort the CDs to produce terms that will cancell due to properties of hTT **)
ToStdFormStep2=%//EHTruncate

(** Let's turn this into a function that we can apply on any component of hΓ2hMatrix **)
ToStdForm[field1_[inds1___],field2_[inds2___]][expr_]:=
Module[{dummy},
dummy=expr//EHTruncate//SortCovDsToDiv[field1]//SortCovDsToDiv[field2]//CollectTensors//EHTruncate;
(************************************************************************************************)
dummy=If[field1===field2,1/2,1]×field1[inds1]VarD[field1[inds1],CDbar][dummy]//CollectTensors;
dummy=FixedPoint[(#//SortCovDsToDiv[field2]//SortCovDsToBox[field2]//CollectTensors //EHTruncate)&,dummy];
(************************************************************************************************)
dummy=If[field1===field2,1/2,1]×field2[inds2]VarD[field2[inds2],CDbar][dummy]//CollectTensors ;
dummy=FixedPoint[(#//SortCovDsToDiv[field1]//SortCovDsToBox[field1]//CollectTensors//EHTruncate)&,dummy];
(************************************************************************************************)
 dummy
]

(** Checking the result in the TT-TT sector **)
hΓ2hMatrix[[1,1]]//ToStdForm[hTT[μ,ν],hTT[α,β]];
ToStdFormStep2-%//CollectTensors

(** Constructing a matrix hΓ2h with all terms in the standard form **)
Monitor[hΓ2hMatrixStdForm=Table[Module[{list1={hTT[μ,ν],ξfield[μ],σfield[],htr[]},list2={hTT[α,β],ξfield[α],σfield[],htr[]}},hΓ2hMatrix[[i,j]]//ToStdForm[list1[[i]],list2[[j]]]],{j,1,4},{i,1,4}];,Row[{i+4(j-1),"/",16}]]

hΓ2hMatrixStdForm

Faddev-Popov ghost sector

(** T+L decomposition of ghosts **)
DefTensor[GhT[-μ],SpaceTime,PrintAs->"(Subscript[c, T])"]
DefTensor[GhL[],SpaceTime,PrintAs->"(Subscript[c, L])"]
DefTensor[GhBarT[-μ],SpaceTime,PrintAs->"(Subscript[(OverBar[c]), T])"]
DefTensor[GhBarL[],SpaceTime,PrintAs->"(Subscript[(OverBar[c]), L])"]

(** Defining properties of (Subscript[c, T]) **)
GhT/:CDbar[α_]@GhT[μ_]/;(α+μ===0):=0 (** transversality condition **)
GhBarT/:CDbar[α_]@GhBarT[μ_]/;(α+μ===0):=0 (** transversality condition **)

(** Performing T+L decomposition and bringing the result into matrix form **)
DefConstantSymbol[{ϵT,ϵL}]
Γghost/.{Gh[μ_]:>ϵT GhT[μ]+ϵL CDbar[μ]@GhL[],GhBar[μ_]:>ϵT GhBarT[μ]+ϵL CDbar[μ]@GhBarL[]}//CollectTensors;
GhBarΓ2GhMatrix=Table[Coefficient[%,x y],{x,{ϵT,ϵL}},{y,{ϵT,ϵL}}];

(** Constructing a matrix (OverBar[c])Γ2c with all terms in the standard form **)
GhBarΓ2GhMatrixStdForm=Table[Module[{list1={GhT[μ],GhL[]},list2={GhBarT[μ],GhBarL[]}},GhBarΓ2GhMatrix[[i,j]]//ToStdForm[list1[[i]],list2[[j]]]],{j,1,2},{i,1,2}];

GhBarΓ2GhMatrixStdForm

From 2-point functions to the final result

Field redefinition and regularisation

Field redefinition

We have used the TT decomposition in both graviton and ghost sector - this gives rise to a Jacobian. At the same time, our fields now have higher order kinetic terms. Via a suitable field redefinition, we get another Jacobian that solves both problems: it exactly cancels the original Jacobian, and gives all terms a standard kinetic term.

The Hessian is now a function of the Laplacian and the constant Ricci scalar - they thus commute and we can treat them as constant symbols (we already have done so with R).

DefConstantSymbol[Laplacianbar,PrintAs->"(OverBar[Δ])"];
ToLaplacian=(#//.CDbar[μ_]@CDbar[ν_]@X__/;(μ===-ν):>-Laplacianbar X)&;

hΓ2hMatrixStdFormSimplified1=Collect[hΓ2hMatrixStdForm//ToLaplacian//CollectTensors[#,SimplifyMethod->FullSimplify]&,Rbarconst]//EHTruncate;
%//MatrixForm

hΓ2hMatrixStdFormSimplified1/.ξfield[μ_]:>(1/(Sqrt[Laplacianbar-(1/dim)Rbarconst]))ξfield[μ]/.σfield[]:>(1/(Sqrt[(Laplacianbar^2)-(1/(dim-1))Laplacianbar Rbarconst]))σfield[];
hΓ2hMatrixStdFormSimplified2=Series[%,{Rbarconst,0,1}]//Normal//CollectTensors[#,SimplifyMethod->(FullSimplify[#,Assumptions->{Laplacianbar>0}]&)]&;
%//MatrixForm

Now all entries in the Hessian are linear in the Laplacian! Let us do the same in the ghost sector.

GhBarΓ2GhMatrixStdForm//ToLaplacian//CollectTensors
GhBarΓ2GhMatrixStdForm1=%/.GhL[]:>(1/(Sqrt[Laplacianbar]))GhL[]/.GhBarL[]:>(1/(Sqrt[Laplacianbar]))GhBarL[]

Actual Hessian

To compute the Hessian, we have to take the second functional derivative again.
Here, there is a subtlety: the functional derivative of a constrained field w.r.t. itself is a projector! We take this into account by hand, using the corresponding heat kernel coefficients and just the scalar part of the Hessian in all sectors.

TTHessian=2hΓ2hMatrixStdFormSimplified2[[1,1]]//RemoveTensors;
ξHessian=2hΓ2hMatrixStdFormSimplified2[[2,2]]//RemoveTensors;
scalarHessian={{2,1},{1,2}}hΓ2hMatrixStdFormSimplified2[[3;;4,3;;4]]//RemoveTensors;

In a careful treatment, we have to take into account the Grassmann nature of the ghosts. We shall play it dirty here, as we can get away with it at this level:

GhTHessian=2GhBarΓ2GhMatrixStdForm1[[1,1]]//RemoveTensors;
GhLHessian=2GhBarΓ2GhMatrixStdForm1[[2,2]]//RemoveTensors;

Regulator

As a final ingredient, we need to specify the regulator. For simplicity, we simply replace the Laplacian by a regularised Laplacian.

DefScalarFunction[Rh];
DefScalarFunction[Rhdot];
DefScalarFunction[Rgh];
DefScalarFunction[Rghdot];

TTRegulator=D[TTHessian,Laplacianbar]Rh[Laplacianbar];
ξRegulator=D[ξHessian,Laplacianbar]Rh[Laplacianbar];
scalarRegulator=D[scalarHessian,Laplacianbar]Rh[Laplacianbar];

GhTRegulator=D[GhTHessian,Laplacianbar]Rgh[Laplacianbar];
GhLRegulator=D[GhLHessian,Laplacianbar]Rgh[Laplacianbar];

To get the RG derivative of the regulator, note that by our construction, the mass dimension of the shape function is 2, and the argument is dimensionless. This implies the following rule:

(*this refers to the dimensionful beta function!*)
DefConstantSymbol[Gdot,PrintAs->"(OverDot[G])"];
DefConstantSymbol[αgaugedot,PrintAs->"(Subscript[(OverDot[α]), gf])"];
DefConstantSymbol[βgaugedot,PrintAs->"(Subscript[(OverDot[β]), gf])"];

ComputeRdot=Simplify[Grad[#,{G,αgauge,βgauge,Rh[Laplacianbar],Rgh[Laplacianbar]}].{Gdot,αgaugedot,βgaugedot,2Rh[Laplacianbar]-2Laplacianbar Rh'[Laplacianbar],2Rgh[Laplacianbar]-2Laplacianbar Rgh'[Laplacianbar]}]&;

TTRegulatordot=ComputeRdot[TTRegulator];
ξRegulatordot=ComputeRdot[ξRegulator];
scalarRegulatordot=ComputeRdot[scalarRegulator];

GhTRegulatordot=ComputeRdot[GhTRegulator];
GhLRegulatordot=ComputeRdot[GhLRegulator];

Inversion and trace

We can finally combine everything and compute the trace, sector by sector. HK coefficients can be found, e.g., in the Universal RG Machine. Let us specify d=4 for convenience.
We use some tricks here. First, as noted before, there are projectors for all tensor fields. This is handled by using appropriate HK coefficients - in a more general case, this is more complicated. Second, we do not explicitly perform the inverse Laplace transform, and introduce a new symbol for the integration variable in the threshold integral - rather, we use Δ to serve a dual role. Thus, all final expressions still have to be integrated over Δ from 0 to ∞.

DefConstantSymbol[etaG,PrintAs->"(Subscript[η, G])"];
GdotToEta={Gdot->etaG G};

DefConstantSymbol[HKs,PrintAs->"𝔰"];

(*TT sector*)
Simplify[(1/(TTHessian+TTRegulator))TTRegulatordot/.GdotToEta];(*G·Rdot*)
(1/2)%(1/(((4π HKs))^2))(5-(5/6)HKs Rbarconst)/.dim->4//Simplify;(*traced HK coefficients*)
(Table[(1/(Gamma[n]))(Laplacianbar^(n-1)),{n,1,Length[#]}].#)&[CoefficientList[%,(1/HKs)][[2;;]]];(*Inverse Laplace transformation*)
TTtrace=Simplify/@(Series[%,{Rbarconst,0,1}]//Normal)(*Expand to linear order in R*)

(*ξ sector*)
Simplify[(1/(ξHessian+ξRegulator))ξRegulatordot/.GdotToEta];(*G·Rdot*)
(1/2)%(1/(((4π HKs))^2))(3+(1/4)HKs Rbarconst)/.dim->4//Simplify;(*traced HK coefficients*)
(Table[(1/(Gamma[n]))(Laplacianbar^(n-1)),{n,1,Length[#]}].#)&[CoefficientList[%,(1/HKs)][[2;;]]];(*Inverse Laplace transformation*)
ξtrace=Simplify/@(Series[%,{Rbarconst,0,1}]//Normal)(*Expand to linear order in R*)

(*scalar sector - this is a 2x2 matrix!*)
Simplify[Tr[Inverse[scalarHessian+scalarRegulator].scalarRegulatordot]/.GdotToEta];
(1/2)%(1/(((4π HKs))^2))(1+(1/6)HKs Rbarconst)/.dim->4//Simplify;(*traced HK coefficients*)
(Table[(1/(Gamma[n]))(Laplacianbar^(n-1)),{n,1,Length[#]}].#)&[CoefficientList[%,(1/HKs)][[2;;]]];(*Inverse Laplace transformation*)
scalartrace=Simplify/@(Series[%,{Rbarconst,0,1}]//Normal)(*Expand to linear order in R*)

(*ghT sector*)
Simplify[(1/(GhTHessian+GhTRegulator))GhTRegulatordot/.GdotToEta];(*G·Rdot*)
-2(1/2)%(1/(((4π HKs))^2))(3+(1/4)HKs Rbarconst)/.dim->4//Simplify;(*traced HK coefficients*)
(Table[(1/(Gamma[n]))(Laplacianbar^(n-1)),{n,1,Length[#]}].#)&[CoefficientList[%,(1/HKs)][[2;;]]];(*Inverse Laplace transformation*)
GhTtrace=Simplify/@(Series[%,{Rbarconst,0,1}]//Normal)(*Expand to linear order in R*)

(*ghL sector*)
Simplify[(1/(GhLHessian+GhLRegulator))GhLRegulatordot/.GdotToEta];(*G·Rdot*)
-2(1/2)%(1/(((4π HKs))^2))(1+(1/6)HKs Rbarconst)/.dim->4//Simplify;(*traced HK coefficients*)
(Table[(1/(Gamma[n]))(Laplacianbar^(n-1)),{n,1,Length[#]}].#)&[CoefficientList[%,(1/HKs)][[2;;]]];(*Inverse Laplace transformation*)
GhLtrace=Simplify/@(Series[%,{Rbarconst,0,1}]//Normal)(*Expand to linear order in R*)

(*summing everything up*)
totaltrace=TTtrace+ξtrace+scalartrace+GhTtrace+GhLtrace;

Phase diagram example

To illustrate our beta functions, we will pick a specific gauge and regulator now and compute the phase diagram, the fixed point and its critical exponents.
We use the fact that in the Landau limit (α->0), we have (OverDot[α])=-η α+(Subscript[η, α]) (α^2) and (OverDot[β])=(Subscript[η, β]) α (see 1707.01397). Notably, in our flow both (Subscript[η, α]) and (Subscript[η, β]) drop out of the equations, so we do not have to compute them. We will furthermore use the Litim regulator.

DefConstantSymbol[ηα];
DefConstantSymbol[ηβ];
smallαscaling={αgaugedot->-etaG αgauge+ηα ((αgauge)^2),βgaugedot->ηβ αgauge};

DefConstantSymbol[Λdot,PrintAs->"(OverDot[Λ])"];
(*dimensionless couplings*)
DefConstantSymbol[g];
DefConstantSymbol[λ];
DefConstantSymbol[gdot,PrintAs->"(OverDot[g])"];
DefConstantSymbol[λdot,PrintAs->"(OverDot[λ])"];

todimensionless={Gdot->gdot-2g,etaG->((gdot-2g)/g),G->g,Λdot->λdot+2λ,Λ->λ};

(*LHS*)
((ΓEH)/(Sqrt[Detmetric[]]))/.RicciScalarCD[]->Rbarconst;
LHSd4=Collect[Grad[%,{G,Λ}].{Gdot,Λdot}/.todimensionless,Rbarconst,Simplify]

(*RHS*)
Normal@Series[totaltrace/.smallαscaling,{αgauge,0,0}];
%/.Rh->(1-#&)/.Rgh->(1-#&);
RHSd4=Collect[Integrate[%,{Laplacianbar,0,1}]/.todimensionless,Rbarconst,Simplify/@Apart[#,λ]&]

(*solve for β functions*)
βfuncs=Solve[Thread[CoefficientList[LHSd4-RHSd4,Rbarconst]==0],{gdot,λdot}][[1]]//Simplify

(*find fixed points for general β*)
{gdot,λdot}/.βfuncs;
FPs=Solve[Thread[%==0],{g,λ}]//Simplify;
stabmat=-Grad[%%,{g,λ}]//Simplify;

(*there are four fixed points in total*)
Length[FPs]

(*the first is the Gaussian fixed point*)
{g,λ}/.FPs[[1]]
Eigenvalues[stabmat/.FPs[[1]]]

(*the second is the real fixed point with a pair of complex conjugate critical exponents*)
{g,λ}/.FPs[[2]];
Plot[%,{βgauge,-3,2},AxesLabel->{"β",""},PlotLegends->{"(SubStar[g])","(SubStar[λ])"}]
Eigenvalues[stabmat]/.FPs[[2]];
Plot[{Re[%[[1]]],Im[%[[1]]],Im[%[[2]]]},{βgauge,-3,2},AxesLabel->{"β",""},PlotLegends->{"Re θ","Im (Subscript[θ, 1])","Im (Subscript[θ, 2])"}]

(*the other two fixed points are complex*)
{g,λ}/.FPs[[3]];
Plot[Im[%[[1]]],{βgauge,-3,2},AxesLabel->{"β",""},PlotLegends->"Im (SubStar[g])"]
{g,λ}/.FPs[[4]];
Plot[Im[%[[1]]],{βgauge,-3,2},AxesLabel->{"β",""},PlotLegends->"Im (SubStar[g])"]

Let us show the phase diagram for β=1:

(*compute separatrix*)
{λdot,gdot}/.βfuncs/.βgauge->1/.g->gtraj[t]/.λ->λtraj[t];
separatrix=NDSolve[Join[Thread[{λtraj'[t],gtraj'[t]}==%],{λtraj[0]==((10^(-5))/(8π)),gtraj[0]==(10^(-5))}],{λtraj,gtraj},{t,-3,10}][[1]];

(*phase diagram for β=1*)
{λdot,gdot}/.βfuncs/.βgauge->1;
Show[StreamPlot[%,{λ,-3/4,1/2},{g,0,3/2},PlotRange->{{-0.78,0.52},{-0.03,1.52}},RegionFunction->Function[{λ,g},72 π (((1-2 λ))^2)+5 g (-7+6 λ)>0],RegionFillingStyle->None,RegionBoundaryStyle->None,FrameLabel->{"λ","g"},RotateLabel->False,Epilog->Inset["β=1",{0.4,1.3}]],ParametricPlot[{λtraj[t],gtraj[t]}/.separatrix,{t,-3,10},PlotStyle->Darker@Green],Graphics[{Red,Disk[{0,0},0.015]}],Graphics[{Purple,Disk[{λ,g}/.FPs[[2]]/.βgauge->1,0.015]}],ContourPlot[72 π (((1-2 λ))^2)+5 g (-7+6 λ)==0,{λ,-3/4,1/2},{g,0,3/2},ContourStyle->{Dashed,Red},PlotPoints->50]]