Chaos emerging in soil failure patterns observed during tillage: Normalized deterministic nonlinear prediction (NDNP) and its application

Soil consists of a very large number of interconnected granular particles with internal friction at the particle interfaces. It has been shown that the soil cutting forces of a tillage tool can be modelled by an autoregressive model, which has parameters related to soil properties such as soil moisture and bulk density.¹ Numerous studies have indicated that simple tools interacting with agricultural soils result in horizontal forces that fluctuate periodically.^2–4 The fast Fourier transform of the tillage data has shown the presence of a clear peak in the frequency domain.^5–7 These phenomena are often observed by farm operators engaged in tillage. Some suggest that the phenomenon is a periodic motion caused by the shape-plane dynamics. However, Glancey et al.⁵ reported that there are significant variations in the periods and amplitudes of individual cycles, even in well-controlled soil bin conditions. In other words, the motion should be recognized as a pseudo-periodic motion consisting of inter and intra dynamics. A clear understanding of this phenomenon is important for soil tillage research and for all areas related to soil–tool interactions.

Stick-slip models are popular methods for representing nonlinear vibrations in mechanical engineering, structural engineering, and geophysics.^8–11 Walker et al.¹² applied complex systems analysis methods to investigate stick-slip dynamics observed in a laboratory fault as a proxy for a geological fault and suggested that stick events are governed by the formation of force chains, whereas slips are due to the collective failure of force chains caused by buckling and elastic unloading. We should note that their suggestion is consistent with observations regarding soil failure in farm tillage operations.

Our objective in this study was to clarify the nonlinear characteristics of soil failure patterns observed in soil tillage operations. We used two steps to reach this goal. First, we validated the existence of a deterministic process. Second, we identified the nonlinear dynamics characterized by an orbital instability. We applied the translation error (Wayland test) and Fourier transform phase randomization surrogate (FPR) to check the determinism in the first step. To reduce the effect of autocorrelations caused by contaminated noise, we normalized the deterministic nonlinear prediction (DNP) using an autocorrelation function (ACF). We applied the proposed method to distinguish the motions from the deterministic and linear stochastic Gaussian processes in four different soil conditions and two soil cutting speeds. We demonstrated that normalized deterministic nonlinear prediction (NDNP) can distinguish the effects of additive noise from nonlinear dynamics.

Field tests were conducted in a Yolo loam field located near the Western Center for Agricultural Equipment, University of California at Davis. The experiment followed a split-plot design where the main treatments were tilled and undisturbed conditions, and the sub-treatments were irrigated/wet and dry soil conditions. The tilled treatment consisted of ripping the soil to a depth of 60 cm and subsequently disking the surface to a depth of 15 cm. In the undisturbed treatment, the land was left fallow during the previous fall and winter seasons. The irrigation treatment consisted of 30 h of sprinkler irrigation until the soil was saturated with water to a depth of 40 cm. In each subplot, an instrumented tine (texture soil compaction index [TCI] sensor, developed at UC Davis by Liu et al.¹³ and Andrade-Sánchez et al.¹⁴) was used to obtain soil cutting force data at two different speeds, which we subsequently denote as FAST (1.1 m/s) and SLOW (0.22 m/s). The experiment was replicated in four blocks. This lead to four distinct soil conditions: tilled dry (TD), tilled wet (TW), untilled dry (UD), and untilled wet (UW). At the time of the field tests, cone index measurements¹⁵ were obtained using a standard tractor-mounted, hydraulically driven, self-recording cone penetrometer that was capable of obtaining cone index data to a depth of 53 cm. Moreover, soil moisture contents and density profiles from 5 to 60 cm were obtained in depth increments of 5 cm using a soil sampling device.

