Computationally efficient Monte Carlo electron transport algorithm for nanostructured thermoelectric material configurations

A plethora of new materials and their alloys have recently been synthesized for a variety of applications, and typically most of them are poly- or nano-crystalline with embedded defects and often a large degree of non-uniformities, as, for example, in thermoelectric (TE) materials.¹ Specifically for TEs, these materials are nanostructured, and they are so on purpose. In order to evaluate their transport properties, we typically employ the semi-classical Boltzmann transport equation (BTE). The solution of the BTE can be evaluated either analytically or numerically.^2–6 Since it is a seven-dimensional integrodifferential equation (six dimensions in the phase space and one in time), its analytical solution is cumbersome and can be solved under very restrictive assumptions.^7–9 A grid-based deterministic numerical method, such as the renewed spherical harmonics approach, is sometimes used, but it requires very powerful computational platforms.^10–15 Another way to solve the BTE is based on stochastic solution methods using Monte Carlo (MC) computational algorithms, which use statistical sampling to solve the BTE numerically and are frequently used for electronic device applications.^16–22 Over the last several years, MC techniques also found extensive use in the fields of charge and energy transport in semiconductor materials.^23–28 MC simulation methods are particularly useful for nanostructured materials, where the carriers encounter a plethora of defects along their transport direction and interact with them, while analytical solutions in these cases are not as accurate. The simplest way to consider theoretically the influence of nanostructuring on carrier transport is to use Matthiessen’s rule to add an additional geometry dependent scattering rate on top of the internal bulk material scattering processes in the analytical BTE.^29,30 However, it is never clear how to determine this geometrical length scale in complex nanostructured materials with multiple and irregularly placed defect types. Using stochastic MC simulations, the transport in nanostructured configurations can be modeled in a real space domain, and proper insight can be gained.^24,31–35

One category of materials, which has benefited from nanostructuring over the last two decades, is thermoelectric materials, which convert the heat energy from a temperature gradient to electrical energy and vice versa.^1,36,37 Their conversion efficiency is controlled by both their charge and heat/energy transport. New-generation TE materials are typically highly nanostructured, with nano-features spanning from macro- to nano-scale (including boundaries, potential barriers, pores, nanoinclusions, atomic defects, second-phasing, etc.).^38–43 The recent improvement in the performance of TE materials originates in most cases from reduced thermal conductivity due to phonon scattering with the boundaries of the nano-defects, and, thus, MC simulations are typically performed for phonon transport in real space in such materials.^{23,31,33,44,45} Nevertheless, studies on electronic transport using MC are also emerging for nanostructured TE materials after it was pointed out that specific designs can also improve the so-called power factor, which directly determines the efficiency of a TE material as well.^46,47

Monte Carlo simulations involve the ray-tracing of particle trajectories rather than the direct solution of partial differential equations. These particles are allowed to move in the domain in both left/right directions under the influence of a driving force [as shown in Fig. 1(a)], and the net flux is computed in a statistical manner. Although this method served well simulations in bulk materials with success over many years, for nanostructured materials large difficulties are encountered, which make simulations computationally extremely expensive and logistically cumbersome. The presence of multiple scattering sites from nanostructuring reduces the flux in the domain at such a degree, which makes it very difficult to gather enough statistics for converged flux results. This is particularly difficult under linear response, where the two opposite going fluxes vary only slightly. Typical simulations in the literature require tens of thousands to even millions of trajectories for adequate results.^33,48–52

In order to address the above issues and enable efficient and accurate MC simulations for nanostructures, we have developed a hybrid MC algorithm which: (i) merges information from analytical BTE solutions with the numerically extracted flux, (ii) considers only a single flux initialized from the left only and injected into the channel where it is ray-traced to either of the contacts, and (iii) does not require the application of a driving force. The method we present provides the same accuracy as common methods, but with a significantly reduced computational cost.

This paper is organized as follows: in Sec. II, we describe the challenges encountered in existing MC algorithms for nanostructures and our method, which provides improved behavior. A stochastic error analysis and convergence is presented in Sec. III. How this method is very efficient and effective for transport in nanostructures is discussed in Sec. IV. Finally, in Sec. V, we conclude.

