Perovskite mobile ion solver

The mobile ion solver is implemented after the work of Calado

\[\label{eq:pdrive} \boldsymbol{J_a} = q \mu_a a_{f} {\nabla E_{v}} - q D_a {\nabla a_{f}}.\]

\[\label{eq:contp} \nabla \boldsymbol{J_a} = - q \frac{\partial a}{\partial t}.\]

Average free carrier mobility

In this model there are two types of electrons (holes), free electrons (holes) and trapped electrons (holes). Free electrons (holes) have a finite mobility of \(\mu_e^0\) (\(\mu_h^0\)) and trapped electrons (holes) can not move at all and have a mobility of zero. To calculate the average mobility we take the ratio of free to trapped carriers and multiply it by the free carrier mobility.:

\[\mu_e(n)=\frac{\mu_e^0 n_{free}}{n_{free}+n_{trap}}\]

Thus if all carriers were free, the average mobility would be \(\mu_e^0\) and if all carriers were trapped the average mobility would be 0. It should be noted that only \(\mu_e^0\) (\(\mu_h^0\)) are used in the model for computation and \(\mu_e(n)\) is an output parameter.

The value of \(\mu_e^0\) (\(\mu_h^0\)) is an input parameter to the model. This can be edited in the electrical parameter editor. The value of \(\mu_e(n)\), and \(\mu_h(p)\) are output parameters from the model. The value of \(\mu_e(n)\), and \(\mu_h(p)\) change as a function of position, within the device, as the number of both free and trapped charge carriers change as a function of position. The values of \(\mu_e(x)\), and \(\mu_h(x)\) can be found in \(mu\_n\_ft.dat\) and \(mu\_p\_ft.dat\) within the \(snapshots\) directory. The spatially averaged value of mobility, as a function of time or voltage can be found in the files \(dynamic\_mue.dat\) or \(dynamic\_muh.dat\) within the dynamic directory.

Were one to try to measure mobility using a technique such as CELIV or ToF, one would expect to get a value closer to \(\mu_e(n)\) or \(\mu_h(p)\) rather than closer to \(\mu_e^0\) or \(\mu_h^0\). It should be noted however, that measuring mobility in disordered materials is a difficult thing to do, and one will get a different experimental value of mobility depending upon which experimental measurement method one uses, furthermore, mobility will change depending upon the charge density profile within the device, and thus upon the applied voltage and light intensity. To better understand this, try for example doing a CELIV simulation, and plotting \(\mu_e(n)\) as a function of time (Voltage). You will see that mobility reduces as the negative voltage ramp is applied, this is because carriers are being sucked out of the device. Then try extracting the mobility from the transient using the CELIV equation for extracting mobility. Firstly, the CELIV equation will give you one value of mobility, which is a simplification of reality as the value really changes during the application of the voltage ramp. Secondly, the value you get from the equation will almost certainly not match either \(\mu_e^0\) or any value of \(\mu_e(n)\). This simply highlights, the difficult of measuring \(a\) value of mobility for a disordered semiconductor and that really when we quote a value of mobility for a disordered material, it really only makes sense to quote a value measured under the conditions a material will be used. For example, for a solar cell, values of \(\mu_e(n)\) and \(\mu_h(n)\), would be most useful to know under 1 Sun at the \(P_{max}\) point on a JV curve.

Configuring the electrical solver

Behind OghmaNano are a series of non-linear solvers that solve the electrical equations in a highly efficient way. These can be configured by going to the electrical tab. There you will see the Drift diffusion button, to the left of that is an arrow. If you click on this it will bring up a window which allows you to configure the "Newton solver". The options are described below.

Related YouTube videos:

How to optimize simulations in OghmaNano so they run faster

• Max Electrical iterations (first step): The maximum number of steps the solver can after it’s cold started onto a new problem. This is usually at 0V in the dark. The solver usually takes more steps on it’s first go.

• Electrical clamp (first step): This is a number by which the maximum newton step is clamped to. 0.1 will make the solver very stable but very slow, 4.0 will make the solver very fast but unstable. A recommended value of 1.0 is suggested for normal problems. If you are solving for high doping or other unusual conditions it can be worth reducing the step. Likewise if you want the solver to be fast and you know the problem is easy set the value to 2.0 or higher. For the first step, I would consider setting this value to be slightly lower than for the subsequent steps.

