Wave Dispersion in One-Dimensional Nonlinear Local Resonance Phononic Crystals with Perturbation Method

: Nonlinear phononic crystals are receiving increasingly greater attention in the ﬁeld of sound absorption and vibration reduction. In this paper, we use the perturbation method to investigate elastic wave propagation in one-dimensional discrete local resonance nonlinear phononic crystals. The nonlinear force on the inner resonator is expressed in the form of a linear part plus a cubic nonlinear ﬂuctuation. By combining Bloch wave theory and the perturbation method, the nonlinear dispersion relation is obtained by a ﬁrst-order approximate analytical solution. The results show that the band’s cut-off frequency is not only affected by the degree of nonlinearity but is closely related to the wave amplitude. In addition, the ﬁnite element method is used for comparison and veriﬁcation. Finally, an application example of a wave ﬁlter is provided based on the nonlinear characteristics.


Introduction
In 2000, Liu [1] put forward the concept of mass-in-mass structure phononic crystals, which are also called local resonance phononic crystals. The local resonance mechanism is different from Bragg scattering, but is determined by the strong resonance characteristics of the inner resonator inside the outer substrate. In a linear periodic structure, the principle of local resonance acoustic metamaterials/phononic crystals is to confine the elastic waves of relevant frequencies to an inner spring-mass system. The starting frequency of the acoustic band gap is around the natural frequency of the resonator [2,3]. The waves attenuated in the frequency range near the resonance frequency form a forbidden band gap. For decades, many experiments and theoretical analyses have been carried out around this new type of material [4,5]. The reasonable design of parameters such as spring stiffness and mass can realize the function of adjusting the band gap interval. Local resonance phononic crystals can produce band gaps at larger wavelengths with smaller lattice sizes. As a result, they have important application prospects in low-frequency sound insulation and vibration reduction.
Although the linear periodic structure has many interesting dynamical characteristics, nonlinear factors cannot be ignored in engineering practice. As the excitation wave amplitude increases, the response of the material/structure usually tends to be nonlinear, and the wave propagation becomes more complicated. By studying the nonlinearity of the periodic structure, we can reveal more peculiar phenomenon compared with the linear case, which includes nonlinear resonance, bifurcation, mixing, or self-trapping. Based on these topics, nonlinear devices have potential new applications, such as frequency conversion and energy harvesting.
In 2009, B. Ling [6] combined phononic crystals and nonlinearity and proposed a new mechanism of acoustic rectification, which realizes the one-way propagation of sound amplitude. With the introduction of new rectification ideas for acoustic diodes, scientific researchers from all over the world have paid close attention to and conducted in-depth research into the acoustic rectification effect. N. Boechler et al. [7] proposed the use of one-dimensional nonlinear chains to achieve the asymmetric control of sound waves. There are many numerical methods for studying nonlinear wave phenomena. For nonlinear continuous systems, there are already some theories [8,9]. However, for discrete systems, Vakakis [10] introduced a continuity assumption and used multi-scale analysis techniques to describe a periodic nonlinear oscillator. The results showed that the boundary frequency of the band is dependent on the amplitude. In addition to the perturbation method, another mapping technique [11] was proposed to determine the band gap of the discrete chain nonlinear harmonic oscillator. The proposed nonlinear transfer matrix contains amplitude and nonlinear parameters.
In this article, a perturbation analysis method is applied to the discrete nonlinear local resonant periodic structure to predict the amplitude-dependent dispersion relation. The displacement and frequency are expressed as a first-order approximation and substituted into the differential equation for solution. Different from the method of expanding the wavenumber in linear processing [12], the dispersion relation in the form of the amplitude and degree of nonlinearity is derived. The simulated transmittance of a finite one-dimensional lattice is compared with the dispersion relation. After introducing the nonlinear spring, the change in the transmittance curve is studied. Finally, an application example based on the proposed nonlinear phononic crystal is provided.

