## Spatial phase synchronisation identified in the orchard

Figure 1. Spatial distributions of the in-phase synchronisation in the orchard.

To explain the spatial synchrony of the orchard ,the indices for phase synchronisations were defined as described in the Method-Statistics section. Figure 1(a) displays the spatially distributed fraction of the in-phase (\({f}_{in}^{i}\left(t\right)\)) for 9,562 individual trees from 2002–2007. To conduct this analysis the orchard was divided into fourteen blocks by rows and columns: (row, column) = (\(2\times 7\)). \({f}_{in}^{i}\) were calculated for each tree *i* within the block to which tree *i* belonged to.\({f}_{in}^{i}\) was high enough to reach 1.0 in the west blocks and was as low as 0.7 in the east blocks. Figure 1(c) shows \({F}_{IN}^{K}\) (*K* = 1, 2, …, 30), which is the average \({f}_{in}^{i}\) for every five columns. \({F}_{IN}^{L}\) (*L* = 1, 2, …, 8) is the south to north directional average of \({f}_{in}^{i}\) for every eight rows. \({F}_{IN}^{K}\)drops from 0.92 to 0.7 from west to east (Fig. 2(b)). Conversely, only a small gradient of \({F}_{IN}^{L}\) was observed along the south-north direction (Fig. 2(c)). These results indicate that the strength of the phase synchronisation significantly decreases from west to east but is trivial in the south-north direction. In other words, a significant phase transition is observed from order to disorder in the west-east direction since perfect synchronisation (order) and desynchronization (disorder) of a population corresponds to \({F}_{IN}^{}\) to 1.0 and 0.5 respectively. We applied *I*(*d*), which is known as Moran’s I, as the index for the spatial correlation29−31 as shown in Fig. 1(d). The six years averaged *I*(d) of the yield data exhibit a 1/3 power-law scaling as plotted in red.

Figure 2 shows the time evolution of the spatial distribution of phase synchronization in the five periods of two successive years: [2002-2003], [2003-2004],[2004-2005],[2005-2006] and [2006-2007] in five rows, respectively. The first and second columns represent the spatial distributions of \({f}_{in}^{i}\) and \({F}_{IN}^{K}\left(t\right)\) vs. *K* plot. In the first period [2002–2003], the strength of phase synchronization and spatial distribution of phase synchronization were weak compared to other periods, which is shown in Fig. 2(a) and Fig. 2(f).

In the second to the third periods;[2003-2004] and [2004-2005], the strength of phase synchronization increased, and the spatial distribution of phase synchronization from west to east became distinct, that is, it was stronger in the west and weaker in the east. In fact, \({F}_{IN}^{1}\left(2004\right)\) and \({F}_{IN}^{30}\left(2004\right)\) were 0.98 and 0.56, indicating perfect phase synchronization and desynchronization, respectively. For the fourth [2005–2006] and fifth [2006–2007] periods, the perfect phase synchronisations dominate; thus, \({f}_{in}^{i}\left(2006\right)\) and \({f}_{in}^{i}\left(2007\right)\) reach 1.0 over the entire orchard, and \({F}_{IN}^{K}\left(2005\right) and {F}_{IN}^{K}\left(2006\right)\) values are almost always 1.0 for the whole orchard ( i.e., *K*=1,2,..,30). It is very important that the perfect synchronization mode lasted for two periods [2005-2006] and [2006-2007]. This is clear experimental evidence that mode-locking can occur in a real orchard.

Figures 1 and 2 show that the phase transition between “order” and “disorder” emerge both spatially and temporally (annually) in the orchard. As described above, the index of the fractions of the in-phase such as \({f}_{in}^{i}\left(t\right),{f}_{in}^{i}, {F}_{IN}^{K}\left(t\right) \text{a}\text{n}\text{d} {F}_{IN}^{K}\) successfully describe detailed information on how the features of phase synchronisations behave both spatially and temporally.