• Desired solver error (first step): This is the desired error, smaller is more accurate and slower. I would generally not accept answers above \(1x10^{-5}\)

• Max Electrical iterations: Maximum number of electrical iterations on all but the first step.

• Electrical clamp: Electrical clamp (first step): This is a number by which the maximum newton step is clamped to. 0.1 will make the solver very stable but very slow, 4.0 will make the solver very fast but unstable. A recommended value of 1.0 is suggested for normal problems. If you are solving for high doping or other unusual conditions it can be worth reducing the step. Likewise if you want the solver to be fast and you know the problem is easy set the value to 2.0 or higher.

• Desired solver error: This is the desired error, smaller is more accurate and slower. I would generally not accept answers above \(1x10^{-5}\)

• Newton solver clever exit: If the solver starts bouncing in the noise then assume we can’t get a better answer and quit.

• Newton minimum iterations: Don’t allow the solver to quit before doing this number of steps. Often the error in the first few steps of the solution can be below "Desired solver error", thus the solver can quit before finding the true answer.

• Solve Kirchhoff’s current law in Newton solver: Solve Kirchhoff’s current law in the main Newton Jacobian.

• Matrix solver: This selects the matrix solver to use.

• Newton solver to use:

  • none: No electrical solver is selected, this is used when only solving optical or thermal problems.

  • newton: The standard 1D Newton solver.

  • newton2D: The standard 2D Newton solver.

  • newtonnorm: The standard 1D Newton solver but with Slotboom normalization. This is handy when solving systems with large difference in density between minority and majority carrier density.

  • poisson2d: A 2D Poisson solver with no drift diffusion equations.

• Complex matrix solver:

• Slotboom T0: Slotboom variable for the newtonnorm solver.

• Slotboom D0: Slotboom variable for the newtonnorm solver.

• Slotboom n0: Slotboom variable for the newtonnorm solver.

• Use newton cache (experimental): Cache large problems to disk - experimental.

• Quit on convergence problem: Quit on convergence problem. Quite often

• Quit on inverted Fermi-level:

• Solver output verbosity:

Solver stability

Avoiding very big and very small numbers

Try opening up MATLAB (Octave if you are on Linux) and typing in the following equation \(((1e-1+1e1)-1e1)/1e-1\). Before pressing enter, try to evaluate it in your head. the \(1e1\) and the \(-1e1\) cancel leaving \(\frac{1e-1}{1e-1}\) which equates to \(1\). Now try replacing the powers to 1 with to the 19, so type in \(((1e-19+1e19)-1e19)/1e-19\), again evaluate this in your head. Again , \(1e19\) and the \(-1e19\) cancel leaving \(\frac{1e-19}{1e-19}\) which equates to \(1\) Now let the computer evaluate the expression. In fact this time the computer does not give you \(1\) but gives you \(0\). Double check that you typed it in correctly... you did so what is happening. Why is the computer giving me an answer which is 100% wrong. The answer is easy, computers have a limited precision. This means that they can only store a limited number of decimal places. On a modern PC it’s about 15 decimal places. After this the computer starts ignoring the numbers. So when we added \((1e-19+1e19)\) the computer could not keep track of the decimal places so it assumed that the answer was exactly \(1.000000000000000e19\) and not \(1.0000000000000000001e19\), then when we subtracted \(-1e19\) from the answer the computer gave us zero instead of \(1e-19\). The \(1e-19\) was lost in the precision.

All computers are affected by this no matter how powerful they are, this has important implications when solving device equations. If you have too big a spread of numbers in your simulation (matrix/Jacobian) the computer won’t be able to solve it easily. So if you have very low values of mobility say \(1e-19\) and very big values say \(1e5\) the computer wills start to have problems solving the electrical problem. There fore generally try to reduce the spread of parameters in you model. This is important when simulating insulators.

Avoid zeros

Zeros are bad because they cause divide by zero errors. So don’t have zero mobilities, carrier cross sections, tail slopes or densities of states. It’s fine to have zero recombination constants though.

Very big steps in the band gap

Big steps in the band gap will produce very small and very large carrier densities - see Avoiding very big and very small numbers above.