Future sixth-generation (6G) low Earth orbit (LEO) satellite systems will utilize large-scale multi-beam phased arrays operating at millimeter-wave (mmWave) frequencies, with embedded power amplifiers (PAs) at each antenna element. However, the severe path loss drives PAs near saturation, and dynamic beamforming leads to nonuniform, time-varying PA conditions, jointly resulting in complex memory nonlinearities that degrade system performance. Such in-band distortions require mitigation by PA linearization, which remains highly constrained in LEO satellite systems. Therefore, accurate modeling and quantification of these nonlinearities are critical for realistic performance evaluation and to guide the development of effective linearization strategies. Current analytical methods focus on single-PA or memoryless-array models, without addressing the frequency-dependent statistical characteristics and spatially varying radiation patterns of nonlinear distortions in multi-beam phased arrays. In this paper, we develop statistical models for the PA output and the received signals in multi-beam phased arrays, explicitly decomposing the distortion structure and deriving closed-form expressions for each distortion component’s statistical properties. Furthermore, we analyze the spatial beam-pattern distortions and derive a closed-form expression for the signal-to-distortion-plus-interference-and-noise ratio (SDINR) under nonlinear distortion. Extensive simulations, conducted measurements, and over-the-air (OTA) experiments validate that the proposed model reduces optimal input back-off estimation error by 2 dB and improves SDINR prediction accuracy by approximately 3 dB. This research establishes a robust analytical foundation for future LEO satellite system design and linearization technology development.