Many books and literature references describe the MC process;^{5,16,18,19,48,53} thus, in this work, we provide only the essential details that help with the discussion of the method we present. In general, to evaluate the transport behavior of a system, a synchronous ensemble of particles [shown in Fig. 1(a)] are randomly initialized in the domain according to the material’s density of states and carrier statistics. Their trajectories are simulated by alternating between free-flight and scattering events according to the intrinsic material scattering rules. The simulation is performed under a driving force, which could be an electric field, or a temperature gradient. Typically, after initialization, a time step $Δt$ is selected, during which a few free-flights followed by scattering events (including self-scattering²) occur. All particles’ trajectories are traced in the channel material. Both left-to-right and right-to-left going particles are simulated, and the transport properties are evaluated from the difference in the fluxes of the opposite-going particles. Customarily, there are two MC methods employed in the literature—the “ensemble” and the “incident-flux” MC methods.^2,35 In the ensemble MC, a synchronous ensemble of particles is simulated and traced simultaneously with periodic boundary conditions for particles that reach the edges of the domain [see Fig. 1(a)]. In the incident-flux (or single-particle) approach, an adequate amount of particles are initialized at the contacts and injected and traced one by one in the domain until they exit from a contact [see Fig. 2(a) further below]. In both cases, the simulation stops when adequate statistics have been collected to have convergence in the conductivity. These techniques have been successfully employed in many different settings, e.g., electronic devices, thermal materials, even nanostructured materials.^24,50,54,55 In this work, we focus on electronic transport under linear response when discussing the MC specifics, however, with appropriate modifications the method we describe could apply to phonon transport as well. The method will specifically be applied here for nanostructured TE materials.

MC challenges encountered in nanostructures: The main difficulty that arises for MC methods for nanostructures is that the presence of nanostructuring features, such as built-in potential barriers, multiple backscattering events due to embedded features of grain boundaries, pores, nanoinclusions, etc. [see Figs. 1(b) and 1(c)], to name a few, reduce the particle flux and our ability to gather adequate statistics. Furthermore, classically, potential barriers limit the number of carriers that participate in transport to the ones that can actually flow over them [see Fig. 1(b)]. Figure 1(d) compares the time taken for simulating $104$ electrons per energy point (injecting them from the left, and ray-tracing them until they exit the domain) vs the material’s channel length. We simulate 100 energy points in total; thus, the overall trajectories simulated are $106$⁠. Four cases are shown: (i) pristine material, (ii) material with grain boundaries forcing reflection of carriers with 50% probability, (iii) a porous material with 40% porosity, and (iv) hierarchical nanostructuring with pores added in addition to grain boundaries, as in a typical situation for thermoelectric materials.⁴¹ As nanostructuring is added in the channel, especially pores in this case, the simulation time required for the same number of particles injected in the channel increases by several orders of magnitude. In Fig. 1(d), the barchart shows the relative time $tsim$ with respect to the pristine material, again indicating orders of magnitude increase in simulation time, promoting the need for more efficient MC algorithms for nanostructures.

In addition, typically TE materials operate under low-field transport conditions, where the driving force is a small potential [see Fig. 1(e)] or small temperature difference at the two ends of the material^56–62 (see  Appendix C for details). This results in only a small variation between the right- and left-going fluxes around the Fermi level, which presents additional difficulties in gathering statistics for the net flux. A larger applied voltage to allow for better statistics [see Fig. 1(f)] can lead to deviation from low-field transport, and other numerical difficulties, such as the need to include inelastic, in addition to elastic scattering processes, as well. In a different case, any carrier energetically positioned below the highest potential energy level (left side) will be reflected back to the right side eventually and only carriers above the source (left side) potential energy will contribute to the flux—this can also lead to difficulties in the treatment of the periodic boundary conditions. This situation becomes increasingly worse for larger applied fields. The electron–phonon scattering processes in TE materials are typically strongly influenced by elastic acoustic phonon scattering processes, and this is what we focus on in this work, although optical/inelastic processes can offer significant design opportunities.^46,47,63,64

Another numerical peculiarity that is encountered in MC simulations for low-energy electrons even in pristine materials, makes the computational difficulty even larger. Noticeably, the low-energy carriers near the band edge have small velocities. Thus, there exists a population of slow moving electrons that end up dominating the computational time (even so in certain cases where the Fermi level is placed higher into the bands).³⁵ To make things worse, in the presence of potential barriers [as shown in Fig. 1(b)], the contribution of those low-energy electrons to transport is insignificant. In addition, at these low energies, the scattering rate of ionized impurity scattering (IIS), an important mechanism especially for TEs, diverges (under low carrier densities).^2,35,65 Thus, this slow moving population of carriers could scatter very strongly and contributes to flux even less and, thus, ultimately requires unnecessarily even more computational resources.

