\section{Details of the calculations for Si}

First, we compute the self-consistent potential and Kohn-Sham states on a
$12\times12\times12$ Monkhorst-Pack {\bf{k}}-point grid using DFT and lattice-dynamical properties with
DFPT\cite{baroni2001phonons} on a $3\times3\times3$ {\bf{q}}-point grid,
as implemented in the Quantum Espresso distribution\cite{giannozzi2009quantum}
using the Perdew-Burke-Ernzerhoff exchange-correlation functional\cite{perdew1996generalized}.
We used a full-relativistic norm-conserving
optimized Vanderbilt pseudopotential\cite{scherpelz2016implementation}.
The unit cell consists of Si in the diamond structure with
an experimental lattice parameter of $5.43~\si{\angstrom}$.
The {\it{e-ph}} matrix elements are first computed
on coarse grids, then they are determined in the significantly finer grids using
both Wannier-Fourier (W-F) interpolation only, through EPW code
and our dual interpolation method, Turbo-EPW (T-EPW).
Maximally localized Wannier functions\cite{marzari1997maximally} for the
wannierization procedure are obtained
from Wannier90.\cite{mostofi2008wannier90} Thus, Bloch-to-Wannier
rotation matrices and then Wannier-to-Bloch diagonalizer matrices are used
to interpolate {\it{el-ph}} matrix elements.

\section{Convergence analysis}

\begin{figure}
\centering
\includegraphics[width=0.48\textwidth]{figS1.png}
\vspace*{6mm}\caption{\label{figureS1}
Convergence analysis of the scattering rate of Si at 300K
due to el-ph coupling in the Fan-Migdal approximation as a
function of the electron energy. The calculation has been
performed on different {\bf{q}} meshes, namely, $(20)^3$ (green crosses),
$(60)^3$ (blue squares), $(100)^3$ (red dots) {\bf{q}}-meshes in the
first W-F interpolation, while keeping fixed the number of irreducible
{\bf{k}} points in the second plane-waves interpolation,
leading to $(100)^3$ {\bf{k}}-mesh.}
\end{figure}

Fig.~\ref{figureS1} shows how the scattering rate aproaches
convergence by increasing the mesh size of W-F interpolation,
from $(20)^3$ to $(100)^{3}$ ${\bf{q}}$ points,
while keeping the number of irreducible points
fixed at $5216$ ${\bf{k}}$ points.
The second interpolation by star functions leads
into a converged grid with $(100)^3$ ${\bf{k}}$ points.
Fig.~\ref{figureS2} shows the difference in
scattering rates calculated by different approaches, namely,
different parameters in the second plane-waves interpolation (PWI),
over equivalent meshes.
The analysis shows that the accompanying error due to the increase in $M$
decreases by enlarging the number of irreducible points, $n_{{\bf{\bar{k}}}}$ (see main text).
Moreover, the difference between scattering rates calculated over $(100)^3$ ${\bf{q}}$/ $1661$ ${\bf{k}}$
points expanded by using $M = 5$ to reach $(100)^3$ ${\bf{k}}$-mesh and
data from $(100)^3$ ${\bf{q}}$/ $256$ ${\bf{k}}$ points expanded by using $M = 30$
to reach $(100)^3$ ${\bf{k}}$-mesh,
is within $\approx \pm 6\%$.
For the remaining, the difference between the data points is about
the same, within $\approx \pm 3\%$.

\begin{figure}
\centering
\includegraphics[width=0.48\textwidth]{figS2.png}
\vspace*{6mm}\caption{\label{figureS2}
Difference between scattering rates computed by using
different parameters in the second plane-waves interpolation,
over equivalent meshes. Green crosses correspond to the difference
between scattering rates calculated over $(100)^3$ {\bf{q}}/ $1661$ {\bf{k}}
points expanded by using $M = 5$ to reach $(100)^3$ {\bf{k}}-mesh and
data from $(100)^3$ {\bf{q}}/ $256$ {\bf{k}} points expanded by using $M = 30$
to reach $(100)^3$ {\bf{k}}-mesh. Blue squares correspond to the
difference between scattering rates calculated over $(100)^3$ {\bf{q}}/
$5216$ {\bf{k}} points expanded by using $M = 10$ to reach $(180)^3$ {\bf{k}}-mesh
and data from $(100)^3$ {\bf{q}}/ $1661$ {\bf{k}} points expanded by
using $M = 30$ to reach $(180)^3$ {\bf{k}}-mesh. Red dots
correspond to the difference between scattering rates calculated over
$(100)^3$ {\bf{q}}/ $5216$ {\bf{k}} points expanded by using $M = 8$
to reach $(160)^3$ {\bf{k}}-mesh and data from $(100)^3$
{\bf{q}}/ $1661$ {\bf{k}} points expanded by using $M = 25$ to reach
$(160)^3$ {\bf{k}}-mesh.}
\end{figure}

