CIGS: stationary with traps and Schottky contacts.

(source code)

Simulating stationary charge transport for CIGS with hole traps and mixed Schottky/Ohmic contact conditions. Assume that SRH recombination only happens within a small regime.

module Ex108_CIGS_WithTraps

using ChargeTransport
using ExtendableGrids
using PyPlot

# function to initialize the grid for a possible extension to other p-i-n devices.
function initialize_pin_grid(refinementfactor, h_ndoping, h_pdoping_left, h_pdoping_trap, h_pdoing_right)
    coord_ndoping = collect(range(0.0, stop = h_ndoping, length = 2 * refinementfactor))
    coord_pdoping_left = collect(range(h_ndoping, stop = (h_ndoping + h_pdoping_left), length = 3 * refinementfactor))
    coord_pdoping_plus = collect(
        range(
            (h_ndoping + h_pdoping_left),
            stop = (h_ndoping + h_pdoping_left + h_pdoping_trap),
            length = refinementfactor
        )
    )
    coord_pdoping_right = collect(
        range(
            (h_ndoping + h_pdoping_left + h_pdoping_trap),
            stop = (h_ndoping + h_pdoping_left + h_pdoping_trap + h_pdoing_right),
            length = 3 * refinementfactor
        )
    )
    coord = glue(coord_ndoping, coord_pdoping_left)
    coord = glue(coord, coord_pdoping_plus)
    coord = glue(coord, coord_pdoping_right)

    return coord
end

you can also use other Plotters, if you add them to the example file you can set verbose also to true to display some solver information