Novel MC method with improved robustness: In this work, we develop a MC method that is more suitable for nanostructured materials, specifically addressing the difficulties described above. We use the incident-flux (single-particle) approach, where the electrons are initialized at the domain boundaries one-by-one and propagate through the domain until they exit at either boundary [propagated to the other side or back-scattered—see Fig. 2(a)].

The first issue we tackle is the large computational time associated with the ray-tracing of low energy electrons. The simulated electrons with their low group velocities end up taking a very large number of free-flight steps intervened by most probably lots of (unnecessary) self-scattering events where the electron does not scatter at the end of the free-flight as is common practice.^2,5 To avoid these redundant occurrences, we consider a mean-free-path (⁠$mfp$⁠) approach, rather than the picking of random free-flight time and the self-scattering approach. We compute the total scattering rate of the particle, and using its band structure velocity, we calculate its mean-free-path. The particle propagates one $mfp$⁠, and then undergoes (enforced) scattering at all times, as depicted in Fig. 2(b). In the case of acoustic phonon scattering in 3D, for example, the mean-free-path is constant in energy [as shown in Fig. 2(b)], and, therefore, electrons from all energies are treated in the same way, with only different free-flight durations. Since in this work we focus on introducing the new method, we only consider elastic, i.e., acoustic phonon scattering process and, thus, a constant $mfp$⁠.

Following the incident-flux method, we use a single-injected flux from the left end of the channel and neglect the injection of flux from the right of the channel. Thus, we refer to this from here on as the “single-flux” method, and we essentially consider only one of the two separate injections of the otherwise “two-flux” method, in which particles are injected from both left and right ends and the flux difference is computed. For these simulations, we consider a two-dimensional domain for numerical simplicity, with two open boundaries at the left/right, whereas the top/bottom boundaries are closed as shown in Fig. 2(a). The simulation procedure is as follows: we initialize and inject electrons in the channel domain only from the left side. We initialize those electrons uniformly in energy, rather than according to their density of states (DOS), as in typical MC methods. We typically use 1000 electrons per energy, with a uniform energy discretization steps of 5 meV. We use these for convenience, since we only consider elastic transport conditions. In the inelastic scattering case, this initialization might need to be performed according to the DOS, but it is beyond the scope of this work (such situations will be addressed in subsequent works). Then, one-by-one, the electrons are ray-traced, by alternating between free-flight events of a $mfp$ distance and intermediate scattering events. Here, scattering only changes the particles’ directions because we consider only elastic processes; thus, their energy stays the same. The upper/lower closed boundaries simply specularly reflect back the electrons in the domain (surface scattering details are not considered and they are out of scope).

To gather the flux statistics, we record the time spend in the domain by those electrons, which propagate all the way from the left to the right end of the domain and exit from there [red line in Fig. 2(a)]. All electrons that are back-scattered to the left do not contribute to the flux and are not considered [yellow line in Fig. 2(a)]. The time an electron spends in the simulation domain until it exits to the right end is referred to as its time-of-flight (⁠$ToF$⁠). The average $ToF$ per particle is then computed as

where $tr(E)$ is the time taken by a particular electron to exit from the right end, and $Nr(E)$ is the number of electrons that make it to the right end. We chose to keep the energy dependence since we only deal with elastic transport conditions. Then, the $⟨ToF(E)⟩$ is used to calculate the flux per simulated electron at each energy as

Using the flux per electron, we can form the overall flux in energy by multiplying by the density of states (DOS), $g(E)$⁠, which essentially is proportional to the transport distribution function (TDF) of the analytical BTE as

This is simply because the product of flux with DOS provides essentially the flow of charge, which is directly related to conductivity in the same way the TDF determines the conductivity.

The proportionality constant $C$ in the equation of the TDF accounts for the super-electron charge that is typically used in MC (the fact that we only simulate a finite number of electrons), and geometrical factors related to the simulation in a finite 2D domain rather than an infinite 3D domain (and connects the conductance to conductivity). Essentially, we use $C$ to map the MC simulated TDF to the TDF that can be derived and extracted from the analytical solution of the BTE, for the case of pristine material alone. We discuss this mapping in detail further below. After having the transport distribution function numerically from MC, and calibrated to the analytical one, we can substitute it in the place of the analytical function $Ξ(E)$⁠, which is given by $τs(E)v(E)2g(E)$⁠, in the usual BTE formulas below.