\section{Calculation of electron self-energy and thermoelectric properties}

The expression for the imaginary part of electronic self-energy due to
{\it{el-ph}} coupling in the
Fan-Migdal approximation can be derived from quantum field theory\cite{giustino2017electron}
and it is expressed as
\begin{equation}
\label{imag}
 \begin{split}
\Sigma^{\prime\prime}_{n,{\bf{k}}}(\omega,T)=\pi \sum_{m,\beta} \int_{BZ} \frac{d{\bf{q}}}{\Omega_{BZ}} |g_{mn,\beta}({\bf{k},{\bf{q}}})|^2 \\
       \times \Bigg[\left[n_{{\bf{q}}\beta}(T) + f_{m{\bf{k}}+{\bf{q}}}\right]\delta(\omega - (\epsilon_{m{\bf{k}}+{\bf{q}}} - \epsilon_F) + \omega_{{\bf{q}}\beta}) \\
        + [n_{{\bf{q}}\beta}(T) + 1 - f_{m{\bf{k}}+{\bf{q}}}]\delta(\omega - (\epsilon_{m{\bf{k}}+{\bf{q}}} - \epsilon_F) - \omega_{{\bf{q}}\beta})\Bigg]~,
\end{split}
\end{equation}
where $n_{{\bf{q}}\beta}(T)$ and $f_{m{\bf{k}}+{\bf{q}}}$ are the Bose-Einstein and the Fermi-Dirac distributions,
$\Omega_{BZ}$ is the BZ volume, $m$ and $n$ are the corresponding electronic states, while $\beta$ represents the
phonon branch, $\epsilon_{m{\bf{k}}+{\bf{q}}}$ are the electronic eigenenergies of the state $m{\bf{k}}+{\bf{q}}$ and
$\omega_{{\bf{q}}\beta}$ are the corresponding eigenfrequencies with wavevector ${\bf{q}}$ and phonon branch $\beta$.
Basically, from the first W-F interpolation we can get
$f({\bf{\bar{k}}}) = \Sigma^{\prime\prime}_{n,{\bf{\bar{k}}}}(\omega,T)$,
over the irreducible points,
which will be interpolated throughout the whole BZ by star functions,
resulting in $\Sigma^{\prime\prime}_{n,{\bf{k^{\prime}}}}(\omega,T)$ over denser
${\bf{k^{\prime}}}$ grids.

$\Sigma^{\prime\prime}$ is directly related to the scattering rate, that is, inversely
proportional to the relaxation time
\begin{equation}
\label{tau}
\frac{1}{\tau_{n,{\bf{k}}}} = 2\Sigma^{\prime\prime}_{n,{\bf{k}}}(\omega = 0,T)~,
\end{equation}
which enters in kinetic transport equations.
Indeed, the kinetic coefficient tensors can be expressed through
\begin{equation}
\label{Lambda}
\Lambda^{(\alpha)}(\mu;T)= e{^2}\int\Xi(\epsilon,\mu,T)(\epsilon - \mu)^{\alpha}\left(-\frac{\partial f^{(0)}(\mu;\epsilon,T)}{\partial \epsilon}\right)d\epsilon~,
\end{equation}
where $\mu$ is the chemical potential and $\Xi(\epsilon,\mu,T)$
is the transport distribution kernel given by
\begin{equation}
\label{Xi}
\Xi(\epsilon,\mu,T) = \int \sum_n{{\bf{v}}_{n,{\bf{k}}}\otimes{\bf{v}}_{n,{\bf{k}}}{\tau}_{n,{\bf{k}}}}(\mu,T)\delta(\epsilon - \epsilon_{n,{\bf{k}}})\frac{d{\bf{k}}}{8\pi^3}~,
\end{equation}
with ${\bf{v}}_{n,{\bf{k}}}$ being the electron velocity.
From both experimental conditions of zero temperature gradient ($\nabla T = 0$)
and zero electric current, the kinetic coefficient tensors
can be identified with the electrical conductivity tensor, $\sigma = \Lambda^{(0)}$,
the Seebeck coefficient tensor, $S = (eT)^{-1}\Lambda^{(1)}/\Lambda^{(0)}$,
and the charge carrier contribution to thermal conductivity tensor,
$\kappa_e = (e^2T)^{-1} \left({\Lambda^{(1)}\cdot{\Lambda^{(0)}}^{-1}}\cdot{\Lambda^{(1)}} - \Lambda^{(2)}\right)$.
Consequently, once we have $\tau_{n,{\bf{k}}}$ from Eq.~(\ref{tau}),
we can compute the electrical conductivity, the Seebeck coefficient, Lorenz function and
the charge carrier contribution to thermal conductivity from the solution of Eq.~(\ref{Xi}) and~Eq.~(\ref{Lambda}).
Note that the both bandstructure and phonon dispersion have also been interpolated by the method
presented in the main text over the same grid as electron self-energy.
Indeed, we have implemented these equations on top of the BoltzTraP code,\cite{madsen2006boltztrap}
from which we can obtain these transport properties
directly from first principles.

