Neutral modes of surface temperature and the optimal ocean thermal forcing for global cooling

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.


INTRODUCTION
Equilibrium climate sensitivity concerning how much the global mean surface temperature increases in response to a doubling of the concentration of carbon dioxide (CO 2 ) has been the central focus of climate research. But two important aspects of climate change response-regionality and nonlinearity-are much less studied and understood. For the first, climate change response is not uniform, exhibiting considerable regional heterogeneity even if the forcing is uniform 1 . What really matters to society and environment is regional climate projection since no one lives under the average climate. But regional downscaling approach often suffers the 'garbage-in-garbage-out' conundrum, especially as global climate models often disagree even on the large-scale features of the future climate 2,3 . Predicting regional climate change response is one of the grand challenges to the climate community 1,4,5 , while confidence in regional climate projection is sorely needed for climate mitigation and adaptation planning. An example in point is the projection of hydrological cycle: only over a small fraction of the globe do climate models show consensus on the sign of precipitation response 3 . The climate community is only beginning to tackle the challenge of regionality from the perspective of dynamical system 6,7 to understand and quantify the systematic relationships between the response and the forcing, each being geographically varying, to improve climate projection and confidence therein. This study stands as a pilot effort towards systematically linking the regional thermal forcing around globe to the global temperature through building linear response relationships and performing neutral mode analysis. As to be demonstrated later, such an effort may potentially be beneficial for exploring geoengineering experiments 8,9 for achieving the optimal cooling effect.
For the second aspect (nonlinearity), climate sensitivity research has been mostly focusing on the one-sided or symmetric response to the external forcing perturbation; many have assumed implicitly that the forcing of doubling CO 2 is small enough to warrant linearity 10,11 . However, if a system is sufficiently nonlinear, as some modeling studies have suggested [12][13][14][15] , one must account for the higher-order response to make accurate prediction. This nonlinearity question inquires into the intrinsic asymmetric response of the climate system to warming vs. cooling, a subject largely left unexplored in the climate change community 16 and will be another focus of this study.
Here, we develop a linear response framework making use of Green's function experiments, which comprise a large number of paired perturbation runs, with each pair being forced by two localized surface heat flux anomalies of the same magnitude but opposite signs, to probe the most excitable temperature patterns and their corresponding optimal forcing. This is achieved not only for the linear response, but also the quadratic nonlinear response 17 . The latter unveils the intrinsic propensity of our climate system towards global cooling, or in other words, our climate has a tendency to cool than to warm.