The electron conductivity is calculated as

and the Seebeck coefficient as

where $q$⁠, $kB$⁠, and $T$ are the electronic charge, Boltzmann constant, and domain temperature (assumed constant at $T=300K$ throughout), respectively. Parameter $Ef$ is the Fermi energy and f is the Fermi–Dirac distribution function. Further, the TE power factor and the electronic thermal conductivity are evaluated, respectively, as

and

The above transport coefficients are computed using MC initially for the pristine material configuration, where the calibration of the constant $C$ takes place. The idea is that once this is calibrated, we can perform MC simulations for complex nanostructured domains without further calibration and benefit from the robustness of the single-flux injection method.

Considering 3D carriers and acoustic phonon scattering limited processes alone (which are elastic), the transport distribution function $Ξ(E)$ calculated analytically from BTE is linear in energy as shown in Figs. 3(a) and 3(b) (left axes, blue lines). In order to form the integrand in Eqs. (4) and (5) above, essentially to account for linear response, we need to multiply the TDF with the differential of the Fermi distribution function with respect to energy (i.e., $∂Ef=∂f∂E$⁠) for the electrical conductivity and with respect to temperature (i.e., $∂Tf=E−EfkBT×∂f∂E$⁠) for the Seebeck coefficient calculations. These products are shown by the red lines in Figs. 3(a) and 3(b) for an arbitrarily chosen Fermi level of $Ef=100$ meV. Using our developed single-flux MC method, we simulate the transport distribution function as well, plotted in Figs. 3(c) and 3(d). The linear trend of the analytical $Ξ(E)$ is also achieved in this case. However, it takes approximately $103$ simulated electrons per energy to gather enough statistics to match the linear trend. $103$ electrons per energy results in an adequate result for the transport coefficients (as we will also show below), although the TDF is still slightly noisy (green lines), while the more or less “noise-less” $105$ electrons case (not shown) does not provide a noticeable conductivity accuracy advantage compared to $104$ (blue lines). For the case of $104$ electrons per energy point, we plot the $Ξ×∂fE$ and $Ξ×∂fT$ products in Figs. 3(c) and 3(d) (red lines) that match the trends of the analytical results as well. In this way, by obtaining the TDF from MC, and multiplying by the differential of the Fermi distribution, we effectively eliminate: (i) the need for two counter-propagating simulation fluxes (now this is captured by the differential functions) and (ii) the need for an application of a driving force, either a voltage difference or a temperature difference, thus avoid the peculiar situation described in Figs. 1(e) and 1(f) where a small enough potential difference window does not provide enough statistics, while a large enough could make the simulation deviate from linear response or from the range of voltages that TE materials utilize (see  Appendix C). Note also that the acquisition of adequate statistics is even more difficult in the case of a temperature gradient for the Seebeck coefficient in common bi-directional flux methods. We do not only need to differentiate between the right- and left-going fluxes but also the ones which flow above and below the Fermi level.

Mapping constant $C$⁠: Here, we provide details on how we obtain the mapping constant $C$ between the simulated MC TDF and the analytical BTE TDF for the pristine material. For this, we calculate the electrical conductivity as a function of the Fermi energy from both approaches (by integrating the TDF times the Fermi derivative in energy) and then find the mapping factor $CEf$ at each Fermi energy, essentially the ratio of the two conductivities. As we can see in the inset of Fig. 4(a), that value is almost constant for all $Ef$⁠, so we take the average of that, i.e., $Cavg$ as the overall final $C$⁠, which gives only one constant factor altogether for the whole range of $Ef$ for mapping the MC conductivity to the analytical one. Since $C$ multiplies the flux, it will be used to obtain the rest of the thermoelectric transport coefficients as a function of $Ef$ as shown in Fig. 4. Clearly, we obtain an excellent match for the electrical conductivity, Seebeck coefficient, and power factor between the MC and analytical BTE as a function of Fermi energy $Ef$⁠, as plotted in Figs. 4(a)–4(c), respectively. The electronic thermal conductivity also matches similarly well (not shown).

