where n D and n e are described in Eqs. ( ) and . Since n e is a function of ϵ F (see Eq. ( )) and n D is a function of the depletion region, the Fermi energy is calculated by solving Eq. ( ) using a standard root-finding procedure in the HADOKEN m-files charge_neutral.m and find_epsilon_F.m The normalize_and_sqrt_m_triangular.m routine calculates d x d ỹ |ψ n ( x, ỹ)| 2 m * ( x, ỹ)/m 0 for each of the n wavefunctions, which is needed as input to find_epsilon_F.m. The depletion region itself is calculated/stored in the variable heaviside_n_D within the HADOKEN code. Specifically, the geometric regions of the NW that satisfy V T ( x, ỹ) ≥ E F become ionized and contribute a positive-valued dopant density, n D , that is taken into account in the charge neutrality condition of Eq. ( ). The procedure for satisfying charge neutrality is carried out during each iterative cycle until self-consistency in the total potential, V T ( x, ỹ), is reached (described further in the following section).