Q-flux Green's function experiments
To circumvent the uncertainty and ambiguity in defining the radiative forcing for climate change 12,16,18 , we perturb a slabcoupled climate model (NCAR CAM5-SOM, see Methods for details) with 99 pairs of ocean q-flux forcings. These overlapping forcing patches are arranged in such a way that they jointly cover all the open ocean grid points (see Supplementary Fig. 1); and the total forcing of all the (positive) patches is approximately a uniform 12 W m −2 over all ocean grid points, except near the peripheries of the ocean basins. The globally integrated forcing exerted in each patch experiment ranges from 0.03 W m −2 near the pole to 0.15 W m −2 near the equator, thus it is at least one order of magnitude smaller than the doubling CO 2 radiative forcing. A positive q-flux anomaly represents a surface energy gain of the slab ocean through processes such as surface radiation or energy convergence induced by ocean circulation, and the same amount of energy is passed to the overlying atmosphere as the slab ocean must be balanced thermodynamically. This allows us to quantitatively control the amount of energy perturbation into the atmosphere. This large suite of CAM5 experiments will be referred to as the q-flux Green's function experiments in the sense that the response becomes a Green's function solution in the limit of the forcing patch approaching a Dirac delta function. See Liu et al. 19,20 for further details about the experiment configuration. The approach herein should be distinguished from the similar sea surface temperature (SST) Green's function experiments performed to quantify the global radiative feedbacks to local SST perturbations 21,22 . The slab-coupled model here has interactive SST and the q-flux perturbation is in energy flux form (just as the radiative perturbation), affording a construction of forcingresponse relationship for global surface temperature.
The linear and nonlinear response Denoting δx + the response to the positive q-flux perturbation and δx − the response to the negative one of the pair, the linear response (or symmetric response as it is the component symmetric about the abscissa in the illustration of Fig. 1) is defined as δx 0 ðδx þ À δx À Þ=2; the nonlinear (asymmetric) response is defined as δx 1 ðδx þ þ δx À Þ=2. As elaborated in Methods and ref. 23 the linear response is governed by the Jacobian (the first derivative) of the dynamical system, thus capturing the component linearly proportional to the forcing perturbation (∝δf). On the other hand, the nonlinear response is forced by the variance of the linear response through the Hessian (second derivative) of the system, therefore representing the quadratic component of the response (i.e., ∝|δf| 2 ). These definitions are delineated in the right side panel in Fig. 1, which shows the sum of the responses to the 12 concatenated q-flux patches along 44°S. Although globally the linear component (gray lines) is overall larger than the nonlinear component (hatched), they can be comparable in magnitude at certain latitudinal range, especially near the southern pole. It is notable that the nonlinear component is a polar amplified global cooling, a feature that can be generalized to other forcing cases. The linear and nonlinear response to forcing from each of the 11 latitudinal bands are agglomerated and presented in Supplementary Fig. 2. The fact that every single one of the 11 cases drives a high-latitude amplified pattern for the linear and nonlinear response is strongly implicative of a mode behavior intrinsic to the climate system.
The most excitable mode of the linear and nonlinear response Given the non-normal dynamics of the climate system or its numerical representation, neutral mode (NM) analysis proves to be a powerful tool for understanding the modal behavior of the system. There is a rich trove of literature establishing neutral modes as the least damped or most excitable modes of a forceddissipative system [24][25][26][27] . For a linearized system governed by dynamical operator L (i.e., the Jacobian of the full dynamics of the system), the neutral modes can be extracted through the singular value decomposition (SVD) of L, with the leading mode corresponding to the smallest singular value and thus the least damping rate. For the climate model used here, surface temperature is governed by complex interactions among ocean, land, atmosphere, and sea ice, thus the Jacobian operator L cannot be derived explicitly. Instead, we resort to Green's function experiments to obtain an empirical representation of L, which is referred to as linear response function (LRF). See Methods for the specific procedure.
To the best of our knowledge, there is no prior study on the neutral modes of the nonlinear response, and the development of which represents a theoretical advance of this study. Through perturbation analysis, we show in Methods that the nonlinear response is driven by the variance of the corresponding linear component. Aided by the 99 pairs of q-flux Green's function experiments, we can construct the LRF relating the nonlinear response to the variance of the corresponding linear response and obtain the neutral modes for the nonlinear response via SVD analysis as well. As the forcings, responses, and the resultant LRFs, are all expressed in a hyperspace spanned by 62 EOFs, this is a 62dimension dynamical system. The leading right singular vector of the LRF gives the neutral mode, representing the most excitable mode of the climate system examined, while the leading left singular vector of the LRF give the patterns of the optimal forcing that is the most efficient in exciting the corresponding neutral mode. This mode decomposition applies to both the linear/ symmetric and nonlinear/asymmetric components of the response.
The leading neutral modes for the linear response (lNM1) and nonlinear response (nNM1) are displayed in Fig. 2a, b, and their corresponding optimal forcing patterns (i.e., the left singular Fig. 1 The zonal mean TS response to the q-flux pair placed along 44°S. As denoted, the red line is the response to positive forcing (δx + ), the blue the response to negative forcing (δx − ). The gray lines indicate the linear/symmetric component of the response (δx 0 ), the sign of which follows that of the forcing. Thence the δx 0 response component to the forcing pair are symmetric about the x-axis in the diagram. The hatched areas are for the nonlinear/asymmetric component (δx 1 ). As illustrated in the side panel, the nonlinear component can also be expressed as the difference of the full response minus the symmetric response and it is negative irrespective of the sign of the forcing.
vectors) are presented in Fig. 2c, d. The leading mode of the linear response imprints strongly on globally integrated warming, with the sign dependent on the sign of forcing, while the leading mode of the nonlinear response is a global cooling irrespective of the sign of forcing, reflecting the quadratic relationship between the nonlinear response and the forcing. An important attribute of the neutral modes is their power in capturing the variance of the response with very few modes. As shown in Supplementary Fig. 3, the linear response in each Green's function experiment mainly projects onto the first two modes; the nonlinear response can be predominantly accounted for by the first mode alone. As such, the leading modes offer a simple explanation for the key nonlinear/ asymmetric features in Fig. 1: In response to a pair of opposite qflux perturbations located at 44°S, the total TS response is the sum Fig. 2 Spatial structures of the leading linear and nonlinear NM. The TS patterns (a, b), the left singular vectors or the optimal forcing patterns (c, d), SLP patterns (e, f), precipitation patterns (g, h), and ice area patterns (i, j) of lNM1 and nNM1, respectively. In a, b, e-j, the black lines in the right side panels display the corresponding zonal mean profiles. In g, h, the red lines indicate the meridional distribution of AHT. Since these are derived from SVD analysis, the amplitude is all nondimensionalized. Note that the sign of the lNM1 corresponds to the sign of the forcing pattern in panel c, while the nNM1 is always a cooling irrespective of the sign of the forcing.
of linear warming (cooling) plus a nonlinear cooling for the positive (negative) forcing case. Because the nonlinear response is always a cooling regardless of the sign of the forcing, the total response tends to skew towards cooling.
Other physical variables associated with a given TS neutral vector (denoted by ΔT n ) can be retrieved from dZ dT ΔT n , where ∂Z ∂T is the sensitivity matrix to TS of the variable in question. See Methods for how ∂Z ∂T is estimated. The spatial patterns of TS, SLP and precipitation (Fig. 2a, e, g) of the lNM1 resemble markedly those associated with the Southern Oscillation patterns in an earlier version of CAM 28 and the Interdecadal Pacific Oscillation pattern in both observation and models 29 . The SLP pattern also bears considerable resemblance to the Pacific-North America pattern 30,31 and the teleconnection pattern associated with the positive zonal index in the Northern Hemispheres 32 , while it projects on both Pacific-South America and Southern Annular Mode (SAM) patterns in the Southern Hemisphere. The leading linear mode is associated with extensive sea ice loss in both polar regions, with an interhemispheric asymmetry that favors the Southern Hemisphere (Fig. 2i).
The first nonlinear NM (nNM1) shows an opposite circulation feature to that of lNM1 over North Pacific, consistent with the opposite SST pattern there. In terms of zonal index, it picks out an out-of-phase relationship between the North Pacific westerly and the North Atlantic one 32 . In Southern Hemisphere, it imprints strongly onto the SAM in both SLP (Fig. 2f) and atmospheric heat transport (AHT) (Fig. 2h, side panel, cf. Fig. 4a in ref. 33 ). Consistent with the conventional wisdom that AHT fluxes from the relatively warmer to the cooler hemisphere 34,35 , there is a strong crossequatorial AHT anomalies going along with nNM1, in conformation with interhemispheric asymmetry in the TS pattern (Fig. 2b, side panel). Specifically, the strong southward AHT near the equator gives rise to a northward shift of the zonal mean ITCZ, as expected from the energetics constraint on the position of the ITCZ [36][37][38] . Both the second linear and nonlinear modes (lNM2 and nNM2) capture the component of the interhemispheric thermal contrast; the corresponding TS, SLP, precipitation, sea ice concentration, and AHT patterns are illustrated in Supplementary  Fig. 4.
An important hallmark of the neutral mode is the bimodality of its PDF 39 . Fig. 3 show the PDFs of the projections of the TS time series onto the leading nonlinear neutral mode from (a) a 900-year CAM5-SOM simulation, (b) concatenated historical simulations of 36 CMIP5 models, and (c) two concatenated 20th century reanalyses (see Methods for the data sets used), respectively.
Indeed, all three PDFs are characterized by a bimodal distribution and a negative skewness, the latter being implicative of nonlinearity. Viewed in a nonlinear dynamical system framework, this perhaps reflects a greater potential barrier for the climate state to take on the negative phase of nNM1 relative to its positive counterpart. As such, the positive phase, which corresponds to a global cooling state, is more accessible in the phase space of the system. The fact that similar skewed bimodality is also found in reanalyses suggests the relevance of the discovered nNM1 here to reality, with an important implication for climate change response to external perturbations, that is, our climate has an intrinsic propensity towards cooling than warming.
Origin of the nonlinearity How does the nonlinear cooling come about in the first place? Several sources of nonlinearity previously identified in climate system (asides from ocean dynamics) can give a cooling effect, to name a few: the nonlinear cooling effect of vertical masking of both clear-sky quantities (such as water vapor) and cloud 22 , the remnant effect of the stabilized small ice-cap instability giving rise to a greater cooling than a warming under opposite forcing perturbations 40 , the negative nonlinear cloud radiative effect at surface attached to the seasonality in the polar climate 41 , plus the nonlinear atmospheric dynamics 39,42 . An parallel investigation has been undertaken through a suite of methodically designed experiments wherein sea ice formation is disabled to form an ice-free climate (see Methods and ref. 23 ). The result that the nonlinearity disappears completely in an ice-free climate (Supplementary Fig. 5) points to the presence of sea ice and the related feedback processes as the culprit for the nonlinearity. More specifically, the nonlinearity is seeded in the asymmetry between the ice-free and ice-covered coupling regimes: under a warming perturbation, sea ice melting exposes the deep mixed layer of the high-latitude ocean to the atmosphere and increases the effective heat capacity of the local system substantially, making it hard to warm further. On the other hand, under a cooling perturbation of the same magnitude, the newly formed ice shields the deep ocean mix-layer from interacting with the atmosphere, thereby reducing the effective heat capacity and making the system easier to cool. The same mechanism has been found operating in the phase delay and amplitude reduction in the annual cycle of the surface temperature in high latitudes under global warming 43 . This seed of asymmetry can then be further amplified by albedo, lapse rate, and evaporative feedbacks. Since these feedback processes are not the central focus of this paper, we refer the readers to ref. 23 for further supporting evidence for the assertions above.
Predicting the response in a pseudo-global cooling experiment In deriving the relationship between the quadratic nonlinear response and the variance of the linear response, an ansatz is posited in Methods that the climate system in question is a weakly nonlinear system so that the perturbation expansion in power series of ε is appropriate. As the strength of the nonlinearity of the system cannot be known a priori, this ansatz and the resultant LRFs must be validated with independent experiments. To this end, we devise a pseudo-global cooling experiment by imposing a negative 2 Wm −2 heat flux into the slab ocean everywhere, mimicking the effect of making the ocean surface more reflective by foaming the sea water 44 or reflecting the sunlight back to space by injecting reflective aerosols in stratosphere, though only over ocean 45,46 . Then, we utilize both the linear and nonlinear LRFs to predict the simulated TS response (see Methods for details), and the skill of the prediction is evident in Fig. 4, vindicating the LRFs constructed from q-flux perturbations of 12 Wm −2 magnitude. The breakdown to linear and nonlinear response and the corresponding predictions are presented in Supplementary Fig. 6. It is also clear that the response projects strongly onto the leading neutral mode as evidenced by the features such as the horseshoe shaped pattern in North Pacific, the cooling in the southeastern equatorial Pacific, and the polar amplified cooling. Remarkably, the patterns in Fig. 4 show great resemblance to the global TS response to the brightening of the ocean surface in the Southern Hemispheric subtropics in the G4Foam experiment in ref. 44 and that to the radiative flux perturbation due to the brightening of the subtropical marine stratocumulus in ref. 47 despite of the different models used in these two studies. Arguably, it is the same dynamical mode as identified in this study that is excited by the ocean surface foaming and cloud brightening.
The skillful prediction bodes promise to more targeted climate modification at least in the context of perfect model. In view of that a negative forcing (in terms of Wm −2 ) gives a much greater global temperature response than a warming one of the same magnitude, the approach of geoengineering to mitigate the GHGinduced warming might be somewhat easier than it would otherwise be without the nonlinearity. As far as regionality is concerned, one could even construct an optimal forcing for a targeted temperature pattern. If one is to cool the global climate through producing a negative pattern of lNM1 (Fig. 2a), its corresponding right singular vector (opposite of Fig. 2c) gives the forcing pattern it is most responsive to. Using Fig. 2c as a guide, the regions dotted with negative values should better be avoided because a cooling source there would even lead to global warming by projecting onto the positive lNM1 pattern.
Interestingly, Hill and Ming 47 found a much larger climate sensitivity to the marine stratocumulus brightening over subtropical south Pacific than subtropical south Atlantic, in agreement with the optimal forcing pattern for lNM1 revealed in Fig. 2c. For nNM1, one can see from Fig. 2d that an energy perturbation near the Maritime Continents or the southeastern coast of Australia might be the best location to excite global cooling as the nonlinear response. In contrast, North Sea may be a poor choice to try solar geoengineering there, as it not only is inefficient in exciting the linear response associated with lNM1, but also tends to excite a negative nNM1 that projects onto a positive global mean temperature.