Dispersion Relation of Linear Mass-in-Mass Lattice Model
For a one-dimensional linear mass-in-mass phononic crystal, the unit cells (see Figure 1a) are connected by spring k 0 periodically at a spacing of L. The mass of the outer substrate is m 0 while the mass inside is m 1 in one unit cell. The internal spring-mass structure is regarded as a resonator with its local resonance frequency ω 1 = √ k 1 /m 1 , where k 1 is the spring stiffness. By introducing harmonic wave propagation into the mass-in-mass system, the equations of motion for the j-th unit cell can be expressed as where u j,α indicates the displacement of mass α in the j-th cell (α = 0, 1).
The harmonic waveform is where B α is the complex wave displacement amplitude, q is the wave number, and ω is the angular frequency. By combining these equations above, the matrix form of the equations of motion can be obtained as The next step is to solve the eigenvalues of the coefficient matrix. The determinant of the coefficient matrix in Equation (4) should be equal to zero. The dispersion relation can be expressed as If the real Bloch wave vector qL is given, the characteristic frequency ω can be solved; this method is called the ω(q) method. While the characteristic frequency ω is given, the Bloch wave vector qL can be expressed in a complex form; this method is called the q(ω) method.
An example is now presented to illustrate the band gap characteristics of a onedimensional periodic mass-in-mass system (see Figure 2) with these two methods. The parameter ratios are m 1 /m 0 = 2 and k 1 /k 0 = 2. The vertical axis of Figure 2 shows the dimensionless frequency ω 1 = ω/ω 1 . Figure 2a contains the results solved by the ω(q) method and the real part solved by the q(ω) method. Meanwhile, Figure 2b contains the imaginary part of the q(ω) method. Obviously, the curves obtained by the two methods are completely overlapped. From the real part of the wave vector Re(qL), we can observe a prominent band gap between the optical and acoustic modes (upper and lower dot-dashed line, respectively). The attenuation characteristics of the band gap can be observed from the variation of the imaginary part of the wave vector Im(qL). There is a sharp peak in the imaginary part of the wave vector, which represents the local resonance frequency (marked in light green at ω 1 = 1 in Figure 2b) of the internal mass-spring system. The starting frequency is close to the local resonance frequency when the local frequency of the internal resonator √ k 1 /m 1 is sufficiently low relative to the natural frequency of the external mass √ k 0 /m 0 .
The ending frequency of the first band gap should satisfy the prior that both the real and imaginary parts of qL are equal to zero, and we have The width of the band gap can be described as the difference between the upper and lower band frequencies.
In Figure 2b, the peak frequency of the imaginary dispersion curve lies on the resonant frequency of the internal mass-spring system. The ending frequency is √ 3, which satisfies the above derivation when m 1 /m 0 = 2.

Nonlinear Dispersion Relations
In this section, we investigate the band gap characteristics of the nonlinear lattice structure shown in Figure 1b with perturbation analysis. The stiffness coefficients of a nonlinear spring include a linear term and nonlinear fluctuation. The weak nonlinear internal force of the spring on the outer mass is expressed as Among them, δ represents the relative displacement of the spring, ε represents a small parameter, and Γ 1 is a parameter that controls the degree of nonlinearity. Introducing nonlinearity in the internal mass-spring resonator, the equations of motion for the j-th unit cell can be expressed as where u j,α (α = 0, 1) denotes the displacement of each mass α in the global coordinate system for the j-th unit cell. Then, the perturbation technique is applied to obtain the dispersion curves corrected up to the first order. Displacements, as well as wave frequency, can be denoted as asymptotic series u j,α = u j,α (0) + εu j,α (1) + o ε 2 (12) The superscripts (0) and (1) represent the linear part and the first-order approximation, respectively. In order to facilitate the calculation, the dimensionless parameters are introduced as τ = ωt, ω 0 = ω/ω 0 , ω 1 = ω/ω 1 , η = k 0 /k 1 , λ = m 1 /m 0 . Then, Equations (10) and (11) turn into Substitute Equations (12) and (13) into Equations (14) and (15) and collect together all the terms of the same hierarchy. Hence, the subproblem of the system yields to Without loss of generality, the displacement of the α-th dimension in the j-th cell can be broken down into a set of conjugate forms as For the mass with subscript α, the position parameter for one built-in resonator lattice system can be expressed as Substituting the displacement format (18) into ε (0) , we have By solving Equation (20), we can obtain the linear dispersion relation as Replacing the dimensionless parameters with original expressions, the updated linear dispersion relation is in the format which is in the same form as Equation (5). Thus, we have successfully obtained the linear dispersion relation from the first-order term of the nonlinear expression.
In order to explore the influence of the weak nonlinearity in the spring parameters on the band gap characteristics, the coefficient equation in front of ε (1) can be written as The form on the left side of these Equations (24) is similar to the expressions in front of ε (0) . In order for the equations to be true, the coefficient before the exponential term on the right side must be zero, i.e., b 1 = 0. Therefore, the first-order correction term can be solved as ; thus, the final dispersion relation with first-order correction can be expressed as where . It needs to be mentioned that the amplitude ratio ξ takes the amplitude of the inner mass as a reference and dominates the trend of the acoustic mode. The relationship between the amplitude ratio and dimensionless wave number is shown in Figure 3. It can be observed that A 0 , which means the vibration of the inner mass is larger than the outer mass. This means the inner mass takes on the main role of energy consumption. This is also the original intention of imposing local resonators on the host structure. The energy is mainly dissipated by the vibration of the internal mass-spring system.
We can also obtain another form of the nonlinear dispersion relation as In this case, χ takes the amplitude of the outer mass as a reference and dominates the trend of optic mode. The displacement ratios ξ and χ are both indirect calculation results from formula (20). The dispersion relations are expressed in two forms (Equations (26) and (27)) for convenience to understand some extreme conditions. For example, when the vibration amplitude of the external mass is close to zero, i.e., |A 0 | → 0 , it is easier to track the changing trend of nonlinear dispersion curves compared with the linear ones by applying formula (26). When taking the amplitude of the external mass as a reference, the ratio χ > 1 is substituted into formula (27) to better understand the trend of the dispersion relation with the introduction of nonlinearity varying with |A 0 | and Γ.

