De Luca
email: p.deluca@vu.nl

| Messori

| Pons

| Faranda

Dynamical Systems Theory Sheds New Light on Compound Climate Extremes in Europe and Eastern North America

Keywords: Compound Extremes, Multi-Hazards, Climate Dynamics, Climate Extremes

We propose a novel approach to the study of compound extremes, grounded in dynamical systems theory. Specifically, we present the co-recurrence ratio (α), which elucidates the dependence structure between variables by quantifying their joint recurrences. This approach is applied to daily climate extremes, derived from the ERA-Interim reanalysis over the 1979-2018 period. The analysis focuses on concurrent (i.e. same-day) wet (total precipitation) and windy (10m wind gusts) extremes in Europe and concurrent cold (2m temperature) extremes in Eastern North America and wet extremes in Europe. Results for wet and windy extremes in Europe, which we use as a test-bed for our methodology, show that α peaks during boreal winter. High α values correspond to wet and windy extremes in north-western Europe, and to large-scale conditions resembling the positive phase of the North Atlantic Oscillation (NAO). This confirms earlier findings which link the positive NAO to a heightened frequency of extra-tropical cyclones impacting north-western Europe. For the Eastern North America-Europe case, α extremes once again reflect concurrent climate extremesin this case cold extremes over North America and wet extremes over Europe. Our analysis provides detailed spatial information on regional hotspots for these compound ex-1 2 DE LUCA ET AL.

treme occurrences, and encapsulates information on their spatial footprint which is typically not included in a conventional co-occurrence analysis. We conclude that α successfully characterises compound extremes by reflecting the evolution of the associated meteorological maps. This approach is entirely general, and may be applied to different types of compound extremes and geographical regions.