The soil cutting forces under the four soil conditions are shown in the time (s) and the spatial (m) domains in Figs. 1 and 2 for FAST (1.1 m/s) and SLOW (0.22 m/s), respectively. For FAST (Fig. 1), a clear cyclic motion is observed in TD (Fig. 1(a)). The discrete Fourier transform (DFT) is shown in a linear scale (Fig. 1(e)) and logarithmic scale (Fig. 1(f)). The broadband components indicating 1/f noise were observed for the four soil conditions. However, the predominant peak and their harmonics are clear for TD, and a weak peak is present for TW. The spatial frequencies (f_S) of TD are located in the range of 0.7–0.9 cycle/m and clear cracks appear repeatedly at 1.1–1.5 m spaces. Figure 2 shows the same results for SLOW. The predominant spatial frequencies of TD-SLOW (0.22 m/s) are located at the same range as TD-FAST (1.1 m/s). However, for UD and UW, there are no predominant peaks but the broadband components indicating 1/f noise appear in the DFT spectrum for FAST and SLOW. This implies that soil failure patterns in UD and UW originate from a linear stochastic Gaussian process. In these measurements, the sampling interval, Δt, was set to 0.01 (s) for FAST and 0.05 (s) for SLOW. These sampling intervals in the time domain correspond to 0.011 (m) in the spatial domain for both cases. As described above, we used eight experimental data sets of soil cutting forces, in two speeds and four soil conditions.

It is useful to compare the eight experimental data sets to time series produced using known processes. To be used as references, the time series should be derived from clear and known processes. We created four sets of time series data for reference. Lorenz63 denotes the Lorenz model with parameters σ = 10, ρ = 28, and β = 8/3,¹⁶ where

We solved Equation (1) using the Runge-Kutta method (RK4) with step size Δt = 0.01. Lorenz63 is the reference for a pure dynamical system.

Lorenz63A and Lorenz63B were derived from the following equation:

Additive noise (system noise) was incorporated using the term αδ, where δ is the Gaussian noise with mean 0 and standard deviation 1, and α is a scalar multiplier. Lorenz 63A and Lorenz 63B correspond to α = 0.3 and α = 0.6, respectively. The noise when α = 0.3 and 0.6 corresponds to 0.8% and 1.6% of the attractor size, defined as |max(x) − min(x)| = 36.6. Lorenz63A and Lorenz63B are the references for noise contaminated dynamical systems.

A data set referred to as “Colored Noise(Brown noise)” (that corresponds to a linear Gaussian process) was generated by the autoregressive model AR(1), as the distributions of x(t + 1) − x(t) for the soil cutting force data show clear normal distributions¹ 

where ε_t is Gaussian noise with mean 0 and standard deviation 1. Colored Noise is the reference for autocorrelative noise. The unit of t is the sampling interval.

Nonlinear time series analysis (NTSA) was developed to distinguish deterministic chaos from stochastic motion.^17–19 The main features of deterministic chaos such as the underlying dynamics and orbital instability can be characterized using nonlinear deterministic prediction and Lyapunov exponents.^20–22 

However, these conventional tools sometimes misidentify colored noise and 1/f noise as deterministic chaos. Theiler^23,24 pointed out the difficulty in positively identifying chaos in experimentally obtained time series data, because of the mixture of chaos, non-chaotic but nonlinear determinism, linear correlations, and noise. Therefore, it is necessary to use different NTSA tools to capture the different features of chaos. The false nearest neighbor (FNN) method determines the embedding dimension. The translation error quantifies phase space continuity to measure determinism in time series data. The Fourier transfer phase randomization surrogate (FPR) identifies determinism using the null-hypothesis of a linear Gaussian process. Deterministic nonlinear prediction (DNP) is a standard technique for identifying the underlying dynamics and orbital instability of a data set (i.e., sensitive dependency on the initial condition) by showing the coincidence of the short-term predictability and long-term un-predictability.

We assume that the dynamics can be written in the form

The possible dynamics are reconstructed from the observed soil cutting forces, x_t, using a time-delayed embedding technique. The state of the reconstructed phase space²⁵ is

where τ is the time lag and m is the embedding dimension. In our analysis, the lag τ = 16 is consistent in all cases.

The time lag τ = 16 corresponds to a quarter of the predominant period observed in TD-FAST and TD-SLOW shown in Figs. 1(e) and 2(e). The τ = 16 is also similar to the value obtained from the Lorenz model by the mutual information.²⁶ 

FNN determines the embedding dimension.²⁷ The state vector X(t) defined in Eq. (5) is rewritten as X(t,m) to include information regarding the embedding dimension m, i.e.,

We set

where n(t,m) is the index of the nearest neighbor vector X(n(t,m),m) of X(t,m).

If the embedding dimension is smaller than the minimum embedding dimension, a(t,m) becomes too large because of false neighbors effects. This approach can be realized by counting the cases where a(t,m) becomes larger than a given threshold TH, i.e.,

where H is the heaviside step function.

