Inquiry into the climate response to external forcing perturbations has been the central interest of climate dynamics. But the understanding of two important aspects of climate change response—nonlinearity and regionality—is still lacking. Here a Green’s function approach is developed to estimate the linear response functions (LRFs) for both the linear and quadratic nonlinear response to ocean thermal forcing in a climate model, whereby the most excitable temperature modes, aka the neutral modes, can be identified for the current Earth climate. The resultant leading mode of the nonlinear response is characterized by a polar-amplified global cooling pattern, unveiling an intrinsic inclination of the modern climate towards cooling. Moreover, optimal forcing patterns are identified to most efficiently excite the corresponding neutral mode patterns. The forcing-response framework developed herein can be utilized to determine the optimal forcing patterns to inform solar geoengineering experiments and to interpret regional climate response and feedback in general.