## The states of “on-year” and “off-year”

In this section, we identify the states of the “on-year” and “off-year” of a population by introducing the phase synchronisation analysis technique (see METHODS section). The phase angle \({\theta }_{i}\left(t\right)\) of each tree for each year calculated from eq.(7) determined the “on” and “off” states for the individual trees. The states of “on-year” and “off-year” are defined by the period of time covered as follows; an on-year is one in which the phase angle \(\theta\) is positive, or the yield \(x\left(\text{t}\right)\) is greater than the average yield of the period \(\stackrel{-}{x}\left(\text{t}\right)\) and vice versa16.

$$\begin{array}{c}\text{o}\text{n}-\text{y}\text{e}\text{a}\text{r}: \text{c}\text{o}\text{s} \theta >0 or x\left(\text{t}\right)>\stackrel{-}{x}\left(\text{t}\right)\\ \text{o}\text{f}\text{f}-\text{y}\text{e}\text{a}\text{r}: \text{c}\text{o}\text{s} \theta \le 0 or x\left(\text{t}\right)\le \stackrel{-}{x}\left(\text{t}\right)\end{array} .$$

1

The vector point for each subpopulation block is given by the mean of all coordinates of the dots i.e. (\(\frac{1}{n}\sum _{i=1}^{n}\text{c}\text{o}\text{s}{\theta }_{i}\left(t\right), \frac{1}{n}\sum _{i=1}^{n}\text{s}\text{i}\text{n}{\theta }_{i}\left(t\right))\), where *n* is the subpopulation size. The amplitude of the vector is the order parameter \(p\left(t\right)\) of the subpopulation.

$$p\left(t\right)=\sqrt{({\frac{1}{n}\sum _{i=1}^{n}\text{c}\text{o}\text{s}{\theta }_{i}\left(t\right))}^{2}+({\frac{1}{n}\sum _{i=1}^{n}\text{s}\text{i}\text{n}{\theta }_{i}\left(t\right))}^{2}}.$$

2

Values *p* = 1 and 0 represent perfect phase synchronisation and desynchronisation, respectively. The “on-year” and “off-year” states of the population are distinguished by the red and black arrows. For example, in 2002 and 2007, all seven blocks are “on-year” states and the states are identical between the seven sub-populations and the whole population. However, in 2003, the state of each subpopulation shifted from a weak “on-year” state to “off-year” states in west to east. The coexistence of the “on-year” and “off-year” states within pistachio orchards has only been known as general qualitative knowledge. With this novel technique, we can quantitatively visualise the coexistence of the “on-year” and “off-year” states of individual trees in both spatial and temporal (annual) domains.

Figure 4 shows the composition of periods in individual trees and the population. Several researchers have reported that alternate bearing of pistachios do not always exhibit a clear two-year cycle but may feature various types of synchronisations both spatially and annually20,36−38. Sakai et al. (2019) proposed a practical approach to determine the composition of periodic components of alternate bearing and masting of plant populations16. For instance, period-two and period-three sequences are defined as ‘ON ⇒ OFF ⇒ ON’ and ‘ON ⇒ OFF ⇒ OFF ⇒ ON’, respectively (see METHODS-*Fraction of period* section). We mathematically determined the states of On and Off in each year for individual trees with \({ON}_{i}\left(t\right)\) and the spatial distributions of \({ON}_{i}\left(t\right)\) for 6 years are displayed in Fig. 4-(a)-(f). Wild tree spaces, generally have a wide-ranged spectrum of periods in which variable periods such as period-two, period-three, period-four, and more, including aperiodic motions, exist. For example, period-three was dominant for *Zerkova serrata* but period-four and period-five were non-trivial16. For tree crops, their synchronisations occurred generally on a two-year cycle (i.e. period-two), which is why “alternate” bearing and/or “biennial” bearing has become common terminology in pomology1. Using the fractions of periods *FR*(*Q*) (see METHODS-*Fraction of period* section), we successfully display the spatial distributions of the compositions of period-two, period-three, period-four and period-five for individual trees of the orchard in Fig. 4 (g)-(j). Figure 4 (g)-(j) indicate the fractions of period-two, period-three, period-four and period-five, which are *FR*(2) = 0.517, *FR*(3) = 0.103, *FR*(4)=0.0244 and *FR*(5)=0.0060, respectively. Obviously,”period-two” is dominant in the orchard. However, from a farm management point of view, it is informative to identify those trees whose cycles are three years, four years or five-years. For example, in the eastern area, the fraction of period-two *FR*(2) is smaller than that in the west; further, the fraction of period-three *FR*(3) in the east is larger than that in the west. For another example, from the 9,562 trees, we can identify 57 trees whose period are five-years and investigate why they behave in this unique manner.

