How to reduce PolyLog function

Question

I want to solve an equation involving Fermi-Dirac integrals for which I am using PolyLog functions. But I get the error:

NSolve::ifun: Inverse functions are being used by NSolve, so some solutions may not be found; use Reduce for complete solution information.

Here is the code:

ClearAll["Global`*"] (*Clear all variables*)
h := 6.63*10^-34; (* J.s, Planck's constant *)
hbar := (663*10^-34)/(   100*2 π); (* J.s, Reduced Planck's constant *)
q := 16/10*10^-19 ;  (* Coulomb, electron charge *)
m0 := 91/10*10^-31; (* kg free electron rest mass *) 
kb := 8.61*10^-5; (*in eV/K *)
gs = 2;
gv = 1;
mc = 0.2*m0;
T = 300;
nc1d  = gs*gv *((2*π*mc*kb*T)/h^2)^0.5;
n1d = 5*10^6;
v = 1;
NSolve[0.5*nc1d*(-PolyLog[1/2, -Exp[x]] - PolyLog[1/2, -Exp[x - (q*v)/(kb*T)]]) == n1d, x];

Also I want to solve $x$ for different $v$ and plot it, how do I go about it?

Answer 1

Here are your parameters with v=1;

ClearAll["Global`*"] (*Clear all variables*)
h := 6.63*10^-34;(*J.s,Planck's constant*)hbar := (663*10^-34)/(100*2 \
\[Pi]);(*J.s,Reduced Planck's constant*)q := 
 16/10*10^-19;(*Coulomb,electron charge*)m0 := 
 91/10*10^-31;(*kg free electron rest mass*)kb := 
 8.61*10^-5;(*in eV/K*)gs = 2;
gv = 1;
mc = 0.2*m0;
T = 300;
nc1d = gs*gv*((2*\[Pi]*mc*kb*T)/h^2)^0.5;
n1d = 5*10^6;
v = 1;

We can have a look at the equation graphically:

Plot[{0.5*
   nc1d*(-PolyLog[1/2, -Exp[x]] - 
     PolyLog[1/2, -Exp[x - (q*v)/(kb*T)]]), n1d}, {x, -26, -25}]

showing that the solution is close to 25. Now one can solve the equation numerically:

Clear[v];
lst = Table[{v, 
   FindRoot[
     0.5*nc1d*(-PolyLog[1/2, -Exp[x]] - 
         PolyLog[1/2, -Exp[x - (q*v)/(kb*T)]]) == n1d, {x, -25}][[1, 
     2]]}, {v, 1, 10, 0.1}]

this yields a list that can be plotted:

ListPlot[lst]

As you see, it is not visible that anything varies with v. It is since the value of

 q/(kb*T)

(*  6.19435*10^-18  *)

is very small. It is up to you to think, what does it all mean.

A general advice: before starting solving such an equation it is good if you can rescale it and eventually, replace parameters (by combinations of your actual ones) such that you do not operate neither with very large, not with very small numbers.

Have fun!

Answer 2

In addition to Henrik Schumacher's and Alexei Boulbitch's answers, we note the following.

Due to the extreme scale set by $q/k_BT\sim 10^{-18}$, the root of your equation varies extremely slowly with v, and you need $v\sim 10^{18}$ to notice any change in the plot. For reasonable values of $v$, the plot is flat. If you increase the value of v to the physical scale of your problem, you get Underflow[].

We therefore measure all constants in natural units, $c=\hbar=m_0=k_B=1$:

ClearAll["Global`*"]
h = 1/(2 \[Pi]);
hbar = 1;
q = .3;
m0 = 1;
kb = 1;
gs = 2;
gv = 1;
mc = 0.2*m0;
T = 1;
nc1d = gs*gv*((2*\[Pi]*mc*kb*T)/h^2)^0.5;
n1d = 5;
Plot[FindRoot[
     0.5*nc1d*(-PolyLog[1/2, -Exp[x]] - PolyLog[1/2, -Exp[x - (q*v)/(kb*T)]]) == n1d, {x, 1}][[1, 2]]
, {v, 0, 20}]

Note: I've set $T=1$ and $n_{1d}=5$ mostly out of laziness. I leave it to you to transform these dimensionful parameters into the correct numerical values.

Answer 3

In addition to what AccidentalFourierTransform has said, it is a good idea to use exact numbers (numbers without the floating point dot). This way, you can use FindRoot to compute solutions of your equation in (almost) arbitrary precision:

ClearAll["Global`*"] (*Clear all variables*)
h = 663 10^-36;(*J.s,Planck's constant*)
hbar = (663 10^-34)/(100 2 \π);(*J.s,Reduced Planck's constant*)
q = 16/10 10^-19;(*Coulomb,electron charge*)
m0 = 91/10 10^-31;(*kg free electron rest mass*)
kb = 861 10^-7;(*in eV/K*)
gs = 2;
gv = 1;
mc = 2/10 m0;
T = 300;
nc1d = gs gv ((2 π mc kb T)/h^2)^(1/2);
n1d = 5 10^6;
v = 1;
eq[v_?NumericQ] := 1/2 nc1d (-PolyLog[1/2, -Exp[x]] - PolyLog[1/2, -Exp[x - (q v)/(kb T)]]) - n1d;
sol[v_?NumericQ] := FindRoot[eq[v], {x, -25}, AccuracyGoal -> 30, WorkingPrecision -> 60]
eq[v] /. sol[v]

{x -> -25.3646655741586561398370381257338714952592650205634307458808}

0.*10^-39

A plot can be obtained this way:

vlist = Range[1, 3000, 10];
xlist = x /. sol /@ vlist;
ListLinePlot[Transpose[{vlist, xlist}]]

Admittedly, this plot is very unspectacular...