Effect of the Degree of Nonlinearity
We use the structural parameters in Figure 2 to calculate the dispersion relation. It is easy to identify from the nonlinear force f = kδ + εΓδ 3 , including cubic term, that when the degree of nonlinearity is greater than zero (Γ > 0), it performs a nonlinear hardening system. In addition, when the degree of nonlinearity is smaller than zero (Γ < 0), it performs a nonlinear softening system. Figure 4 describes the variation of the dispersion relation with the degree of nonlinearity. The flat blue arrow indicates the directional change of the starting frequency with Γ. When Γ > 0 (see Figure 4a), both the optical mode and acoustic mode curves move upward. The relative curves are all above linear conditions (black dashed line). The rising amplitude grows with the increase in the wave number, resulting in a narrowing of the band gap. Conversely, when the degree of nonlinearity is less than zero (Γ < 0, see Figure 4b), the curves are below linear conditions. The bound frequency would decrease, which leads to a wider band gap interval.
It is very interesting that due to the introduction of nonlinearity in the lattice system, a band gap related to the nonlinear parameter Γ is formed. For the case of Γ > 0, the width of the band gap gradually decreases with the increase of Γ. If Γ is large enough, this band gap may disappear when the acoustic mode is stretching up to the optic mode. In this situation, the system is considered strongly nonlinear. For the other case of Γ < 0, the starting frequency (see the red thin arrow) of the band gap is changeable due to the non-monotonic increase in the acoustic branch of the dispersion curve. If the curve does not conform to monotonicity, the starting frequency is located at the point where the slope of the acoustic mode is zero. To conclude, a hard built-in nonlinear resonator (Γ > 0) ought to be shrinking the band interval between the acoustic and optic modes, while a soft built-in nonlinear resonator (Γ < 0) does the opposite.

Effect of the Wave Amplitude
The introduction of nonlinearity into the resonator not only makes the band gap adjustable and achievable but also brings new amplitude-related wave propagation characteristics. The influence of wave amplitude on the nonlinear dispersion relation has a similar conclusion to that of the degree of nonlinearity. It can be observed from the expression of the nonlinear dispersion curve (see Equations (26) or (27)) that the first-order fluctuation of first-order dispersion frequency is proportional to the square of the amplitude.
The position of the band gap is not only related to the degree of nonlinearity but also affected by the wave amplitude, as shown in Figure 5. The flat orange arrow indicates the direction of change of the starting frequency with the amplitude of outer mass |A 0 |. For the hardening nonlinear lattice structure (see Figure 5a), the starting frequency of the band gap increases as the amplitude of the wave increases, resulting in a decrease in the first band gap interval. For the softening nonlinear lattice structure (see Figure 5b), the starting frequency shows a downward trend until the amplitude increases to a certain limit. It is when the acoustic mode curve no longer maintains monotonicity that the starting frequency is located at the point where the slope of the acoustic brunch is zero. Thereby, we have proposed a one-dimensional nonlinear phononic crystal with a band gap characteristic related to the vibration amplitude.
In summary, both the degree of nonlinearity and wave amplitude may influence the distribution of the dispersion curves. A hard nonlinear resonator will shrink the band gap. A softening nonlinear system is able to achieve the aim of extending band intervals within a limitation due to the non-monotonicity of the acoustic branch in the dispersion relation. By carefully arranging the relevant parameters, an adjustable band gap is available with a nonlinear phononic crystal.

