A common question has been about the built-in potential
temperature mapping. It is not simple to reverse-engineer
the basis for this mapping, so the explanation is given
here. The original development was by R. V. H. Booth,
however this contained an error, which was corrected
by C. C. McAndrew to give the final form below.

Note that below we will map from temperature Tnom to temperature T.
This appears to be a "1-way" mapping. However physically for any
temperature mapping for any parameter P(T) we should have

P(T1->T2->T1)=P(T1)
P(T1->T3)=P(T1->T2->T3)

so the mappings need to be invertible and transitive. This
is the reason for the final expression, it is meant as a general
mapping from one temperature to another, not as a mapping from
Tnom to other temperatures.

The standard SPICE built-in potential mapping with temperature is:

VJ(T)=VJ(Tnom)*(T/Tnom)-(k*T/q)*ln(exp(-EA*(1-T/Tnom)/(k*T/q))*(T/Tnom)^3) (1a)
     =VJ(Tnom)*(T/Tnom)-(k*T/q)*ln(exp(-(EA*k/q)*(1/T-1/Tnom))*(T/Tnom)^3) (1b)
     =VJ(Tnom)*(T/Tnom)+EA*(1-T/Tnom)-3*(k*T/q)*ln(T/Tnom)                 (1c)

where Tnom is the reference temperature (in K), at which
the initial built-in potential VJ(Tnom) (in V) is known, T is
the temperature (in K) at which the built-in potential VJ(T)
is desired to be known, k is Boltzmann constant (J/K), q is
the magnitude of the electronic charge (Coulomb), and EA is
the activation energy (in V) used to characterize the junction.
In some cases this may be taken as the bandgap voltage, which can
be modeled as a function of temperature also. Here it is taken to
be constant.

The above expression is based on the standard formula that

VJ(T)=(k*T/q)*ln(NA*ND/ni(T)^2) (2)

where NA and ND are the acceptor and donor concentrations on either
side of an assumed abrupt p-n junction, and ni(T) is the temperature
dependent intrinsic concentration. From this expression

VJ(T)=(k*Tnom/q)*(T/Tnom)*ln((NA*ND/ni(Tnom)^2)*(ni(Tnom)^2/ni(T)^2)) (3a)
     =(T/Tnom)*VJ(Tnom)-2*(k*T/q)*ln(ni(T)/ni(Tnom))                  (3b)

The temperature dependence of the intrinsic carrier concentration is

ni(T)=ni(Tnom)*exp(-(0.5*EA*k/q)*(1/T-1/Tnom))*(T/Tnom)^(3/2)

which, substituted into eqn. (3b) gives eqn. (1c).

This model is flawed in that at sufficiently high temperature it
gives a negative value for the built-in potential VJ, which is
physically incorrect. Empirical, piecewise models have been formulated
to avoid this, by asymptotically limiting VJ to zero for high
temperatures. But this is not physical. The error is that ni(T)^2 in
eqn. (2) exceeds NA*ND.

The root cause of the error is that at high temperatures the semiconductor
on either side of the junction becomes intrinsic, and a more correct
expression for the built-in potential is

VJ(T)=(k*T/q)*ln(ppo*nno/ni(T)^2) (4)

where nno is the mobile electron concentration (in equilibrium, with no
bias applied) on the n-type side of the junction and ppo is the
mobile hole concentration (in equilibrium, with no bias applied)
on the p-type side of the junction. At low temperatures ppo is very
nearly NA and nno is very nearly ND, and so eqn. (2) is a good
approximation. At high temperatures, when the doped semiconductor
material becomes intrinsic, then both ppo and nno are very nearly ni,
therefore VJ in eqn. (4) approaches zero, which is physically correct.

The error in the formulation is the use of NA and ND rather than ppo and nno.
This can be fixed in the following manner. For the n-type side of the
junction we have the equilibrium and charge neutrality conditions

n*p=ni^2 (5)
n=p+ND   (6)

and solving this for the equilibrium mobile electron concentration gives

nno=ND*0.5*(1+sqrt(1+4*(ni/ND)^2)) (7)

and a similar analysis for the p-type side gives

ppo=NA*0.5*(1+sqrt(1+4*(ni/NA)^2)) (8)

So making one simplification on the product of the 1+sqrt terms gives

VJ=(k*T/q)*ln((NA*ND/ni^2)*0.25*(1+sqrt(1+4*(ni^2/(NA*ND))))^2) (9)

Use the representation

Pj=(k*T/q)*ln(NA*ND/ni^2) (10)

so that

NA*ND/ni^2=exp(Pj/(k*T/q)) (11)

then substituting this into the sqrt part of eqn. (9) gives

VJ=Pj+2*(k*T/q)*log(0.5*(1+sqrt(1+4*exp(-Pj/(k*T/q))))) (12)

The standard SPICE model gives the temperature mapping for Pj,
as it is based on (k*T/q)*ln(NA*ND/ni^2). Eqn. (12) then gives
the mapping from the temperature dependent calculation of Pj
to the built-in potential VJ. For a SPICE model for the complete
temperature dependence of VJ, i.e. for a model of VJ(T) given
VJ(Tnom), we therefore need to

1. determine Pj from VJ
2. apply the (k*T/q)*ln(NA*ND/ni^2) to Pj (see eqn. (1))
3. apply the mapping eqn. (12) from Pj to VJ at the new temperature

So the final expression needed is for the first of the above steps.
It requires inverting eqn. (9), and this is a bit tricky.
First, let X=NA*ND/ni^2 and E=exp(VJ/(k*T/q)), then eqn. (9)
can be rewritten as

4*E=X*(1+sqrt(1+4/X))^2
   =2*X+4+2*sqrt(X^2+4*X)

=> 2*E-2-X=sqrt(X^2+4*X)
=> 4*E^2+4+X^2-8*E-4*E*X+4*X=X^2+4*X
=> X=E-2+1/E
    =(sqrt(E)-1/(sqrt(E))^2

this gives

Pj=(k*T/q)*ln(X)
  =2*(k*T/q)*ln(exp(0.5*VJ/(k*T/q))-exp(-0.5*VJ/(k*T/q)))

which is the desired transformation from VJ to Pj.

The explicit steps in mapping VJ(Tnom) to VJ(T) are therefore:

Pj(Tnom)=2*(k*Tnom/q)*ln(exp(0.5*VJ(Tnom)/(k*Tnom/q))-exp(-0.5*VJ(Tnom)/(k*Tnom/q)))
Pj(T)   =Pj(Tnom)*(T/Tnom)+EA*(1-T/Tnom)-3*(k*T/q)*ln(T/Tnom)
VJ(T)   =Pj(T)+2*(k*T/q)*log(0.5*(1+sqrt(1+4*exp(-Pj(T)/(k*T/q)))))

which, with a slight difference in notation, is exactly the psibi
function used in VBIC. (The main difference in notation is that the
arguments provided are P (which here is VJ), Vtv=k*T/q, and rT=T/Tnom).

This mapping is physically based, has only one approximation in the
formulation, is completely invertible, transitive, and smooth, and
asymptotically approaches zero as T approaches infinity, which
is physically correct.