\textbf{How the self-interaction mechanism affects zonal flow drive and convergence of turbulent transport simulations with system size}

We use gyrokinetic flux-tube simulations to report a decrease in the shearing rate of ExB zonal flows with increasing system size measured by 1/$\rho $*$=$a/$\rho_{\mathrm{i}}$, where a is the tokamak minor radius and $\rho_{\mathrm{i}}$ is the ion Larmor radius. This is done in practice by decreasing k$_{\mathrm{y,min}}\rho_{\mathrm{i\thinspace }}$(\textasciitilde $\rho $*), where k$_{\mathrm{y,min}}$is the minimum wavenumber along the direction y, bi-normal to the magnetic field. The corresponding gyro-Bohm normalised heat and particle fluxes also increase with decreasing k$_{\mathrm{y,min}}$. We find that this results from the non-adiabatic passing electron dynamics and corresponding fine structures at mode rational surfaces associated to each k$_{\mathrm{y\thinspace \thinspace }}$mode. The related strong self-interaction mechanism disrupts resonant 3-wave interactions involving the zonal modes. As a consequence, the different k$_{\mathrm{y}}$ contributions to Reynolds Stress driving the zonal flows tend to get decorrelated, which results in the shearing rate level developing a statistical dependence on k$_{\mathrm{y,min}}$. In adiabatic electron simulations, the scaling is not as severe, owing to a weaker self-interaction mechanism at play.