Issue with boundary conditions in depletion approximations

[phoebe-p]

Describe the bug There is an error in the depletion approximation which means that if the surface recombination velocity (SRV) corresponding to the rear surface is increased, the EQE (and current) actually increase and can eventually be > 1. This appears to be due to a sign error on the boundary conditions.

To Reproduce See file attached which (under current Solcore version 5.9.1) produces the plot below (issue noted and example created by Zubair Abdullah-Vetter at UNSW).

Expected behavior Increasing SRV should never lead to increased EQE/current. Certainly EQE should not be > 1 (i.e. current exceeds possible current through photogeneration, IQE > 1).

Additional context The relevant equations for the boundary conditions from Jenny Nelson's book are 6.18 (front surface) and 6.23 (rear surface). Because the electron and hole current go in opposite directions, the signs are flipped, but this is not reflected in depletion_approximation (in both the wavelength-dependent and integrated versions of the code):

    if side == 'top':
        def bc(ya, yb):
            left = ya[1] - S / D * (ya[0] - y0)
            right = yb[0]
            return np.array([left, right])
    else:
        def bc(ya, yb):
            left = ya[0]
            right = yb[1] - S / D * (yb[0] - y0)
            return np.array([left, right])

This is also incorrect in the Green's function implementation, where the same boundary conditions is also applied to both carrier types. GaAs_EQE_BL_looped.txt

[phoebe-p]

@vathomass hello! See above - a user identified an issue with the depletion approximation solver in Solcore. Interestingly, the behaviour is exactly the same whether I use the solve_bvp method or the Green's function solver you implemented (see plot attached mode with 'da_mode' : 'green'), leading me to think there is a similar issue with the boundary condition enforcement in both methods. I have (I believe) fixed the issue in the solve_bvp method (see #253), but I don't understand the details of your Green's function code so I am not sure what to change (in the solve_bvp method it was just changing the sign in front of the (S/D) term to a + for the bottom half of the junction).

[phoebe-p]

@dalonsoa (not expecting you to fix this, just tagging you in case you have any thoughts or some more knowledge of how the Green’a function works)

[vathomass]

Hi. Firstly, the fact that the error is independent of the solving method is expected. During the implementation I have derived the Green's function method from (and also tested against) the solve_bvp equations. That is to say that, in theory, both methods solve the same boundary value problem.

I had a quick look at changes in #253, where I also see several stylistic/PEP8 changes, so I have might missed something. I see that both boundary conditions have changed (around l. 382 named as leftand right). The left boundary condition is also updated in the "bottom" case. Is this correct?

Moreover, I do not really understand what you mean when you say

(in both the wavelength-dependent and integrated versions of the code)

in the Additional context section.

My first thought is that changing the two boundary conditions is not too hard in the Green's function code. I suggest something like (around l. 606)

    # if L too low in comparison to junction width, avoid nan's
    if xbL > 1.0e2:
        if side == "top":
            cadd = 0.
            fun = partial(_conv_exp_top, xa=xa, xb=xb, g=g, L=L, phoD=ph_over_D)
        else:
            cadd = 0.
            fun = partial(_conv_exp_bottom, xa=xa, g=g, L=L, phoD=ph_over_D)
        cp = 1.0
    else:
        if side == "top":
            cp = -np.cosh(xbL) - crvel * np.sinh(xbL)
            cadd = 0.
            fun = partial(
                _conv_green_top, xa=xa, xb=xb, g=g, L=L, phoD=ph_over_D, crvel=crvel
            )
        else:
            cp = np.cosh(xbL) + crvel * np.sinh(xbL)
            cadd = 0.
            fun = partial(
                _conv_green_bottom, xb=xb, g=g, L=L, phoD=ph_over_D, crvel=-crvel
            )

Some explanation, you can skip if you want: the Greens function solves for the excess of the carriers, so setting right = yb[0] - y0 corresponds to cadd=0. Changing the sign of the S parameter, corresponds to changing the sign of the crvel variable which also propagates the sign change to the convolution kernels when passed as argument.

Disclaimer: I have neither implemented nor tested any of these changes myself, so I maybe I have missed something. The simplest test is to compare against the 'solve_bvp` results. If you have problems with it, please let me know, but I will need some time to check the issue myself.

[phoebe-p]

Thank you for the fast response @vathomass! First of all, for the "both the wavelength-dependent and integrated versions of the code" point, this is only relevant for the solve_bvp method (there are two functions, one of which is wavelength-resolved to calculate EQE and one where the whole current is calculated at once in case you only want the current-voltage behaviour).

Sorry for the stylistic changes making it hard to see what changed; compared to the code snippet in my original post, I now have:

    if side == 'top':
        def bc(ya, yb):
            left = ya[1] - S / D * (ya[0] - y0)
            right = yb[0] - y0
            return np.array([left, right])
    else:
        def bc(ya, yb):
            left = ya[0] - y0
            right = yb[1] + S / D * (yb[0] - y0)
            return np.array([left, right])

so the sign has changed on the right boundary condition in the bottom ('else') case, and you are correct that I also changed the other boundary condition. I did not mention this in the original issue, but looking at some textbooks I think that setting n = n0 / p = p0 at the edge of the depletion region is the correct boundary condition (so in the Greens function code that would be excess carriers = 0 at that boundary), rather than setting n = 0 / p = 0. However, in practice for solar cells this doesn't really make a difference since n0 and p0 are usually very small compared to the photogenerated carrier concentrations.

Thank you for your (very fast) input, I wasn't completely sure which parts of the Green's function code were setting the boundary conditions but I think it makes sense to me now. I will test it against the solve_bvp method, and I have also just written up the analytical solution for Beer-Lambert law absorption in Python so I can also compare both against that to make sure it is all correct.

[phoebe-p]

Thanks again @vathomass, I tried your changes as suggested and am getting good agreement with the corrected solve_bvp method and analytical solutions for Beer-Lambert absorption. Will be fixed when #253 is merged.

[vathomass]

Hi. You are welcome. Glad to know that the issue was resolved. Cheers