Why code is running and not giving any output for NSolve?

Question

I am trying to solve the given equation in the code and want to plot 'N1' with 'P1'. But the code is not giving me any output and still showing running. What is the problem with it? If anyone can figure it out will be appreciated.

dela = -1.5;(*Subscript[\[CapitalDelta], L] ; Laser-Lower polariton \
detuning*)
g1 = 0.0003; (* Subscript[g, 0] ; vacuum optomechanical strength *)
del0 = -0.006; (* \[Delta] ;  exciton-photon detuning*)
ome = 3.014; (* \[CapitalOmega] ; Rabi splitting *)
omem = 0.125; (* Subscript[\[CapitalOmega], m] ; mechanical \
resonator's frequency *)
kexc = 0.002;  (* Subscript[\[Kappa], exc  ;  ]exciton decay rate*)
kk = .2;   (* \[Kappa]; cavity decay rate *)
gma = 0.00001;  (* \[CapitalGamma] ; phonon dcay rate *)
nth = 0;
kal = (kk + kexc)/2 - del0*(kk - kexc)/(2*ome);
g0 = g1*(1 - del0/ome)/2;
kb = 0.1;
eq1 = -I*(P1/Sqrt[2]*(1 - del0/(2*ome)))*(Conjugate[a1] - a1) - 
I*(P1*g1/(2*Sqrt[2]*ome))*(alab - Conjugate[alab] + 
   Conjugate[betab] - betab) + kal*N1 + 
kb*(alab*a1 + Conjugate[alab]*Conjugate[a1] + 
   betab*Conjugate[a1] + 
   Conjugate[betab]*
    a1 + (N1 - 2*(Abs[a1])^2)*(b1 + Conjugate[b1])) == 0;
eq2 = I*dela*alab + I*omem*alab - 
I*g0*((2*N1 - Abs[a1]^2 - Np + 2*2*(Abs[b1])^2 + b1^2)*
    Conjugate[a1] - (Conjugate[betab] + 2*alab)*b1 - 
   alab*Conjugate[b1]) + I*(P1/Sqrt[2]*(1 - del0/(2*ome)))*b1 - 
I*(P1*g1/(2*Sqrt[2]*ome))*(N1 - Np + (Conjugate[a1])^2 - 
   b1^2) + (kal + gma)*alab/2 + 
kb/2*(Conjugate[betab]*b1 + 
   alab*Conjugate[b1] + (Np - 2*(Abs[b1])^2)*Conjugate[a1] + 
   2*alab*b1 - b1^2*Conjugate[a1]) == 0;
eq3 = -I*g0*((N1 - 2*(Abs[a1])^2)*(Conjugate[b1] - b1) + 
   Conjugate[alab]*Conjugate[a1] - alab*a1 + Conjugate[betab]*a1 -
    betab*Conjugate[a1]) - 
I*(P1*g1/(2*Sqrt[2]*ome))*(Conjugate[betab] - alab + 
   Conjugate[alab] - betab) + gma*Np - gma*nth == 0;
eq4 = -I*dela*betab + I*omem*betab - 
I*g0*((2*N1 - (Abs[a1])^2 + Np - 2*(Abs[b1])^2)*
    a1 - (2*betab - a1*b1 + Conjugate[alab])*b1 + 
   betab*Conjugate[b1]) - I*(P1/Sqrt[2]*(1 - del0/(2*ome)))*b1 - 
I*(P1*g1/(2*Sqrt[2]*ome))*(N1 + Np + a1^2 + b1^2 + 1) + (kal + 
   gma)*betab/2 + 
kb/2*((2*betab - a1*b1 + Conjugate[alab])*
    b1 + (Np - 2*(Abs[b1])^2)*a1 + betab*Conjugate[b1]) == 0;
eq5 = -I*dela*a1 - I*g0*(Conjugate[alab] + betab) - 
I*(P1/Sqrt[2]*(1 - del0/(2*ome))) - 
I*(P1*g1/(2*Sqrt[2]*ome))*(b1 + Conjugate[b1]) + kal*a1/2 + 
kb*(Conjugate[alab] + betab)/2 == 0;
eq6 = I*omem*b1 - I*g0*N1 - 
I*(P1*g1/(2*Sqrt[2]*ome))*(a1 + Conjugate[a1]) + gma*b1/2 == 0;
sol[P1_?NumericQ] := FindRoot[{eq1, eq2, eq3, eq4, eq5, eq6}, {N1, alab, Np, betab, a1,
 b1}];
Plot[{Evaluate[N1] /. sol[P1]}, {P1, 0, 20}]

Answer 1

I had thought that the skeleton of the code I showed above in my comment would be enough for you to be able to fill in the details. Apparently that wasn't true. So here is all the code.

dela=-1.5;g1=0.0003;del0=-0.006;ome=3.014;omem=0.125;kexc=0.002;kk=.2;gma=0.00001;nth=0;
kal=(kk+kexc)/2-del0*(kk-kexc)/(2*ome);g0=g1*(1-del0/ome)/2;kb=0.1;

expr1=-I*(P1/Sqrt[2]*(1-del0/(2*ome)))*(Conjugate[a1]-a1)-I*(P1*g1/(2*Sqrt[2]*ome))*(alab-
Conjugate[alab]+Conjugate[betab]-betab)+kal*N1+kb*(alab*a1+Conjugate[alab]*Conjugate[a1]+
betab*Conjugate[a1]+Conjugate[betab]*a1+(N1-2*(Abs[a1])^2)*(b1+Conjugate[b1]));

expr2=I*dela*alab+I*omem*alab-I*g0*((2*N1-Abs[a1]^2-Np+2*2*(Abs[b1])^2+b1^2)*Conjugate[a1]-
(Conjugate[betab]+2*alab)*b1-alab*Conjugate[b1])+I*(P1/Sqrt[2]*(1-del0/(2*ome)))*b1-I*(P1*
g1/(2*Sqrt[2]*ome))*(N1-Np+(Conjugate[a1])^2-b1^2)+(kal+gma)*alab/2+kb/2*(Conjugate[betab]*
b1+alab*Conjugate[b1]+(Np-2*(Abs[b1])^2)*Conjugate[a1]+2*alab*b1-b1^2*Conjugate[a1]);

expr3=-I*g0*((N1-2*(Abs[a1])^2)*(Conjugate[b1]-b1)+Conjugate[alab]*Conjugate[a1]-alab*a1+
Conjugate[betab]*a1-betab*Conjugate[a1])-I*(P1*g1/(2*Sqrt[2]*ome))*(Conjugate[betab]-
alab+Conjugate[alab]-betab)+gma*Np-gma*nth;

expr4=-I*dela*betab+I*omem*betab-I*g0*((2*N1-(Abs[a1])^2+Np-2*(Abs[b1])^2)*a1-(2*betab-a1*
b1+Conjugate[alab])*b1+betab*Conjugate[b1])-I*(P1/Sqrt[2]*(1-del0/(2*ome)))*b1-I*(P1*g1/(2*
Sqrt[2]*ome))*(N1+Np+a1^2+b1^2+1)+(kal+gma)*betab/2+kb/2*((2*betab-a1*b1+Conjugate[alab])*
b1+(Np-2*(Abs[b1])^2)*a1+betab*Conjugate[b1]);

expr5=-I*dela*a1-I*g0*(Conjugate[alab]+betab)-I*(P1/Sqrt[2]*(1-del0/(2*ome)))-I*(P1*g1/(2*
Sqrt[2]*ome))*(b1+Conjugate[b1])+kal*a1/2+kb*(Conjugate[alab]+betab)/2;

expr6=I*omem*b1-I*g0*N1-I*(P1*g1/(2*Sqrt[2]*ome))*(a1+Conjugate[a1])+gma*b1/2;

tbl=Table[{P1,N1}/.NMinimize[Abs[expr1]+Abs[expr2]+Abs[expr3]+Abs[expr4]+Abs[expr5]+Abs[expr6],
{N1,alab,Np,betab,a1,b1}][[2]],{P1,0,20,1}];

ListPlot[tbl,Joined->True]

Now, unless I've made a mistake, I have made no change to any of your six equations other than removing the right hand side so that I can NMinimize the sum of the Abs of each of those.

To be very clear, trying to do NMinimize[Abs[expression]] is not the same as Solve[expression==0]. The first one of those tries to get close to expression==0.

So, as I wrote above, you should also look carefully at

Table[NMinimize[Abs[expr1]+Abs[expr2]+Abs[expr3]+Abs[expr4]+Abs[expr5]+Abs[expr6],
{N1,alab,Np,betab,a1,b1}],{P1,0,20,1}]

and see how well it did. The first item in each of those 21 results will be the minimum of the sum that it was able to find. You might like all those first items to be exactly zero, but they won't be. The question is: Are the sums good enough to give you an idea what your plot looks like.

There is another issue. Plot might do several hundred or even several thousand iterations, finding the value of the N1 at each of those points. This only does 21 points. So it can easily be that there is behavior between any pair of those points which is not displayed in the resulting ListPlot. Changing the step size for P1 might help you get a little better idea what the resulting graph will be like, but without taking hundreds or a thousand times longer to display a result.

There is another issue. You almost certainly intended that alab,Np,betab,a1,b1 could have complex values, not just real values. If you read the documentation for NMinimize and you click on the orange Details and Options you will find: By default, all variables are assumed to be real. Perhaps you can explore changing that default or looking at other functions that might let you find better minima if you search among complex values for alab,Np,betab,a1,b1 and still be able to rapidly get an estimate of what your graph looks like.

All this was ONLY intended to rapidly give you some idea what your graph might look like.

Hopefully this is enough to show you how to use this cheap reasonably fast trick that I intended.

If you can quickly determine whether this is good enough for your purposes then that was the only goal. If it isn't good enough then perhaps letting Solve run for a few hours or days or months might be enough to get you the original Plot that you desired.

Answer 2

Rationalize all of your constants so that you can specify WorkingPrecision later to use arbitrary-precision.

dela = -3/2;(*Subscript[Δ,L];Laser-Lower polariton detuning*)
g1 = 3*^-4;(*Subscript[g,0];vacuum optomechanical strength*)
del0 = -6*^-3;(*δ;exciton-photon detuning*)
ome = 3014*^-3;(*Ω;Rabi splitting*)
omem = 1/8;(*Subscript[Ω,m];mechanical resonator's frequency*)
kexc = 1/500;(*Subscript[κ,exc;]exciton decay rate*)
kk = 1/5;(*κ;cavity decay rate*)
gma = 10^-5;(*Γ;phonon dcay rate*)
nth = 0;
kal = (kk + kexc)/2 - del0*(kk - kexc)/(2*ome);
g0 = g1*(1 - del0/ome)/2;
kb = 1/10;
eq1 = -I*(P1/Sqrt[2]*(1 - del0/(2*ome)))*
     (Conjugate[a1] - a1) -
    I*(P1*g1/(2*Sqrt[2]*ome))*
     (alab - Conjugate[alab] + Conjugate[betab] - betab) +
    kal*N1 + kb*(alab*a1 + Conjugate[alab]*Conjugate[a1] +
       betab*Conjugate[a1] + Conjugate[betab]*a1 +
       (N1 - 2*(Abs[a1])^2)*(b1 + Conjugate[b1])) == 0;
eq2 = I*dela*alab + I*omem*alab - 
    I*g0*((2*N1 - Abs[a1]^2 - Np + 2*2*(Abs[b1])^2 + b1^2)*
        Conjugate[a1] - (Conjugate[betab] + 2*alab)*b1 - alab*Conjugate[b1]) +
     I*(P1/Sqrt[2]*(1 - del0/(2*ome)))*b1 - 
    I*(P1*g1/(2*Sqrt[2]*ome))*(N1 - Np + (Conjugate[a1])^2 - b1^2) + (kal + 
       gma)*alab/2 + 
    kb/2*(Conjugate[betab]*b1 + 
       alab*Conjugate[b1] + (Np - 2*(Abs[b1])^2)*Conjugate[a1] + 2*alab*b1 - 
       b1^2*Conjugate[a1]) == 0;
eq3 = -I*g0*((N1 - 2*(Abs[a1])^2)*(Conjugate[b1] - b1) + 
       Conjugate[alab]*Conjugate[a1] - alab*a1 + Conjugate[betab]*a1 - 
       betab*Conjugate[a1]) - 
    I*(P1*g1/(2*Sqrt[2]*ome))*(Conjugate[betab] - alab + Conjugate[alab] - 
       betab) + gma*Np - gma*nth == 0;
eq4 = -I*dela*betab + I*omem*betab - 
    I*g0*((2*N1 - (Abs[a1])^2 + Np - 2*(Abs[b1])^2)*
        a1 - (2*betab - a1*b1 + Conjugate[alab])*b1 + betab*Conjugate[b1]) - 
    I*(P1/Sqrt[2]*(1 - del0/(2*ome)))*b1 - 
    I*(P1*g1/(2*Sqrt[2]*ome))*(N1 + Np + a1^2 + b1^2 + 1) + (kal + gma)*
     betab/2 + 
    kb/2*((2*betab - a1*b1 + Conjugate[alab])*b1 + (Np - 2*(Abs[b1])^2)*a1 + 
       betab*Conjugate[b1]) == 0;
eq5 = -I*dela*a1 - I*g0*(Conjugate[alab] + betab) - 
    I*(P1/Sqrt[2]*(1 - del0/(2*ome))) - 
    I*(P1*g1/(2*Sqrt[2]*ome))*(b1 + Conjugate[b1]) + kal*a1/2 + 
    kb*(Conjugate[alab] + betab)/2 == 0;
eq6 = I*omem*b1 - I*g0*N1 - I*(P1*g1/(2*Sqrt[2]*ome))*(a1 + Conjugate[a1]) + 
    gma*b1/2 == 0;

eqns = {eq1, eq2, eq3, eq4, eq5, eq6};

Using a modified version of what Bill suggested for NMinimize along with using arbitrary-precision (specified WorkingPrecision).

data = Table[{N1, P1} /. 
    NMinimize[Total[Abs[First[#]]^2 & /@ eqns], {N1, alab, Np, betab, a1, b1},
       WorkingPrecision -> 15][[2]], {P1, 0, 20, 1/2}];

ListLinePlot[data]

Note that in your revised code you used the wrong syntax for FindRoot. See the documentation. You must supply initial estimates for each variable.