The results from the novel approach applied to the yield data of the orchard uncover unique and important features of the phase synchronisations in the orchards:

(1) The gradient of the strength of in-phase synchronisation from west to east.

(2) The phase transition from order (perfect in-phase synchronisation) to disorder (desynchronisation).

(3) The occurrence of perfect synchronisation in the whole orchard (mode locking).

To infer from these findings and explain the observed phase synchronisation, we developed a model. The structure of the model was a network of oscillators, which were the RBMs for all the trees; a common noise (external force) \({e}_{C}\)(t) was identically imposed on all the trees in each year *t*. We employed diffusive coupling as the direct coupling of the network and considered the coupling timing in the growth stage of a plant. The detail of the model is given in the METHODS- *model development* *section*)

## Essential factors in the spatial synchronisation of the orchard

Figure 5 shows that the common noise, local coupling, and gradient of the cropping coefficient are essential factors. Given the combinations of the essential parameters (ε, eC,α,β), the best fit common noise (external force) *e*C(t) and initial values of *S*i(1) (i = 1, 2, …, 9,562) for 25-year periods are determined as follows. The randam initial values of *S*i(1) for 9,562 trees are given by homogenious random number in (LT-P0, LT). The model runs until *t* becomes 500 and the best fit 25 years period in terms of \({F}_{\text{I}\mathbb{N}}^{K}\) (*K*=1,2,…,30) is subsequently selected. The spatial correlation *I*(d) diagrams and the spatial distribution of \({f}_{\text{i}\text{n}}^{\text{i}}\) maps are exhibited in the left column and the right column, respectively.

According to the above three findings, the results, and previous work, we assume that the common noise (external environmental force) (*e*C), gradient of cropping coefficient (*m* ), and direct coupling (\(\epsilon\)) are the three essential factors in the model. The range of coupling is set r=13 m so that one tree coupled with neighbouring twelve trees at most. Three characteristics of the spatial correlation are observed in the orchard: (i) high short-range spatial correlation, (ii) long-range spatial correlation with 1/3 power-law scaling, and (iii) wide range variation of *I*(*d*) on the time (year) domain.

In this section, we validate the effects of the three parameters on phase synchronisation by considering the characteristics of the spatial correlation.