However, it is obvious that FNN(TH) is substantially affected by TH. To avoid the effect of TH, Cao defined a new measure²⁸ 

where $E(m)=1N−mτ∑t=1N−mτa(t,m)$⁠. E1(m) stops changing when m becomes greater than the appropriate embedding dimension. To check the determinism, he also proposed a quantity

where $E*(m)=1N−mτ∑t=1N−mτ|x(t+mτ)−x(n(t,m)+mτ)|$⁠. For random data that are independent of the past, E2(m) should be equal to 1 for any m. Consequently, there must exist some m at which E2(m) is not equal to 1.

We applied Cao's method to the eight experimental datasets and two reference datasets (i.e., Lorenz63 and colored noise). At first, we checked the determinism using E2 in Fig. 3. Lorenz63 is a completely deterministic, so we expect E2(1) of Lorenz 63 to be 0.1 (≠1 at m = 1). The E2(1) values for TD-FAST, TD-SLOW, and TW-SLOW were 0.74, 0.75, and 0.75 (≠1 at m = 1), respectively. We can interpret these values as satisfying the criteria of determinism.²⁸ This information allows us to use FNN to determine the dimension m using TD-FAST, TD-SLOW, and TW-SLOW. When the dimension m is larger than 5, the percentage of FNN is smaller than 0.5% (Fig. 4(a)). All the E1 asymptotically approach to 1.0, with m reaching 7 (Fig. 3(b)). Based on the results of Fig. 3, we decided that the embedding dimension, m = 6.

Note that the estimated dimension m = 6 satisfies the necessary condition²² D_C < 2log10(N), that is, the embedding dimension D_C is smaller than 7.40 when N = 5000.

Wayland²⁹ proposed a method for measuring determinism in time series using the “phase space continuity” observed in the time series. Let x(t) be an arbitrary point of a reconstructed attractor and y(t, j) denote the jth nearest neighbors of x(t). y(t + p, j) is the projection after p evolution steps of y(t, j). Here, y (t, 0) is defined as x (t) itself.

If the data are deterministic, we can assume that there will be a parallel vector field and the M translation vectors can be denoted by v(t + p,j) = y(t + p,j) − y(t,j).

The translation error is defined as

where $‖.‖$ is the Euclidian distance and

We set E_trans (t,p) to the median of $etrans(t,p)$⁠, t = 1,2, …, N.

Figure 4 shows E_trans(p) in a logarithmic scale (the horizontal axis is scaled by log₂(p)). FAST (1.1 m/s) and SLOW (0.22 m/s) are expressed by solid and dotted lines, respectively. There is no significant difference between FAST and SLOW for the four soil conditions (i.e., TD, TW, UD, and UW). With an increase of p, E_trans(p) decreases gradually except in the case of Lorenz63. E_trans(p) of Lorenz63 is independent of p and remains a constant with small values in the range of (0.00725–0.00633). These results imply that the slopes of all of the lines except Lorenz63 depend on additive and/or imposed noise. The values of E_trans(p) for TD-FAST and TD-SLOW are in the range between Lorenz63A and Lorenz63B. This means that phase space continuity of the reconstructed phase space of TD-FAST and TD-SLOW is similar to that of Lorenz63A and Lorenz63B. These results suggest that TD-FAST and TD-SLOW are additive-noise dynamical systems. However, it is difficult to distinguish colored noise from UD and UW for both FAST and SLOW. This implies that UD and UW mostly arise from linear Gaussian processes. The trends of TW are between TD and the two untilled conditions (i.e., UD and UW) for both FAST and SLOW.

The surrogate algorithm applied to the results of Wayland test can be utilized for determinism check using a null hypothesis test. Several different surrogate algorithms have been developed for an identical null hypothesis, e.g., random shuffle surrogate (RSS), Fourier transform phase randomization surrogate (FPR), attractor trajectory surrogate (ATS), and pseudo periodic surrogate (PPS).³⁰ According to the Fourier spectrum in Section II B, linear stochastic Gaussian processes must be a part of all the data because there is a clear broad 1/f component in Figs. 1 and 2. In this paper, we applied the FPR surrogate with a null hypothesis of a linear Gaussian process.

Figure 5 shows E_trans(p) of 100 surrogates as histograms and the original E_trans(p) as red lines for each data set (eight experimental data). The evolution step was set to p = 16, which corresponds to the lag of the embedding τ = 16. For TD-FAST and TD-SLOW, the original E_trans(16) are obviously out of the range of the E_trans(16) corresponding to 100 surrogates, as shown in Figs. 5(a) and 5(e).