We used the overall integrated conductivity to extract the constant $C$⁠, but it is interesting to compare how the transport distribution function extracted from the MC flux method and afterward multiplied by the computed $C$⁠, compares to the analytical TDF from BTE. Figure 4(d) shows this comparison, indicating excellent agreement with the analytical TDF (essentially the black line for the analytical curve resides directly below the blue line for the MC flux TDF). Even the MC lines extracted using a rather small number of simulated carriers, resulting in noisy $Ξ(E)$⁠, also reside around the analytical line. Finally, notice the increasing variations in the flux with energy in MC, more noticeable for the low numbers of simulated electrons per energy. The oscillations in $Ξ(E)$ are stronger at higher energies, but as a matter of fact these oscillations, which are just statistical variations, appear uniformly in energy in the calculations of the $ToF(E)$ and flux $F(E)$⁠. They are of similar relative amplitude compared to the mean values of $ToF(E)$ [or $F(E)$] in energy. For the $Ξ(E)$⁠, however, at low energies: (i) these oscillations have a smaller absolute amplitude since the flux is lower there, (ii) but also they are scaled by the lower density of states, compared to the high energy part of $Ξ(E)$⁠, which further reduces their amplitude. We explain this in more detail in  Appendix A.

The method we present requires only a single-flux injection from one contact into the channel, does not require the application of a driving force, and bypasses the time consumed by the low energy electron ray-tracing process, all of which make it computationally very efficient. It is tested at this point for elastic scattering processes alone and further work will be required to make it compatible with inelastic processes as well. To further stress the error suppression and computational effectiveness of this single-flux method we present, in Figs. 5(a), 5(c), and 5(e), we plot the conductivity, Seebeck coefficient, and power factor with error bars vs the number of simulated electrons. We show results for 10 independent simulation runs for each data point, and the error bars are a measure of the variation within these ten runs. In each simulation run we use ten up to $105$ electrons per energy point in increments of an order of magnitude. As mentioned earlier, we use 100 energy points in total. A single Fermi energy (⁠$Ef=100$ meV) is used in the calculations. In Figs. 5(b), 5(d), and 5(f), we plot the percentage standard deviation (% error) of the corresponding transport coefficients vs the simulated electrons per energy point (error bar sizes). However, we normalize these values to the average value of their corresponding simulation run and present the percentage error. For both the conductivity and the Seebeck coefficient, it appears that $103$ electrons per energy point (i.e., $105$ in total) already significantly reduce the simulation variation to a few percentage points while the transport coefficients converge to their final value. The power factor, since it composes of quantities which are inversely proportional (conductivity and Seebeck), is, in general, a less-fluctuating quantity and requires only $102$ electrons per energy to reach very close to its convergent value. Note that these simulations are executed within minutes (⁠$≈$2–3 min). This makes this method much more computationally efficient compared to common MC works.^24,30,33 For further comparison details, see  Appendix B.

We now demonstrate the effectiveness of the proposed single-flux method for nanostructured materials. We simulate the geometries shown in Fig. 6. First, we populate the two-dimensional domain with noncrystalline grain boundaries, as shown in Fig. 6(a), using the Voronoi algorithm,^66,67 as implemented in Matlab. The algorithm takes seeding points and the domain dimension as inputs to create the grain boundaries. The average size of the grain is calculated by taking the mean of the length of all lines joining all nearest neighboring grain seeds by the Delaunay triangulation algorithm.⁶⁶ During transport, when an electron encounters these grain boundaries, a random decision is made whether the electron will reflect or transmit through, as depicted in Fig. 6(d). Many detailed models exist to describe the transmission and reflection probabilities, which can be momentum, energy, and angle of incidence dependent as well.^25,33,45,68 However, since here our focus is to demonstrate our new algorithm and not the details of boundary scattering, we allow for 50% transmission and 50% reflection probability, independent of the carriers’ energy and direction.

For the second class of nanostructures, we create randomly distributed nanopores in the domain (with pore sizes $5<d<20$ nm) as shown in Fig. 6(b). The average porosity in the domain depends on the number of pores and their sizes. The shape of the pores can be circular, elliptical (oval), with the assumption of specular or diffusive scattering, i.e., with sharp or irregular rough edges in a realistic scenario, respectively [as shown in Fig. 6(d)]. Upon scattering of a particle on a pore, the angle of incidence is computed, and the direction of reflection is determined according to the nature of the pore boundary (diffusive or specular). In this case, for simplifying the scattering process, we use specular reflection at the pores, which means that the angle of incidence is equal to the angle of reflection. We then combine these two features and make the domain hierarchically more complex as shown in Fig. 6(c).