function main(; n = 3, Plotter = PyPlot, plotting = false, verbose = "", test = false, AdditionalTrapSpecies = false)

    if plotting
        Plotter.close("all")
    end
    ################################################################################
    if test == false
        println("Set up grid and regions")
    end
    ################################################################################

    # region numbers
    regionDonor = 1                           # n doped region
    regionAcceptorLeft = 2                           # p doped region
    regionAcceptorTrap = 3                           # p doped region with trap
    regionAcceptorRight = 4                           # p doped region
    regions = [regionDonor, regionAcceptorLeft, regionAcceptorTrap, regionAcceptorRight]
    numberOfRegions = length(regions)

    # boundary region numbers
    bregionDonor = 1
    bregionAcceptor = 2
    bregionDALeft = 3
    bregionALeftATrap = 4
    bregionATrapARight = 5

    # grid
    refinementfactor = 2^(n - 1)
    h_ndoping = 0.5 * μm
    h_pdoping_left = 1.0 * μm
    h_pdoping_trap = 0.1 * μm
    h_pdoing_right = 1.0 * μm
    w_device = 0.5 * μm  # width of device
    z_device = 1.0e-4 * cm  # depth of device
    h_total = h_ndoping + h_pdoping_left + h_pdoping_trap + h_pdoing_right
    coord = initialize_pin_grid(
        refinementfactor,
        h_ndoping,
        h_pdoping_left,
        h_pdoping_trap,
        h_pdoing_right
    )

    grid = simplexgrid(coord)

    # set different regions in grid
    cellmask!(grid, [0.0 * μm], [h_ndoping], regionDonor) # n doped
    cellmask!(grid, [h_ndoping], [h_ndoping + h_pdoping_left], regionAcceptorLeft) # p doped
    cellmask!(grid, [h_ndoping + h_pdoping_left], [h_ndoping + h_pdoping_left + h_pdoping_trap], regionAcceptorTrap) # p doped with traps
    cellmask!(grid, [h_ndoping + h_pdoping_left + h_pdoping_trap], [h_total], regionAcceptorRight) # p doped

    bfacemask!(grid, [h_ndoping], [h_ndoping], bregionDALeft, tol = 1.0e-18)
    bfacemask!(grid, [h_ndoping + h_pdoping_left], [h_ndoping + h_pdoping_left], bregionALeftATrap, tol = 1.0e-18)
    bfacemask!(grid, [h_ndoping + h_pdoping_left + h_pdoping_trap], [h_ndoping + h_pdoping_left + h_pdoping_trap], bregionATrapARight, tol = 1.0e-18)

    if plotting
        gridplot(grid, Plotter = Plotter, legend = :lt)
        Plotter.title("Grid")
    end

    if test == false
        println("*** done\n")
    end

    ################################################################################
    if test == false
        println("Define physical parameters and model")
    end
    ################################################################################

    iphin = 1 # index electron quasi Fermi potential
    iphip = 2 # index hole quasi Fermi potential
    if AdditionalTrapSpecies
        iphit = 3 # index trap quasi Fermi potential
        numberOfCarriers = 3 # electrons, holes and traps
    else
        numberOfCarriers = 2 # electrons and holes
    end

    # physical data
    T = 300.0 * K

    # band edge energies
    Ec_ZnO = 3.4 * eV
    Ev_ZnO = 0.0 * eV

    Ec_CIGS = 3.4 * eV
    Ev_CIGS = 2.3 * eV

    EC = [Ec_ZnO, Ec_CIGS, Ec_CIGS, Ec_CIGS]
    EV = [Ev_ZnO, Ev_CIGS, Ev_CIGS, Ev_CIGS]

    # effective densities of states
    Nc = 4.35195989587969e17 / (cm^3)
    Nv = 9.139615903601645e18 / (cm^3)

    NC = [Nc, Nc, Nc, Nc]
    NV = [Nv, Nv, Nv, Nv]

    # mobilities
    mun_ZnO = 100 * (cm^2) / (V * s)
    mup_ZnO = 25 * (cm^2) / (V * s)
    mun_CIGS = 100.0 * (cm^2) / (V * s)
    mup_CIGS = 25 * (cm^2) / (V * s)

    μn = [mun_ZnO, mun_CIGS, mun_CIGS, mun_CIGS]
    μp = [mup_ZnO, mup_CIGS, mup_CIGS, mup_CIGS]

    # relative dielectric permittivity
    εr_ZnO = 9 * 1.0
    εr_CIGS = 13.6 * 1.0

    ε = [εr_ZnO, εr_CIGS, εr_CIGS, εr_CIGS]

    # trap information
    zt = 1 # hole traps
    Et = 2.8 * eV
    ET = [0.0, 0.0, Et, 0.0]
    Nt = 1.0e18 / (cm^3)
    NT = [0, 0, Nt, 0]

    mu_t = 0 * (cm^2) / (V * s)
    μt = [0.0, 0.0, mu_t, 0.0]

    # recombination information parameters
    ni_ZnO = sqrt(Nc * Nv) * exp(-(Ec_ZnO - Ev_ZnO) / (2 * kB * T))   # intrinsic concentration
    n0_ZnO = Nc * Boltzmann((Et - Ec_ZnO) / (kB * T))                   # Boltzmann equilibrium concentration
    p0_ZnO = ni_ZnO^2 / n0_ZnO                                        # Boltzmann equilibrium concentration
    ni_CIGS = sqrt(Nc * Nv) * exp(-(Ec_CIGS - Ev_CIGS) / (2 * kB * T)) # intrinsic concentration
    n0_CIGS = Nc * Boltzmann((Et - Ec_CIGS) / (kB * T))                  # Boltzmann equilibrium concentration
    p0_CIGS = ni_CIGS^2 / n0_CIGS                                      # Boltzmann equilibrium concentration

    p0 = [p0_ZnO, p0_CIGS, p0_CIGS, p0_CIGS]
    n0 = [n0_ZnO, n0_CIGS, n0_CIGS, n0_CIGS]