DISCUSSION
Treating the climate system or its model representation as a highdimensional, weakly nonlinear dynamical system and probing it in a large set of q-flux Green's function experiments allow us to extract the most excitable modes-neutral modes-from the associated LRF for both linear and nonlinear TS response. Although the approach used here for estimating the LRFs is empirical, the resultant neutral modes reflect the deterministic physical and dynamical constraints built into the climate model. The leading nonlinear mode is a pattern of global cooling, implying an intrinsic tendency of skewing towards cooling of the climate system, with an intriguing corollary that our current climate is easier to cool than to warm given the amount of energy perturbation at one's disposal. Given the fact that the leading neutral mode is the most excitable mode and it projects most strongly onto global mean temperature, its corresponding left singular vector-representing the optimal forcing pattern-would be the most efficient forcing in producing a global warming/ cooling. The forcing map presented in Fig. 2c may potentially be of some value in the search for the optimal geographic locations for implementing the radiative cooling for geoengineering purpose.
Admittedly, we have a long way to go before having a map like that for cooling the climate in reality. Firstly, the prescribed q-flux perturbations to ocean slab are too idealized, perturbations at top of the atmosphere such as solar insolation or reflection by stratospheric aerosols may be more reflective of the geoengineering practice, plus the complexity due to the drift of the reflective substance in the ocean or stratosphere. Secondly, the q-flux forcing patch used here is too large in size, and thence the effective resolution of the forcing array is very limited and far from being able to resolve the true effective degrees of freedom of the forcing sensitivity patterns. Thirdly, the slab model used here does not allows the feedback from ocean dynamics, which can modify the temperature and atmospheric circulation patterns in a fundamental way 48,49 . It remains to be seen how the ocean dynamical feedback would modify the spatial structures of the  Supplementary Fig. 6. Regions where the response passes the Monte Carlo-based significance test at 99% level are marked with cross symbols. The ratio of the variance of the predicted signal to that of the simulated signal is also labeled on the top of the side panel of (b).
leading modes and the corresponding optimal forcings. Fourthly, we are not yet certain that the mode behavior in the model examined here represents the mode in reality. Further confidence may be built in this regard from inter-model comparison effort: if multiple climate models agree on the patterns of the neutral modes, considerable confidence can be assigned to these modes.
Notwithstanding all above, we believe the linear response framework developed here represents a methodological advance, it can be equally applied to more sophisticated models and/or other climate models. If the consensus on the modes can indeed be established across climate models one day, these modes can be used not only to guide more targeted implementation of solar geoengineering, but also as fingerprints for detecting and attributing the robust regional patterns of the climate change response in the past and future alike.