In Fig. 7, we show the electrical conductivity for these structures as a function of Fermi energy $Ef$⁠. The uncertainly related to the error bars shown is a result of the MC statistical variation of simulating the same structure five times. The conductivity decreases when we nanostructure the domain as compared to the pristine material (blue line). There is an average drop of (i) $30%$ in the conductivity of the channel which includes the grain boundary features, (ii) $≈40%$ in the case of nanoporous features, and (iii) a drastic decrease of $≈80%$ when both features are combined. These reduction levels are expected, considering that we have set a 15 nm mean-free-path for the pristine matrix, and the average grain size and distance between pores is of the order of 30 nm. We also find a good agreement in the MC results as compared to Matthiessen’s rule (depicted by the red-dotted line) and computed by adding the resistivity obtained at each step (as a measure of the scattering rates) as

where $σph+GB$ is the conductivity due to phonon and grain boundary scattering and $σph+pores$ is due to phonons and pore scattering.

We now examine the error in the stochastic calculation of transport coefficients vs the required number of simulated electrons for convergent results, as well as the variability due to different simulation runs (here, we have executed five MC runs per data point). As shown in Figs. 8(a), 8(c), and 8(e), the calculated quantities are well converged for $103$ electrons per energy grid point and with almost negligible errors [shown in Figs. 8(b), 8(d), and 8(f)] in the variability of different simulation runs. This is an observation very similar to the pristine channel, indicating that the accuracy of the method does not degrade with nanostructuring geometry complexity. Beyond $103$ electrons per energy grid point the variation per MC run is almost stagnant, so we used the $103$ carrier case for the data in Fig. 7 earlier. However, more carriers can reduce the error more at the expense of more computational resources. A total of $105$ ray-racing trajectories is significantly smaller compared to typical works using standard MC algorithms for such structures, which require millions of trajectories. Note that these simulations with $105$ total trajectories take approximately 20 h on a single CPU to be executed.

Experimental transport studies in nanoporous materials overwhelmingly consider mostly heat transport. We were able to identify a study with measured electronic transport data in nanoporous Si.⁶⁹ The paper compared the mobility of pristine Si with 60% porosity and reported $∼4$ times mobility reduction for n-type samples and between $∼3$–7.5 reduction for p-type samples. Despite the many uncertainties of the experiment, we have simulated a channel with 60% porosity (same as in the experiment) and computed the reduction in electrical conduction compared to the pristine channel to be $∼3.5$⁠, which is similar to what was measured in the experiment.

Finally, we would like to mention a couple of things that were omitted in our implementation, which are typically present in transport codes. We have considered only acoustic phonon scattering, but including more scattering mechanisms, such as by optical phonons, ionized impurities, and defect scattering, would be treated in the same way as in any existing Monte Carlo method,^2,5,48 but a mean-free-path will be derived rather than a time-of-flight. This does not alter the performance advantage of our method, which is based on eliminating the subtraction of two very similar fluxes. We also note that in all of our simulations we used a flat conduction band profile at 0 eV. In reality, nanostructured materials have a varying conduction band profile in the channel, as a result of a varying electrostatic potential due to charging effects caused by the nanoinclusions and the grain boundaries. This is typically captured in most transport simulators by self-consistent coupling of carrier statistics (for the local charge density) with the Poisson equation (for the potential). We have not considered a varying potential in the simulation domain, but this is a separate calculation that can be performed independently prior to the MC ray-tracing part. In that case, as the particles move in the channel, their different attributes, e.g., their velocity, scattering rates, mean-free-paths, etc., would change as dictated by their energy relative to the underlying band profile. With regard to time-dependent simulations, this incident-flux method might not be the most relevant one, to the best of our knowledge.

In summary, we have presented a novel Monte Carlo ray-tracing algorithm with computational efficiency of at least an order of magnitude compared to existing algorithms.^24,25,33 Our new method is a hybrid between the analytical Boltzmann transport and stochastic Monte Carlo and considers many deviations from the common Monte Carlo transport algorithms, which make it specifically efficient for low-field transport in nanostructured materials. For example, in our algorithm, we do not consider the timing griding as it involves several self-scattering events, which especially for the low-energy electrons, it is a computationally expensive process with little influence on the results. Instead, we employ a mean-free-path approach. The method uses significantly less number of ray-tracing electrons compared to standard Monte Carlo algorithms for similar problems, avoids the statistically challenging subtraction of two opposite going fluxes, and avoids the application of driving external forces all together. We demonstrated the algorithm’s efficiency and strength for accurate simulations in large nanostructured domains with multiple defects. The method is convenient for studying electronic transport in highly complex nanostructured materials, especially in the field of thermoelectrics. It is tested for elastic transport and scattering conditions but can be extended to include inelastic processes, as well as go beyond electronic, to phonon transport (see  Appendix D for a summary of how this can be implemented). We believe the new method provides an efficient and user friendly algorithm, which will enable the proper study of highly nanostructured materials under low-field steady-state conditions.