Thus, we reject the null hypothesis of a linear Gaussian process for TD-FAST and TD-SLOW with 1% significance level and TD-SLOW with 5% significance level. This result is strong evidence of determinism in TD-FAST and TD-SLOW. The null hypothesis was not rejected in TW-FAST, UD-FAST, UD-SLOW, UW-FAST, and UW-SLOW (Figs. 5(b)–5(d), 5(g), and 5(h)). In arriving at these conclusions, we conducted the conventional hypothesis testing for large scale and the Monte Carlo hypothesis testing for determining the significance levels.³¹ 

Deterministic nonlinear prediction (DNP) is a popular and useful tool for identifying deterministic chaos.^32–38 We applied the prediction method initially proposed by Lorenz. Let X(t) be the predictee to find the M of nearest neighbors of X(t) denoted by y(t,1), y(t,2), …, y(t,M) in the reconstructed phase space. Here, the ten nearest neighbors of X(t) are selected from {X(l),l = 1,2, …, N; t − 4τ < l < t + 4τ}. The direction of the ith nearest neighbor in p-step forward is expressed by v(t + p,i): = y(t + p,i) − y(t,i). The p-step forward prediction $X⌢(t+p)$ is

where $ΔYAV(t+p)=1M∑i=1Mv(t+p,i)$⁠.

The prediction performance is usually evaluated using the correlation coefficient (CC). We define DNP as

where

If there is deterministic chaos, the underlying dynamical system could be used to predict the short-term state, but long-term prediction is impossible because of the orbital instability. The underlying dynamical system and orbital instability are two major features and necessary conditions of deterministic chaos.

To see the weakness of DNP with respect to colored noise, we plot the DNP of Lorenz63, Lorenz63A, and Colored Noise in Fig. 6(a). It is obvious that Lorenz63 is predictable in the short-term and un-predictable in the long-term. Lorenz63A and Lorenz63B are similar. The DNP(p) of Colored Noise superficially shows short term predictability with a very high CC (almost 1.0) in the short term. However, it is obvious that the colored noise arises from a stochastic process not a deterministic process. This contradiction is the typical example of the weakness of DNP when it is applied to noise-contaminated data. For instance, Fig. 6(b) shows that the autocorrelation (ACF) of Colored Noise is much higher than that of Lorenz63 and the other sets. Therefore, it is reasonable to assume that the short term trend of DNP(p) that is observed in Colored Noise expresses its high autocorrelation.

To resolve this shortcoming, we propose normalizing DNP by considering the autocorrelation. Here, we define the normalized deterministic nonlinear prediction $NDNPCC(p)$ as

where

To validate the performance of the proposed NDNP, we applied it to the four reference data sets in Fig. 6(c). The NDNP_CC(p) values of Lorenz63, Lorenz63A, and Lorenz63B are always larger than 1.0 and increase monotonically. It only decreases monotonically for Colored Noise and is always smaller than 1.0. This result demonstrates the ability of NDNP_CC(p) to distinguish nonlinear dynamics (Lorenz63, Lorenz63A, and Lorenz63B) from colored noise (i.e., a linear Gaussian process). NDNP_CC(p) > 1 is equivalent to DNP(p) > ACF(p), so this condition means that the component of time series fluctuations attributable to the deterministic dynamical system is larger than the effect of the autocorrelation.

Based on the discussion presented above, we applied DNP, ACF, and NDNP to soil cutting force data from the FAST and SLOW experiments. The results are also shown in Fig. 6. All the DNP values had high CC in the short-term and low CC in the long term. In particular, the DNP values for UW have the highest CC (Figs. 6(d) and 6(g)). The ACF values of UW and UD were also higher than the others, and some cyclic motions can be observed in TD and TW (Figs. 6(e) and 6(h)).

It is very impressive that NDNP(p) is larger than 1.0 for TD-FAST and TD-SLOW. As shown in Figs. 6(a)–6(c), NDNP(p) > 1.0 in the short term is a criterion for underlying dynamics, so we believe that TD-FAST and TD-SLOW have significant underlying dynamics. However, similar to Colored Noise (Fig. 6(c)), NDNP(p) is smaller than 1.0 in five cases (UD-FAST, UD-SLOW, UW-FAST, UD-SLOW, and TW-FAST). For these five cases, the effect of the autocorrelations is stronger than the underlying dynamics, even if they exist. It is of interest that NDNP(p) of TW-SLOW is very close to 1.0 in the short term. We discuss this point in Section IV, in combination with the Fourier transform surrogate results in Fig. 5(f).

