## 2 Theoretical Model
The propagation of the wave packet is obtained by applying the time-evolution operator  
\[\Psi(x,y,t+\Delta t)=e^{-\frac{i}{k}H\Delta t}\Psi(x,y,t), \tag{3}\]  
where the Hamiltonian \(H\), assumed to be time-independent, is the one for low energy electrons in graphene in the presence of a potential landscape \(V(x,y)\), [39; 40]  
\[H=v_{F}\left(\vec{\sigma}\cdot\vec{p}\right)+V(x,y)\mathbf{I}, \tag{4}\]  
with \(\mathbf{I}\) as the \(2\times 2\) identity matrix and wave functions written as pseudo-spinors \(\Psi=(\Psi_{A}\ \Psi_{B})^{T}\), so that \(|\Psi_{A}|^{2}\) (\(|\Psi_{B}|^{2}\)) is the probability of finding the electron in the sub-lattice \(A\) (\(B\)) of graphene. We separate the potential and kinetic energy terms of the time-evolution operator through the split-operator technique [55; 56],  
\[e^{-\frac{i}{k}H\Delta t}=e^{-\frac{i}{k}V(x,y)\Delta t}e^{-\frac{i}{k}v_{F} (\vec{p}\cdot\vec{p})\Delta t}e^{-\frac{i}{2k}V(x,y)\Delta t}, \tag{5}\]  
where terms of order higher than \(O(\Delta t^{3})\) are neglected as an approximation, by taking a small time step \(\Delta t\). The advantage of this factorization is that it allows us to perform multiplications in real and reciprocal spaces separately, thereby avoiding writing the momentum operator as a derivative, just by making a Fourier transform of the functions and writing \(\vec{p}=\hbar\vec{k}\). Furthermore, the exponential of Pauli matrices can be re-written exactly as matrices, which simplifies the calculations even more [46; 48]. We performed wave-packet propagation with a time step as small as \(\Delta t=0.1\) fs and kept track of the probabilities of finding the electron before (\(P_{1}\)), inside (\(P_{2}\)), and after (\(P_{3}\)) the scattering region. The large-\(t\) limit of the probability \(P_{3}\) was then identified as the transmission probability through the perturbed region. Notice that performing a Fourier transform while solving Equation (5) is similar to applying periodic boundary conditions to the system. Thus, in order to avoid spurious reduction (enhancement) of the transmission (reflection) probability \(P_{3}\) (\(P_{1}\)) when the wave packet reaches the upper edge of the numerical grid, we considered a very long grid in the propagation direction. Specifically, the graphene ribbon contained 1024 \(\times\) 65,536 grid points regularly spaced by 0.125 nm in both \(x\) and \(y\) directions. Fourier components of the wave packet were computed on discrete wave-vector nodes regularly distributed in reciprocal space with 0.049 nm\({}^{-1}\) and 0.00077 nm\({}^{-1}\) step sizes in the transverse and longitudinal directions, respectively. Nitrogen doping of graphene usually leads to randomly distributed substitutional sites [49; 57; 58; 59]. Accordingly, doping was simulated by a random distribution of scattering centers within an area of \(A=128\) nm \(\times\) 40 nm around the \(y=0\) axis of the ribbon, as sketched in Figure 1. The N atoms, which transfer 0.35-0.45 electron to the graphene host [10; 60], are positively ionized. As revealed by _ab-initio_ calculations [30; 61], the perturbation produced on the \(\pi\) electrons extends somewhat around each impurity, due to Coulomb attraction. Each dopant atom is thus represented by a potential well of depth \(|U|\) and width \(\sigma\)[61]:  
\[v(\rho)=-|U|\exp\left(-\frac{\rho^{2}}{2b^{2}}\right), \tag{6}\]  
where \(\rho\) is the distance to the atom. The parameters to be used for well-separated nitrogen atoms are \(b=0.15\) nm and \(U=-4.0\) eV [19; 61]. In order to avoid any overlapping of the potential wells generated by the dopants, the selection of the scattering centers was constrained by the requirement that the distance between two nitrogen centers remains greater than 2 nm (for dopant concentrations 0.1 at.%-0.4 at.%) or 1.8 nm (for dopant concentration 0.5 at.%).