CAM5-SOM model and q-flux Green's function experiments
The model used in this study is the atmospheric component (Community Atmospheric Model version 5, CAM5) of Community Earth System Model version 1.1 coupled to a slab representing a thermodynamic ocean mixedlayer, the Community Land Model version 4 (CLM4), and the Community Ice CodE (CICE). Vertically, CAM5 is spaced in a hybrid pressure-σ coordinate with 30 levels; the horizontal resolution of CAM5 is 2.5°l ongitude × 1.9°latitude. The horizontal resolution of CICE and SOM is at nominal 1°, telescoped meridionally to~0.3°at the equator. The model will be referred to as CAM5-SOM, whose thermal coupling allows the SST to be calculated interactively as a prognostic variable.
Prior to perturbation runs, a long control must be established. This is achieved by forcing the CAM5-SOM with a repeating seasonal cycle of qflux estimated from a 900-year coupled CESM1.1 preindustrial control simulation, together with greenhouse gases, aerosols, and solar insolation at preindustrial level. The mixed-layer depth is also estimated based on the coupled CESM1.1 control simulation; it is geographically varying, but time independent. To facilitate the Green's function approach, we perturb the CAM5-SOM with an array of 99 q-flux anomaly patches, as illustrated in Supplementary Fig. 1. For each patch, both a positive and a negative forcing is added to the control field of q-flux and the simulations branch out from the 51st year of the control run and run for 40 years. Only the last 20 years are used for analysis. Following Barsugli and Sardeshmukh 50 , the q-flux patch is specified functionally to be a cosine hump ± Qcos 2 π 2 ϕÀϕ k ϕ w cos 2 π 2 λÀλk λw ; for Àϕ w ϕ À ϕ k ϕ w and Àλ w λ À λ k λ w 0; otherwise ( (1) where Q is the peak value of the hump, set to be 12 Wm −2 , (ϕ k , λ k ) denotes the latitude and longitude of the center of the patch, ϕ w = 12°and λ w = 30°are the half-widths of the rectangular patch in zonal and meridional directions, respectively, they also indicate the zonal and meridional distances between the adjacent patches. Making use of the forcing pair, we can define the linear (or symmetric) component and nonlinear (or asymmetric) component of the response in variable x as x þ Àx À 2 and x þ þx À 2 for each patch, with the superscript +(−) denoting the response to the positive (negative) forcing.
An extra pair of experiments with a forcing that is equivalent to the sum of all the 99 patches but reduced amplitude (−2 Wm −2 ) is also conducted as a test case for evaluating the skill of the prediction by the LRFs. More details of these runs are listed in Supplementary Table 1.

