In this work, we present a comprehensive Kinetic Monte Carlo (KMC) modeling based statistical framework to evaluate the device-to-device variation of thin-film HfO2 ferroelectric FET (FeFET). We conclude that the closing of the memory window in a FeFET array with device scaling can be attributed to: 1) limited number of domains; 2) variation among domains; 3) intrinsic stochasticity of individual domain switching. To enable further scaling of FeFET, co-optimization approaches from material, process, and device operation to control variation are proposed: i) increase the number of domains through material/process optimization (e.g. decrease of deposition temperature, etc.); ii) improve the uniformity of domains (e.g. minimizing the domain size variation and defect distribution, etc.); iii) increase the pulse amplitude/width to ensure deterministic switching of individual domains.