Transmittance with Linear Spring Constant
The finite element analysis method (FEM) is a recent, effective, and convenient way to verify structural characteristics [13][14][15]. The large-scale application of finite element software has provided researchers with great convenience in terms of accurate calculation and time-saving. We are using the commercial software COMSOL Multiphysics to mimic the wave transformation in a phononic crystal structure.
The multi-body dynamics module in COMSOL Multiphysics can simulate systems containing various types of joints, as well as flexible and rigid bodies. For example, it can design flexible parts using nonlinear materials. Using the multi-body dynamics module of two-dimensional drawing, the mass-in-mass structure can be drawn more intuitively. The schematic diagram of the mass-in-mass phononic crystal using the multi-body dynamics module in COMSOL Multiphysics is shown in Figure 6. The outer mass is represented by a concentric circle with an outer diameter of 0.4 m and inner diameter of 0.32 m. The internal mass is a circle with a radius of 0.2 m. The thickness of the two-dimensional geometry is set as 1 m. The lattice constant is also set as 1 m for the convenience of calculation. The attached material properties to the inner and outer mass are shown in the description in Figure 6. These parameters are carefully chosen to produce an example of an acoustic band interval. Since these two masses are rigid bodies, their mimic density is easy to calculate based on the given mass quality and volume. Our purpose is to solve the transmittance of waves of different frequencies after passing through the lattice structure. To ensure calculation accuracy, at least 10 units are used in the calculation.
Next, two spring models are created. Using the spring damper module in multi-body dynamics, a single spring can be formed by connecting the source and target points. In total, 10 internal springs and 9 external springs are created due to the limited periodicity in model construction. The first spring model is set as a constant, and the second is expressed as a function of the elongation. The former can be applied to solve the linear transmittance, and the latter is used to calculate the nonlinear one.
Probes are added to the input and output points of the designed multi-body structure, and the built-in expression of each probe is set as the displacement in the direction of periodic extension. In the multi-body dynamics module of the software, the mass is treated as a rigid body, so there is only translational displacement. A spring is an element that connects two rigid bodies, and its elastic force is affected by the relative distance between the two rigid bodies. Because we only consider the one-dimensional vibration mode, the vibration tendency of the spring and mass is restricted to the direction arranged along the period (x-axis). The displacement deviate from the periodic direction is not considered. Thus, the spring-mass multi-body dynamics system constructed by FEM vividly imitates the one-dimensional mass-in-mass phononic crystal structure with finite period. It can be regarded that there are real springs between the inner and outer masses, and between unit cells. The springs and masses here are not necessarily concrete materials in real life but virtual finite element models that provide them physical parameters. We use displacement excitation. We specify the unidirectional displacement along x-axis at the input point (see the input probe in Figure 6) as 1 m and conduct the frequency domain analysis. The frequency-domain analysis calculates the steady-state solution of the system at different frequencies. The considered frequency range in this simulation is 50-450 Hz, which is the approximate frequency range estimated from the analytical solution, including the band gap interval (see Figure 7a). The expression of the transmittance is Tr = 20 log 10 ( d out d in ) (28) Among them, d out is the displacement collected from the output probe after the frequency domain analysis and d in is the initial input displacement excitation of the input probe. Then, the transmittance curve can be drawn, as shown in Figure 7b.
In Figure 7a, the calculated first band interval (126~276 Hz) between the acoustic mode and optic mode is marked in light blue shadow, and the second band interval (above 400 Hz) is marked in light red. The frequency interval where the transmittance calculated by COMSOL is less than zero corresponds to the forbidden band interval obtained from the analytical solution. It needs to be mentioned that the frequency of the transmittance peak corresponds to the position of the peak of the imaginary part of the theoretical solution (159 Hz). This indicates that the transmittance curve and imaginary part of the dispersion relation express both the speed and tendency of the wave propagating in the phononic crystal. They are consistent in revealing the principle of sound absorption.
In addition, we plot displacement cloud images with frequencies of 80 Hz and 350 Hz outside the sound insulation interval and 150 Hz and 270 Hz inside the interval (see Figure 8). The first and last cloud images show that each cell is in a state of vibration, and there is no obvious energy loss when the frequency is not in the forbidden band. At the frequency 270 Hz near the boundary of the sound insulation interval, the input energy is gradually consumed after about five cells, and the energy transferred to the end of the chain reaches zero. Therefore, the lattice chain structure achieves the effect of energy absorption in the sound insulation frequency range. Typically, at the frequency 150 Hz near the peak of the sound insulation interval, the input energy is basically consumed by the local resonator of the first unit cell, and the energy transferred to the second cell is almost zero. This phenomenon illustrates the local resonance mechanism. It effectively explains why the absorption peak point corresponds to the local resonance frequency of the unit cell. The rationality of the local resonance method in terms of energy absorption has been perfectly embodied. The above results show that using the multi-body dynamics module to simulate the lattice chain structure is accessible. The calculation of the transmittance can also obtain the band gap characteristics of the one-dimensional phononic crystals instead of the dispersion relation, with the accuracy of the calculation being guaranteed. In addition, the wave propagation mechanism can be explained with cloud images from simulation.