To reinforce the rationality of NDNP as a new tool for NTSA, we examined the behavior of state vectors in the reconstructed phase space. The 10 nearest neighbors, y(t,j), for an arbitrary selected predictee X(t) are defined in Eq. (6). Figure 7 shows the six-dimensional reconstructed phase space projected onto a phase plane whose coordinates are the first and second largest components of the increment vector defined as ΔX(t + 1): = X(t + 1) − X(t). The 10 green arrows denoted by v(t + 1,j): = y(t + 1,j) − y(t,j) indicate their directions. The average of the 10 nearest neighbors and their one step evolutions is denoted by Y(t) and Y(t + 1), respectively. The increment vectors of the predictee ΔX(t + 1) and the average of the 10 nearest neighbors ΔY_AV(t + 1) are expressed in thick blue and red arrows, respectively. $X⌢(t+1)$ denotes the prediction of X(t + 1) and is defined as

Equation (16) is equivalent to Eq. (13) when p = 1.

For TD-FAST and Lorenz63A (Figs. 7(a) and 7(c)), the 10 green vectors v(t + 1,i) and ΔX_AV(t) move in the same direction as ΔX(t + 1) so that the prediction $X⌢(t+1)$ is close to its predictee X(t + 1). In these cases, DNP works as expected. On the contrary, for UW-FAST and Colored Noise (Figs. 7(b) and 7(d)), the 10 green vectors v(t + 1,i) move randomly and ΔX_AV(t + 1) is not related to ΔX(t + 1). Consequently, the prediction $X⌢(t+1)$ remains around the initial position of the predictee X(t). These two subplots explain why the short-term correlation coefficients of the DNP and ACF values for UW and UD are as high as Colored Noise. We also calculated the largest Lyapunov exponent (LLE), as it is one of the most popular measures for quantifying the orbital instability of a dynamical system. An LLE of 0.0098 is a reasonable value for Lorenz63. However, the other LLE values for Lorenz63A, Lorenz63B, and eight data are much larger and similar to the LLE's of Colored Noise. This outcome is consistent with the suggestion from many previous researchers^24,39 that chaos artifacts can be misled by autocorrelations.

The results shown in Figs. 6 and 7 also reveal the fundamental mechanism that leads conventional NTSA tools such as DNP and LLE to overestimate features of deterministic chaos when applied to noise-contaminated motion. This insight provides us a theoretical basis for adapting NDNP to resolve DNP's shortcoming.

We comprehensively applied NTSA to the soil cutting force data and the outcomes from the four methods are listed in Table I. We were most interested in the determinism and underlying dynamics of the data. The determinism was investigated using three methods: the translation error function (TEF), Fourier transform phase randomization surrogates (FPR), and Cao's E2. The TEF identified significant phase space continuity for TD-FAST, TD-SLOW, TW-FAST, and TW-SLOW. The FPR surrogate clearly rejected the null hypothesis of a linear Gaussian process for TD-FAST and TD-SLOW, with a 1% significance level. Coa's E2 indicated that TD-FAST, TD-SLOW, and TW-SLOW were not random data. For the second step, we proposed NDNP, which divides DNP by ACF. The NDNP values demonstrated that the short-term predictability detected in TD-FAST and TD-SLOW originated from their inherent dynamical systems.

All four outcomes were positive for TD and not positive for UD and UW (Table I). These results are clear, so we can say that TD was produced from a deterministic process and UD and UW are from stochastic processes. However, the results are not straightforward for TW. For TW-FAST, only E_trans(p) was positive, and for TW-SLOW only E_trans(p) and Cao's E2 were positive. FPR surrogate rejected the null hypothesis with a 5% significance level. Additionally, NDNP(p) is not larger than 1.0 but is sufficiently close to 1.0 to be counted as neutral. Therefore, we may say that, for TW-SLOW, the existence of determinism is recognized, but the noise effect is still significant.

