Energy Filtering in Doping Modulated Nanoengineered Thermoelectric Materials: A Monte Carlo Simulation Approach

Using Monte Carlo electronic transport simulations, coupled self-consistently with the Poisson equation for electrostatics, we explore the thermoelectric power factor of nanoengineered materials. These materials consist of alternating highly doped and intrinsic regions on the scale of several nanometers. This structure enables the creation of potential wells and barriers, implementing a mechanism for filtering carrier energy. Our study demonstrates that by carefully designing the nanostructure, we can significantly enhance its thermoelectric power factor compared to the original pristine material. Importantly, these enhancements stem not only from the energy filtering effect that boosts the Seebeck coefficient but also from the utilization of high-energy carriers within the wells and intrinsic barrier regions to maintain relatively high electronic conductivity. These findings can offer guidance for the design and optimization of new-generation thermoelectric materials through improvements in the power factor.


I. INTRODUCTION
In the past several years, large efforts have been devoted to the field of thermoelectric (TE) materials, resulting in significant progress in their performance [1][2][3][4].
The performance of TE materials is generally quantified by the figure of merit ZT , which is a measure of the ability of a material to convert heat into electricity.
The improvement in the figure of merit primarily arises from the significant reduction in lattice thermal conductivity observed in nanostructured materials, reaching values approaching the amorphous limit of κ l = 0.1 − 2 W/mK and even lower in certain cases [1,6,[18][19][20][21][22]. For instance, amorphous silicon typically has a thermal conductivity between 0.1 and 2 W/m-K, with common values around 0.7 to 1.5 W/mK for thin films [23].Similarly, amorphous germanium has a thermal conductivity between 0.1 and 1 W/mK, depending on the conditions and sample preparation [24].Hierarchical nanostructuring, a highly successful approach for reducing thermal conductivity [6,7], involves the introduction of embedded atomic defects, nanoinclusions, and grain boundaries, leading to extensive phonon scattering and consequently lower thermal conductivity.These structural distortions scatter phonons across various wavelengths, effectively reducing phonon transport across the entire spectrum [25].Nonetheless, despite achieving ultralow thermal conductivities, sometimes falling well below the amorphous limit, further enhancements to ZT are expected to stem primarily from improvements in the power factor.
The challenge in enhancing the power factor (P F ) is attributed to the adverse interdependence between electronic conductivity and the Seebeck coefficient via the carrier density, which tends to keep the P F low.However, a promising direction to address this challenge lies in energy filtering techniques, achieved through the incorporation of energy barriers.These barriers effectively block low-energy carriers while permitting the flow of high-energy carriers [7,[26][27][28][29][30][31][32][33][34].Consequently, the Seebeck coefficient (which is a measure of the energy of current flow) increases.Various methods for implementing energy filtering exist, including the use of superlattices [26,35], poly/nanocrystalline materials (as illustrated in Figure 1a, where barriers are formed on the grain boundaries) [36,37], and materials featuring dislocation loops [38], among others.Indeed, recent efforts focusing on employing energy filtering for the grain/grainboundary systems across several materials have shown improved thermoelectric performance [31][32][33].
In this computational study, we investigate the TE performance of material systems in which the wells and energy filtering potential barriers are electrostatically formed by alternating heavily doped and intrinsic regions, as depicted in Figure 1b.This investigation is inspired by our previous work on such energy filtering systems, as outlined in Ref. [39], utilizing simple series resistance models, as well as experimental studies on nanocrystalline materials (Figure 1a).These studies revealed that ultra-high P F s can be achieved by forming energy barriers at grain boundaries, surpassing five times the optimal pristine material P F [36,37].In Ref. [39], we presented the concept of achieving such PF improve- (a) A nanocrystalline material, having notably high thermoelectric power factors through high doping and the formation of potential barriers along grain boundaries.(b) A nanoengineered controlled version of (a), featuring regularly spaced highly doped square regions separated by undoped regions.(c) A schematic 1D cross-section of the conduction band Ec profile for the nanostructure depicted in (b), derived from the self-consistent solution of the Poisson equation.LW and LB denote the widths of the well and barrier, respectively, and EF represents the Fermi level.Consequently, the energy filtering features allow electrons at higher energies to propagate while hindering electrons at lower energies, depicted with green and pink arrow respectively.ments in detail.However, the simple series resistance analytical model we used, by its nature, includes many uncertainties, simplifications, and assumptions with regard to the treatment of electronic and thermoelectric transport.The purpose of our current study is to reevaluate this design concept using a full-scale advanced Monte Carlo simulation formalism, which relaxes many of the approximations of the simplistic series resistance model and provides more confidence in the proposed design, a capability that we did not have at the time of the original concept paper.Other than transport specifics, one important addition in this work is that we employ the self-consistent solution of the Poisson equation to obtain the accurate potential barriers and their shapes for a specific underlying doping profile.
We employ Monte Carlo (MC) simulations to solve the semi-classical Boltzmann transport equation (BTE) in a 2D domain, based on an advanced method and code that we have developed specifically for efficient transport computation in nanostructures [40].The basic structure we consider is shown in Figure 1b, with blue domains representing highly doped regions and gray domains representing undoped regions.This approach is then coupled self-consistently with the 2D Poisson equation to capture the electrostatics of the domain and the shape of the energy band profile, E c , that consists of electrostatic potential barriers and wells, as depicted in Figure 1c through a 1D schematic example.
Our computational framework accounts for phonon scattering, as well as phonon plus ionized impurity scattering (IIS), with the latter being particularly relevant in potential well regions.We then present a novel design for nanostructured materials capable of significant P F enhancements.This design relies on several key concepts or "ingredients", whose contributions are comprehensively described in Refs.[36,37], and involves the following: (i) energy filtering in the presence of potential barriers, (ii) reducing the well length to effectively transfer the Seebeck effect of the barrier into most of the well region (i.e., allowing carriers to propagate at high energies), (iii) high doping to utilize carriers with high velocities and position the Fermi level near the top of the potential barrier, and (iv) an undoped barrier region to mitigate the reduction in electrical conductivity introduced by the barriers (referred to as the "clean filtering" approach [39]).
The paper is structured as follows: In Section II, we introduce our approach and simulation method along with the key features of the design parameters.Section III describes the P F performance by investigating the improvements as a consequence of the design.We summarize our findings and conclude in Section IV.

II. METHOD AND APPROACH
In our approach, we introduce a well-barrier structure within a defined domain by selectively doping the domain in a periodic manner, as depicted by the dark-colored square regions in Figure 1b.This doping variation induces an electrostatic potential energy variation, which is computed using the Poisson equation, where ϕ is the potential, ρ is the charge density, and ϵ 0 (ϵ r ) is the vacuum (relative) permittivity of the medium.The charge density in the doped semiconductor is obtained by where g(E) is the density of states and f (E) is the Fermi-Dirac distribution function.We then solve the Poisson equation self-consistently across the entire 2D domain (with dimensions of 110 nm × 90 nm), subject to Neumann boundary conditions at the boundaries.This computational process yields the charge density distribution and determines the profile of the conduction band energy E c .Subsequently, the E c profile serves as an input for the Monte Carlo ray-tracing simulation to extract the electron flux.The term "ray tracing" is well-known in fields like particle physics and computer graphics for Monte Carlo simulations.But, similarly, it has been adopted in the electronic transport community, primarily in the past, to trace the pathways for electrons in transistor devices in the presence of potential variations, boundaries, and scattering events [41,42].Here, we have adopted this approach to investigate how charge carriers move through the nanostructured domains in real space and use the term "ray-tracing" as a reference to keeping track of the electronic trajectories and the time electrons spend in the domain.
Our ray-tracing method, recently developed and detailed in Ref. [40], is specifically tailored for intricate nanostructured materials and optimized for complex thermoelectric material geometries, ensuring computational efficiency.While comprehensive methodological details are presented in Ref. [40], we provide a brief overview for completeness, as depicted in the flowchart in Figure 2.
Our method adopts a single-particle incident flux, where electrons are initialized sequentially at the left boundary of the domain and proceed either to the opposite side or undergo backscattering.Instead of randomly selecting free-flight times and implementing selfscattering as in previous Monte Carlo approaches [43], we use a mean-free-path (mfp, λ) approach [44], where electrons propagate an mfp, followed by definite scattering after each mfp is completed.For electron-phonon scat-tering, we assume an mfp of λ ph = 30 nm, and for the bandstructure, we use a parabolic band effective mass m * = 0.9m 0 .This choice leads to a low carrier density mobility of n-type Si of roughly 1500 cm 2 /V.s (here, we consider only elastic phonon scattering for simplicity).We emphasize, though, that it did not attempt to simulate specific materials, as this concept is material agnostic.We also calculate the IIS rate using the Brooks-Herring model combined with the strong screening scattering model in each discretization cell within the simulation domain, which aligns better with experimental mobility measurements, particularly for pristine Si [36,39].Using the band structure velocity, we determine the mfp (λ IIS ) due to IIS and then employ Matthiessen's rule to compute the total mean free path.An electron in a simulation domain cell advances one total mfp, λ total , before undergoing enforced scattering.Notably, because of the IIS process, λ total varies as a function of energy (although λ ph is generally energy-independent, at least for acoustic phonon scattering).
In a 2D domain featuring periodically arranged doped square sections ( see Figure 1b), the spatial profile of E c across the domain is obtained from the solution of Equation ( 1), with a typical example shown in Figure 3b.For the calculation of the electron flux, we inject electrons at all energies from the left side of the domain and all possible directions.Electrons undergo free flight and scattering events, and we ray-trace their paths until they exit the domain from the right side.We record the time an electron spent in the domain and express this as its time of flight (ToF).Electrons backscattered to the left are excluded from flux calculations as they do not contribute.The inverse of this quantity ensembled from many trajectories provides the flux [43,45].We compute the average ToF per particle and utilize it to determine the flux per simulated electron at each energy: Multiplying the flux per electron by the density of states, g(E), allows us to estimate the transport function as a function of energy.
where the constant C in the TDF equation accounts for two things: the super-electron charge used in Monte Carlo simulations, because we only simulate a limited number of electrons, and the geometric factors due to simulating in a finite 2D domain instead of an infinite 3D domain.This constant is determined to map the conductivity computed from Monte Carlo to the bulk material conductivity that we can compute analytically (see [40] for details).This relationship is analogous to the transport distribution function (TDF) Ξ(E) in the BTE.Essentially, the product of flux with the density of states represents the flow of charge, which is directly correlated with the conductivity, much like how the TDF influences conductivity [40].We calculate the conductivity and the Seebeck coefficient as where q is the electron charge and k B is the Boltzmann constant, while we set T = 300 K in all simulations.

III. POWER FACTOR IN A STRUCTURE WITH DOPED/INTRINSIC REGIONS
We begin our investigation by simulating the TE properties of a 2D nanostructure with dimensions L = 110 nm (length) and W = 90 nm (width), with periodically doped and undoped regions, as depicted in Figure 3a.The Fermi energy is set to E F = 0, and we consider three sizes of the doped regions: L W = 10, 15, and 22 nm, corresponding to barrier widths of L B ≈ 19, 13, and 6 nm, respectively.We begin by solving the Poisson equation self-consistently for a doping density of N D = 3.5 × 10 26 m −3 to determine the profile of the conduction band energy (E c ), as shown in Figure 3b for L W = 22 nm.We deliberately select a high density in accordance with previous studies [36,37,39].In Figure 3c, we present a crosssectional view of the 2D E c profile at 45 nm across the domain for all three sizes of the doped regions, i.e., passing through the middle of the doped wells.Increasing the well size shrinks the barrier region and reduces the barrier height (compare the purple and blue lines in Figure 3c).Increasing the doping density leads to deeper wells (and lower barriers) and shifts the top of the barriers toward the Fermi level; see Figure 3d for the band profile (and Figure 3e for the carrier density in the L W = 22 nm scenario).Thus, the barrier height and its distance from the Fermi level, which determines the Seebeck coefficient and energy filtering, can be controlled by the geometrical features in the domain (doped/intrinsic region sizes) and the doping level.
For comparison, we first calculate the TE coefficients of the pristine material for the two cases: (i) pristine material subjected only to electron-phonon (e-ph) scattering (representing the ultimate performance of the material) and (ii) doped pristine material undergoing e-ph and IIS scattering, reflecting the realistically achievable perfor-mance of the material.It is worth noting that, for IIS, both the scattering rate and the screening length depend on the doping concentration [43].In Figure 4, we present the conductivities, Seebeck coefficients, and P F s as functions of charge density for these cases.For the pristine material in the absence of IIS (solid blue lines), the P F exhibits a peak at high densities (around 10 25 −10 26 /m 3 ), exceeding 2× the P F in the presence of IIS (dashed blue lines).This is mainly due to the suppression of conductivity, which is primarily due to the stronger influence of IIS compared to phonon-limited scattering, despite the slight increase in the Seebeck coefficient.
When the nanostructure is periodically doped, electrons traverse through alternating potential wells and potential barriers; see Figure 3c.In this context, we simulate three systems with different well lengths, as indicated in the legend of Figure 3a.For each system, we vary the doping density within the range N D = 3 × 10 24 − 3 × 10 26 /m 3 .We compare the P F results of these structures with the realistically achieved performance of the pristine material (dashed blue line), as, in practice, doping facilitates high carrier densities despite reduced mobility due to IIS.We also compare the results of the filtering structures with the solid blue line, the ultimate performance achievable by our fictitious material.This comparison covers scenarios where gating or modulation doping enables high carrier densities in the absence of dopants in the lattice.
In the filtering systems, electrons in the wells are not filtered and traverse over the barriers and possess higher energies, indicating a larger Seebeck coefficient.However, the formation of barriers simultaneously decreases conductivity.This is observed in Figure 4 for L W = 10 and 15 nm.Despite the lower conductivity compared to the pristine case, the Seebeck coefficient is large, leading to a slight increase in the P F at high densities (around an order of magnitude higher than the densities at which the P F peaks in the pristine case).As the well size increases to L W = 22 nm, the barriers become narrower (L B ≈ 6 nm), and their height decreases, which mitigates the resistance they introduce.At such higher densities (as shown in Figure 3c-e), the resistance is further mitigated by the use of higher velocity carriers coming from the wells, resulting in a significant increase in conductivity.Note that the Seebeck coefficient decreases as the barrier height reduces, but it still remains considerably higher compared to the pristine material case.This increase in the Seebeck coefficient, in combination with the relative robustness of the electrical conductivity, leads to high P F improvements.In fact, the P F reaches approximately 33 mW m −1 K −2 for L W = 22 nm at a charge density of 7 × 10 26 m −3 .In Figure 4c, we indicate with black circle the four cases for which we simulate the band profile and carrier density in Figure 3d-e.Across all cases, the P F is enhanced compared to the pristine material, with the most significant benefits observed at higher well doping and narrower/reduced barrier height cases.It is important to note that our simulations re-veal a monotonic increase in P F .This increase will persist as long as the top of the barrier remains above E F .Our final simulation point for high density lies at the borderline, as shown in Figure 3d.Beyond this point, we anticipate a reduction in P F , although we encounter convergence issues in our simulations for such ultra-high doping densities.We also indicate some findings of the PF in Figure 4c, marked with black stars.These findings align with recent experimental results showing similar power factor enhancements in nanocrystalline materials [36][37][38]46].
Typically, the introduction of potential barriers in material systems reduces electrical conductivity, as carriers encounter exponentially increased difficulty in traversing over the barriers.However, in the case where L W = 22 nm and for the highest density considered (right-most purple line in Figure 4), the electrical conductivity is relatively high, reaching up to half of what the pristine material can offer (Figure 4a, comparing the purple and bluedashed lines at the highest density region).There are a couple of reasons why the conductivity is not substantially degraded despite the barrier regions posing higher resistance: (i) the barrier height, in that case, is relatively small (as seen with the purple line in Figure 3c), and (ii) the barriers are undoped, signifying that carriers entering those regions having high mobility, compensating for the lower velocities of electrons near the band edge.
To provide more insight, we illustrate the average flux of a single electron as a function of energy, computed via Monte Carlo, in Figure 5.This quantity represents the inverse of the average time of flight taken for an electron to travel from the left to the right contacts of the material system.It considers the different scattering mfps and carrier bandstructure velocities across different regions of the material.We take as reference energy the position of the band profile in the potential well, i.e., E C,min = −200 meV, and show two cases.The flux for a pristine highly doped material (blue line in Figure 5) and the flux for the well-barrier system yield the highest reported P F (rightmost point on the purple line in Figure 4c).The flux for the well-barrier structure (red line in Figure 5) starts at the energy of the barrier height, which is a few meV above E F = 0 eV in this case.A key observation is a sharp increase in flux after encountering the barrier, reaching approximately 70 percent of the value for the pristine material within a few meV.This observation explains the limited degradation of conductivity by potential barriers in this specific case, as observed in Figure 4a.

IV. SUMMARY AND CONCLUSIONS
In this study, we investigated the thermoelectric properties of nanoengineered materials featuring highly doped regions periodically separated from undoped regions.Employing Monte Carlo simulations for electronic transport, coupled self-consistently with the Poisson equation to capture the electrostatics of the domain, we demon- strated the potential of these structured materials to achieve thermoelectric power factors up to five times higher than the optimal power factor of pristine materials.These findings, derived from advanced numerical simulations and software, confirm predictions made by much simpler models [39] and are in line with recent experimental results showcasing such power factor enhancements in nanocrystalline materials [36,37,46].Our results can motivate the significance of further research into nanostructured thermoelectric materials aiming to achieve ultra-high power factors through efficient utilization of energy filtering mechanisms.
FIG. 1.(a) A nanocrystalline material, having notably high thermoelectric power factors through high doping and the formation of potential barriers along grain boundaries.(b) A nanoengineered controlled version of (a), featuring regularly spaced highly doped square regions separated by undoped regions.(c) A schematic 1D cross-section of the conduction band Ec profile for the nanostructure depicted in (b), derived from the self-consistent solution of the Poisson equation.LW and LB denote the widths of the well and barrier, respectively, and EF represents the Fermi level.Consequently, the energy filtering features allow electrons at higher energies to propagate while hindering electrons at lower energies, depicted with green and pink arrow respectively.

FIG. 3 .
FIG. 3. (a) The 2D nanostructure with periodically doped regions (square regions) and dopant-depleted regions simulated.Three sizes for the well regions are considered: LW = 10, 15, and 22 nm.(b) Conduction band Ec profile for LW = 22 nm as calculated from the self-consistent solution of the Poisson equation.(c) A 1D cross-sectional view of the Ec profile along the domain at 45 nm.(d,e) The behavior of the Ec and carrier density for the LW = 22 nm case channel with increasing doping density ND.

FIG. 4 .FIG. 5 .
FIG. 4. (a) Conductivity, (b) Seebeck coefficient, and (c) P F versus charge density for (i) pristine material with e-ph (solid blue line) and e-ph plus ionized impurity scattering (blue dashed line) and (ii) for the periodically doped nanostructures of three different well sizes as indicated.Notably, for a relatively large well size (LW = 22 nm), the P F shows a significant increase at high densities.In (c), the data points encircled correspond to the charge and potential profiles depicted in Figure3d-e.The gray star marked data show recent experimental measurements from[36][37][38]46], which show similar power factor enhancements in nanocrystalline materials.
[40]2.Flowchart that describes the calculation method for electron transport in the domain.It utilizes the selfconsistent solution of the Poisson equation and a Monte Carlo ray-tracing algorithm, as elaborated in Ref.[40].