Legacy effects of past thinnings modulate drought stress reactions at present

ABSTRACT Thinning experiments were primarily established for analysing how treatment variants determine the amount and quality of wood production. Given climate change, they may also explain how silvicultural treatment determined drought resistance. Especially for tree species cultivated in regions beyond their natural range, silvicultural treatment may help mitigate drought stress. Here, we used the 47-years-old combined spacing-thinning trial Fürstenfeldbruck 612 and metrics for quantifying the trees’ recent and past growth to test if the information of tree treatment and development in the past significantly improved the prediction of their growth at present and if spacing and density regulation, kind of thinning, and temporal sequence of thinning significantly co-determined tree and stand growth during drought. Based on linear models, we revealed the following ecological legacy effects: (i) information of tree treatment and development in the past significantly improved the prediction of their growth at present, and (ii) higher densities, past thinnings from below, and low variations in thinning strength were beneficial for the tree and stand growth during drought. Thus, the prevailing repeated strong thinnings from above for promoting a selected collective of future crop trees may be questioned because of climate change.


Introduction
Compared with the past, the present generation of trees and stands in Central Europe developed with less competition and growing more open due to stand density reductions by natural disturbances (Sommerfeld et al. 2018;Senf and Seidl 2021) or stronger thinnings (Mund et al. 2002;Bosela et al. 2016). Therefore, they have a different legacy when exposed to stress such as drought or pest infestations. Whereas in the past, thinning experiments were mainly focused on the amount and quality of the produced wood (Zeide 2002(Zeide , 2001, there is currently an increasing interest in how to make stands more resistant to drought to ensure their survival by spacing and thinning (Lagergren and Lindroth 2004;Kohler et al. 2010;Giuggiola et al. 2013;Elkin et al. 2015). At the beginning of classical sustainable forestry, stands were mostly kept dense (Fernow 1913). This was primarily a backlash to forest devastation by over-exploitation in the industrial age (Wiedemann 1923;Kreutzer 1972;Knapp et al. 2015). The maximisation of mass production (Yaffee 1999;Zeide 2001) and the protection of the internal climate of the stands (Blume et al. 2010) were the further main arguments for keeping stands on high-density levels by leaving them unthinned or just moderate thinning from below (Zeide 2001;Pretzsch 2020).
During the ensuing period, stands were increasing stronger thinned mainly to promote the size growth and quality of a selected number of final crop trees (Schädelin 1942). The extensification of forest management by reduction of harvesting costs (Soman et al. 2019), the forest infrastructure provision by forest roads and log trails (Weintraub and Navon 1976), advanced harvesting techniques and machines (Sirén and Aaltio 2003;Eberhard and Hasenauer 2021), and the improvement of individual tree stability (Cremer et al. 1982;Gardiner et al. 1997) further promoted the tendency towards lower stand densities. Since extreme climatic events such as droughts are predicted to become more frequent, more severe and longer-lasting (Intergovernmental Panel on Climate Change (IPCC), 2021), the focus of forest management and research is increasingly shifting to tackling these upcoming challenges. Amongst others, recent studies recommend stand density reductions to reduce drought stress and growth losses or tree mortality as a consequence (Martín-Benito et al. 2008;Sohn et al. 2016b;Aldea et al. 2017).
Thinning is put to the test again after this long history and turnaround (Zeide 2001(Zeide , 2002. Empirical findings raise the question of whether repeated abrupt strong thinnings may reduce internal and external tree stability (Brüchert et al. 2000;Marchi et al. 2017). The stress gradient hypotheses provide the theoretical background for a trade-off between growth and stress resistance (Callaway and Walker 1997;Holmgren et al. 1997) and suggest that an acceleration of growth by competition release, fertilisation, or nitrogen deposition may increase the growth rate but also the susceptibility to biotic damages.
Intense thinning may be detrimental to the stand climate (Thom et al. 2020) and promote forest floor vegetation's growth and water consumption (Hedwall et al. 2013;Trentini et al. 2017). Drought stress damage can be most severe on the more large trees with wide annual rings and a low basic density (Persson 1994).
The beneficial shading effect on smaller trees in years with low rainfall and strong solarisation may be reduced by stand density reduction . The inter-individual shading may be interrupted, especially in the case of smaller crowns with low self-shading effects; sunburn and bark beetle attack may be promoted (Nicolai 1986;Jactel et al. 2009;Björkman et al. 2015;Corcobado et al. 2020;Hlásny et al. 2021).
The interest in the relationship between past and future growth originates from the desire to mitigate drought, and other climate change-induced stresses through silvicultural measures (Ogle et al. 2015;Camarero et al. 2018). For example, initially wide spacing or repeated heavy thinning leads to large xylem vessels, which improve water conduction and growth under normal conditions and may improve drought resistance (Martín-Benito et al. 2008;Sohn et al. 2016b;Aldea et al. 2017). However, this may also increase the risk of embolism during drought. Conversely, other studies have found opposite effects, where density reduction had positive effects in the short term, followed by adverse effects in the long term (D'Amato et al. 2013;Calama et al. 2019;Steckel et al. 2020). These discrepancies may be clarified when analysing the past tree-ring structure. They could help re-evaluate and improve stand density regulation, as Sohn et al. (2016b) recommended.
Knowledge of the legacy effects of stand establishment, spacing, and past thinnings on drought stress reactions is highly relevant for improving stand management and developing silvicultural guidelines under climate change. It is well known that current growth is strongly determined by the present environmental conditions, competitive status, size and crown structure (e.g. Forrester 2019; Forrester et al. 2018;Grote et al. 2016;Pretzsch 2021a). Thus, our overarching question was how tree growth under drought stress is modulated by different spacing, kind, strength, and frequency of thinning in the past. Here, we were especially interested in the co-determination of tree growth by the past treatment, represented by metrics such as variation of growth and competitive status in the past, kind of thinning, the strength of competition release in the past, and temporal sequence of thinning.
We based our study on the combined spacing and thinning experiment in Norway spruce (Picea abies [L.] KARST.) Fürstenfeldbruck 612 (FFB 612). This experiment comprises 21 plots with a wide range of initial spacing, variants of density reduction, and kinds of thinning. The stands are stocked on sites representing excellent growing conditions in standard years but growth limitation by water in dry years. The stands were established in 1974 and surveyed seven times by 2020. Trees were 51 years old at the last survey, and their growth represents a period with several intense drought events (1976,2003,2015,2018). At the first survey in 1991, there were 4614 trees on the plots, and at the last survey in 2020, the experiment comprised 928 trees. The tree and stand data of this experiment was appropriate for scrutiny of the following hypotheses: H1: Stand growth during drought is independent of the silvicultural treatment of the stands in the past.
H2: Metrics characterizing the establishment, spacing, and silvicultural treatment of the trees and their growth in the past significantly improve the prediction of their growth at present. H3: Spacing and density regulation in the past co-determine tree growth during drought.
H4: Kind of thinning (from above vs from below) in the past codetermine tree growth during drought.
H5: Temporal sequence of thinning (regularity or irregularity) in the past co-determine tree growth during drought.

Location
The combined spacing-thinning experiment Fürstenfeldbruck 612 is part of an international series of experiments initiated by the International Union of Forest Research Organizations (IUFRO, working group 1.05.05) in 1967 to analyse the effects of different silvicultural treatments on the tree and stand growth of N. spruce under different site conditions. Here, we explain the location and the essentials of the experimental design of FFB 612 that we used for this study; similar experiments were established in several European countries. For a more detailed description, see Abetz (1981).
FFB 612 lies in an N. spruce plantation that was 51 years old at the last survey in spring 2021. The stand is located 40 km northeast of Munich (48°14' N, 11°05' E) in South Bavaria, Germany, at 550 m a.s.l. The mean annual temperature is 7.5°C, the annual precipitation is 825 mm, and the respective characteristics in the growing season are 15.2°C and 380 mm. The stand belongs to the ecoregion "13 Schwäbisch-Bayerischen Schotterplatten-und Altmoränenlandschaft" and is located in the ecodistrict "13.5 Landsberger Altmoräne" (Walentowski 2004). Stands are stocked on a luvisol that developed from a loess layer above a base moraine of the Würm glaciation. The natural vegetation would be a Galio-Fagetum, i.e. European beech dominated forest ecosystems with Norway spruce interspersed only occasionally (Walentowski 2004).
The site has experienced pronounced environmental changes since 1970 ( Figure 1). The mean temperature increased by approximately 2°C, whereas the annual precipitation sum only slightly decreased from 1970 to 2020 (DWD Climate Data Center 2021a, 2021b). Both trends result in an unfavourable progression of soil water availability, as indicated by the soil moisture index (Figure 1c; Zink et al. 2016Zink et al. ). 2003Zink et al. , 2015Zink et al. , and 2018 were exceptionally warm and dry years (Figure 1, broken vertical lines).

Stand description
The experiment was established in 1974 with 4-year old plants (2/2 transplanted) on an area of 3.36 ha with altogether 21 plots (900 m 2 each) assigned to two blocks ( Figure 2) on similar sites and in close neighbourhoods. The experimental factors are spacing and thinning. Each of the two blocks represents plots with different initial spacing ( Figure 2, framed with thin lines) and plots with the same spacing of 2.5m × 1.6m and 2500 trees ha −1 , but different thinnings (0, Z1-Z4). The initial spacing varies between 1m × 1m and 5m × 5m, i.e. a tree number of 400−10.000 ha −1 . No treatment of the plots took place between 1974 and 1991. Then the factor thinning varied between unthinned and variants Z1-Z4 (Table 1). The thinning on the plots was based on set curves of tree number depending on top height and a crop tree thinning from above, thinning from below or a crop tree thinning combined with thinning from below. The curves arrive at 500 trees per hectare at a top height of 27.5 m, and the velocity of density reduction differed as follows 0 < Z1 < Z2 < Z3 < Z4 (see Table 1). The thinning on the spacing trials was similar to that on the thinning trials (see Table 1).

Measurements
At the first survey after the growth period in 1991, all tree coordinates were measured as, in some cases, there were minor deviations of the tree positions from the spacing scheme due to, e.g. stumps or stones. Since its establishment in 1974, the stands have been inventoried seven times until now. The first time in spring 1992, followed by surveys in 1996, 2001, 2006, 2012, 2017 (all autumn) and 2021 (spring). Inventories comprised the measurement of all stem diameters (d 1.3 in cm) by tape measurement. On each plot, at least 30 trees distributed over the stem diameter range were selected for measuring tree height (h in m) and height to the crown base (hcb in m) using a Haglöf Vertex (Haglöf, Långsele, Sweden). The survey in 2001 comprised measurement of eight crown radii per tree cardinal directions (N, NE, E, SE, S, SW, W, NW) by crown mirror according to standards described by Pretzsch (2009) and Röhle and Huber (1985). Trees were logged only in years with surveys, except for trees that died naturally due to wind, bark beetle or competition. Those trees were removed due to pest control. Because plots were only treated in survey years, it was possible to observe detailed changes in stand composition.
Soil moisture Index (SMI): The Soil Moisture Index (SMI) is an indicator to quantify agricultural droughts. The SMI is simulated by the hydrologic model mHM using current weather data. Model mHM consists of a digital elevation model, geological map, soil map and land cover and leaf area information derived from satellite data (see Samaniego et al. (2013) for further explanation of the model). The model provides 30-day soil moisture indices SMI Marx 2017). The range of the values is SMI = 0−1. A value of SMI < 0.2 is considered a drought. The index is primarily used to provide information on the current moisture status of the soil. Applied mainly in climate research, the index is used for national drought monitoring in Germany (Samaniego et al. 2013). Furthermore, the SMI has been successfully related to the tree-ring series in dendroecological studies to identify drought events (Schwarz et al. 2020;Uhl et al. 2021).

Tree level evaluations
By modelling the relationship between height, stem diameter, and age, we calculated the height of each tree. To estimate the individual tree height (h) depending on the stem diameter (d) and tree age (age), we parameterised the model below using all available measurements of tree heights, stem diameters, and tree ages.
All regression coefficients and the model as a whole were significant, at least at the level of p , 0.001 (n = 4372, R 2 = 0.91). For the model parameters, see Supplementary Table 1.
From the crown projection measurements on the 21 plots, we used the eight radii to calculate the mean crown radius cr = (r 2 1 + r 2 2 + . . . + r 2 8 )/8 and crown diameter cd = 2 × cr. To estimate the individual crown diameter (cd) depending on stem diameter, we parameterised the model below using all available measurements of tree crown diameter and stem diameters.
We derived the allometric relationship cd = exp(−0.48 + 0.56 × ln(d)). It reflects the quadratic mean crown diameter (cd in m) concerning the stem Figure 2. Experimental design of the combined spacing-thinning experiment Fürstenfeldbruck 612 (FFB 612) with 21 plots represented in two blocks (left and right). Both blocks include plots that vary in initial spacing (1.0m × 1.0m, 1.25m × 2.0m, 1.6m × 2.5m, 2.5m × 4.0m, 5.0m × 5.0m) and are framed by thin lines. They also have plots with equal spacing (2.5 m ×1.6 m) but different thinning regimes (0, Z1-Z4) framed by bold lines. See Table 1 for details of the different spacing and thinning variants. Table 1. Description of the spacing and thinning variants on the 21 plots of the combined spacing-thinning experiment Fürstenfeldbruck 612 (FFB 612). The plots in rows 1-5 differ in spacing and thinning; the plots in rows 6-10 have the same initial spacing but differ in thinning. Thinning on the plots in rows 1-5 and 7-10 applied tree number reduction following set curves for tree number-top height (h O ) relationships and c = crop tree thinning, b = thinning from below, or c, b = crop tree thinning combined with thinning from below. Plots in row 6 remained unthinned. Thinning of plots in row 10 followed defined tree volume removal. diameter (d in cm). It was used for adding missing crown measurements.
For analysing the effects of the competition, competition release by thinning, kind and frequency of thinning, and structure on the stem diameter growth of each tree, we constructed a circle with a radius sr = 0.25 × h last survey around its standpoint; with h last survey being the height of the central tree at the last survey. All trees within the constructed circle were used to quantify local competition, competition release by thinning, kind of thinning, and diversity in structure. In the constructed circles, there were, on average, 9-12 trees and at least 5-6 most impactful neighbours (Prodan 1968) in the last survey in 2020. This study and preliminary studies (Pretzsch and Schütze 2021;Pretzsch 2022aPretzsch , 2022b showed that this circle size resulted in the highest correlations between stand structure characteristics and growth. To quantify the local density and competitive status of each tree, we used the Hegyi-index (Hegyi 1974), ci, representing a distance dependent type of competition index as well as the local stand density index (Reineke 1933), sdi, as a distance independent competition index. Note that sdi refers to the stand density at the individual tree level, whereas SDI was the density calculated at the stand level (see section 2.4.2 Stand level evaluations).
For quantifying local sdi, all trees within the search radius, sr, except the central tree were used to calculate the local density n on circle area a. N = 10.000/a × n was the respective tree number upscaled to one hectare. For the n trees, we calculated the quadratic mean stem diameter d q ; based on N and d q we then calculated the local density sdi = N × (25/d q ) around each tree. The local sdi was calculated using the species-specific allometric exponents of = −1.664 derived by Pretzsch and Biber (2005). The used exponent = −1.664 deviated from the species-overarching exponent of −1.605, as proposed by Reineke (1933). Note that this exponent α was derived on unthinned and A-grade plots of long-term experiments in South Germany, located in the same area as FFB 612.
In contrast to the ci by Hegyi, the local sdi does not weigh the competition by neighbours by their distance to the central tree. However, the sdi is easier to interpret as it represents the number of trees per hectare with an index diameter of 25 cm. In contrast, the index by Hegyi provides a dimensionless relative index.
The ci and sdi values were calculated with and without the removed trees for all circles and surveys; the relationships reported in the result section were based on the ci and sdi values of the remaining stand at the end of each survey period.
Before calculating the local sdi values and Hegyi-index for neighbourhood analysis, we established a toroidal shift of the plot to all eight directions of the plot periphery for edge bias compensation (Radtke and Burkhart 1998;Pommerening and Stoyan 2006;Pretzsch 2009); see Supplementary Figure 1. This was necessary as, in several cases, the search radius reached even beyond the buffer zone around the individual plots ( Figure 2). The plot edges were located in the middle distance between the planting rows. By the toroidal shift, we extended the same planting pattern and distances in all eight directions and avoided any overdensity as it can be produced by other techniques such as reflection via a reflecting line through the edge trees; and a random arrangement of interior trees around the plot (Radtke and Burkhart 1998).

Metrics for tree level
(a) stem and crown size, and competitive status d, h, hd, cl, cr, cd, cv: Current stem diameter, d, tree height, h, slenderness ratio, hd, crown length, cl, crown ratio, cr, crown diameter, cd, and crown volume, cv, are classical predictors for the tree growth (Table 2).
ci before , ci after , Dci, sdi before , sdi after , Dsdi: The Hegyi-index before and after thinning, ci before and ci after , and the competition release Dci = ci before − ci after before a growth period in question are also frequently used predictors. The sdi before , sdi after , and Dsdi analogously describe the individual trees competition and competition release. In the cases of ci before , ci after , sdi before , and sdi after values of zero mean open grown conditions, and the higher the ci or sdi, the stronger the competition. strudiv: The coefficient of variation of the stem diameters in the vicinity of each tree, strudiv, was calculated for characterisation of the structural composition of the competitors. The calculation was based on the trees within the search radius also used for calculating the ci and local sdi; the respective central tree was excluded from calculating the coefficient of variation of the stem diameters. The higher strudiv, the higher the structural diversity of the stand within the vicinity of the respective central tree.
(b) current growth rate and metrics of the tree development in the past id, id mean , id cvar , id rg : This part lists the periodical annual diameter growth in the last survey period, id, the mean stem diameter growth over all previous survey periods, id mean , and the coefficient of variation, id cvar , of the diameter increment over all previous periods since stand establishment. The relative range of annual diameter increment was calculated as id rg = (id max − id min )/id mean , with id max and id min as the maximum and minimum diameter increment, respectively, and id mean as the mean diameter growth over all previous surveys.
(c) metrics for the treatment of the trees This set of metrics describes the trees' establishment and treatment history.
sdi trend : Thinnings may be temporally staggered. e.g. first strong and later moderate or vice versa. To characterise this sequence of thinning strength, we divided the considered time period by two. i.e. in the first and second half of the observation period, p 1 and p 2 . We then calculated the mean sdi for the second period SDI p2 and divided it by the overall mean of the sdi in the whole observation time, sdi p1,p2 . The ratio sdi trend = sdi p2 /sdi p1,p2 is sdi trend = 1 if the stand density is similar in both periods (e.g. moderate thinning according to (Schober 1988a(Schober , 1988b) throughout the whole rotation). A value of sdi trend . 1 indicates an increasing stand density, i.e, a staggered thinning first with decreasing density reductions later (e.g. staggered thinning by Assmann and Franz (1965)). A value of sdi trend , 1 indicates an increasing stand density management with high density first and stronger density reductions in the later phase (e.g. increment thinning according to Schädelin (1942)).
ci cvar , sdi cvar : For calculating the respective coefficient of variation, we used all ci and sdi indices before and after thinning to calculate the standard deviation and mean of the competition. Low ci cvar and sdi cvar values indicate a continuous development of competition without abrupt change, whereas high values indicate substantial variation due to repeated strong density reductions.
ci mean , sdi mean : They reflect the mean level of density; low stand densities around the individual trees result in low means, high densities in high ci mean and sdi mean values. ci rg , sdi rg : Analogously to the relative range of the stem annual diameter increment, id rg , these two values represent cr cr = (r 2 1 + r 2 2 + . . . + r 2 8 )/8 crown diameter cd cr × 2 crown projection area cpa cr 2 × p crown ratio cratio cl/h, relative length of the crown crown volume cv cv = cl × cpa search radius for neighbourhood analysis sr sr = 0.25 × h last survey for analysing local competition before thinning ci before ; sdi before local competition in the circle before thinning local competition before thinning ci after ; sdi after local competition in the circle after thinning competition release by thinning Dci; Dsdi competition removed by thinning kind of thinning kindth kindth = d removed /d in the vicinity of each tree diversity in structure strudiv coefficient of variation of stem diameters in the vicinity of each tree (b) current growth rate and metrics of the tree development in the past periodic annual mean stem diameter growth id periodical diameter increment/period length mean stem diameter growth id mean id over all previous survey periods coefficient of variation of stem diameter increment id cvar id cvar = sd(id)/id mean maximum stem diameter increment id max id max = max(id) minimum stem diameter increment id min id min = min(id) the relative range of diameter increment id rg id rg = (id max − id min )/id mean (c) metrics for the treatment of the trees initial growing area initfl initial growing area of each tree mean competition ci mean ; sdi mean ci before ci after over all previous survey periods; same for sdi values coefficient of variation of competition ci cvar ; sdi cvar ci cvar = sd(ci before , ci after )/ci mean ; same for sdi values maximum competition ci max ; sdi max ci max = max(ci before , ci after ); same for sdi values minimum competition ci min ; sdi min ci min = min(ci before , ci after ); same for sdi values the relative range of competition ci rg ; sdi rg ci rg = (ci max − ci min )/ci mean ; same for sdi values mean competition release by thinning the range of competition variation; the more substantial the variation of the thinning, the higher the two metrics. Dci cvar , Dsdi cvar : These metrics indicate the variation in thinning strength; the higher the values, the stronger the variation of density reductions by the repeated thinnings.
kindth: To characterise the kind of thinning, we calculated the ratio between the mean stem diameter of the trees removed by thinning and the stem diameter of the central tree in question. We used all removed trees within the respective search radius to calculate the mean stem diameter of the removed trees. High kindth values (≥ 1) indicate thinning from above. Low values (, 1) thinning from below. kindth = 1 means that the removed trees had the same size as the central tree, indicating a thinning from above.
ci rg , sdi rg : Both values reflect the relative range of the density reductions by thinning. We calculated the highest and lowest values (ci max , ci min ) of density reduction within the tree development since stand establishment and divided the range by the mean ci value over the whole observation period, ci rg = (ci max − ci min )/ci mean . Analogously, we calculated sdi rg . Low values of ci rg and sdi rg indicate low peaks and a regular sequence of thinning, high values abrupt rises and peaks and irregularity of thinning interventions; occasionally, strong thinnings after longer pauses of interventions result in high values of ci rg and sdi rg .

Stand level evaluation
To overview the 21 plots, we also derived common dendrometric stand characteristics (Table 3). The stand level characteristics were derived from the successive inventories of the tree diameters, tree heights, and records of the removed trees. We used standard evaluation methods according to the DESER-norm recommended by the German Association of Forest Research Institutes (in German "Deutscher Verband Forstlicher Forschungsanstalten"; Johann 1993; Biber 2013). For estimating the merchantable stem volume in dependence on tree diameter, tree height and form factor, we used the approach by Franz et al. (1973) with the stem form equations and coefficients published by Pretzsch (2002, p. 170, Table 7.3). The results encompassed all common stand characteristics, such as the quadratic mean diameter and height of the mean and dominant trees, the stand basal area, standing stem volume, stand stem volume growth, and total yield; all variables were calculated for the total, remaining and removal stand. The Stand Density index, SDI, was calculated analogously to the local sdi at the tree level (see section 2.4.1 Metrics for tree level); whereas at the tree level, the central tree was excluded, the calculation of the SDI at the stand level included all trees. The stand data was used for introducing the stand and for analysing how past silvicultural treatments affect stand growth under drought.

Stand level
To test how drought periods have affected stand volume growth, we have applied the following simple linear mixed model. (3) For further testing, if density regulation or temporal sequence of thinning in the past co-determine stand volume growth during drought periods (H I), we extended Model 3 with covariates that characterise the silvicultural treatment of the stands in the past (see Table 2). This procedure resulted in the following models: ln (IV ik ) = a 0 + a 1 × ln (V remain,ik ) + a 2 × SMI k + a 3 × SDI mean, ik + a 4 × SMI k × SDI mean,ik + b i + e ik (4) ln (IV ik ) = a 0 + a 1 × ln (V remain,ik ) + a 2 × SMI k + a 3 In case of the Models 3-5 the independent variables were the standing stem volume of the remaining stand at the beginning of the survey period, V remain , soil moisture index, SMI, mean stand density index over all previous survey periods, SDI mean , and the coefficient of variation of density reduction by thinning over all previous survey periods, DSDI cvar . The indexes i and k represent the k th observation of the i th plot.
The fixed effects were covered by the parameters a 0 -a n . With the random effect b i N(0, t 2 ) we cover the correlation between the single observations on the plot level. With 1 ik we denote the independently and identically distributed errors.

Individual tree level
To test if metrics characterising the establishment, spacing, and silvicultural treatment of the trees and their growth in the past significantly improve the prediction of their growth at present (H II), we applied linear models to the periodic annual mean stem increment data of the last survey period (2018-2020). The fixed effect variables, such as stem diameter, tree state variables, metrics for silvicultural treatment, and metrics for the trees' past growth (see Table 2), represented the influence of the trees' present and past characteristics on its growth; the parameters covered the fixed effects by a 0 -a n . Model selection from the extensive models was carried out with a principal component analysis (PCA, see supplement Figure 2) and further supported by testing all possible mathematical models using all combinations of variables by the Akaike information criterion (AIC; Akaike 1981).
We used the following model as a baseline.
ln(id 2018−2020,i ) = a 0 + a 1 × ln (d i ) It represents the effect of the initial tree diameter (d) in spring 2018 only on the mean annual growth (id) in the 2018-2020 period.
We have extended this baseline model with tree state variables in a second step.
ln (id 2018−2020,i This model represents the effect of initial size and tree state variables on growth id during the last survey period. In a further step, we have added metrics for silvicultural treatment to model 7, which resulted in the following model. ln (id 2018−2020,i ) =a 0 + a 1 × ln (d i ) + a 2 × sdi after, i + a 3 × sdi trend, i + a 4 × Dsdi cvar, i + a 5 × sdi mean,i + a 6 × kindth mean,i + a 7 This model represents the effect of initial size, tree state variables, and metrics for silvicultural treatment on growth during the last survey period. In the last step, we also included metrics of the tree development in the past. The best model that took all metrics into account was the following model.
In addition to current tree characteristics and silvicultural treatment, this model considered metrics of tree development in the past. In case of the Models 6-9, the independent variables were the individual tree diameter, d, local stand density index after thinning, sdi after , coefficient of variation of the stem diameters in the vicinity of each tree, strudiv, competition release by tree removal, Dsdi, sequence of thinning strength, sdi trend , mean competition index over all previous survey periods, sdi mean , coefficient of variation of competition release by thinning over all previous survey periods, Dsdi cvar , mean kind of thinning over all previous survey periods, kindth mean , and the coefficient of variation of stem diameter increment over all previous survey periods, id cvar . In equations (6)-(9), the index i represents the i th tree.
For testing, if spacing and density regulation, kind of thinning, or temporal sequence of thinning in the past co-determine tree growth during drought periods (H III -H V), we applied linear mixed models as follows. (12) We chose double-logarithmic relationships in all cases as they appeared biologically more plausible than linear relationships and resulted in higher R 2 values. In case of the Models 10-12, the fixed effects cover the main effects of individual tree diameter, d, mean competition index over all previous survey periods, sdi mean , mean kind of thinning over all previous survey periods characterised by "thinning from below" (values < 1) and "thinning from above" (values ≥ 1), kindth mean , coefficient of variation of competition release by thinning over all previous survey periods, Dsdi cvar , soil moisture index, SMI, and all their two-way interactions. When fitting the model, non-significant interactions were removed, and the model was re-fitted. Still, if the interaction was significant, the contributing main effects were kept in the model even when not significant, following a protocol suggested by Zuur et al. (2009). The indexes i and k represent the k th observation of the i th tree. The fixed effects were covered by the parameters a 0 -a n . With the random effect b i N(0, t 2 ) we cover the correlation between the single observations on the tree level. To not eliminate the effects of spacing and density regulation, kind of thinning, or temporal sequence of thinning in the past, we did not implement random effects at the plot level. Random effects at the plot level would have eliminated the effects on tree growth which we wanted to analyse. With 1 ik we denote the independently and identically distributed errors. All fitted models were subject to the usual visual residual diagnostics. For all models, the residuals were plotted against the fitted values. In no case the plots suggested a violation of variance homogeneity. Likewise, the normality of errors was verified by making normal q-q plots of the residuals.
All data processing, including the toroidal shift and analyses, was conducted using the statistical software R

Stand characteristics
The dendrometric stand characteristics (Table 3) show that tree sizes differed considerably on the different plots, although all trees had the same age (e.g. quadratic mean stem diameter 26.16-50.39 cm). As a result of the different initial spacing and thinning frequency and intensity, the tree number, basal area, or standing volume of plots without any silvicultural density regulation (see Table 1) had considerably higher values than the most sparsely stocked plots (e.g. 340.54 vs 908.41 m 3 ha −1 ). The average mean annual volume increment since stand establishment, MAIV, on the plots was 17.14 m 3 ha −1 year −1 , but varied between plots from 9.17-21.27 m 3 ha −1 year −1 .
The periodic mean annual stand volume growth plotted over the calendar year revealed a clear trend (Figure 3, a-e). The lower the initial density in the initial stand, the lower the mean annual stand volume growth. Irrespective of the initial density, however, all plots had in common that they showed a drop in volume increment in the survey period 2001-2006, which included the extremely dry year 2003 ( Figure 3e). In all plots, growth peaked in the second to last recording period (2012-2017; age 43-48 years) and showed a downward trend in the last survey period (2017-2021), which was also characterised by several dry years (see Figure 1).
Tree characteristics Table 4 lists three categories of variables that characterise (a) the trees' size and competitive status at present, (b) growth rate, and (c) neighbourhood and treatment in the past that were analysed regarding their effect on the tree growth at present. It reflects the wide range of tree state variables and stem diameter increment in the 21 plots. At 51 years, the tallest trees achieve 63.3 cm stem diameter, 34 m tree height, 22.4 m crown length, and 6.39 m crown diameter. In contrast, the smallest trees of the same age achieve 10.2 cm stem diameter, 18.7 m tree height, 1.0 m crown length, and 2.3 m crown diameter.
All other variables and metrics, such as competition status, current growth rates and metrics on the development of the trees in the past, and metrics for the treatment of the trees in the past, also showed apparent differences between the individual trees.
More than 16,000 stem records from the repeated surveys and their variation in size, tree development, and competitive state provide a solid database for the subsequent evaluation of the tree growth depending on trees' metrics for silvicultural treatment and past growth. The periodic mean annual stem diameter growth plotted over age revealed a clear trend (Figure 4, a-e). N. spruce from wide initial spacing showed a unimodal course of growth with a peak at a tree age of 27 years (Figure 4, d and e). Trees from plots with an initial spacing of 1.6m × 2.5m tended to show a degressive trend in growth (Figure 4, c). The trees had high growth rates at the beginning and lower growth rates at the end. In contrast, N. spruces from plots with an initial spacing of 1.25m × 2.0m and 1.0m × 1.0m showed a relatively constant growth (Figure 4, a and b). At least the average growth did  Table 4. Overview of the tree characteristics (n = 928) at the last survey of FFB 612 in 2020. Shown are the mean, standard deviation, minimum, and maximum values of (a) stem and crown size and competitive status, (b) current growth rate and metrics of the tree development in the past, and (c) metrics for the treatment of the trees (explanation of variables and metrics see  Figure 1), was evident in all cases. Additionally to the course of the mean stem diameter growth depending on age (Figure 3), Supplement Figure 2 shows all individual trees' growth versus age and calendar year. The shown growth trends may be co-determined by the plots' stand structure and silvicultural treatment. To avoid respective biases, we considered those effects in our analyses.

Co-determination of silvicultural treatment in the past on stand growth during drought (H I)
The statistical characteristics of the stand volume growth prediction with and without the inclusion of information of stand silvicultural treatments in the past (Models 3-5) are shown in Table 5. All three models showed a significant positive effect of the soil moisture index (SMI) on stand volume growth (Figure 5a). This means that stand growth was significantly lower in dry periods than in years with sufficient water supply. The mean stand density over all previous survey periods (SDI mean ) had a significant positive effect on stand growth (Table 5, Model 4). The higher the density, the higher was the stand volume growth. It is interesting to note that the interaction between SMI and SDI mean showed a significantly negative effect. However, although we found a significant positive effect of higher stand densities on stand growth during dry periods, the effect size is small (Figure 5b). Consideration of the coefficient of variation of density reduction by past thinnings (DSDI cvar ) on stand volume growth (Table 5, Model 5) showed a significant negative effect of the main effect DSDI cvar . Of high importance, however, is our finding that the interaction between SMI and DSDI cvar was significantly positive (Table 5, Model 5). This indicates that high variations in thinning intensity had only a minor effect on stand volume growth in years with sufficient water availability. In dry periods, however, high variations in thinning intensity in the past had a strong negative influence on stand volume growth ( Figure  5c). The high values of the Bayesian factor (BF) from the model comparisons with Model 3, which did not include any information on past silvicultural treatment, in the denominator, suggest strong evidence that the silvicultural treatment of the stands in the past had a strong influence on the growth of the investigated stands during drought periods.

Contribution of the past to estimation of current tree growth (H II)
The variable correlation plot (Supplement Figure 2) shows that the variables describing the tree dimension (e.g. d, cpa, cl) described in section 2.4.1 (Table 2) are all highly positively correlated with each other. The mean diameter increment and the initial stand area are also positively correlated with these tree dimension variables. On opposite sides, and thus negatively correlated with the tree dimension variables, were variables assessing the competitive status of the trees (e.g. sdi before , sdi after ). The second group of highly positively related variables were those describing the temporal sequence of thinning (sdi cvar , Dsdi cvar , sdi rg Dsdi rg , etc.). These variables, which describe the steering of the trees, were positively correlated with the within-tree variables id rg or id cvar . A strong negative correlation existed between id rg or id cvar with the current diameter increment id. The statistical characteristics of the growth prediction with and without the inclusion of information about the trees' past are demonstrated using Models 6 and 9 (Tables 6 and 7). For the characteristics of Models 7 and 8, the details are presented in Supplementary Tables 2 and 3. Supplementary  Table 4 summarises the stepwise model improvement for all four models, starting with baseline Model 6 and progressing to Model 9. It is shown that all three model characteristics, AIC, R 2 , and RMSE, can be improved by several percentages compared with the baseline model. The high Bayesian factor (BF) values from the model comparisons with the baseline Model 6 in the denominator showed that we had found clear evidence that considering metrics that account for a tree's past silvicultural treatment and past growth helps to improve the prediction of stem diameter growth further.
Co-determination of silvicultural treatment in the past on tree growth during drought

Spacing and density regulation (H III)
Our results on the influence of the soil moisture index, SMI, and spacing and density regulation in the past on the periodic annual mean stem growth, id, showed a clear trend  (Table 8). We found significant positive main effects of stem diameter, d, and SMI. While the mean local stand density index over all previous survey periods, sdi mean , had no significant (p = 0.0919) influence on id, there was a significant negative interaction of SMI with sdi mean . Regardless of the diameter of the trees, this led to reduced growth during dry periods. It was striking that the growth reductions of large trees were more pronounced than those of smaller trees (Figure 6a). The significant negative interaction of SMI with sdi mean leads to a positive effect of higher sdi mean values on id in dry periods, whereas the opposite was true in periods with sufficient water supply (Figure 6b).

Kind of thinning (H IV)
The results obtained by fitting the regression model from Equation 10 are listed in Table 9. All fixed effect parameters were significant with p < 0.0001, indicating apparent effects of soil moisture index, SMI, and kind of thinning, kindth mean , on the periodic annual mean stem growth, id. We found positive main effects of stem diameter, d and SMI, a negative main effect of kindth mean , and a positive effect of the interaction between SMI and kindth mean . Hence, a trend towards thinning from below positively affected growth in dry periods, whereas thinning from above had a slightly positive impact on growth in periods with a good water supply (Figure 6c).

Temporal sequence of thinning (H V)
The relationship between the periodic annual mean tree diameter increment, id, the soil moisture index, SMI, and the temporal sequence of thinning according to model 12 is visualised in Figure 6d. All fixed effect parameters were significant with p < 0.0001, indicating apparent effects of SMI and the coefficient of variation of competition release by thinning, Dsdi cvar , on id (Table 10). We found positive main effects of stem diameter, d, and SMI, and a negative main effect of Dsdi cvar . This negative main effect of Dsdi cvar was somewhat mitigated by the positive effect of the interaction between SMI and Dsdi cvar , resulting in a low impact of Dsdi cvar in periods with sufficient water supply. However, in dry periods, higher Dsdi cvar values had a negative effect on id (Figure 6d).

The results in the context of other studies
Our study revealed significant differences in mean tree growth by soil moisture index. In addition, tree growth responses to drought stress were significantly influenced by the kind of thinning, the temporal sequence of thinning, and local stand density to a lesser degree. Several studies have already examined the effect of thinning on drought responses aiming to increase resilience or resistance (e.g. Kohler et al. 2010;D'Amato et al. 2013;Sohn et al. 2016b). In accordance with our findings, McDowell et al. (2006) showed that ponderosa pine (Pinus ponderosa Douglas ex C. Lawson) trees growing in open stands in Arizona, USA, were bigger and exhibited a more significant relative growth decline under drought than trees growing at a Table 5.
Model 4: Model 5: slower rate in high competition environments under high stand densities. They concluded that trees growing in highly competitive environments were less likely to benefit from wet or average growing conditions and less vulnerable under drought conditions. Furthermore, extreme reductions in stand density have been reported to counteract the generally positive effects of reduced stand density on drought resistance and resilience over time (D'Amato et al. 2013). This was attributed to significantly greater tree sizes attained within the lower-density stands through stand development, resulting in higher water demand than smaller trees. Steckel et al. (2020) showed that the resistance and resilience to drought events of ponderosa pine significantly increased with increasing stand density and that the growth sensitivity peaked under maximum stand density. Martínez-Vilalta et al. (2012) observed that fast-growing trees were more affected by drought for Scots pine (Pinus sylvestris L.), suggesting the benefits of denser stands with slower growth rates. Low stand densities may also negatively affect the stand climate. Increased wind speeds and intensified direct solar radiation may reduce soil humidity and higher evapotranspiration (Aussenac 2000;Lagergren et al. 2008;Primicia et al. 2013). A widely open canopy could further promote the herbaceous understorey (Hedwall et al. 2013) and thus, increase competition for water. Furthermore, low stand densities may also negate the positive shading effects for smaller trees (Pretzsch 2018) However, contradictory results have also been reported in the case of Scots pine (Pinus sylvestris L.). Particularly for trees growing at dry sites, low stand densities seemed to be beneficial (Giuggiola et al. 2013;Olivar et al. 2014;Manrique-Alba et al. 2020;Steckel et al. 2020) as it helped to reduce the climatic sensitivity of the remaining trees (Martín-Benito et al. 2010;Magruder et al. 2013). Sohn et al. (2016a) found that growth recovery following drought was highest after the first thinning intervention and in recently and heavily thinned stand for Scots pine stands in Germany. However, with time since the last thinning, this effect decreased and could even become negative compared to unthinned stands. In a meta-analyses, Sohn et al. (2016b) also found that benefits for post-drought growth decrease with time since the last thinning. This fits with our findings that increased variation in thinning intensity had a negative effect on both stand growth and individual tree growth, especially during dry periods.

Underlying mechanisms
We statistically revealed new relationships between the growth of trees and their silvicultural treatment in the past. Figure 5. Effect of (a) soil moisture index, SMI, (b) mean stand density of all previous survey periods, SDI mean , and (c) coefficient of variation of density reduction by thinning, DSDI cvar , on the current stand volume growth, IV. Lines were generated by fitting the linear mixed effect models from Equations (3)-(5). In case of the graphs (b) and (c), the respective variables SDI mean and DSDI cvar were varied, whereas the volume of the remaining stand, V, was kept at the mean value. We did not fathom the underlying mechanisms but can hypothesise how the silvicultural treatment and thereby generated stand and tree structure may determine tree growth and behaviour under drought stress at a later age.
In principle, any long-term effects of silvicultural treatment (e.g. tending, pruning, thinning) on tree growth (e.g. stem diameter growth) may operate via the environmental conditions modified by the treatment and/or modification of the tree structure and functioning. In terms of growing conditions, silvicultural measures may lower the stand density, improve the stand-internal temperature, accelerate the humus decomposition, nutrient turnover, and thus the growth of the remaining trees, especially on moist sites in temperate areas (Aussenac 2000). On dry sites, in Mediterranean regions, density reduction may even increase the drying of the soil, the raw humus accumulation, and deteriorate the nutrient turnover and growth (Calama et al. 2019). The long-term effect of silvicultural treatment on tree structure and functioning may operate by, among others, the crown length and width, the tree ring width, sapwood width, or cell embolism. All these traits may be immediately affected by silvicultural measures but may also determine the growth at a later age, i.e. cause legacy effects on growth. Our findings that in addition to tree size and state variables, metrics of the local growing conditions, the past silvicultural treatment, and tree growth can considerably improve the estimation of present tree growth (Supplementary Table 4) suggest that the legacy effect is based on both the local growing conditions and the internal tree structure and morphology, both shaped by silvicultural treatment. Table 8. Results of fitting linear mixed effect model of periodic annual mean tree diameter increment, id, as a function of tree diameter at the beginning of the survey period, d, mean local stand density index over all previous survey periods, sdi mean , and mean soil moisture index during the growth period, SMI. The underlying Model 10 was ln (id ik ) = a 0 + a 1 × ln (d ik ) + a 2 × sdi mean, ik + a 3 × SMI k + a 4 × SMI k × sdi mean, ik + b i + e ik .  Figure 6. Effect of (a) diameter at breast height, dbh, (b) mean local stand density index over all previous survey periods, sdi mean , (c) kind of thinning, kindth, and (d) the coefficient of variation of competition release by thinning, Dsdi cvar on the current mean diameter increment ∼ soil moisture index relationship. Lines were generated by fitting the linear mixed effect models from Equations (10)-(12). In each of the graphs, the respective variables dbh, sdi mean , kindth, and Dsdi cvar were varied, whereas the other variables were mean centred.
Effect of treatment and resulting stand structure: Among the causality mechanisms invoked to explain the adverse impact of reduced stand density on drought responses, increased evaporation and transpiration will likely prove detrimental following a reduction in stand density (Aussenac 2000;Lagergren et al. 2008;Brooks and Mitchell 2011). An explanation of the observed beneficial drought responses of N. spruce under higher stand densities may be found in facilitative interactions (e.g. wind shelter and a reduction in radiation) that may have overridden intraspecific competition (Fajardo and McIntire 2011;Calama et al. 2019;Owen 2019). In our case, similar mechanisms may have decreased N. spruce growth responsiveness to the climate in dense compared to more open locale densities and improved growth responses under episodic drought. This result is somewhat surprising, as intraspecific competition is generally expected to be fiercer than interspecific competition (Fajardo and McIntire 2011;del Río et al. 2014;Pretzsch 2020). Pretzsch (2021a) was able to show a strong effect of the trees' past on the diameter growth at present. In his study of an approximately 200-year-old European beech (Fagus sylvatica L.) stand, he found that trees with lower interannual variations of growth in the past had significantly higher growth rates at present than trees with higher interannual variations of growth in the past. We can confirm this result with our study, especially for growth during dry periods. However, we were also able to show that the variation in thinning intensity (e.g. Dsdi cvar ) and the kind of thinning served the same principal component as the variation in diameter growth (e.g. id rg ; cf. Supplement Figure 2). "External" and "internal" indicators are therefore linked.
Effects of tree structure and morphology: Past conditions may affect present growth via epigenetic (Rico et al. 2014;Bose et al. 2020), transcriptional, proteomic, and metabolic changes (Fleta-Soriano and Munné-Bosch 2016), and changes in the plant structure (Netzer et al. 2019). Trees in seasonal forests have a special property: their annual rings and crowns structurally store their past growth rhythm and can therefore be used to quantify their past using appropriate metrics (Pretzsch 2021a). Thus, tree-ring and crown patterns represent a structural memory embedded in the stem and crown (Backhaus et al. 2014;Ogle et al. 2015). They may affect the trees' functioning and growth, e.g. via light interception.
Strong adverse effects of the inter-annual fluctuations represented by id cvar and id rg , may be caused by a higher demand for photosynthates for repair and a depletion of the tree's reserve pool in the past, resulting in low resistance. High id cvar values may indicate repeated embolisms and disturbance of the internal conduction system, which increases susceptibility and predisposition to drought stress. Netzer et al. (2019) addressed the ambivalent effect of early high growth rates on growth and drought resistance. Trees with larger xylem vessels have a greater risk of embolism formation and lower hydraulic conductivity during drought. Thus, the larger the vessel size and water conductivity, the higher the embolism risk; however, the better the xylem transport system is protected against the formation of embolisms, the less efficient it is in the water pipeline (Pretzsch 2021b).

Methodological considerations
Although the experiment FFB 612 is suitable for the investigations we made in this study, the experimental design for analysing effects on the stand level could have been better Table 9. Results of fitting linear mixed effect model of periodic annual mean tree diameter increment, id, as a function of tree diameter at the beginning of the survey period, d, mean kind of thinning over all previous survey periods, kindth mean , and mean soil moisture index during the growth period, SMI. The underlying Model 11 was ln (id ik ) = a 0 + a 1 × ln (d ik ) + a 2 × kindth mean, ik + a 3 × SMI k + a 4 × SMI k × kindth mean, ik + b i + e ik . replicated. At the individual tree level, however, the study design is well replicated by our approach of describing the competitive situation of the trees via tree-individual circles. Despite the selected search radius of 0.25 x tree height to describe the competition situation of the trees showed the highest correlation to growth in our study as well as in other studies (Pretzsch and Schütze 2021;Pretzsch 2022aPretzsch , 2022b, it should be taken into account that the results could change with the selection of other search radii. Since stand density may affect stem forms (Jacobs et al. 2020), further studies should also examine tree volume growth of the individual trees in addition to tree diameter growth. However, to reconstruct the tree volume growth in detail, trees would have to be removed for such analyses to examine stem slices along the stem axis.
To test the hypothesis on memory effects at the tree level in more detail, it would have been helpful to record crown characteristics on the experimental plots repeatedly. Due to the allometric relationships between stem diameter and crown growth, the development of tree sizes such as tree height, branch diameter and branch length can follow a similar pattern as stem growth. Thus, if the accumulated trunk and crown diameters are equal, this can lead to a different internal trunk and crown structure due to different growth patterns. Such patterns may represent a structural memory anchored in the stem and crown that can influence tree function and growth (Pretzsch 2021a). Future studies should examine such patterns by analysing stem growth and addressing crown structure in detail, for example, by TLidar.
About our analyses on the contribution of the past to the estimation of current tree growth (H II), it should be taken into account that the last recording period we examined only included two years and that this recording period included the year 2018, which was a drought year in Europe.
Our findings only apply to N. spruce trees that grew on a relatively productive site up to a medium age. Even though  was also able to prove a strong effect of the trees' past on the diameter growth of approximately 200year-old European beech (Fagus sylvatica L.) trees at another site, we propose to include a site gradient and extend the design to additional tree species to verify our initial findings in further studies.

Implications for modelling
Current tree growth may be co-determined by the tree and stand development of the longer past. Past tree development may strongly correlate with growth via the trees' internal structure and morphology (Pretzsch 2021b). For example, the growth reactions during drought may also strongly depend on the trees' past. Even in stands of the same age and stocked by trees with similar tree diameters, the drought reaction may be much stronger in stands with strong crop tree thinning in the past than in stands with a moderate and slow opening in the youth. To avoid such a bias due to past development, the ecological memory effect can be considered through metrics such as those we have introduced in this study. Growth models could be improved by considering such metrics that reflect the memory effect or the treatment history at the tree level. Models at the stand level can exemplarily test the relevance of the model variables at the tree level. However, the treelevel results need to be put on a broader basis to derive general conclusions on the influence of tree-level model variables on stand dynamics by forest growth simulators.
Our study examined diameter increment (id) at tree level and volume increment (IV) at the stand level. We did not investigate biomass or carbon increments. To conclude biomass or carbon increments, allometric relationships between id or IV and biomass or wood density would have to be considered. However, it might be that these allometric relationships vary from treatment variant to treatment variant and further increase or decrease the differences between the variants.

Consequences for silvicultural treatment
Our results showed that the timing and strength of stand density regulation and kind of thinning could strongly codetermine tree and stand growth, especially under drought stress. The introduction addressed the change of silvicultural prescriptions in Central Europe towards stronger and more abrupt density reductions by heavy thinning from above. This means that the presently common silvicultural approaches may be detrimental to growth, especially under drought. Given climate change, the prevailing repeated strong thinnings from above for promoting a selected collective of future crop trees may be called into question. We found a positive effect of thinning from below, i.e. a positive effect on growth when future crop trees are released by removing small instead of large neighbours.
Furthermore, the growth during drought decreased with the variation of the strength of the competition release. However, tree growth increased with the variation of the strength of competition release during periods with sufficient water supply. The response pattern was similar at the tree and stand level; the analyses at the stand level confirmed that a low variation in thinning intensity positively affected stand growth during dry periods. However, trees may develop uniformly in continuously dense stands. In the event of humans opening up the stands or natural disturbances, all trees may respond similarly (Pretzsch 2021a). On the contrary, two or multi-layered forests contain a mixture of trees with different past development. They may have a higher capacity to respond to silvicultural inferences and disturbances, enabling a stable stand growth on a high level.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have influenced the work reported in this paper.