The four soil conditions were prepared using irrigation and tillage, which change two major soil parameters: moisture and bulk density. Our results suggest that these two parameters are control parameters that determine the soil failure pattern as they influence more fundamental soil failure properties such as cohesion and internal angle of friction. Additionally, considering the above discussion regarding TW-FAST and TW-SLOW, the soil cutting speed is also a control parameter. For speeds in the range of 1.1–0.22 m/s, the influence of the speed on soil failure dynamics is smaller than that of moisture and bulk density. Beyond this range, the speed and therefore strain rate or visco-elastic effects may have more of an effect on the soil failure dynamics.

Particle based modeling such as DEM (discrete element modeling) and GD (granular dynamics)^40–43 is a very powerful technique for granular materials, because it simulates mechanical interactions between granular micro scale particles and their collective interactions with tools and machines in a macro scale. It is challenging to develop a reliable method for determining the mechanical parameters of each particle so that the model can represent macro level behavior such as different soil failure patterns (e.g., share plane failure, brittle failure, flow failure, bending failure, and tensile failure⁴⁴). Macro soil physical properties such as moisture, bulk density, and soil type must change along with the micro mechanical parameters of particles. The method proposed in this paper is useful for comparing DEM and/or GD models and field experiments and can be used to create possible evaluation functions for parameterization. Walker et al.¹² suggested that stick-slip dynamics are governed by the continual evolution of a so-called force chain, i.e., self-organized columnar structures comprising of particles, which carry the majority of the load in the system. Changing the physical micro parameters of particles should affect collective behavior, so we must crucially determine the sensitive physical micro parameters of granular materials. Our proposed approach for soil should be also applicable to other granular dynamics.

Soil tillage tools are classified into vibratory tillage and non-vibratory tillage. In vibratory tillage, tools are oscillated by tractor engine power.⁴⁵ Sakai⁴⁶ experimentally indicated that oscillating tools can reduce the draft force and total energy consumption for certain oscillating parameters (frequency, amplitude, and tractor velocity). The relation between these oscillation parameters and the soil's natural frequency is interesting. In non-vibratory tillage, tractors tow the tools at a constant speed. The tine type tillage tool oscillates due to self-excitation. The stiffness and shape of the spring tine are determined empirically in terms of energy consumption and soil pulverization. The relation between spring tine parameters and soil's natural spatial frequency should also be investigated and used to optimally design the tillage tools. Each tillage operation has a critical depth. Tillage operations that occur above this critical depth create a crescent failure zone.⁴⁷ When the tillage depth is deeper than the critical depth, the soil failure pattern at the cutting-edge changes from crescent failure to plastic flow failure in the region that is below the critical depth. This suggests that a combination of the soil's physical properties and tool (machine) geometries determines the state of the soil granular dynamics and can be viewed as a type of control parameter.

Soil is a typical granular material. The soil failure process observed during tillage is assumed to be a mixture of dynamical processes and stochastic processes (e.g., a linear Gaussian process). In this paper, we focused on detecting determinism and orbital instability in soil cutting force data, because these two features are the most important properties of deterministic chaos. Soil conditions depend on changes in soil moisture content and soil bulk density, so we prepared four different conditions based on these factors. Only TD (tilled dry condition) showed significant determinism through the translation error function and Fourier transfer phase randomization surrogate (FPR). DNP generally overestimated the underlying dynamics of the deterministic chaos because of autocorrelations in the data sets. To overcome this, we proposed NDNP to reduce the effect of autocorrelative noise. We demonstrated that NDNP can clarify whether the short-term predictability expressed by DNP is due to autocorrelative noise or an underlying dynamical system. We suggest that NDNP(p) > 1 is a new necessary condition of deterministic chaos and is especially useful for noise contaminated dynamics. Our proposed approach could also be used to investigate the collective failure patterns of granular material in terms of nonlinear dynamics.

Kenshi Sakai thanks the Japanese Ministry of Education, Culture, Sports, Science and Technology (MEXT) long-term study abroad program at University of California, Davis in 2001–2002 for initiating this collaborative project with Shrini Upadhyaya and Pedro Andrade-Sanchez. Professor William J. Chancellor provided us with insightful comments, suggestions, and encouragements. Kenshi Sakai is grateful to Professor Richard Goodwin for his invaluable suggestions during the course of this study. Nina Sviridova thanks the Tokyo University of Agriculture and Technology intensive study abroad program and MEXT for committing to this project. This work was supported by JSPS Grant-in-Aid Nos. 25660204 and 15H04572.