Process-disabling experiments
Driven by the hypothesis that the large nonlinearity might be rooted in the sea ice component of the CAM5-SOM climate system, we perform an extra set of experiments. First we disable the formation of sea ice by letting seawater to be supercooled so that sea ice cannot form (NI configuration). In the ice-disabled configuration, we further perturb the q-flux along 55°S with a pair of positive and negative anomalies of the same latitudinal structure as defined in Eq. (1) and of peak value Q = 12 Wm −2 . Through this set of experiments, we can examine the effect of the sea ice on the response nonlinearity. Further details about these experiments are provided in Supplementary Table 1. The comparison of nonlinear responses to q-flux forcing at 55°S in the ice-disabled and ice-permitting configurations are reported in Supplementary Fig. 5.

LRFs for linear and nonlinear response
Leveraging the paired q-flux Green's function experiments described above, we adopt a perturbation expansion approach to the governing dynamical system to formulate the LRFs for both the linear and nonlinear component of the response. Supposing the perturbations are with respect to a reference state x (represent a climate variable, say, surface temperature) governed by a nonlinear system: N x ð Þ ¼ f; the perturbation responses to a pair of small and opposite forcings δf + and δf − (δf − = −δf + ) can be expressed as follows via Taylor expansion: where L ∂N ∂x x is the Jacobian of the operator N and H is the Hessian of N, expressed as ∂N ∂x and respectively. To reduce dimension, x variable is expressed in a hyperspace spanned by n EOFs (x 2 R n ) and thus L is a matrix in R m n space (m is the dimension of physical space, or number of grid points in a numerical model) and H is a third-order tensor with frontal-slice symmetry in R n n m space.
To the extent that the forcing perturbations are sufficiently small, we assume the linear response is larger than the higher order nonlinear response and the full response may expanded in the powers of a small parameter ε: for the Oð1Þ order linear response and the O ε 1 ð Þ order nonlinear response, respectively. We hence refer to δx 0 as the linear response, as it is governed by the linear dynamics of the system and proportional to the forcing, and refer to εδx 1+ as O ε 1 ð Þ nonlinear response, which is independent of the sign of the forcing, in the context of this perturbation approach.
Taking difference between Eqs. (2) and (3) and dividing it by 2 yields the dynamical system governing the O ε 0 ð Þ linear response, wherein the cross term between δx 0 and δx 1 at the O ε 1 ð Þ order has been neglected. For the q-flux Green's function experiments, we have 99 pairs (indexed by k) in total and all the linear responses satisfy Eq. (6). Stacking them together, we have where δX 0 2 R n k is the stacked O ε 0 ð Þ linear response, F 2 R m k is the forcing matrix containing the forcing vectors as its columns, L 2 R m n is the linear response function (LRF) in matrix form. As Eq. (7) implies, δX 0 is the linear response directly forced by the prescribed q-flux perturbations. To estimate L, we in practice estimate the Green's function matrix G following the procedure of Barsugli and Sardeshmukh 50 . The column vector of G is computed as through normalized linear response by the area integrated q-flux within the patch (k). Note that subscript j here represents the grids within the kth patch. The Green's function-based LRF can be then computed as the pseudo-inverse of G.
The singular value decomposition (SVD) factors L into orthogonal left and right vectors via where U, V are orthogonal matrices and ∑ is the diagonal matrix of the singular values. Substituting it back into Eq. (6), we can see that a response δx 0 to an arbitrary forcing δf + can be expressed as where the angle bracket represents inner product. As such, the solution δx 0 can be retrieved by a linear combination of the v vectors (NMs), each weighted by 1 σi times the projection of the forcing on the corresponding u vector. If the hu i ; δf þ i factor is the same for all singular modes, since the leading mode v 1 is weighted by the inverse of the smallest singular value σ 1 , it would be manifested with the largest amplitude as the most excitable mode (or neutral mode) of the system. Whereas, for a given v i vector, a forcing that has the same spatial structure as the corresponding u i vector can excite the largest variance. Therefore, u i vector gives the optimal forcing pattern for vector v i .
A simple manipulation [(2) + (3)]/2 and ignoring the terms of Oðε 2 Þ order lead to the equation linking the symmetric/linear response to the asymmetric/nonlinear response Note here H 2 R n n m and is a tensor of frontal symmetry and the depth dimension m is independent of the dimensions of the frontal slices. Thus for each m, the corresponding frontal slice is diagonalizable: where the columns in P contain the unit eigenvectors of H at a given depth m; D is a diagonal matrix containing n non-zero eigenvalues. Defining a new vector δy 0 P À1 δx 0 , (11) becomes εLδx 1 % Àδy 0T Dδy 0 ; with the variance of δy 0 being equal to the variance of δx 0 (i.e., σ y 0 ¼ σ x 0 ) as required by P T P = I. An important insight gained from this exercise is that the O ε 1 ð Þ order nonlinear response is forced by the weighted variance of the linear response, with the weight being provided by D.
Grouping all the m grid points together, we have the following relation in matrix form (with the dimensionality expressed explicitly): where D (m,1) is a vector containing elements d m ¼ P n ðe n Þ m δy 2 n =σ 2 y 0 , with e n being the eigenvalues of the corresponding frontal slice of H and δy n the component of δy 0 . Multiplying (13) with the pseudo-inverse of L ðm;nÞ , we yield δx 1 ðn;1Þ % À σ 2 x 0 ε L À1 ðn;mÞ D ðm;1Þ ; (14) or δx 1 ðn;1Þ % ÀM n;1 ð Þ σ 2 x 0 ; (15) showing that the ε 1 -order nonlinear response is forced by the variance of the O 1 ð Þ order linear response. Since the weight factor d m is dependent on δx 0 , M also varies from case to case, reflecting its sensitivity to the location of the forcing perturbation. Stacking all the 99 q-flux forcing cases (indexed by k) together, we yield finally δX 1 % ÀMΣ; where δX 1 2 R n k is the collection of k number of ε 1 -order nonlinear response, Σ 2 R k k is a diagonal matrix whose non-zero diagonal elements are σ 2 x 0 . Finally, M À1 2 R k n gives the LRF of δx 1 in response to the variance of the corresponding linear response.
Since the response is expressed in an orthogonal hyperspace spanned by 62 EOFs, M 2 R 62 99 . The retainment of 62 EOFs is based on the North criterion 51 for well-resolved EOFs by the 900-year control simulation. Making use of the 99 pairs q-flux experiments, M can be estimated by simply computing its columns m ¼ δx 1 =σ 2 x 0 for each k. The LRF of the nonlinear response δx 0 can then be estimated by the pseudoinversion of M. Further, as done with the LRF for the linear response, SVD can also be performed on M −1 to extract the leading dynamical modes of the nonlinear response and we refer to them as nonlinear neutral modes.
The resultant LRFs can be utilized to predict the TS response to arbitrary q-flux perturbations. For a given forcing δf of arbitrary pattern, the linear response can be estimated via L À1 δf. The variance of the predicted linear response is then substituted into the diagonal of ∑ in Eq. (14) to predict the O ε 1 ð Þ order nonlinear response. An example of predicting the response to a pseudo-global cooling forcing is presented in Fig. 4 and Supplementary  Fig. 6, exhibiting promising skills in both pattern and magnitude, despite the fact that the LRFs are built from forcing perturbations of a much larger amplitude.
The associated patterns with the neutral modes in other variables (δZ), such as SLP, precipitation etc., can be obtained by the following procedure. First, we estimate the sensitivity matrix as dZ dT % δZðδTÞ À1 where δZ is a matrix comprised of vectors (or columns) representing the response to each of the forcing case, and δT the response in TS, each being expressed in the corresponding EOF basis. The patterns of the nth mode in Z variable is then where ΔT n denotes the TS pattern of the nth neutral mode. The same procedure can be applied to both linear and nonlinear modes.

DATA AVAILABILITY
The NOAA-CIRES twentieth-century reanalysis data version 2 (20CR) used in this research are available at http://www.esrl.noaa.gov/psd/data/. The ECMWF's first atmospheric reanalysis of the 20th century (ERA-20C) are available at https://apps. ecmwf.int/datasets/. The model data simulated with CAM5-SOM analyzed in the current study are available from the corresponding authors upon request.