Transmittance with Nonlinear Spring Coefficient
In this subsection, we use this simulation method to calculate the nonlinear lattice system. Since the spring force can be expressed as a function of the elongation part in the spring damper module, we introduce nonlinear force f = k 1 δ + εΓ 1 δ 3 into the inner spring.
The parameters referred to in Section 4.1 are taken into consideration. Figure 9a shows the analytical dispersion curves of different wave amplitudes |A 0 | = 1, 2.5, 5. The starting frequencies of each band gap are 129 Hz, 145 Hz, and 203 Hz separately. We apply the displacement of the input probe d in to calculate the nonlinear dispersion relation. The linear band gap interval or sound insulation area is marked in blue with a transparency of 0.3. The nonlinear band gap intervals are marked in green with a transparency of 0.3 individually. If the intervals overlap, the color transparency increases. Obviously, for a linear system, the distribution of the transmittance curve will not be affected by the initial displacement of the input probe with a sound insulation range of 126~276 Hz. In the nonlinear simulation, we choose three sets of variables (d in = 1, 2.5, 5 m) and obtain the transmittance, as shown in Figure 9b. The invariants include the degree of nonlinearity (Γ = 1) and a small nonlinear parameter (ε = 0.005). It turns out that the change of the initial input displacement perturbation does not affect the distribution of the sound insulation interval (144~276 Hz) of the nonlinear system either. As a result, the simulation method cannot reflect the influence of wave amplitude on the transmittance of the nonlinear system. Consequently, the input displacement d in in the simulation cannot be equal to the wave amplitude |A 0 | in the analytical solution. Therefore, it is not possible to simulate nonlinear systems under different wave amplitudes. Although the finite element simulation method has disadvantages, it is still evident that the sound insulation interval narrows after the system changes from linear to hardening nonlinear. The trend of the band gap change caused by nonlinearity coincides with the analytical results.
Researchers somehow find it challenging to realize the nonlinear device with accurate parameters in practice. On the one hand, the small parameter ε can usually be obtained by fitting the force-displacement curve in an experiment. It has no real physical meaning but evaluates how weak the nonlinearity is. The smaller it is, the weaker the nonlinearity. By contrast, another variable parameter Γ has physical meaning. It is a parameter that characterizes stiffness in the nonlinear form, and the dimension is the ratio of the force to the cubic of the displacement.
On the other hand, the real material is defined by the elastic modulus, while the tensile parameter in the simulation is characterized by the spring constant. Although it is easy to control the mass magnitude for each material, it is hard to define the effective spring constant. Scholars have made many attempts in parameter equivalence [16] and parameter identification [17] research. For example, a classic mass-in-mass structure is already provided by Liu [1]. The composite material is composed of a solid core material with a relatively high density and an elastic soft material coating. The core material corresponds to m 1 in our structure. The coating silicone rubber can be regarded as k 1 , with epoxy as the hard matrix material, which we can refer to as m 0 and k 0 . The equivalent of the nonlinear spring model usually arranges two vertical springs on the basis of the original linear spring-mass structure [9].