Numerical experiments with the developed model were conducted to investigate the effect of the essential factors in terms of *I*(d). Here, we use *I*(*d*), known as Moran’s I, as the index of the spatial correlation39, the map of \({f}_{\text{i}\text{n}}^{\text{i}}\) and \({F}_{IN}^{K}\) to compare with the real data shown in Fig. 1-(d),(b) and (a), respectivelly. Four different values of the set \({(\epsilon ,e}_{C},\alpha ,\beta )\) were tested (Fig. 5). *m* and \(\epsilon\) are the control parameters of the network model and \({e}_{C}\) is the common noise (environmental external force) that generates the phase synchronisation. *α* and *β* are the spatial gradients of *m* in west-east and south-north directions, respectively. Based on a preliminary parameter study (see S-1), we estimated the values for the essential factors as *e*C = 0.2, *ε *= 0.1 and *m* = *R*0 + *α*x/*L*WE + *β*y/*L*WE (*R*0 =1.1,* α *= 0.6 and *β *= −0.1) i.e. \({(\epsilon ,e}_{C},\alpha ,\beta )=(0.1, 0.2, 0.6, -0.1)\) for Fig. 5(a)-(c). In this case, the spatial correlations (*I*(*d*)) of the model satisfies the three characteristics of the orchard as follows; (i) for short-range spatial correlation, *I*(*d*) is up to 0.64 at *d* = 5.2 m, (ii) for long-range spatial correlation, *I*(*d*) is larger than 0.37 even at *d* = 100 m with 1/3 power row scaling, and (iii) *I*(*d*) fluctuates widely over the same range as actual data (Fig. 1-(d)). The map of \({f}_{\text{i}\text{n}}^{\text{i}}\) (Fig. 5(b)) shows good agreement with that of the real data shown in Fig. 1-(a), and the actual and model plots of \({F}_{IN}^{K}\) plots (Fig. 5(c)) also shows good agreement. In the absence of coupling (second row), the short range spatial correlation (Fig. 5(d)) is much smaller than that of the presence of coupling case (Fig. 5(a)), while the map of \({f}_{\text{i}\text{n}}^{\text{i}}\) (Fig. 5(e)) and the actual and model of \({F}_{IN}^{K}\) plots (Fig. 5(f) show some agreement with the real data in the western area, but not the eastern area. Without common noise case (third row), in Fig. 5(g), the short-range correlation is as high as 0.3 because of the local direct coupling (*ε* = 0.1). However, all the *I*(*d*)s rapidly drop to zero when *d* = 20 m, and there is no long-range spatial correlation, no power row scaling, and no yearly fluctuation in *I*(*d*). The map of \({f}_{\text{i}\text{n}}^{\text{i}}\)and the actual and model of \({F}_{IN}^{K}\) plots indicate that there is no significant phase synchronization over the entire orchard (Fig. 5(h) and (i)). These results suggest that common noise is indispensable as an external force and the synchronization observed in the orchard (Fig. 5) should be a type of common noise induce synchronisation.

As the cropping efficiency gradient is given by *m* = *R*0 + α*x/L*WE + *βy/L*SN, the absence of *m*-gradient is implemented by setting α = *β *= 0 ,therefore, *m* is 1.1 over the orchard. In this case of without the spatial gradient of the cropping efficiency (Fig. 5(j)), very strong short-range spatial correlation, *I*(*d*) = 0.98 exists at *d* = 5.2 m; however, there is no long-range spatial correlation. Figure 5-(k) and (l) exhibits the spatial phase synchronization with \({f}_{\text{i}\text{n}}^{\text{i}}\) =1 over the entire orchard and for 25 years. The results of Figs. 5(j)-(l) suggest that even with the presence of both the common noise (\({e}_{C}\)) and coupling (\(\epsilon\)), the gradient of *m* is still necessary to explain the features of long-range spatial correlations and spatial distribution of phase synchronisation in the orchard as observed in Fig. 1.

## Numerical interpretation of the experimental results: How endogenous dynamics and exogenous force generate variable types of phase synchronisations

Using the model with the estimated parameters ; \({(\epsilon ,e}_{C},\alpha ,\beta )=(0.1, 0.2, 0.6, -0.1)\) demonstrated in Fig. 5-(a) and (b), Fig. 6 exhibits how endogenous dynamics and exogenous force generate variable types of phase synchronisations. Figure 6-(a) shows the spatial distribution of \({F}_{IN}^{K}\left(\text{t}\right)\) of the model. The phase transitions between order and disorder in both space and time (annual) domains were the substantial nature of the dynamics and important behaviour in terms of phase synchronisation.

Resource switching and weather cues act together to cause synchronisation such as masting and alternate bearing (Lyles et al. 2009)7. The results described in Fig. 6 confirm this knowledge mathematically. The combination of endogenous dynamics and exogenous forces cause the perfect phase synchronisations. Figure 6(a) displays the \({F}_{IN}^{K}\left(\text{t}\right)\) for 25 years. An interesting finding of this analysis is that significantly strong phase synchronisations generally occur when the common noise (external force) is high; however, this is not always the case. Furthermore, almost complete desynchronisation occurs even when a significantly large common noise (external force) is present. The perfect phase synchronisations occur when *t* = 12, 13, and 14 with \({F}_{IN}^{}\left(\text{t}\right)\ge 0.997\) (i.e. \({F}_{IN}^{}\left(12\right)=0.9975,{F}_{IN}^{}\left(13\right)=1.0,\text{a}\text{n}\text{d} {F}_{IN}^{}\left(14\right)=1.0\)). The almost complete desynchronisations occur when *t* = 4 and 7 with \({F}_{IN}^{}\left(\text{t}\right)\le 0.58\).

Three consecutive perfect in-phase synchronisations occur when *t* = 12, 13, and 14; *t* = 13 is an “off-year” as the year (*t*=12) is the second largest “on-year” within 25 years; thus, the production (*C*f (12)) was zero, because the excess of the threshold (*L*t=100) of the total amount of *S*(t), *P*S(t), and *CE*(t) becomes the flowering cost *C*f(t) which is proportional to the production *C*a(t). Consequently, the resource remained in the bodies of the tree, meaning that *S* (13) nearly obtained the almost maximum value for all trees. Additionally, that year (*t*=13) had the best condition as eC(13)=0.4 (*CE*(13)=4) within 25 years, which allowed all the trees to change their states into the maximum on-year state with the perfect in-phase synchronisation (in fact \({F}_{IN}^{}\left(13\right)=1.0)\). As a result, the smallest amount of resources remained in the bodies of the trees meaning that *S* (14) obtained the minimum value within 25 years (see Fig. 6(c) at *t*=14). The external force was not large (*CE*(14)=2) with no excess (i.e. *C*f(14)=0) for all trees so that the perfect synchronisation (\({F}_{IN}^{}\left(13\right)=1.0)\)still remained in “off-year” state. Thus, three consecutive perfect phase synchronisations in other words, this is the mechanism of the three years consecutive perfect phase synchronisations (i.e. three years mode locking)). If *CE*(14) is at least greater than 6, *t*=14 will be a weak "on-year" state with partial desynchronization and this mode lock will terminate.

In the opposite case, *t*=6 was an “on-year” with a weak synchronisation mode (see Fig. 6(a)), *t*=7 should be succeeded by “off-year” if *CE*(7) was a normal value. However, *CE*(7) was large enough (i.e. *CE*(7) =3.5 is the second largest within 25 years,) to keep very weak “on-year” and the model turned to the almost complete desynchronization (\({F}_{IN}^{}\left(7\right)=0.568\)). It is interesting to note that the simulation results show that such a large external force causes desynchronization instead of synchronization.

Owing to the historically imposed exogenous forces (accumulated footprints of exogenous forces) on the endogenous dynamics, a mode-locking emerged. We note that two types of “in-phase mode-locking” exist, one which occurred in the “off-year” and another which occurred in the “on-year”. The first one occurred at *t* = 12 and 14, and the second one occurred at *t* = 13. Also two types of “desynchronization mode” exist, one which occurred in the “off-year” and another which occurred in the “on-year”, occurred at t=4 and t=7, respectively. The perfect in-phase synchronisation (order) and perfect desynchronisation (disorder) were the two extremes. Most of the modes of phase synchronisation appeared either in between or were mixtures of both extremes e.g. Fig. 6.(a) at t=20 and 23. The mode of the mixture may be considered as a “chimaera like state”, which was first experimentally captured in optical-coupled map lattices in nonlinear physics40. Here, we detected the “chimera like state” in the field experiments of the orchard in 2004-2005 period as shown in Fig. 2(h).