set the lifetime value high in all other regions, such that SRH recombination can be neglected there

    SRH_LifeTime = [1.0e100, 1.0e100, 1.0e-3 * ns, 1.0e100]

    Auger = 1.0e-29 * cm^6 / s
    Radiative = 1.0e-10 * cm^3 / s

    # Schottky contact information
    An = 4 * pi * q * mₑ * kB^2 / Planck_constant^3
    Ap = 4 * pi * q * mₑ * kB^2 / Planck_constant^3
    vn = An * T^2 / (q * Nc)
    vp = Ap * T^2 / (q * Nv)
    barrier = 0.7 * eV

    # doping information
    Nd = 1.0e18 / (cm^3)
    Na = 5.5e15 / (cm^3)

    # we will impose this applied voltage on one boundary
    voltageAcceptor = 1.0 * V

    if test == false
        println("*** done\n")
    end
    ################################################################################
    if test == false
        println("Define System and fill in information about model")
    end
    ################################################################################

    # initialize Data instance and fill in data
    data = Data(grid, numberOfCarriers)
    data.modelType = Stationary # R = Rn = Rp, since the model type is stationary
    if AdditionalTrapSpecies
        data.F = [FermiDiracOneHalfTeSCA, FermiDiracOneHalfTeSCA, FermiDiracMinusOne]
    else
        data.F .= FermiDiracOneHalfTeSCA
    end

    data.bulkRecombination = set_bulk_recombination(;
        iphin = iphin, iphip = iphip,
        bulk_recomb_Auger = true,
        bulk_recomb_radiative = true,
        bulk_recomb_SRH = true
    )

    data.boundaryType[bregionAcceptor] = SchottkyContact
    data.boundaryType[bregionDonor] = OhmicContact
    data.fluxApproximation .= ExcessChemicalPotential

    if AdditionalTrapSpecies
        # Here, we enable the traps and parse the respective index and the regions where the trap is defined.
        enable_trap_carrier!(; data = data, trapCarrier = iphit, regions = [regionAcceptorTrap])
    else
        # pass trap data in stationary setting since there is no separate trap species
        add_trap_density_Poisson!(data = data, zt = zt, Nt = NT)
    end

    if test == false
        println("*** done\n")
    end

    ################################################################################
    if test == false
        println("Define Params and fill in physical parameters")
    end
    ################################################################################

    # physical parameters
    params = Params(grid[NumCellRegions], grid[NumBFaceRegions], numberOfCarriers)
    params.temperature = T
    params.UT = (kB * params.temperature) / q
    params.chargeNumbers[iphin] = -1
    params.chargeNumbers[iphip] = 1
    if AdditionalTrapSpecies
        params.chargeNumbers[iphit] = zt
    end

    for ireg in 1:numberOfRegions           # interior region data

        params.dielectricConstant[ireg] = ε[ireg] * ε0

        # effective DOS, band-edge energy and mobilities
        params.densityOfStates[iphin, ireg] = NC[ireg]
        params.densityOfStates[iphip, ireg] = NV[ireg]
        params.bandEdgeEnergy[iphin, ireg] = EC[ireg]
        params.bandEdgeEnergy[iphip, ireg] = EV[ireg]
        params.mobility[iphin, ireg] = μn[ireg]
        params.mobility[iphip, ireg] = μp[ireg]

        if AdditionalTrapSpecies
            params.densityOfStates[iphit, ireg] = NT[ireg]
            params.bandEdgeEnergy[iphit, ireg] = ET[ireg]
            params.mobility[iphit, ireg] = μt[ireg]
        end

        # recombination parameters
        params.recombinationRadiative[ireg] = Radiative
        params.recombinationSRHLifetime[iphin, ireg] = SRH_LifeTime[ireg]
        params.recombinationSRHLifetime[iphip, ireg] = SRH_LifeTime[ireg]
        params.recombinationSRHTrapDensity[iphin, ireg] = n0[ireg]
        params.recombinationSRHTrapDensity[iphip, ireg] = p0[ireg]
        params.recombinationAuger[iphin, ireg] = Auger
        params.recombinationAuger[iphip, ireg] = Auger

    end

    # doping -- since we do not set any doping for the traps it is automatically zero
    params.doping[iphin, regionDonor] = Nd
    params.doping[iphip, regionAcceptorLeft] = Na
    params.doping[iphip, regionAcceptorTrap] = Na
    params.doping[iphip, regionAcceptorRight] = Na

    # values for the schottky contacts
    params.SchottkyBarrier[bregionAcceptor] = barrier
    params.bVelocity[iphin, bregionAcceptor] = vn
    params.bVelocity[iphip, bregionAcceptor] = vp

    data.params = params
    ctsys = System(grid, data, unknown_storage = :sparse)

    if test == false
        show_params(ctsys)
        println("*** done\n")
    end

    ################################################################################
    if test == false
        println("Define control parameters for Solver")
    end
    ################################################################################

    control = SolverControl()
    control.verbose = verbose
    control.tol_round = 1.0e-7
    control.damp_initial = 0.5
    control.damp_growth = 1.2
    control.maxiters = 30
    control.max_round = 3

    if test == false
        println("*** done\n")
    end

    ################################################################################
    if test == false
        println("Compute solution in thermodynamic equilibrium")
    end
    ################################################################################

    # solve thermodynamic equilibrium and update initial guess
    solution = equilibrium_solve!(ctsys, control = control)
    inival = solution

    if plotting
        label_solution, label_density, label_energy = set_plotting_labels(data)

        if AdditionalTrapSpecies
            # add labels for traps
            label_energy[1, iphit] = "\$E_{\\tau}-q\\psi\$"; label_energy[2, iphit] = "\$ - q \\varphi_{\\tau}\$"
            label_density[iphit] = "\$n_{\\tau}\$";        label_solution[iphit] = "\$ \\varphi_{\\tau}\$"
        end

        # ##### set legend for plotting routines #####
        Plotter.figure()
        plot_energies(Plotter, ctsys, solution, "Equilibrium", label_energy)
        Plotter.figure()
        plot_densities(Plotter, ctsys, solution, "Equilibrium", label_density)
        Plotter.figure()
        plot_solution(Plotter, ctsys, solution, "Equilibrium", label_solution)
    end

    if test == false
        println("*** done\n")
    end

    ################################################################################
    if test == false
        println("Stationary bias loop")
    end
    ################################################################################

    endVoltage = voltageAcceptor       # final bias value
    biasValues = collect(range(0, stop = endVoltage, length = 52))

    IV = zeros(0)
    chargeDensities = zeros(0)

    for i in eachindex(biasValues)

        Δu = biasValues[i] # bias

        # Apply new voltage: set non equilibrium boundary conditions
        set_contact!(ctsys, bregionAcceptor, Δu = Δu)

        if test == false
            println("bias: Δu = $(Δu) V")
        end

        # solve time step problems with timestep Δt
        solution = solve(ctsys, inival = inival, control = control)
        inival = solution

        # save IV data
        current = get_current_val(ctsys, solution)
        push!(IV, w_device * z_device * current)

        # store charge density in donor region (ZnO)
        push!(chargeDensities, charge_density(ctsys, solution)[regionDonor])


    end # bias loop

    # compute static capacitance: check this is correctly computed
    staticCapacitance = diff(chargeDensities) ./ diff(biasValues)

    # plot solution and IV curve
    if plotting
        Plotter.figure()
        plot_energies(Plotter, ctsys, solution, "bias \$\\Delta u\$ = $(endVoltage) V", label_energy)
        Plotter.figure()
        plot_densities(Plotter, ctsys, solution, "bias \$\\Delta u\$ = $(endVoltage) V", label_density)
        Plotter.figure()
        plot_solution(Plotter, ctsys, solution, "bias \$\\Delta u\$ = $(endVoltage) V", label_solution)
        Plotter.figure()
        plot_IV(Plotter, biasValues, IV, "bias \$\\Delta u\$ = $(biasValues[end]) V", plotGridpoints = true)
        Plotter.figure()
        plot_IV(Plotter, biasValues, chargeDensities, "bias \$\\Delta u\$ = $(biasValues[end]) V", plotGridpoints = true)
        Plotter.title("Charge density in donor region")
        Plotter.ylabel("Charge density [C]")
        Plotter.figure()
        plot_IV(Plotter, biasValues, staticCapacitance, "bias \$\\Delta u\$ = $(biasValues[end]) V", plotGridpoints = true)
        Plotter.title("Static capacitance in donor region")
        Plotter.ylabel("Static capacitance [C/V]")

    end

    if test == false
        println("*** done\n")
    end

    testval = sum(filter(!isnan, solution)) / length(solution) # when using sparse storage, we get NaN values in solution
    return testval

end #  main

function test()
    testval = 1.484831264268335
    testvalAdditionalSpecies = 1.1334257649339574

    return main(test = true, AdditionalTrapSpecies = false) ≈ testval && main(test = true, AdditionalTrapSpecies = true) ≈ testvalAdditionalSpecies
end

if test == false
    println("This message should show when this module has successfully recompiled.")
end


end # module

This page was generated using Literate.jl.