\section{Scattering by ionized impurities}

For the calculation of thermoelectric transport properties of
$n$- and $p$-type Si polycrystals, we have also considered the scattering
by ionized impurities.
Such scattering has been treated theoretically by
Brooks and Herring (B-H)\cite{brooks1955theory,chattopadhyay1981electron}
by considering a screened Coulomb potential, the
Born approximation for the evaluation of transition probabilities
and neglecting perturbation effects of the
impurities on the electron energy levels and wave functions.
In the B-H theory the electron is scattered independently
by dilute concentrations of ionized centers randomly distributed in the semiconductor.

The per-unit-time transition probability for the scattering of charge carriers
by ionized impurities can be given in the plane-wave approximation as
\begin{multline}
\label{U}
W({\bf{k}}|{\bf{{k}^{\prime}}}) = \frac{2\pi}{\hbar}\frac{N_i}{V} \\
\left\vert{\int U({\bf{r}}) \exp{\left[i({\bf{k}}-{\bf{{k}^{\prime}}};{\bf{r}})\right]}d{\bf{r}}}\right\vert^2 \delta(\epsilon({\bf{{k}^{\prime}}})-\epsilon({\bf{k}}))~,
\end{multline}
where $U({\bf{r}})$ is the scattering potential and $N_i$ is the ionized impurity concentration.

A long-range Coulomb field,
$U(\mat{r}) = e\phi(r) = \pm e^2/\zeta r$, with potential $\phi$ at a point $r$ of the crystal
is created by the presence of positive (donor) or negative (acceptor)
impurity ions, within a medium with dieletric constant $\zeta$.
The straightforward application of this field in Eq.~(\ref{U}) leads
into a logarithmic divergence, and hence, a screened Coulomb potential has to be considered.
From the B-H theory the potential can be expressed in a more rigorous form as
$\phi(r) = \pm e/{\zeta r} \left(\exp{\left(-{r}/{r_0}\right)}\right)$,
where $r_0$ is the radius of ion field screening defined by
\begin{equation}
\label{r0}
r{_0}^{-2}({\bf{k}}) = \frac{4\pi{e^2}}{\zeta{_0}}\int{-\frac{\partial f_0}{\partial \epsilon({\bf{k}})}}g(\epsilon)d\epsilon~,
\end{equation}
where $f^{(0)}(\epsilon)$ is the equilibrium electron
distribution function, $\zeta_0$ is the static dielectric constant,
and $g(\epsilon)$ is the density of states,
calculated numerically on an energy grid
with spacing $d\epsilon$ sampled over $N_k$ ${\bf{k}}$-points
\begin{equation}
\label{dos2}
g(\epsilon) = \int \sum_n \delta(\epsilon - \epsilon_{n,{\bf{k}}}) \frac{d{\bf{k}}}{8\pi^3} = \frac{1}{\Omega N_k}\sum_{n,{\bf{k}}} \frac{\delta(\epsilon - \epsilon_{n,{\bf{k}}})}{d\epsilon}~,
\end{equation}
where $\Omega$ is the volume of the unit cell.

