We investigate stable recovery guarantees for phase retrieval under two realistic and challenging noise models: the Poisson model and the heavy-tailed model. Our analysis covers both nonconvex least squares (NCVX-LS) and convex least squares (CVX-LS) estimators. For the Poisson model, we demonstrate that in the high-energy regime where the true signal $pmb{x}$ exceeds a certain energy threshold, both estimators achieve a signal-independent, minimax optimal error rate $\mathcal{O}(\sqrt{\frac{n}{m}})$, with $n$ denoting the signal dimension and $m$ the number of sampling vectors. In contrast, in the low-energy regime, the NCVX-LS estimator attains an error rate of $\mathcal{O}(\|\pmb{x}\|^{1/4}_2\cdot(\frac{n}{m})^{1/4})$, which decreases as the energy of signal $\pmb{x}$ diminishes and remains nearly optimal with respect to the oversampling ratio. This demonstrates a signal-energy-adaptive behavior in the Poisson setting. For the heavy-tailed model with noise having a finite $q$-th moment ($q>2$), both estimators attain the minimax optimal error rate $\mathcal{O}( \frac{\| ξ\|_{L_q}}{\| \pmb{x} \|_2} \cdot \sqrt{\frac{n}{m}} )$ in the high-energy regime, while the NCVX-LS estimator further achieves the minimax optimal rate $\mathcal{O}( \sqrt{\|ξ\|_{L_q}}\cdot (\frac{n}{m})^{1/4} )$ in the low-energy regime. Our analysis builds on two key ideas: the use of multiplier inequalities to handle noise that may exhibit dependence on the sampling vectors, and a novel interpretation of Poisson noise as sub-exponential in the high-energy regime yet heavy-tailed in the low-energy regime. These insights form the foundation of a unified analytical framework, which we further apply to a range of related problems, including sparse phase retrieval, low-rank PSD matrix recovery, and random blind deconvolution.