This work received funding from the European Union’s Horizon 2020 research and innovation program under Grant Agreement No. 863222 (UncorrelaTEd).

The authors have no conflicts to disclose.

Pankaj Priyadarshi: Conceptualization (equal); Data curation (equal); Methodology (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Neophytos Neophytou: Conceptualization (equal); Investigation (equal); Methodology (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

As depicted in the main text, the transport distribution function $Ξ(E)$ calculated using the developed Monte Carlo ray-tracing algorithm has more fluctuations at higher energies. In order to understand this, we plot in Figs. 9(a) and 9(b) the energy resolved time-of-flight $ToF(E)$ and its inverse, the flux $F(E)$ for the cases of $103$ (blue lines) and $104$ (red lines) simulated electrons per energy grid point. The $ToF(E)$ follows a decreasing trend due to the higher velocities of electrons at higher energies, while $F(E)$ increases. Note that we find that the average number of electrons exited to the right is almost constant with energy in our elastic simulations. We observe statistical fluctuations around the local mean value of the $ToF(E)$—note that we plot this in the logarithmic scale. The inset plot in Fig. 9(a) shows the percentage variation of the simulated values with respect to the mean value (i.e., percentage variation of the ratio of local standard deviation to the local $ToF$ mean value). That percentage is low, but importantly is also rather constant in energy (as expected). The flux $F(E)$ in Fig. 9(b) is plotted in linear scale, as it will enter the $Ξ(E)$ evaluation. The inverse of high $ToF(E)$ fluctuations at low energies appears now reduced compared to those at higher energies, although as we showed above, they are relatively of the same magnitude compared to their local mean value in energy. When it comes to form the $Ξ(E)$⁠, and we multiply by the DOS, $g(E)$⁠, the low-energy fluctuations are further scaled due to the lower DOS at lower energies, but magnified at higher energies due to the larger DOS, as depicted in Fig. 9(c). For an increased number of electrons, the fluctuations become smaller compared to their average local value in energy as shown by the red lines. Thus, it appears that at higher energies, we have the presence of more noise.

In our Monte Carlo method, for a $1000×500nm2$ two-dimensional channel size, under extreme nanostructuring, the transport results converge when using approximately $103$ ray-tracing particles per energy. In the case of pristine material, convergent results with reduced error can be achieved for even less particles (see Fig. 8). By using 100 energy points, adequate results are, thus, obtained with $105$ trajectories in total in the nanostructured channels. The simulations we present in the paper, the length of channel and the excessive nanostructuring, are reasons why only a small percentage of particles make it from the left to the right and contribute in building the transport distribution function. In fact, only 2%–3% make it from the left to the right, while the majority is back-scattered to the left contact; thus, a large portion of the particles do not trace the full length of the channel either, with a reduced computational cost.

To estimate the computational benefits of this single-flux method compared to other similar Monte Carlo works, we provide below some information about the number of ray-tracing trajectories typically employed. We are not aware of Monte Carlo electronic transport works in the nanostructures we consider and the domain sizes we consider, thus, we compare first with our own prior “two-flux” method in phonon transport Monte Carlo works (the ray-tracing part of the codes can still be compared). In Ref. 33, on 2D phonon transport in similar nanostructured domains, $2.5×106$ particles were required under the “two-flux” incident flux method. Our work by Wolf et al.,^24,25 used more than $106$ ray-tracing particles for a 3D $1000nm$ channel length again using the “two-flux MC” incident-flux method. These are typical numbers for many MC (phonon transport) works on micrometer length scale domains, which are more than an order of magnitude larger compared to the single-flux method we present here.

The situation is different if the channel is shorter, where more particles cross the channel and, thus, less overall need to be eventually simulated. In Fig. 10(a), we show the transport distribution function and conductivity from simulations of a pristine material of a $≈200nm$ size channel length. Although the TDF still has some noise when we use only $102$ particles per energy, the error in the conductivity [see blue line in Fig. 10(b)] is small compared to the $103$ particle case (red line), while the $5×102$ case (orange line) shows only minimal error. Indicatively, these simulations for the pristine channel either for the $200$ or the $1000nm$ channel cases are concluded within seconds to minutes.

In common MC methods, a driving force is applied across the contacts of the simulation domain, either a voltage $ΔV$ or temperature difference $ΔT$⁠. In the case of TE materials, these driving forces are small, dictating linear response within Boltzmann transport. In a typical single TE leg simulation analysis,^58,59 an n-type leg of length 10 mm develops roughly a 10 mV open circuit voltage which corresponds to 1 $μ$V per $μ$m voltage drop for the channels we simulate. We also mention here a work on a Bi$2$Te$3$ TE device,⁶⁰ for which the typical open circuit voltage is measured to be around 2.0 mV.^61,62 Such small voltages are very hard to be simulated numerically using two-flux incident-flux or ensemble MC methods, in which the bi-directional flux difference needs to be statistically accurately resolved. For an indication of the level of this difficulty, we simulate fluxes from the left and the right of a $1000nm$ channel independently. For this, we place the Fermi level of the right contact lower compared to that of the left contact by $ΔV$⁠, essentially reducing the inflow from the right. We then obtain the conductivity by subtracting the left/right-going fluxes that were injected from the two contacts. Note that this is just an indicative scenario, since we do not apply the potential drop in the channel. For example, as $ΔV$ increases, we should resume to the single-flux method. We compare the results with our developed single-flux method. We plot in Fig. 11(a), the conductivity as a function of the number of simulated electrons for the single-flux method (red line) and for the two-flux method under different $ΔV$⁠. We also include the stochastic error bar after executing 10 simulations per data point. In Fig. 11(b), we plot the actual error for the different cases with $ΔV$ as indicated in the legend of the figure. The single-flux method (red line) indicates the lowest error, which gets reduced linearly as the number of simulated particles increases. The case that uses $ΔV=5mV$ has the largest error, more than an order of magnitude compared to the single-flux line at the same number of simulated particles. Another way to see this (horizontally in the figure) is that 2 orders of magnitude more particles need to be simulated to reach a similar convergence error compared to the single-flux method. We need to apply a larger potential difference, e.g., $ΔV=100mV$ as it can be seen from Fig. 11(b) to reduce the error at a certain electron count number, which takes us closer to the single-flux method. However, any of these values used here are much larger compared to what thermoelectrics experience in reality.

Here, we provide some details on how the algorithm we present can be materialized to treat phonon transport.^7,8 We start from the expression of the temperature derivative of the Bose–Einstein (BE) distribution (instead of the Fermi–Dirac distribution), which will provide the driving force,

By expressing the specific heat in terms of the BE distribution derivative as well, we reach

A general expression for the thermal conductivity under Boltzmann transport assumptions is given by

Thus, by substituting for the specific heat, we get

Note that this expression is very similar to the one we have for the electronic conductivity, with the addition of the phonon energy term $ℏω$ and the temperature derivative of the BE distribution, rather than the energy derivative of the Fermi–Dirac distribution. Also, note that the common Callaway model can be derived from Eq. (D3) as well, following:

Above we use the phonon (acoustic) density of states as $g(ω)=ω2vs2$⁠. Now, we consider that the scattering rate for the most severe intrinsic phonon–phonon scattering event, the three-phonon Umklapp scattering, typically follows $τ−1∼ω2$⁠. We also consider that the phonon velocity (for acoustic phonons) is constant with $ω$⁠. So, the frequency dependence of the phonon TDF then becomes constant with $ω$ as $Ξph(ω)=τ(ω)vs2(ω)g(ω)∼C$⁠. The flux from a typical ray-tracing for a given group of particles reflects the product of mean-free-path (⁠$mfp$⁠) and the velocity $∼[τ(ω)vs(ω)]vs(ω)∼C/ω2$ (note the $mfp$ is not constant in this case with $ω$⁠). This is the flux simulated in the case of phonons, having a different energy dependence compared to the electronic case (⁠$∼E$⁠). Multiplying this by $g(ω)(ℏω)$ this time, since the thermal conductivity is the flow of energy, rather than by only $g(.)$ as in the case of electrons, we end up with again a linear dependence, $∼ω$⁠, for the argument of thermal conductivity integral. Then, the constant $C$ is used to map the analytical BTE model to the ray-tracing algorithm, in the same way as described for electrons in the main text. As in the case of electrons, all intrinsic material transport behavior is lumped into the simulated flux and $C$⁠. From there on, we can focus on the effect of geometrical scattering on phonons. Thus, the formalism is very similar to the case of electrons.