Within the relaxation time approximation for the
Boltzmann transport equations, the relaxation time for the scattering of the
charge carriers by ionized impurities can be expressed as
\begin{equation}
\label{tau22}
\tau_{imp}({\bf{k}}) = \frac{\hbar\zeta{_0}{^2}}{2{\pi}{e^4}{N_i}F_{imp}({\bf{k}})}{\bf{k}}^2 \left\vert\frac{\partial \epsilon({\bf{k}})}{\partial {\bf{k}}}\right\vert
\end{equation}
where
\begin{equation}
F_{imp}({\bf{k}}) = ln(1+\eta) - \frac{\eta}{1+\eta}~,
\end{equation}
is the screening function with $\eta = (2{\bf{k}}r_0)^2$.
Here, we interpolated $\left\vert\frac{\partial \epsilon({\bf{k}})}{\partial {\bf{k}}}\right\vert$
within Eq.~(\ref{tau22}) by using the plane-waves interpolation (see the main text) and
used Mathiessen's rule to consider both the scattering by ionized impurities and phonons.

\section{Calculation of phonon-assisted optical absorption}

To calculate the phonon-assisted absorption coefficient,
we use the Fermi's golden rule expression\cite{parravicini1975g,noffsinger2012}:
\begin{eqnarray}
\label{imag}
\alpha(\omega) = 2 \frac{4\pi^2 e^2}{\omega c n_r(\omega)}\frac{1}{\Omega} \frac{1}{N_{{\bf{k}}} N_{{\bf{q}}}} \sum_{\beta i j {\bf{k}} {\bf{q}}}|{\mat{\lambda}} \cdot ({\bf{S_1}} + {\bf{S_2}})|^2 \\
\times P \delta(\epsilon_{j,{\bf{k}}+{\bf{q}}} - \epsilon_{i,{\bf{k}}} - \hbar\omega - \pm \hbar\omega_{\beta {\bf{q}}}) \nonumber ~,
\end{eqnarray}

with

\begin{equation}
{\bf{S_1}}({\bf{k}} {\bf{q}}) = \sum_m \frac{{\bf{v}}_{im}({\bf{k}}) g_{m j,\beta}({\bf{k}},{\bf{q}})}{\epsilon_{m,{\bf{k}}} - \epsilon_{i {\bf{k}}} - \hbar\omega + i\Gamma_{m,{\bf{k}}}}~,
\end{equation}

\begin{equation}
{\bf{S_2}}({\bf{k}} {\bf{q}}) = \sum_m \frac{g_{i m,\beta}({\bf{k}},{\bf{q}}) {\bf{v}}_{mj}({\bf{k}} + {\bf{q}})}{\epsilon_{m,{\bf{k}} + {\bf{q}}} - \epsilon_{i {\bf{k}}} \pm \hbar\omega_{\beta {\bf{q}}} + i\Gamma_{m,{\bf{k}}+{\bf{q}}}} ~,
\end{equation}

and

\begin{equation}
P = (n_{\beta {\bf{q}}} + \frac{1}{2} \pm \frac{1}{2})(f_{i {\bf{k}}} - f_{j, {\bf{k}}+{\bf{q}}})~.
\end{equation}

In these equations, $\omega$ is the photon frequency,
$c$ is the speed of light, $n_r$ is the refractive index
(for silicon we used $n_r = 3.4$), and $\mat{\lambda}$ is
the photon polarization. ${\bf{S_1}}$ and ${\bf{S_2}}$ are
the two possible ways for the indirect absorption process, while
$P$ is related to the carrier and phonon statistics.
For calculations of Si, the DFT band gap and all
conduction bands have been shifted up by $0.7$\si{\eV}
to simulate experimental gap.
Our calculations have been carried out over $(60)^3${\bf{q}}/$(40)^3${\bf{k}}
meshes. It took $\approx 36$ minutes to perform the calculation by using
our dual interpolation method on $8$
CPU cores at the Odyssey cluster (Harvard University).

\end{appendices}

\setcounter{secnumdepth}{0}
\section{Acknowledgments}
The authors thank the Harvard FAS Research Computing facility and the Brazilian
CCJDR-IFGW-UNICAMP for computational resources.
A.S.C. and A.A. gratefully acknowledge financial support from the
Brazilian agency FAPESP under Grants No.2015/26434-2,
No.2016/23891-6, No.2017/26105-4, No.2018/01274-0 and No.2019/26088-8.
A.S.C. also acknowledges the kind hospitality of SEAS-Harvard University.

\bibliographystyle{apsrev4-1}