Here, we propose a novel approach grounded in dynamical systems theory which overcomes some of these limitations, and aims to provide a complementary view to more traditional analyses issuing from statistics and climate dynamics (e.g. [START_REF] Martius | A global quantification of compound precipitation and wind extremes[END_REF][START_REF] Waliser | Extreme winds and precipitation during landfall of atmospheric rivers[END_REF]. We specifically propose an objective measure of the co-recurrence of extremes in different climate variables, which can provide both temporal and spatial information and can be linked to the underlying dynamical properties of the climate system, as well as explicitly accounting for its underlying chaotic nature. This is complemented by two other metrics characterising the evolution of large-scale climate fields [START_REF] Lucarini | Extremes and recurrence in dynamical systems[END_REF]Faranda et al., 2017b;[START_REF] Messori | A dynamical systems approach to studying midlatitude weather extremes[END_REF]. The approach is very flexible, and can in principle be applied to any number of variables, geographical regions and types of dataset. Within the text, we adapt the vocabulary used in the literature to fit the novel approach we propose. We will be referring to compound extremes when considering bivariate metrics calculated from dynamical systems theory and to concurrent extremes when discussing the joint occurrences (i.e. same-day) of extreme values in climate variables.

We begin by providing a description of our data (Section 2) and methodology (Sections 3,4). Next, we illustrate the physical interpretation of the different dynamical systems metrics and their seasonality (Section 5. 1.1), and present an application to concurrent wet and windy extremes in Europe (Section 5.1.2), along with concurrent cold extremes in Eastern North America and wet extremes in Europe (Section 5.2). We thus address concurrent extremes in both single and separate geographical regions. Our choice is motivated by the significant scientific and media attention that both wet and windy weather in Europe (e.g. [START_REF] Pinto | Serial clustering of extratropical cyclones over the North Atlantic and Europe under recent and future climate conditions[END_REF][START_REF] Van Oldenborgh | Drivers of the 2013/14 winter floods in the UK[END_REF][START_REF] Huntingford | Potential influences on the United Kingdom's floods of winter 2013/14[END_REF][START_REF] Wild | Was the Extreme Storm Season in Winter 2013/14 Over the North Atlantic and the United Kingdom Triggered by Changes in the West Pacific Warm Pool?[END_REF] and cold spells in North America (e.g. [START_REF] Trenary | Was the cold eastern US winter of 2014 due tro Increased Variability?[END_REF][START_REF] Lee | Compounding effects of warm sea surface temperature and reduced sea ice on the extreme circulation over the extratropical North Pacific and North America during the 2013-2014 boreal winter[END_REF][START_REF] Messori | On cold spells in North America and storminess in western Europe[END_REF][START_REF] Harnik | The circumglobal north american wave pattern and its relation to cold events in eastern north america[END_REF] have attracted in recent years. The analysis in Section 5.1 primarily aims to relate the information provided by our metrics to known features of the atmospheric variability over the Euro-Atlantic region, and can be viewed as a test of the robustness of the metrics we propose. Section 5.2 instead illustrates how our approach can be used to gain novel insights into concurrent extremes in two different regions. We conclude by discussing our results in the context of the literature on these classes of extremes (Section 6).

| DATA

We use ECMWF's ERA-Interim reanalysis [START_REF] Dee | The ERA-Interim reanalysis: configuration and performance of the data assimilation system[END_REF], at a horizontal resolution of 0.75 • , from 1979 to 2018. ERA-Interim outperforms other reanalysis products, such as MERRA, NCEP-NCAR, ERA-40, CFSR and GLDAS, in reproducing temperature and wind-speed observations. Precipitation is notoriously problematic to capture on a relatively coarse horizontal grid, but ERA-Interim nonetheless performs comparatively well in this respect [START_REF] Decker | Evaluation of the reanalysis products from gsfc, ncep, and ecmwf using flux tower observations[END_REF]. The analysis was conducted on daily data at 12 and 00 UTC (12h step) of total precipitation (mm) and 10m wind gusts (ms -1 ). To obtain a unique daily value, these were then summed and averaged, respectively. We also used 2m temperature (K ) at 00, 06, 12 and 18 UTC, which were averaged as daily observations. From now on we will refer to these variables simply as precipitation, wind and temperature. We consider two domains: Europe (18

• W-51 • E, 30 • N-75 • N) and Eastern North America (99 • W-75 • W, 30 • N-51 • N).
Daily data over 1979-2018 for the North Atlantic Oscillation Index (NAOI) [START_REF] Barnston | Classification, Seasonality and Persistence of Low-Frequency Atmospheric Circulation Patterns[END_REF], the dominant mode of climate variability within the North Atlantic region, were downloaded from the NOAA -Climate Prediction Center website (available here). The analyses were also repeated by making use of daily maximum wind gust observations (ms -1 ), i.e. the daily maximum of the 12 and 00 UTC values for the 1979-2018 period. The results are qualitatively very similar compared to daily wind means (see Supporting Information).

| DYNAMICAL SYSTEMS METRICS

| Qualitative interpretation of the dynamical systems metrics

We use three dynamical systems metrics: the local dimension (d ), the local persistence (θ -1 ), and the local co-recurrence ratio (α). d and θ -1 are calculated for both univariate and bivariate cases, whereas α always refers to two different variables. All three metrics are local in phase-space, and instantaneous in time. They can thus be computed for any state ζ on the underlying attractor -a state in our case being a latitude-longitude map of one or more atmospheric variables at a given time. d and θ -1 have recently been applied to a range of different climate variables over different geographical domains and were found to successfully reflect large-scale features of atmospheric motions (Faranda et al., 2017a(Faranda et al., ,b, 2019a(Faranda et al., , 2020;;[START_REF] Messori | A dynamical systems approach to studying midlatitude weather extremes[END_REF][START_REF] Rodrigues | Dynamical Properties of the North Atlantic Atmospheric Circulation in the Past 150 Years in CMIP5 Models and the 20CRv2c Reanalysis[END_REF]Hochman et al., 2019a,b).

The local dimension d (ζ) describes the evolution of the system around ζ, and can be interpreted as a proxy for the number of degrees of freedom active around the state of interest. The local persistence θ -1 (ζ) measures the mean residence time of the system around such state. A high (low) persistence indicates that the system's evolution will slowly (rapidly) lead to a dynamically different state. The higher the value of θ -1 , the more likely it is that the preceding and future states will resemble ζ over comparatively long timescales (Faranda et al., 2017b;[START_REF] Messori | A dynamical systems approach to studying midlatitude weather extremes[END_REF]Hochman et al., 2019a). The two dynamical systems metrics can be computed for both individual variables (e.g. sea-level pressure, Faranda et al., 2017b) and multiple variables jointly (Faranda et al., 2020). Given two different variables, we can define a state on the Poincaré section jointly spanned by x and y as ζ = {ζ x , ζ y }. We will use the notation d x or θ -1

x to refer to the monovariate metrics computed for the atmospheric variable x and the notation d x,y or θ -1

x,y to refer to the bivariate metrics jointly computed on the atmospheric variables x and y .

Along with d and θ -1 , we adopt a phase-space local measure of dependence between variables first introduced in Faranda et al. (2020). We term this the co-recurrence ratio (α).

| Derivation of the dynamical systems metrics

We view the atmosphere as a dynamical system with an observed phase-space trajectory x (t ), and consider the properties of such system near a state of interest ζ. To compute our dynamical systems metrics, we first define our observable, via logarithmic returns as:

g (x (t ), ζ) = -log[dist(x (t ), ζ)] (1)
Where d i st is the Euclidean Norm. This observable discriminates parts of the trajectory x (t ) that are close to ζ (i.e. when the distance approaches 0) from those that are further away. We then define s(q, ζ) a high q -th quantile of the time series g (x (t ), ζ), and we define g (x (t ), ζ) > s(q, ζ) an exceedance u(ζ) = g (x (t ), ζ)s(q, ζ). According to the Freitas-Freitas-Todd theorem [START_REF] Freitas | Hitting time statistics and extreme value theory[END_REF], later modified in [START_REF] Lucarini | Universal Behaviour of Extreme Value Statistics for Selected Observables of Dynamical Systems[END_REF], the cumulative probability distribution F (u, ζ) then converges to the exponential member of the Generalized Pareto Distribution:

F (u, ζ) exp -ϑ(ζ) u(ζ) σ(ζ) (2)
The parameters ϑ and σ depend on the chosen point ζ. ϑ is the so-called extremal index [START_REF] Moloney | An overview of the extremal index[END_REF], and we estimate it here using the approach of [START_REF] Süveges | Likelihood estimation of the extremal index[END_REF]. The local dimension d is obtained as: d (ζ) = 1/σ(ζ), while the local persistence is given by: θ -1 = ∆t /ϑ, where ∆t is the timestep of the data used and the local persistence is in the same units as the timestep. In our case, trivially ∆t = 1 day.

The bivariate d and θ -1 metrics are derived following an analogous procedure (Faranda et al., 2020) (3)

Where . represents the average root mean square norm of a vector's coordinates. One can then obtain the codimension d x,y analogously to the derivation of the univariate d (for ease of notation we hereafter drop the dependence on ζ). The co-dimension has the following property:

min(d x , d y ) ≤ d x,y ≤ d x + d y (4)
The case d x,y = d x + d y implies no coupling bewteen x and y . If instead the variables are deterministically coupled (i.e. x is a function of y or viceversa), d x,y = min(d x , d y ).

Following the definition of joint logarithmic returns given in Eq. 3, one can also define the local co-persistence θ -1 x,y . This is defined as a weighted average of θ -1

x and θ -1 y , where the weights depend on the size of the hyper-ball around ζ in the Poincaré sections x and y (Faranda et al., 2020;[START_REF] Abadi | Dynamical counterexamples regarding the Extremal Index and the mean of the limiting cluster size distribution[END_REF].

Along with d and θ -1 , we adopt a phase-space local measure of dependence termed co-recurrence ratio (Faranda et al., 2020). For two variables x (t ) and y (t ), the co-recurrence ratio 0

≤ α(ζ) ≤ 1 of a state ζ = {ζ x , ζ y } is defined as a ratio: α(ζ) = ν[g (x (t )) > s x (q ) |g (y (t )) > s y (q )] ν[g (x (t )) > s x (q )] (5)
Where s x (q ) and s y (q ) are high q -th quantiles (or thresholds) of the univariate logarithmic returns g (x (t )) and g (y (t )), and ν[-] represents the number of events that satisfy condition [-].

| STATISTICAL ANALYSIS

| Definition of anomalies and extremes

We define anomalies relative to a daily climatology. For example, the climatological temperature at a given gridpoint for the 19 t h February is the mean of all 19 t h February temperatures at that location over the 40 years considered here.

Precipitation and wind extremes are defined as daily values > 99 t h quantile, whereas cold temperature extremes are daily observations < 1 st quantile of the respective distributions at each gridpoint. The selection was performed on absolute values, rather than anomalies, in the spirit of an impacts-based perspective. Indeed, a very strong windstorm is more likely to cause damage than an out-of-season (and thus extreme in terms of anomalies) but moderate (in terms of absolute wind-speed values) windstorm, although out-of-season extremes are associated with their own specific sets of risks. The same analysis was also performed for extremes based on anomalies of the climate variables and the two approaches show good spatio-temporal agreement (see Supporting Information).

The anomalies of precipitation, wind and temperature are stratified according to α daily extremes (> 99 t h quantile), thus informing on the behaviour of the climate variables under strong dynamical coupling conditions. We also consider, for each gridpoint, the composite α mean anomalies associated with univariate and concurrent extremes of precipitation, wind and temperature. A list of extreme occurrences at each gridpoint was compiled, and then the corresponding daily α anomalies were composited, thus providing a different composite α anomaly value at each gridpoint. As an example, consider precipitation and wind extremes for a gridpoint within the European region. First, we identified days when the two extremes (co)occurred at that location, and then we computed the α mean anomalies based on those days (note that the α daily values are the same for both types of extremes, since it is computed using both variables simultaneously). For the case of concurrent cold and wet extremes in Eastern North America and Europe respectively, we adopted a similar procedure. We first identified dates on which at least one gridpoint in each domain displays an extreme observation. We next identified all gridpoints displaying extreme values on each selected date. Finally, we calculated the α mean anomalies for each selected gridpoint over the two regions. Note that, again, α daily values are the same for both variables, regardless of the fact that they are taken over different geographical domains.

| Statistical testing

We evaluate the significance of the anomalies in climate variables associated with α extremes using the Mann-Whitney test [START_REF] Mann | On a test of whether one of two random variables is stochastically larger than the other[END_REF]. This consists of a nonparametric procedure that, given two samples x and y of size n x and n y , tests the null hypothesis that the samples are drawn from populations with the same median. We perform the test in each grid point, in such a way that y is the sub-sample selected using extreme values of α, and x is the remaining of the time series in that grid point. We choose the Mann-Whitney test since the data considered in this analysis, especially wind and precipitation, are strongly non-Gaussian and thus do not meet one of the assumptions of the Student's t test on the sample means. Nevertheless, an exact application of the test would require independently and identically distributed samples, while time series of atmospheric variables are expected to exhibit at least significant autocorrelation, or even nonlinear dependence in time. In general, the presence of serial dependence in the data reduces the power of the test, i.e. the probability to not reject the null hypothesis when the median difference is significant is higher than in the case of samples meeting all the test assumptions (e.g. [START_REF] Yue | The influence of serial correlation on the mann-whitney test for detecting a shift in median[END_REF]. If there are realistic expectations about the dependence structure in the samples, the distribution of the test statistics for the dependent case can be approximated using Monte Carlo sampling. However, this is not the present case, and we decide to apply the regular Mann-Whitney test, aware of the possible loss of power. Indeed, we believe that this issue should not excessively affect the results, since the test is performed over a large number of grid points, and our expectation is to observe regions of significant effects. While the total number of significant points in one of such regions may be reduced due to loss of power, we expect it to be still possible to draw robust general conclusions.

The procedure consists of the following steps: i) the two samples are merged and ranked; ii) the sum of the ranks is computed separately for data belonging to sample x and y , and denoted R x and R y respectively; iii) the Mann-Whitney U test statistics is obtained as U x = R xn x (n x + 1)/2. In case of ties, the rank in the midpoint between the two closest non-tied ranks is used. Under the null hypothesis, U follows a known tabulated distribution, and for samples large enough (n 20) converges to a Normal distribution. Notice that U can be computed indifferently using sample x or y . We state the alternative hypothesis based on prior knowledge that dynamical extremes should correspond to large positive precipitation and wind anomalies, and to large negative temperature anomalies. Then, the alternative hypothesis is that the median of the sample selected with α is larger than the median of the rest of the time series for precipitation and wind, and smaller for temperature.

A similar procedure was also applied to check whether the α mean anomaly values associated with univariate extreme climate variables were significant. In particular, we perform again the Mann-Whitney test, but this time the two samples x and y consists of the percentiles of each value of α in the distribution of the merged samples of α mean anomalies. The alternative hypothesis is now that the values of α corresponding to climate extremes are -in mediansignificantly larger than the ones not corresponding to climate extremes. The strategy for concurrent extremes is similar to the univariate case, merging the time-series of different climate variables by dates (i.e. precipitation-wind in Europe and cold temperatures-precipitation over Eastern North America and Europe), so that only extreme observations co-occurring on the same days are kept, along with the associated α daily observations. We then repeat the procedure as per the univariate case above.

The test is conducted at the standard level a 0 = 0.05. However, the level requires a correction, due to the fact that, for each event, we perform as many tests as the number of gridpoints, n g r . This increases the probability of erroneously rejecting the null hypothesis, leading to false positives. On the other hand, in the case of concurrent extremes, a number n N A of tests is skipped due to missing data. To overcome this issue, we apply a Bonferroni correction [START_REF] Bonferroni | Teoria statistica delle classi e calcolo delle probabilita[END_REF], which consists of using the original level divided (or corrected) by the effective number of tests, a = a 0 /(n g rn N A ).

The robustness of the α anomaly means observed during univariate and concurrent extremes was also assessed by performing a sign test. This quantifies the fraction of individual extreme observations at a given gridpoint which have the same sign anomaly as the overall composite. For example, a sign test value of 66% at a location with a positive α composite anomaly means that 66% of the days used to create the composite display a positive α daily anomaly, while 34% display a negative anomaly. Clearly, the higher the % of members sharing same-sign anomalies, the more robust the results are.

| CLIMATE EXTREMES IN EUROPE AND NORTH AMERICA

| Wet and windy extremes in Europe

| Dynamical climatology and compound extremes

In Figure 1 we show both univariate and bivariate dynamical metrics for precipitation and wind in Europe. α peaks during boreal winter (Dec-Feb) with values in excess of ∼0.17 and reaches its lowest values during spring (∼0.075) and summer (∼0.10, Figure 1a). The α winter peak follows the intuition that during these months the passage of extra-tropical cyclones (ETCs) over western-continental Europe and the British Isles (BI) can lead to widespread flooding episodes and comparatively frequent concurrent wet and windy episodes, whose spatial footprints resemble one another [START_REF] De Luca | Extreme multi-basin flooding linked with extratropical cyclones[END_REF][START_REF] Fink | The European storm Kyrill in January 2007: synoptic evolution, meteorological impacts and some considerations with respect to climate change[END_REF][START_REF] Lavers | Winter floods in Britain are connected to atmospheric rivers[END_REF][START_REF] Matthews | A cyclone climatology of the British-Irish Isles 1871-2012[END_REF]. On the other hand, the lower α values observed during spring and summer may reflect precipitation events which are mainly driven by small-scale convective processes (e.g. [START_REF] Pieri | Sensitivity of Precipitation Statistics to Resolution, Microphysics, and Convective Parameterization: A Case Study with the High-Resolution WRF Climate Model over Europe[END_REF] and therefore have a weaker link to both the synoptic-scale and larger-scale circulation [START_REF] Bisselink | Precipitation Recycling: Moisture Sources over Europe using ERA-40 Data[END_REF].

During α extremes (i.e. strong dynamical coupling), negative d anomalies are observed for both the compound and univariate cases (Figure 1b), with the former and wind displaying the more numerous negative anomalies. An almost opposite pattern of positive values is found for θ -1 anomalies (Figure 1c), with precipitation again being the variable showing the weaker deviation from climatology. Low d and high θ -1 are associated with predictable configurations (Faranda et al., 2017b;[START_REF] Messori | A dynamical systems approach to studying midlatitude weather extremes[END_REF][START_REF] Scher | Predicting weather forecast uncertainty with machine learning[END_REF], and we find here that this matches a strong coupling between the large-scale configurations of the two variables we are analysing.

We next relate α to the NAO index (NAOI) [START_REF] Barnston | Classification, Seasonality and Persistence of Low-Frequency Atmospheric Circulation Patterns[END_REF]. The NAOI daily values during α extremes (Figure 2a) are overwhelmingly positive (mean=0.85, orange dashed line), and well above about what may be obtained by chance (violet dot-dashed line). When the NAO is in a positive phase, Northern/Western Europe experiences enhanced storminess from the passage of ETCs, which bring concurrent flooding and windstorms [START_REF] De Luca | Extreme multi-basin flooding linked with extratropical cyclones[END_REF]Lee et al., 2018;[START_REF] Nobre | The role of climate variability in extreme floods in Europe[END_REF][START_REF] Huntingford | Potential influences on the United Kingdom's floods of winter 2013/14[END_REF]. This agrees with the strong negative and positive anomalies respectively in d and θ -1 , discussed above. Indeed, the positive NAO has been previously identified as a low-d highθ -1 configuration (Faranda et al., 2017b). Moreover, the vast majority of α extremes occur during late autumn and winter, with no observations in May and summer (June-August, Figure 2b), further confirming the plausibility of the link between α extremes and winter ETCs.

Very similar results are obtained when repeating the above analysis on daily wind maxima (Figures S1-S2), instead of wind means, computed over two daily observations. For the maxima, α daily climatological means are lower (Figure S1a), likely reflecting the intrinsically noisier nature of the variable.

| Concurrent extremes

Here, we consider the link between extremes in α computed on the temperature and wind fields and climate extremes.

During α extremes, positive anomalies in both the precipitation and wind maps are observed in Northern/Western Europe, over a broad region spanning the BI, the North Sea and the surrounding coastal areas. Negative anomalies are instead observed in Southern Europe and the Mediterranean (Figure 3a-b). This pattern resembles that associated with the NAO (Figure S3), further strengthening our interpretation that α extremes are able to capture the impacts of ETCs in Northern/Western Europe.

The number of concurrent wet and windy extremes peaks during late boreal autumn and winter (Figure 3c). This closely reflects the timing of the α extremes (Figure 2b), proving additional evidence for a close link between compound dynamical and concurrent climate extremes.

Finally, we test the inverse link, namely whether large α values can be recovered when conditioning on extremes in the climate variables. The α anomaly means, conditioned on monovriate precipitation or wind extremes, are highest in Northern/Western Europe, although for wind significant α anomalies span virtually the whole continent (Figure 3d-e).

The signal for precipitation is somewhat patchy, but there is an overall good spatial correspondence with the regions displaying large anomalies in the climate variables themselves (Figure 3a-b). We also note that the α anomaly means related to wind extremes are higher than those for precipitation. This is not surprising, since precipitation has a noisy and fractal structure, which is reflected in the dynamical systems metrics (Faranda et al., 2017a). α anomaly means calculated for concurrent wet and windy extremes mostly reflect the patterns observed for the univariate extreme cases (Figure 3f). The largest values are again distributed over the BI, Northern/Western Europe and the Nordic Seas.

The α anomaly means for the concurrent extremes show even higher values than for wind alone. We interpret this as a direct result of the α metric reflecting the dynamical coupling of the two variables. The sign test (stippling in Figure 3d-f) largely mirrors the magnitude and spatial pattern of the α anomaly means. The same analyses as shown in Figure 3 and 3c-f were respectively repeated by using daily wind maxima observations (Figure S4) and a different definition of extremes (see Section 4.1 and Figure S5). The results are consistent between the three approaches.

| Cold and wet extremes in Eastern North America and Europe

In this section, we consider temperature over Eastern North America and precipitation over Europe (Figure 4). α is computed on these two variables, each in its respective domain. α extremes are associated with negative temperature anomalies across Eastern North America, with the largest anomalies situated in the north-eastern part of the domain (Figure 4a). Precipitation anomalies are instead mostly positive, and are strongest over Iberia, Western France and Iran (Figure 4b). The vast majority of the concurrent cold and wet extremes are observed during the boreal winter months (Dec-Feb). No co-occurrences are observed from April to October (Figure 4c). North American domain and Iberia/Western France (Figure 4f). These reflect very closely the mean anomaly patterns in Figures 4a-b. Similar to what we found for concurrent wet and windy extremes in Europe (Figure 3f), these α anomaly means for the concurrent extremes display larger values than for the univariate extremes. The monovariate α anomaly means mostly fail the sign test, while the anomalies for the concurrent extremes display extensive sign agreement mainly over Europe (stippling in Figure 4d-f). We interpret these differences as indicating that the α computed on the North American temperature and European precipitation reflect more closely the concurrent extremes in these variables than the monovariate extremes in the individual domains. As for Figure 3c-f, similar results are retrieved for a different definition of extremes (Figure S6).

| DISCUSSION AND CONCLUSIONS

In this study we applied a new approach issued from dynamical systems theory [START_REF] Lucarini | Extremes and recurrence in dynamical systems[END_REF]Faranda et al., 2017b;[START_REF] Messori | A dynamical systems approach to studying midlatitude weather extremes[END_REF]Faranda et al., 2020) to the study of compound and concurrent climate extremes. We use the terms compound extremes when considering bivariate metrics calculated from dynamical systems theory and concurrent extremes when discussing the joint occurrences (i.e. same-day) of extreme values in the climate variables. We specifically focused on wet and windy extremes in Europe and cold and wet extremes respectively in Eastern North America and Europe. We defined a coupling parameter α (the co-recurrence ratio) that measures the joint Poincaré recurrences of specific configurations in two variables. Almost by definition, α is sensitive to the chosen geographical domain(s). Indeed, the dynamical coupling between the variables of interest can change significantly when including or excluding data. It is therefore important to limit the analysis only to regions of interests, whatever the criterion for "interesting" may be. Two more dynamical systems metrics were also computed for both uni-and bivariate cases, namely the local dimension (d ) and the local persistence (θ -1 ). These inform on the number of degrees of freedom active around a given state of the system (here the atmosphere in the Euro-Atlantic and Eastern North American sectors, as represented by the temperature, wind gust and precipitation maps) and on the mean residence time of the system around such state, respectively.

To verify whether our approach can provide new insights into the nature of concurrent extremes, we compare our results for the wet and windy extremes in Europe against those obtained from a Canonical Correlation Analysis (CCA) [START_REF] Härdle | Canonical Correlation Analysis[END_REF], which identifies the linear combinations of two variables with maximum correlation (Figures S7-S9). The spatio-temporal patterns issuing from the CCA largely resemble the ones found in Figures 1a and2-3, lending credence to the robustness of our results. Specifically, the CCA correlation coefficient (cc) peaks during winter, albeit with a more marked seasonal pattern (Figure S7) compared to α (Figure 1a). The cc extreme values (> 99 t h quantile) match positive NAO phases, mainly occurring during winter months (Figure S8), in very close agreement with Figure 2. Lastly, similar results are also obtained when replicating Figure 3, with the cc as a measure of dependence (Figure S9). The main difference is that in Figure S9a-b,e the signal is weaker than in Figure 4a-b,e,

suggesting that α provides a more nuanced picture of the extremes. Additional reasons motivating the use of the co-recurrence metric we propose here over a CCA approach are that α is linked to other dynamical properties of the phase space underlying the atmospheric motions, such as local dimension and persistence. It therefore indirectly provides a rich set of information on the system, including on its intrinsic predictability (Faranda et al., 2019b;[START_REF] Messori | A dynamical systems approach to studying midlatitude weather extremes[END_REF][START_REF] Scher | Predicting weather forecast uncertainty with machine learning[END_REF]. One may, for example, further investigate the predictability implications of the fact that concurrent precipitation-wind extremes over Europe display a low d and high θ -1 and compare this to the empirical results from ensemble reforecast data. Moreover, the definition of α can be easily extended to a multivariate case with more than two variables, whereas the CCA framework requires more complex adaptations (such as partial CCA) when working with more than two sets of variables.

The same approach could in the future be applied to a suite of renalysis products, to evaluating large ensembles of historical Atmosphere-Ocean General Circulation Models runs (e.g. [START_REF] Taylor | An Overview of CMIP5 and the Experiment Design[END_REF][START_REF] Eyring | Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization[END_REF], or to projecting future changes in compound extremes under a warming climate. Moreover, the dynamical systems approach used here can be implemented for any geographical region and variable (as long as the latter shows a chaotic behaviour).

A diverse set of compound and concurrent events can therefore be tested, and global hot-spots of compound extremes identified. There is also the hope that this methodology may help investigating outstanding scientific questions; for example the role of Arctic Amplification in driving northern hemisphere mid-latitude climate extremes (e.g. [START_REF] Cohen | Recent Arctic amplification and extreme mid-latitude weather[END_REF][START_REF] Barnes | Revisiting the evidence linking Arctic amplification to extreme weather in midlatitudes[END_REF][START_REF] Coumou | The influence of Arctic amplification on mid-latitude summer circulation[END_REF][START_REF] Screen | Amplified mid-latitude planetary waves favour particular regional weather extremes[END_REF][START_REF] Mori | Robust arctic sea-ice influence on the frequent eurasian cold winters in past decades[END_REF][START_REF] Ye | Two leading modes of wintertime atmospheric circulation drive the recent warm arctic-cold eurasia temperature pattern[END_REF].

The main findings of the study can be summarised as follows:

• The co-recurrence ratio α captures the seasonal cycle in extreme event occurrences for the domains (Eastern North America, Europe) and variables (precipitation, 10m wind gusts, 2m temperature) considered here.

• α extremes, here termed compound extremes, match same-day extremes in the climate variables, here termed concurrent extremes. High α values indicate a strong coupling between the different climate variables, which intuitively should be the case for spatially and/or temporally co-occurring extremes.

• For the European domain, the close correspondence between extremes in α and in the climate variables is mediated by the strong projection of α extremes on the positive phase of the NAO.

• For the concurrent North American and European extremes, α highlights regional extreme event hotspots. These are regions characterised by similar footprints of the extreme events across the different co-occurrence instances.

These findings highlight the ability of dynamical systems metrics to capture the spatio-temporal variability of the atmosphere. In Europe, our results mostly confirm the well-known link between a positive NAO and stormy weather over the north-eastern part of the continent, leading to widespread flooding linked with severe winds [START_REF] Fink | The European storm Kyrill in January 2007: synoptic evolution, meteorological impacts and some considerations with respect to climate change[END_REF][START_REF] Lavers | Winter floods in Britain are connected to atmospheric rivers[END_REF][START_REF] Huntingford | Potential influences on the United Kingdom's floods of winter 2013/14[END_REF][START_REF] Matthews | A cyclone climatology of the British-Irish Isles 1871-2012[END_REF][START_REF] De Luca | Extreme multi-basin flooding linked with extratropical cyclones[END_REF]. The physical relationship between cold temperatures in Eastern North America and wet extremes in Europe is not as obvious. The notably wet and stormy winter of 2013/2014 in parts of Europe and the British Isles (e.g. [START_REF] Huntingford | Potential influences on the United Kingdom's floods of winter 2013/14[END_REF][START_REF] Van Oldenborgh | Drivers of the 2013/14 winter floods in the UK[END_REF][START_REF] Wild | Was the Extreme Storm Season in Winter 2013/14 Over the North Atlantic and the United Kingdom Triggered by Changes in the West Pacific Warm Pool?[END_REF] co-occurred with a very cold winter in the Eastern United States [START_REF] Trenary | Was the cold eastern US winter of 2014 due tro Increased Variability?[END_REF][START_REF] Lee | Compounding effects of warm sea surface temperature and reduced sea ice on the extreme circulation over the extratropical North Pacific and North America during the 2013-2014 boreal winter[END_REF]. However, the longer-term climatology of concurrent winter cold spells in North America and wet and stormy events over western Europe has received comparatively less attention. [START_REF] Messori | On cold spells in North America and storminess in western Europe[END_REF] explicitly addressed the question of possible physical links between the two sets of extremes, and concluded that cold spells in North America are a systematic, but not necessary, precursor to stormy conditions in western Europe. The question of whether different instances of co-occurring extremes shared similar spatial footprints -of high relevance for investigating their predictability -was left largely unanswered. Our analysis highlights south-western Québec, south-eastern Ontario and part of the western European coastline from Portugal to France as hotspot regions for concurrent extremes with spatially similar characteristics over time. The European region roughly matches that highlighted in Messori et al.

(2016) (see their Fig. 4d). The North American region is, however, shifted to the north-east of the region analysed in a number of studies dealing with the winter of 2013/2014 and, more generally, eastern North American cold spells (e.g. [START_REF] Messori | On cold spells in North America and storminess in western Europe[END_REF][START_REF] Harnik | The circumglobal north american wave pattern and its relation to cold events in eastern north america[END_REF][START_REF] Overland | Potential arctic connections to eastern north american cold winters[END_REF][START_REF] Trenary | Was the cold eastern US winter of 2014 due tro Increased Variability?[END_REF]. We reiterate that our analysis specifically points to regions characterised by concurrent extremes with similar spatial footprints across the different co-recurrence instances, which may explain the discrepancy with some of the results in the literature. Finally, we note that the interpretation of the dynamical systems analysis may pose some challenges -as is for example the case of the α anomalies associated with precipitation considered in Sect. 5.2.

The spatial information on the extreme events intrinsic to our approach can thus provide new insights into multihazards (or compound events) research. Our dynamical systems perspective is flexible, links directly to concepts of predictability and coupling and may be useful in both practical (e.g. numerical forecasting of extremes) and more theoretical (e.g. understanding the atmospheric drivers of extremes) contexts.

AU T H O R C O N T R I B U T I O N

GM and PDL conceived and designed the study. PDL performed the analyses, created the figures and wrote the first manuscript draft. GM and DF provided the MATLAB code to compute the dynamical systems metrics and contributed to the methods. FP designed and performed the statistical tests. All the authors contributed to preparing the manuscript.

  For a given state of interest ζ = {ζ x , ζ y }, we have 0 ≤ α(ζ) ≤ 1. When α(ζ) = 0, no recurrences of ζ = {ζ x , ζ y } are observed in the phase-space when we observe a recurrence of ζ x . Intuitively, when α(ζ) = 1, all recurrences of ζ x correspond to recurrences of ζ = {ζ x , ζ y }, and viceversa. In other words, α quantifies the co-recurrences of two (or more) variables within a hyper-ball around ζ.

  Figure 4e differs from Figure 3e because of the different variables and domains over which α is computed here. Lastly, α anomaly means observed during concurrent cold and wet extremes show statistical significance in the north-eastern

F

  Daily climatological means of the precipitation-wind co-recurrence ratio (α). (b) Cumulative distribution functions (CDFs) of local dimension (d ) anomalies associated with α daily extremes (> 99 t h quantile); (c) as (b) but for local persistence (θ -1 ). In (b)-(c) the metrics are calculated for precipitation (brown), wind (yellow) and jointly for precipitation-wind (blue). Blue shaded areas in (a) represent one standard deviation from the mean. The data cover 1979-2018 over Europe. U R E 2 (a) Daily NAOI and (b) monthly counts for α daily extremes (> 99 t h quantile). In (a) the black solid line represents NAOI=0; the orange dashed line the mean NAOI value for α extremes; and the violet dot-dashed line the NAOI mean 97.5 t h quantile obtained by random sampling (n=10000, no replacement) from the entire NAOI daily time-series. The data cover 1979-2018 over Europe.

  Precipitation and (b) wind anomaly means associated with α extremes (> 99 t h quantile). (c) Monthly counts of concurrent (i.e. same-day) precipitation-wind extremes. α anomaly mean values observed during (d) precipitation, (e) wind and (f) concurrent precipitation-wind extremes (> 99 t h quantile). In (a)-(b) stippling represents statistically significant values. In (d)-(f) only statistically significant values are plotted and stippling represent gridpoints with >66% of daily α positive anomalies (sign test). Statistical significance is assessed via the Mann-Whitney test. As Figure 3 but for temperature in Eastern North America and precipitation in Europe. In (c)-(d) and (f) temperature extremes are daily values < 1 st quantile.

  . Given two trajectories x (t ) and y (t ) -in our analysis two different atmospheric variables -and a state ζ = {ζ x , ζ y }, we can define d x (ζ) and d y (ζ). These are the dimensions of the Poincaré sections defined by x and y around ζ, with respect to the chosen distance metric. We can further define joint logarithmic returns of ζ as:

	g (x (t ), y (t )) = -log dist	x (t ) x	,	ζ x x	2	+ dist	y (t ) y	,	y ζ y	2	1 2

Funding information

PDL was funded by CENTA NERC (Grant No. NE/L002493/1), Loughborough University and Vrije Universiteit Amsterdam. GM was partly supported by the Swedish research council Vetenskapsrådet (Grant. No. 2016-03724).

C O N FL I C T O F I N T E R E S T

The authors declare that they have no conflict of interest.


H J Van Zyl

A A Ilemobade

J E Van Zyl

Y Le Gat

Statistical characterisation and estimation of non-domestic water demand

Estimation of annual average water demand figures is critical for the design and evaluation of water distribution systems. This study evaluated the metered water consumption of more than 67,000 nondomestic consumers in six categories from cities and towns in South Africa. It was found that lognormal distributions provide good descriptions of the annual average daily demand (AADD) distribution in each category. The land use categories Business Commercial, Industrial, Agricultural holdings and Sports & Parks displayed similar median AADDs of between 1.5 and 1.7 kl/property/day. Educational properties used substantially more water (4.7 kl/property/day), while Government & Institutional properties used substantially less water (0.7 kl/property/day). A step-wise regression analyses showed that property size has the greatest impact on water demand for most categories. Finally, a novel statistically based method is proposed for estimating the average AADD of a given number of properties based on an acceptable risk of non-exceedance.

Introduction

Accurate estimation of annual average daily demand (AADD) is important for the design of municipal water distribution systems. AADD values normally form the basis for water supply component sizing including pipelines, storage reservoirs, pump stations and water treatment plants. AADD estimates are also important for the sizing of wastewater distribution and treatment systems.

It is useful to distinguish between domestic and non-domestic consumers in the analysis of water demand. Domestic demand is driven by human needs and behaviour, and can be estimated from historic data and factors such as development level, property size, income and climate.

In comparison, non-domestic demand is highly variable and is primarily driven by the type of activity practiced on the property. For example, a drinks factory and storage depot may have similar property sizes, turnovers and staff numbers, but will consume vastly different quantities of water.

Non-domestic consumers include industrial, commercial and institutional (ICI) users, parks and agricultural holdings (large properties often used for small scale agriculture).

While studies have found relationships between non-domestic demand and parameters such as property size, developed area, economics, water rates and water restrictions [START_REF] Mccuen | Forecasting urban water use: Commercial establishments[END_REF][START_REF] Kollar | The impact of industrial water use on public water supplies[END_REF][START_REF] Stephenson | Water demand patterns in Gauteng[END_REF][START_REF] Van Vuuren | Herevaluering van die Bestaande Riglyne vir Stedelike en Industriele Watervoorsiening gebaseer op Gemete Waterverbruike Fase 1: Pretoria Voorsieningsgebied [Re-evaluation of the existing guidelines for urban and industrial water supply based on measured water demands[END_REF], CSIR 2003[START_REF] Zhou | Water use in China's domestic, industrial and agricultural sectors: An empirical analysis[END_REF][START_REF] Reynaud | An econometric estimation of industrial water demand in France[END_REF]), these models require information at individual property level. They are less suitable for new areas where the types of non-domestic users that will be established are unknown, or existing areas where the types of developments may change over time.

In recent years, municipal data management tools have become common in South African municipalities, giving researchers access to a large quantity of metered consumption data (GLS 2007, Jacobs and[START_REF] Jacobs | Evaluating a tool to increase information processing capacity for consumer water meter data[END_REF]. This data formed the basis for a number of studies on water demand, which identified influencing parameters, investigated typical demand ranges and proposed new design guidelines for the estimation of domestic demands (e.g. [START_REF] Stephenson | Water demand patterns in Gauteng[END_REF][START_REF] Van Vuuren | Herevaluering van die Bestaande Riglyne vir Stedelike en Industriele Watervoorsiening gebaseer op Gemete Waterverbruike Fase 1: Pretoria Voorsieningsgebied [Re-evaluation of the existing guidelines for urban and industrial water supply based on measured water demands[END_REF][START_REF] Jacobs | Estimating residential water demand in Southern Africa[END_REF][START_REF] Husselmann | Effect of property size and income on residential water demand[END_REF][START_REF] Van Zyl | An improved area-based guideline for domestic water demand estimation in South Africa[END_REF][START_REF] Griffioen | Modelling water demand on a suburb level[END_REF], but did not consider non-domestic consumption.

In this study, metered consumption data of more than 67,000 non-domestic properties were linked to various databases and analysed in six categories. The aims of the study were to:

• Estimate the statistical properties and distributions of the data in each category. • Identify the main factors that influence non-domestic water demand in each category. • Develop a method to estimate the average consumption of a given number of non-domestic properties for a chosen level of confidence.

While detailed guidelines for the water demand of different non-domestic uses are often given in table form, these guidelines are only useful if the types and properties of these users The database was then used as the basis for analysing non-domestic demand characteristics and influencing factors.

Non-domestic demand characteristics

In the first phase of the study, the statistical characteristics of the non-domestic demand categories were calculated. It was found that log-normal distributions fit the data remarkably well as shown on the distribution density and cumulative distribution functions in Figures 1 and2 respectively. Although the lognormal distributions didn't fit the data perfectly (the Kolmogorov Smirnov test rejected normality in the data), they clearly provide reasonable descriptions of the data and are mostly conservative for larger demands.

Table 2 summarises the mean, median and lognormal distribution parameters of the different categories.

Due to the long high-end tail of lognormal distributions, the mean demands are considerably larger than the median demands in each category. The median demands for the different categories are compared in Figure 3. The figure shows that Business Commercial, Industrial, Agricultural holdings and Sports & Parks have remarkably similar median demands of between 1.5 and 1.7 kl/property/day.

Figure 2 confirms the similarity of the median values, and further shows that the distributions of Sports & Parks, Industrial and Business Commercial properties are also similar. While Agricultural holdings have a similar median demand as the other three categories, it has significantly lower variability.

The two outermost demand distributions are Educational on the higher side, and Government & Institutional users on the lower side, with median demands of 4.71 and 0.67 kl/property/ day respectively.

To remove the impact of property size on the demand distributions, the AADD of each property was divided by its size to obtain the unit demand (measured in L/m 2 /day). The cumulative distributions are shown in Figure 4 and the mean and median values in Figure 5. are known. However, this is not the case for new developments and even in existing developments the types of non-domestic consumers are not fixed, but change over time. Thus this study approached non-domestic consumption using a statistical approach, rather than at individual property level.

It is hoped that the result of this study will encourage similar research in other countries, allowing for comparison of non-domestic demand parameters and improved estimation of non-domestic demands.

Methodology

The data-set used in this study is the same as that used by [START_REF] Van Zyl | An improved area-based guideline for domestic water demand estimation in South Africa[END_REF] for a study on domestic water demand. It included water meter readings from forty-eight municipal treasury databases totalling more than 2.5 million records in four metropolitan areas (Johannesburg, Tshwane, Ekurhuleni and Cape Town) and 151 other cities and towns. The data is described in detail in [START_REF] Van Zyl | An improved area-based guideline for domestic water demand estimation in South Africa[END_REF].

Approximately 438,400 non-domestic consumer records were extracted from the database, and subjected to a number of verification tests to exclude anomalies. The remaining 67,291 records were then separated into six land use categories, namely:

• Business Commercial. • Industrial. • Government & Institutional.
• Education, including schools, universities and other educational buildings and their grounds. • Agricultural holdings, which are large properties serviced by municipal water distribution systems. Agricultural holdings may be used for small scale farming, although alternative sources of water would often be used for this purpose. These properties are also called 'small farms' or 'residential farms' . • Sports & Parks.

Information on possible factors influencing non-domestic demand was collected from various sources and linked to the demand data to allow correlations to be analysed. This information was obtained from municipal treasury databases (property area, property value and improvements, monthly meter readings), the South African Weather Service (rainfall, average daily minimum and maximum temperatures) and a manuscript on the surface water resources of South Africa [START_REF] Midgley | Surface water resources of South Africa[END_REF].

The AADD for each record in the database was estimated from twelve consecutive months of monthly consumption values. This means that seasonal demand variations did not influence the values, and errors due to differences in reading dates or erroneous reading would be small. The meter readings exclude demand from other sources, such as boreholes and grey water.

To ensure the integrity of the data, three data cleaning phases were implemented. In the first phase, records with inconsistent readings or dates (e.g. due to meter replacements or clock-overs) and records with less than 12 consecutive months of data were excluded. In the second phase, records or properties identified as unmetered, vacant, pre-paid, duplicate (including more than one type of land use) or not classified as one of the target non-domestic categories were excluded. In the third phase, stands with

Factors influencing non-domestic demand

To analyse the various factors influencing non-domestic water demand, a stepwise multiple regression analysis was performed on the data, using the independent variables listed in the Methodology section, and including the natural logarithms of the property area and property value parameters. Table 3 shows the number of statistically significant parameters identified, as well as the three most important parameters for each category.

While Government & Institutional users had the lowest consumption per property (Figure 3), Figure 5 shows that they had the highest median demand per unit area. The median demand per unit area of educational users ranks much lower compared to its demand per property, most likely due to the large stands with school playgrounds, parking and sports facilities that are common at these institutions.

An interesting result is that the demand per unit area of Agricultural holdings is by far the lowest of all the non-domestic categories. pointed out early by [START_REF] Watson | On quantiles of sums[END_REF], who investigated the relationship between the quantiles of a sum of two independent random variables and those of its components.

The problem was recently addressed by [START_REF] Klass | Uniformly accurate quantile bounds for sums of arbitrary independent random variables[END_REF], but the method seems to only provide loose bounds not practically usable for the purposes of this study. Extensive work is reported by Nadarajh (2008) in the case where the distribution of the sum is parametrically defined, but the results related to the log-normal distribution cannot be applied to the consumption data of this study, for which the lognormal fit appears to be skewed by the weight of the upper tail of the empirical distributions.

Thus it was decided to adopt a heuristic approach based on discarding the highest values of the tail of the distribution, and then fitting an approximate model to the remaining cumulative distribution. In practical terms this means that a design based Property size emerged as the most important factor in all categories except for Business Commercial and Agricultural holdings. For Business Commercial, property size was the second most important factor (after property value), and for Agricultural holdings, the third most important factor (after property value and temperature).

Modelling the mean demand of k stands

The final aim of the study was to develop a model to predict the mean demand of k stands of a category at a given level of confidence. The theoretical question of finding a high quantile for a given probability level (95% was used in this study) of a sum of a given number of random variables seems to have received little attention in the statistical literature in the case where the distribution of the random variables is unspecified. This fact was The demand threshold and model parameters determined for the three demand categories are summarised in Table 4, and the functions and data are shown in Figure 6 (for a certainty of 95% that they will not be exceeded).

Proposed design guideline

The models developed above provide a useful design guideline for Business Commercial, Industrial and Agricultural holdings developments where stands are grouped and the details of consumption are not know in advance.

The design mean demand for developments, based on a 95% certainty of non-exceedance, may be estimated using the equations:

Business commercial:

Industrial:

Agricultural holdings:

Where q k is the design AADD in kl/property/day and k the number of properties.

The design models, based on Monte Carlo sampling, are given in Figure 6. The models provide a good description of the demand behaviour, except for Commercial areas with less than six stands, where it is advisable to use the data points themselves.

As an example, consider a new industrial area with 50 stands. From Figure 6, the design demand can be estimated as 7.7 kl/ property/day, or a total AADD of 385 kl/day for the 50 stands. This is significantly higher than the category mean demand of (2)

q k = 3.806 + 12.534 √ k (3) q k = 5.179 + 17.853. √ k (4) q k = 2.276 + 4.674 √ k
on these values will not cater for properties with exceptionally high demands.

The model that was used has the form:

Where q k is the design demand per property, μ the expected value of the distribution, q 1 the chosen percentile of the data, and k the number of properties.

The choice of the analytical form of the equation was guided by the following considerations:

• Compliance with the limit condition q k = q 1 for k = 1, i.e.

when only one property is selected, the design value should be at the chosen percentile of the original data. • A decrease of the design value proportional to 1/ √ k, by analogy to the decrease in the standard deviation.

The threshold value was determined by trading off two objectives: the first to discard as few data points as possible, and the second to find a suitable fit between the model and data. To achieve this, a data-set was first truncated at a given demand threshold, and the 95 percentile of the demand then determined as a function of the number of properties through Monte Carlo sampling. The model parameters giving the best fit to the 95-percentile curve were then estimated.

This model was only applied to the Business Commercial, Industrial and Agricultural holdings categories since, unlike the (1) q k = + q 1 -√ k 5.2 kl/property/day, to accommodate the likelihood of a higher concentration of high-demand stands in the area.

Conclusion

This study analysed six categories of non-domestic demand in South Africa based on metered consumption data for a large number of stands. Basic statistical descriptors and distributions are reported and compared. The lognormal distribution was found to provide a good description of the demand distribution for all categories. A stepwise regression analysis was used to analyse various influencing factors, and property size was found to be the most important factor.

From the results obtained, guidelines for estimating water demand based on a 95% certainty of non-exceedance were developed for Business Commercial, Industrial and Agricultural property areas. It is hoped that this study will encourage comparative studies in other countries leading to a better understanding of the range of demand parameters and improved design guidelines for non-domestic water demand.

Figure 1 .

 1 Figure 1. Demand frequency distributions for the non-domestic categories a) business Commercial, b) Industrial, c) Government & Institutional, d) educational, e) agricultural holdings and f) Sports & Parks. In each graph the X-axis represents the natural logarithm of the aaDD, and the Y-axis the number of properties.

Figure 2 .

 2 Figure 2. Cumulative distribution functions with lognormal distributions for different non-domestic categories. Solid lines show the data and dotted lines the lognormal distributions.

Figure 3 .

 3 Figure 3. Median demands of non-domestic categories.

Figure 4 .

 4 Figure 4. Cumulative distributions of non-domestic demand per unit area.

Figure 5 .

 5 Figure 5. Mean and median demands per unit area for non-domestic categories.

Table 1 .

 1 ranges of property areas and values included in this study.

	Land use category	Property area (m 2 )	Property value (Euro)
	business Commercial	20-50,000	1400-7,000,000
	Industrial	20-50,000	1400-7,000,000
	Government & Institutional	100-50,000	1400-1,400,000
	educational	100-75,000	1400-1,400,000
	agricultural holdings	> 5000	1400-700,000
	Sport and Parks	500-100,000	1400-1,400,000

Table 2 .

 2 Statistical descriptors and lognormal distribution parameters of the different non-domestic demand categories.

					Lognormal distribution parameters
			Mean AADD	Median AADD	Mean Ln	Standard deviation
	Category	No. of data points	(kl/property/day)	(kl/property/day)	(kl/property/day)	Ln(kl/property/day)
	business Commercial	30,850	6.83	1.52	0.4392	1.592
	Industrial	10,851	13.99	1.72	0.5345	1.754
	Government & Institutional	12,731	5.07	0.67	-0.2634	1.399
	educational	2146	13.29	4.71	1.4251	1.714
	agricultural holdings	9612	3.56	1.56	0.3631	1.241
	Sports & Parks	1101	9.05	1.67	0.4044	2.033

Table 3 .

 3 Most important factors influencing non-domestic demand.

	Category	No. of data points	No. of significant parameters	Most important factors	Cumulative adjusted R 2
	business Commercial	24,810	8	ln(property value)	0.230
				ln(property area)	0.283
				location	0.303
	Industrial	8000	8	ln(property area)	0.313
				Property value	0.330
				location	0.345
	Government & Institutional	8259	8	ln(property area)	0.336
				ln(property value)	0.373
				Property value	0.387
	education	1447	5	ln(Property area)	0.408
				Development level	0.422
				Property value	0.437
	agricultural holdings	8895	6	ln(property value)	0.096
				ave minimum temperature	0.113
				Property area	0.118
	Sports & Parks	367	4	ln(Property area)	0.137
				ave minimum temperature	0.176
				Development level	0.191

Table 4 .

 4 Model parameters for the 95 percentile curves.

		Demand
		threshold	Fraction
	Category	(kl/prop-erty/day)	of data excluded	q 1 (kl/prop-erty/day)	μ(kl/prop-erty/day)
	business Com-	50	1.7%	16.340	3.806
	mercial
	Industrial	80	1.8%	23.032	5.179
	agricultural	20	1.3%	6.950	2.276
	holdings

Figure 6. Proposed design demand guidelines for business Commercial, Industrial and agricultural holdings properties.

Acknowledgement

The authors would like to acknowledge the water demand data supplied by GLS Consulting.

Funding

Financial support for the project was provided by the Water Research Commission (Project K5/1525).

Disclosure statement

No potential conflict of interest was reported by the authors.


Ulrich D Jentschura

Jean Zinn-Justin

Calculation of the Characteristic Functions of Anharmonic Oscillators

Keywords: Singular perturbations, turning point theory, WKB methods, Semiclassical techniques including WKB and Maslov methods, General quantum mechanics and problems of quantization

The energy levels of quantum systems are determined by quantization conditions. For one-dimensional anharmonic oscillators, one can transform the Schrödinger equation into a Riccati form, i.e., in terms of the logarithmic derivative of the wave function. A perturbative expansion of the logarithmic derivative of the wave function can easily be obtained. The Bohr-Sommerfeld quantization condition can be expressed in terms of a contour integral around the poles of the logarithmic derivative. Its functional form is B m (E, g) = n + 1 2 , where B is a characteristic function of the anharmonic oscillator of degree m, E is the resonance energy, and g is the coupling constant. A recursive scheme can be devised which facilitates the evaluation of higher-order Wentzel-Kramers-Brioullin (WKB) approximants. The WKB expansion of the logarithmic derivative of the wave function has a cut in the tunneling region. The contour integral about the tunneling region yields the instanton action plus corrections, summarized in a second characteristic function A m (E, g). The evaluation of A m (E, g) by the method of asymptotic matching is discussed for the case of the cubic oscillator of degree m = 3.

MSC numbers 34E20, 81Q20, 81S99

Introduction

The cubic anharmonic oscillator, defined by the Hamiltonian,

h 3 (g) = - 1 2
∂ 2 ∂q 2 + 1 2 q 2 + √ g q 3 , (1)
is a paradigmatic example of a quantum mechanical problem which gives rise to complex resonance energies. A particle initially trapped in the region q ≈ 0 may tunnel through the classically forbidden well and escape toward q → -∞. For positive coupling g > 0, the energies can be determined numerically by the method of complex scaling [START_REF] Balslev | Spectral properties of many-body Schrödinger operators with dilation-analytic interactions[END_REF][START_REF] Yaris | Resonance calculations for arbitrary potentials[END_REF]. One scales the coordinate as q → q e i θ for the cubic oscillator, which results in the Hamiltonian

H c (θ) = e -2iθ - 1 2
∂ 2 ∂q 2 + 1 2
q 2 e 4iθ + √ g q 3 e 5iθ .

(2)

If one diagonalizes this complex operator in the basis of harmonic oscillator wavefunctions {φ n (q)} Nmax n=0 , for large enough N max , one can numerically determine the (complex) resonance energies of the original cubic Hamiltonian [START_REF] Alvarez | Coupling-constant behavior of the resonances of the cubic anharmonic oscillator[END_REF]. As discussed in Refs. [START_REF] Alvarez | Coupling-constant behavior of the resonances of the cubic anharmonic oscillator[END_REF][START_REF] Yaris | Resonance calculations for arbitrary potentials[END_REF], these resonance energies are independent of θ, provided we choose θ sufficiently large so that the rotated branch passes the position of the resonance under investigation. The numerical variation of the resonance energies of the cubic oscillator, ǫ (3) n (g) = Re ǫ (3) n (g) + i Im ǫ (3) n (g) ,

as a function of g under a suitable increase of N max can be used to estimate the numerical uncertainty of the numerical results for the complex resonance energies [START_REF] Jentschura | Structure, time propagation and dissipative terms for resonances[END_REF]. Throughout this article, we keep the superscript "3" of the resonance energy ǫ

n (g) of the cubic oscillator in order to remain consistent with the notation of a brief Letter (see Ref. [START_REF] Jentschura | Unified Treatment of Even and Odd Anharmonic Oscillators of Arbitrary Degree[END_REF]), in which part of the results presented here have already been indicated.

For reasons which will explained below, we associate the resonance energy (with a negative imaginary part) with a value of the coupling parameter g that has an infinitesimal positive imaginary part, whereas a coupling with an infinitesimal negative imaginary part is associated with an antiresonance, that has a positive imaginary part. The resonance energies of the cubic potential as a function of g have a branch cut along the positive real axis, and fulfill the dispersion relation [START_REF] Bender | Large-order perturbation theory for a non-Hermitian PT -symmetric Hamiltonian[END_REF][START_REF] Bender | Numerical evidence that the perturbation expansion for a non-Hermitian PTsymmetric Hamiltonian is Stieltjes[END_REF],

ǫ (M) n (g) = n + 1 2 + g π ∞ 0 ds Im ǫ (M) n (s + i 0) s (s -g) . ( 4
)
We note that the spectrum of resonances is invariant under the transformation √ g → -√ g in Eq. ( 1), which can be compensated by a parity transformation q → -q (the latter leaves the spectrum manifestly invariant).

For weak positive coupling g = |g| → 0 (with Im g > 0), the imaginary part of the resonance energies is determined by a nonanalytic factor, which reads [START_REF] Jentschura | Structure, time propagation and dissipative terms for resonances[END_REF] Im ǫ

(3) 0 (g + i 0) ∼ -

2 3n n! √ π g -n-1/2 exp - 2 15 g for g = |g| → 0 , Im g > 0 . (5)
For large positive coupling g = |g| → ∞ + i 0, the first three resonance energies possess the following asymptotics [see Eq. ( 17) and Table II of Ref. [START_REF] Jentschura | Structure, time propagation and dissipative terms for resonances[END_REF]], ǫ

0 (g + i 0) ∼ g 1/5 (0.617 160 050 -i 0.448 393 023)

for g = |g| → ∞ + i 0 , ǫ (3)
1 (g + i 0) ∼ g 1/5 (2.193 309 731 -i 1.593 532 797)

for g = |g| → ∞ + i 0 , ǫ (3)
2 (g + i 0) ∼ g 1/5 (4.036 380 020 -i 2.932 601 744)

for g = |g| → ∞ + i 0 . (6)
All of the above results have a complex argument of -π/5, as they should. For intermediate g, a large number of resonance energies can be determined numerically by a suitable basis-set method and be used in order to construct a complex-scaled time propagation algorithm for a quantum mechanical wave packet (see Sec. III of Ref. [START_REF] Jentschura | Structure, time propagation and dissipative terms for resonances[END_REF]).

A second available method for the calculation of resonance energies consists in the resummation of the divergent perturbation theory in complex directions of the parameters. The perturbative coefficients ǫ

n,K in the perturbation series

ǫ (3) n (g) ∼ ∞ K=0 ǫ (3) n,K g K , g → 0 , (7)
which describe the nth resonance energy of the cubic oscillator, diverge factorially for large perturbation theory order K. Indeed, for large K, the leading asymptotics read (see, e.g., Ref. [START_REF] Jentschura | Unified Treatment of Even and Odd Anharmonic Oscillators of Arbitrary Degree[END_REF]), ǫ

n,K ∼ -

2 2n-1/2-K 15 n+K+1/2 π 3/2 n! Γ n + K + 1 2 , K → ∞ . (8)
With the exception of the coefficient of order zero, which reads ǫ

(3) n,0 = n+ 1/2, all ǫ (3)
n,K are negative for K ≥ 1. For g > 0, the series (7) is formally not Borel summable. However, numerical evidence [START_REF] Caliceti | From useful algorithms for slowly convergent series to physical predictions based on divergent perturbative expansions[END_REF][START_REF] Franceschini | Complex energies from real perturbation series for the LoSurdo-Stark effect in Hydrogen by Borel-Padé approximants[END_REF][START_REF] Jentschura | Resummation of nonalternating divergent perturbative expansions[END_REF][START_REF] Jentschura | Structure, time propagation and dissipative terms for resonances[END_REF] suggests that the Borel-Padé summation method with the Laplace-Borel integration being carried out along complex contours, as given by Eqs. (196)-(198) of Ref. [START_REF] Caliceti | From useful algorithms for slowly convergent series to physical predictions based on divergent perturbative expansions[END_REF], provides useful numerical approximants to the resonance energies. The method has been put on rigorous mathematical grounds recently [START_REF] Caliceti | Distributional Borel summability of odd anharmonic oscillators[END_REF] in the framework of distributional Borel summability [START_REF] Caliceti | The Distributional Borel Summability and the Large Coupling φ 4 Lattice Fields[END_REF][START_REF] Caliceti | Stark Resonances: Asymptotics and Distributional Borel Sum[END_REF]. When the integration contour C +1 in the conventions of Ref. [START_REF] Caliceti | From useful algorithms for slowly convergent series to physical predictions based on divergent perturbative expansions[END_REF] is employed, a negative sign is obtained for the imaginary part of the resonance energy. The accuracy obtainable using complex Borel summation of the weak-coupling expansion [START_REF] Caliceti | Perturbation theory of odd anharmonic oscillators[END_REF] is restricted by numerical oscillations in the values of the complex Borel-Padé transforms as the order of the Padé transformation is increased. These oscillations cannot be overcome even if extended-precision arithmetic is used in intermediate steps of the calculation. These oscillations have been described in Refs. [START_REF] Jentschura | Structure, time propagation and dissipative terms for resonances[END_REF][START_REF] Surzhykov | Quantum dot potentials: Symanzik scaling, resurgent expansions and quantum dynamics[END_REF] for the complex Borel-Padé resummation method. Indeed, it has been observed in [START_REF] Jentschura | Structure, time propagation and dissipative terms for resonances[END_REF] that, at a relatively moderate coupling g = 0.6, the ground-state energy of the cubic oscillator as determined by resummation cannot be calculated to better accuracy than

ǫ (3) 0 (g = 0.6) = 0.554(1) -0.351(6) i (9)
by Borel resummation. The oscillations, which are described by the numerical uncertainty in the last digit reported in result of Eq. ( 9), cannot be overcome when the (Borel) transformation order is increased and appear to represent a fundamental limit of the convergence of resummed weak-coupling perturbation theory in the case of a moderate and large (modulus of the) coupling parameter g. A more accurate result obtained by complex scaling (2) for the same coupling parameter is ǫ

The applicability of the weak-coupling expansion for the cubic oscillator resonance energies thus is found to be severely limited in the domain of moderate and large coupling, even if the perturbative expansion is a posteriori enhanced by a summation method acting in the complex plane. Likewise, the leading-order analytic result (5) cannot be applied to the domain of moderate and large coupling g. In this article, we explore the question whether one can find higher-order analytic formulas for the corrections to the leading result [START_REF] Bender | Numerical evidence that the perturbation expansion for a non-Hermitian PTsymmetric Hamiltonian is Stieltjes[END_REF] for the decay width of the ground-state resonance energy (and, possibly, an arbitrary excited resonance energy) of the cubic oscillator, and if yes, how large the correction terms are.

Indeed, in Ref. [START_REF] Jentschura | Unified Treatment of Even and Odd Anharmonic Oscillators of Arbitrary Degree[END_REF], an analysis has been presented which allows us to calculate higher-order correction terms to the resonance energies for even and odd anharmonic oscillators of arbitrary degree, and to formulate modified Bohr-Sommerfeld quantization conditions for anharmonic oscillators such as the cubic one. Here, we provide a few more details of the analysis that leads to the results of Ref. [START_REF] Jentschura | Unified Treatment of Even and Odd Anharmonic Oscillators of Arbitrary Degree[END_REF], with a special emphasis on the cubic potential. To this end, we review in Sec. 2 the calculation of higher-order perturbative (in g) approximants to the wave function for the cubic potential. In Sec. 3, the calculation of higher-order corrections to the WKB approximation of the wave function is reviewed. The modified quantization condition pertaining to the cubic oscillator is described in Sec. 4, with a few new results for higher-order terms reserved for Sec. 5. Finally, conclusions are drawn in Sec. 6.

Perturbative expansion

The Schrödinger equation corresponding to the Hamiltonian given in Eq. ( 1) reads

- 1 2
∂ 2 ∂q 2 + 1 2 q 2 + √ g q 3 ϕ = E ϕ . ( 11
)
If we transform the coordinate according to q → q/ √ g, we obtain

- g 2 2 ∂ ∂q 2 + V(q) ϕ = g E ϕ , V(q) = 1 2 q 2 + q 3 . ( 12
)
We may transform to the Riccati equation by setting

ϕ ′ (q) ϕ(q) = - s(q) g , ϕ ′′ (q) ϕ(q) = s 2 (q) g 2 - s ′ (q) g , (13)
The Riccati form of the Schrödinger equation then reads,

g s ′ (q) -s 2 (q) + U 2 (q) = 0 , U(q) = 2 [ V(q) -g E ] . (14)
The zeroth-order term is

s(q) ≈ s 0 (q) = u(q) = 2 V(q) . (15)
The following recursion yields higher-order approximants (in g) for the logarithmic derivative of the wave function,

s(q) = ∞ K=0 g K s K (q) , s 0 (q) = u(q) , s 1 (q) = u ′ (q) -2E 2u(q) , s K (q) = 1 2 u(q) s ′ K-1 (q) - K-1 l=1 s K-l (q) s l (q) , K ≥ 2. ( 16
)
We can divide s(q) into a symmetric and an antisymmetric component under the simultaneous interchange g → -g, E → -E,

s(q) = s + (q) + s -(q) = s + (q, E, g) + s -(q, E, g) , s ± (q, -E, -g) = ±s ± (q, E, g) . ( 17
)
From the following system of equations,

g s ′ -(q) -s + (q) 2 -s -(q) 2 + U 2 (q) = 0 , (18a)
g s ′ + (q) -2 s + (q) s -(q) = 0 , (18b)
we can eliminate the antisymmetric component by setting

s -(q) = g 2 s ′ + (q) s + (q) . ( 19
)
The Bohr-Sommerfeld quantization condition reads

- 1 2πi g C dz s + (z) = n + 1 2 , ( 20
)
where C is a contour that encloses the zeros of the function s + (z) in the anticlockwise (mathematically positive) direction, and n is the quantum number of the level. The function B 3 (E, g), defined by

B 3 (E, g) ≡ - 1 2πi g C dz s + (z) , (21)
[we recall the implicit dependence of s + (z) on g and E] provides for a universal means of determining the perturbative expansion for an arbitrary excited level of the cubic potential, by means of the perturbative quantization condition (

B 3 (E, G) = n + 1 2 . A
) 22

WKB expansion and contour integral

The WKB expansion is an expansion in g at g E fixed, i.e., it implies large values of E. We use this expansion in order to derive, roughly speaking, higher-order corrections to the instanton action which can later be used in order to derive higher-order corrections to resonance energies E of levels whose principal quantum number is not large (this is based on a rather fortunate enhancement of the applicability regions of certain expansions, as explained below). The instanton configuration for the cubic potential has been discussed in Ref. [START_REF] Jentschura | Generalized nonanalytic expansions, PT -symmetry and largeorder formulas for odd anharmonic oscillators[END_REF], with a graphical representation being given in Fig. 3(b) of [START_REF] Jentschura | Generalized nonanalytic expansions, PT -symmetry and largeorder formulas for odd anharmonic oscillators[END_REF]. In order to achieve an instanton configuration pertaining to positive values of the coordinates, we scale the coordinate as q → -q/ √ g. We start from

- 1 2 ∂ 2 ∂q 2 + 1 2 q 2 + g q 3 φ = E φ (23)
and transform according to q → -q/ √ g, which gives

- g 2 2 ∂ ∂q 2 + W(q) φ = g E φ , W(q) = 1 2 q 2 -q 3 . ( 24
)
We note that W(q) differs from V(q) in the sign of the cubic term. Setting

φ ′ φ = - S(q) g , (25)
where the logarithmic derivative is denoted S(q) instead of s(q), we obtain

g S ′ (q) -S 2 (q) + T 2 (q) = 0 , T (q) = 2 [ W(q) -g E ] . (26)
The zeroth-order term is different as compared to the perturbative expansion,

S(q) ≈ S 0 (q) = 2 [ W(q) -g E ] . (27)
A recursion in g can be defined,

S(q) = ∞ K=0 g K S K (q) , S 0 (q) = T (q) , S 1 (q) = S ′ 0 (q) 2 S 0 (q) , S K (q) = 1 2 S 0 (q) S ′ K-1 (q) - K-1 l=1 S K-l (q) S l (q) . ( 28
)
Like the perturbative expansion, the WKB expansion of the logarithmic derivative of the wave function can be calculated recursively in ascending powers of g, and we can also distinguish between even and odd components under the symmetry operator g → -g, E → -E, S(q) = S + (q) + S -(q) = S + (q, E, g) + S -(q, E, g) , S ± (q, -E, -g) = ± S ± (q, E, g) .

The contour integral about the cut of the function S + (q) leads to a more complex structure as compared to the right-hand side of Eq. ( 20). It reads

1 g C dz S + (q) = A 3 (E, g) + 1 2 ln(2π) -ln Γ 1 2 -B 3 (E, g) + B 3 (E, g) ln - g 8 . (30)
Here, C is a contour that encloses the cut of the WKB expansion of the logarithmic derivative of the wave function in the clockwise sense. We thus have to calculate the integral

T = 1 g C dq S + (q) ≈ 1 g C dq S 0 (q) + g 2 S 2 (q) + g 4 S 4 (q) = T 0 + T 2 + T 4 , (31)
Here, the WKB approximants S 0 (q), S 2 (q) and S 4 (q) are defined in Eq. ( 28), and the integrals T 0 , T 2 and T 4 correspond to the contour integrals of the WKB approximants of the respective order. The zeroth-order term is

T 0 = 1 g C dq S 0 (q) (32)
which involves the zeroth-order WKB approximant S 0 (q) = 2W(q) -gE = q 2 -2q 3 -gE .

(33)

If we define the square root function to have its branch cut along the positive real axis, then the above expression has its branch cut (in the limit g E → 0) from q = 0 to q = 1 2 . Directly above the cut, the value of the square root is positive, while directly below, it is negative. As the contour encircles the cut in the clockwise direction, we have an integration interval from zero to 1 2 above the real axis, while below the real axis, we go from 1 2 to zero. With this definition, the contour integral is

T 0 = 1 g C dq q 2 -2q 3 -gE = 2 g Re 1/2 0 dq (q 2 -2q 3 -gE) 1/2 . ( 34
)
For g E > 0, the cut of the WKB approximant does not extend along the full interval (0, 1 2 ) any more. However, by specifying the real part of the final result, we pick up only those terms which are due to the cut.

The method of asymptotic matching (see Ch. 7.4 of Ref. [START_REF] Bender | Advanced Mathematical Methods for Scientists and Engineers[END_REF]) can be used in order to evaluate (34). It involves an overlapping parameter ǫ, that is also used in Lamb shift calculations [START_REF] Feynman | Space-Time Approach to Quantum Electrodynamics[END_REF][START_REF] Jentschura | Higher-order binding corrections to the Lamb shift of 2P states[END_REF][START_REF] Pachucki | Higher-Order Binding Corrections to the Lamb Shift[END_REF] in order to separate the low-energy from the high-energy contribution to the bound-electron self-energy. The overlapping parameter fulfills

0 < g E ≪ ǫ . ( 35
)
We then separate the integration interval into three parts, the first of which extends from (0, ǫ), the second, from (ǫ, 1 2 -ǫ), and the third complementing interval finally reads ( 1 2 -ǫ, 1 2 ). The method of evaluation depends on the particular integration interval. For the interval (0, ǫ), we have to be careful to avoid divergences for very small q = 0. We have to keep the expression q 2 -2g E within S 0 (q) in unexpanded form. However, we can expand the expression q 2 (1 -2q) -2g E = q 2 -2g E -2q 3 in the q 3 term, then do the q integration, subsequently expand the result in g E, and finally in ǫ up to order ǫ 0 , keeping all divergent terms [START_REF] Feynman | Space-Time Approach to Quantum Electrodynamics[END_REF][START_REF] Jentschura | Higher-order binding corrections to the Lamb shift of 2P states[END_REF][START_REF] Pachucki | Higher-Order Binding Corrections to the Lamb Shift[END_REF]. This leads to the result (

I 1 = Re ǫ 0 dq (q 2 -2q 3 -gE) 1/2 = (gE) - 1
) 36
For the interval (ǫ, 1 2 -ǫ), we can expand in g E to any power, because all divergences near q = 0 or q = 1/2 are cut off by the overlapping ǫ parameter. The result reads

I 2 = Re 1/2-ǫ ǫ dq (q 2 -2q 3 -gE)
15015 128ǫ 2 - 45045 64ǫ - 2 √ 2 ǫ 5/2 - 140 √ 2 ǫ 3/2 - 1120 √ 2 ǫ 1/2 + 255255 256 ln ǫ 2 4 + O(g E) 5 . ( 37
)
The result from the interval I 3 compensates the divergences of fractional order in ǫ from I 2 ,

I 3 = Re 1/2 1/2-ǫ dq (q 2 -2q 3 -gE) 1/2 = (gE) 2 2 √ 2 ǫ 1/2 + (gE) 3 4 √ 2 3ǫ 3/2 + 40 √ 2 ǫ 1/2 + (gE) 4 2 √ 2 ǫ 5/2 + 140 √ 2 3ǫ 3/2 + 1120 √ 2 ǫ 1/2 + O(g E) 5 . ( 38
)
In the result for T 0 ,

T 0 = 2 g (I 1 + I 2 + I 3 ) , (39)
the overlapping parameter ǫ cancels, and we obtain (

T 0 = 2 15 g + -E + E
) 43
The perturbative counterterms read

P = 1 2 ln(2π) -ln Γ 1 2 -B 3 (E, g) + B 3 (E, g) ln - g 8 = 7 2880 E 3 - 1 24 E + -E + E ln gE 8 + g - 7
) 44
where the result of B 3 (E, g) from Eq. ( 22) has been used, as well as the asymptotic expansion

ln Γ ( 1 2 + z) = z {ln(z) -1} + 1 2 ln(2π) - 1 24 z + 7 2880 z 3 + . . . . ( 45
)
Note that the logarithm ln(gE/8) is obtained as the combination of B 3 (E, g) ln(-g/8) ≈ E ln(-g/8) and from the asymptotic expansion of the expression -ln Γ ( 1 2 -B 3 (E, g)) ≈ E ln(-E). Finally, the instanton function A 3 (E, g) for the cubic potential is found to read

A 3 (E, g) = T -P = 2
This result for the instanton function is the basis for the calculation of higher-order corrections to the decay width of the resonance energies, as will be explained in the following.

Modified quantization condition

For g < 0, the interaction Hamiltonian √ g q 3 = ±i | √ g| q 3 is PT -symmetric, and its spectrum is purely real. In view of Eqs. ( 20) and ( 21), we then have the purely perturbative quantization condition

B 3 (E, g) = n + 1 2 ⇔ 1 Γ 1 2 -B 3 (E, g) = 0 , g < 0 . (48)
Let us now study how this condition needs to be modified in case instantons exist, i.e., for g > 0. We start from the unperturbed harmonic oscillator

H = - 1 2 ∂ 2 ∂q 2 + 1 2 q 2 . (49)
In this form, we naturally identify

1 Γ ( 1 2 -E) = det(H -E) (50)
as the spectral determinant. The effect of the instanton is to add an infinitesimal imaginary part to the energy of the bound state,

E ≈ E ≡ n + 1 2 + i Im E . (51)
Indeed, the definition of resonance energies via complex scaling and the proof of their existence can be found in Ref. [START_REF] Caliceti | Perturbation theory of odd anharmonic oscillators[END_REF]. Expanding in the nonperturbatively small (in g) imaginary part of the energy, we obtain

1 Γ ( 1 2 -E) ≈ -(-1) n+1/2 n! Im E . (52)
For the cubic oscillator, we approximate the right-hand side of Eq. (52) as

Im E ≈ Im ǫ (M) n (g) ≈ - 2 3n n! √ π g -n-1/2 exp - 2 15 g ≈ - 1 n! √ 8π 8 g B3(E,g) exp (-A 3 (E, g)) . (53)
We now generalize the left-hand side of Eq. ( 52) as 1/Γ ( 1 2 -B 3 (E, g)) and approximate the right-hand side of Eq. ( 52) by the right-hand side of Eq. (53), and obtain

1 Γ 1 2 -B 3 (E, g) = 1 √ 8π - 8 g B3(E,g) exp[-A 3 (E, g)] , g > 0 , (54)
which is our conjecture for the resonance energies of the cubic anharmonic oscillator in the "unstable" regime of positive coupling parameter g > 0. When the logarithm of both sides of Eq. ( 54) is taken, we recover the structure of the terms on the right-hand side of Eq. (30).

In Sec. 3, we have calculated the function A 3 (E, g) in the WKB limit g → 0, with g E fixed. In order to derive the result (46) for the function A 3 (E, g), we had to expand the function B 3 (E, g) in the limit of large E, which implies [according to Eq. ( 45)] the approximation -ln Γ ( 1 2 -B 3 (E, g)) ≈ E ln Γ (-E). A priori, this expansion is not applicable when E ≈ n + 1 2 + δ, where δ summarizes the perturbative and nonperturbative corrections (in g). We have calculated the instanton function A 3 (E, g) in an "unphysical" domain of large quantum numbers n (large E), "away" from the first few resonance energies. However, the result (46) is written in terms of an expansion whose terms actually decrease for typical values of E for the first few resonance energies, where E is close to a small half-integer. We can thus make an-in some sense fortunate-observation: While the expansion (45) for the logarithm of the Γ function is not applicable to the physical domain of E = n + 1 2 + δ, the obtained expansion for the A 3 (E, g) function is applicable in both the "unphysical" WKB domain of g E fixed, and g → 0, thus E → ∞, as well as the physical domain of E ≈ n + 1 2 , and g → 0. We can thus use the modified quantization condition (54) in the regime of low principal quantum numbers, which is of prime interest. This unified result, which gives immediate and systematic access to the higher-order corrections, should be compared to other approaches (Refs. [START_REF] Hamprecht | Tunneling amplitudes from perturbation expansions[END_REF][START_REF] Karrlein | Precise variational tunneling rates for anharmonic oscillator with g < 0[END_REF][START_REF] Voros | Airy function-exact WKB results for potentials of odd degree[END_REF][START_REF] Voros | Exact resolution method for general 1D polynomial Schröodinger equation[END_REF], which also strive to go beyond the simple leading-order results for the decay rates of the resonances. One approach [START_REF] Hamprecht | Tunneling amplitudes from perturbation expansions[END_REF][START_REF] Karrlein | Precise variational tunneling rates for anharmonic oscillator with g < 0[END_REF] uses variationally improved perturbation theory which leads to a numerical improvement over the simple summation of the first few perturbative terms in the strong-coupling domain, but does not lead to a systematic prescription for calculation of higher-order coefficients in the "resurgent" expansion. The other approach [START_REF] Voros | Airy function-exact WKB results for potentials of odd degree[END_REF][START_REF] Voros | Exact resolution method for general 1D polynomial Schröodinger equation[END_REF] is based on a set of quantization conditions derived as variations of the Bohr-Sommerfeld condition.

Results

Based on the results presented in Eqs. ( 22), ( 46) and (54), it is easy to derive higher-order corrections to the imaginary part of the first few resonance energies of the cubic potential, by entering into the condition (54) with an ansatz Im ǫ (3) n (g) = -

2 3n n! √ π g -n-1/2 exp - 2 15 g 1 + a n g + b n g 2 + . . . , for g > 0 , (55)
and equating coefficients. The higher-order coefficients, for the lowest four states, read

Im ǫ (3) 0 (g) = - exp -2 15
Im ǫ

(3)

3 (g) = - 256 e -2/(15 g) 3 √ π g 7/2 1 - 3913 16 g + 8807869 512 g 2 - 1571666861 40960 g 3 - 3214761534593 2621440 g 4 + O(g 5 ) (56d)
For the coefficient of relative order g, a general result (for any n) has been indicated in Eq. ( 22) of Ref. [START_REF] Alvarez | Coupling-constant behavior of the resonances of the cubic anharmonic oscillator[END_REF].

For n = 2, Eq. ( 22) of Ref. [START_REF] Alvarez | Coupling-constant behavior of the resonances of the cubic anharmonic oscillator[END_REF] gives a correction term of -2041g/16 instead of -2101g/16 for the state with n = 2, and a correction term of -10543g/48 instead of -3913g/16 for the state with n = 3. With the help of the dispersion relation (4), these coefficients can be related to the corrections to the leading factorial asymptotics (8) of the perturbative coefficients of relative order 1/K. We have checked our calculation against numerical values of the perturbative coefficients of up to 60th order, which can easily be determined on the basis of the relation

B 3 (E, g) = n + 1/2.

Conclusions

The calculation of higher-order corrections to the decay widths of quantum states in unstable potentials is a challenging problem. Analytic approximations are useful in that regard, if they contain a sufficient number of correction terms which also permit an estimate of the truncation error. As shown in Secs. 4 and 5, the calculation of higher-order corrections to the decay widths of resonances of the cubic anharmonic oscillator, is based on two characteristic functions, namely, the perturbative function B 3 (E, g) = E + O(g) defined in Eq. ( 21) and evaluated in Eq. ( 22) (57)

and on the instanton function A 3 (E, g) = 2 15 g + O(g) defined in Eq. ( 30) and evaluated in Eq. ( 46). ( 58)

The first of these functions is given by an evaluation of a contour integral about the poles of the perturbative approximation to the logarithmic derivative of the wave function (see Sec. 2). The second function, A 3 (E, g), is determined by a contour integral about the cut of the WKB approximation of the wave function (see Sec. 3).

The modified Bohr-Sommerfeld quantization condition (4) then gives us access to the higher-order corrections to the imaginary part of the resonance energies. These are given, for the first few resonances of the cubic potential, in Sec. 5. The full solution of the modified quantization condition has the structure of a generalized nonanalytic expansion, or "resurgent" expansion,

ǫ (3) n (g) = ∞ J=0   i n! √ 8π 8 g n+ 1 2 exp - 2 15 g   J Lmax L=0 ln L - 8 g ∞ K=0 Ξ (3,n) J,L,K g K , g > 0 , (59)
where L max = max(0, J -1) .

The perturbative coefficients ǫ 1,0,1 multiply the correction of relative order g to the one-instanton (J = 1) effect and are given in Eq. (56) for n = 0, 1, 2, 3. They grow quite drastically in magnitude with the principal quantum number. In the notation of Eq. (55), we have for the state with n = 3 a coefficient of a 3 = -244.563. This means that the correction term of order g halves the total result for the decay width of the third excited state of the cubic potential already at a minuscule coupling parameter of g ≈ 0.002. It is quite interesting to see that numerically approximations to such a fundamental physical quantity as the resonance energies of the cubic anharmonic oscillator require a considerable effort in their evaluation. Here, our aim was to carefully explain the calculational approach necessary in order to gain access to the higher-order corrections, which are numerically significant.

  = 0.6) = 0.554 053 519 -0.351 401 778 i .

  K coefficients. The coefficients entering Eq. (56) are the Ξ

Acknowledgments

The authors gratefully acknowledge support from the National Science Foundation (Grant PHY-8555454) and from the Missouri Research Board.