Application Example Based on Nonlinear Dispersion Relations
For a nonlinear hardening system, although there is no advantage in band gap widening, there are still some application values, such as a filter or a unidirectional conduction device [7,18,19]. By introducing parameters that may produce an acoustic stop band, wave filters can be realized based on the amplitude-dependent characteristics of the nonlinear phononic crystal. Taking the parameters described in Figure 9 as an example, we assume that the input signal frequencies of the nonlinear phononic crystal structure include 140 Hz and 170 Hz. Here are three typical cases.
(1) When the adjustable signal amplitude is |A 0 | = 1, the dispersion relation can be observed as red solid lines in Figure 9a. Waves with both 140 Hz and 170 Hz frequencies are forbidden to pass through because they happen to be within the band gap. As shown in Figure 10a, the waves of these two frequencies are totally filtered. (2) When the signal amplitude is |A 0 | = 2.5 (see blue chain lines in Figure 9a), waves with a frequency of 140 Hz are located in the acoustic mode of the dispersion curve, and waves with a frequency of 170 Hz are in the place of the band gap interval. Thus, waves with a frequency of 140 Hz can pass, and those with a frequency of 170 Hz are prohibited. The input waves and filtered frequencies are shown in Figure 10b.
(3) When the wave amplitude increases from |A 0 | = 2.5 to |A 0 | = 5 (violet dotted lines in Figure 9a), the forbidden band is narrowed. Both waves with 140 Hz and 170 Hz frequencies can pass because they are all seated in the acoustic brunch of the dispersion relation. As a result, the nonlinear phononic crystal loses its filtering function in this case. In Figure 10c, the output signal still includes these two frequencies.
Therefore, we have designed a nonlinear phononic crystal structure that can selectively allow the required signals to pass by adjusting the input wave amplitude. It promotes a new idea for designing a functional acoustic wave filter that can selectively achieve the unidirectional transmission of sound waves. It can also provide potential applications in vibration suppression and noise control of mechanical equipment.
The aim of changing the wave amplitude can be achieved through the introduction of a variable cross-section rod. Let the axial direction be x and the cross-sectional area be S(x), then Assuming that the wave propagation form in a variable cross-section rod is u = ue −ikx e iωt , it can be solved that u = e Obviously, the change of wave amplitude from one side of the variable cross-section rod to the other is related to the ratio of the cross-sectional width β 0 . A reasonable choice of β 0 can adjust the amplitude of elastic waves.

Conclusions
In this paper, we focus on the dynamic behavior of a one-dimensional local resonance phononic crystal system containing cubic nonlinear inner springs. A perturbation analysis method is applied to calculate the dispersion relation of this discrete nonlinear periodic structure. By introducing nonlinear parameters, the nonlinear dispersion relation is fluctuated based on the linear one. The boundary frequencies of the forbidden band can be adjusted by changing both the degree of nonlinearity and wave amplitude. When the degree of nonlinearity is greater than zero, the system exhibits hardening nonlinearity, and the band gap width shows a decreasing trend. On the contrary, the system presents softening nonlinearity when the degree of nonlinearity is less than zero and the band gap width is widening. Therefore, using nonlinear characteristics, we can design a phononic crystal with an adjustable band gap that is dependent theoretically on the amplitude or the degree of nonlinearity. Then, a simulation method is applied to calculate the transmittance through finite unit cells with both a linear and nonlinear local resonator. The local resonance is well explained with the displacement cloud image inside the forbidden band interval. It turns out that the nonlinearity would also make the sound insulation area in the transmittance curve move. This simulation result contributes to the analytical solution for the nonlinear phononic crystal. In addition, an amplitude-dependent, one-way nonlinear filter is realized. As a result, this nonlinear local resonance phononic crystal not only provides a new design for adjustable energy band structure but also inspires the design of acoustic components.
Author Contributions: The contribution of the paper is as follows. Y.C.: conceptualization, methodology, software, writing-original draft preparation; G.L.: data curation, writing-reviewing and editing; R.S.: visualization, investigation; G.C.: supervision. All authors have read and agreed to the published version of the manuscript. Data Availability Statement: All data included in this study are available upon request by contact with the corresponding author.