Transient and Steady Sliding Friction of Elastomers: Impact of Vertical Lift

The sliding friction of elastomers has been investigated numerically and theoretically for the line contact between a cylindrical rigid indenter and a “frictionless” Kelvin–Voigt foundation. The onset of sliding under an abrupt increase in the drive velocity has been simulated with different boundary conditions of the rigid indenter. When the rigid indenter is not allowed to move in any direction, just an abrupt change in the friction force appears, which is not accompanied by any transient processes. However, when the rigid indenter is able to move in the vertical direction, the transient sliding friction including three different time constants appears, resembling the typical transition from the static friction to the kinetic friction, in spite of no static friction considered in the simulation. The aforementioned drastic difference is caused by the “vertical lift” of the rigid indenter, which is induced by the damping of the Kelvin–Voigt foundation. In addition, the vertical lift strongly affects the characteristics of the steady sliding friction, which is explained well by using the critical velocity determined from the asymptotes in the master curve of friction coefficient.


INTRODUCTION
It has been well-known that sliding contact of elastomers involves various types of dynamics due to friction. Physical phenomena such as wave propagations (Schallamach, 1971;Barquins, 1985;Rubinstein et al., 2004;Maegawa and Nakano, 2010), wear pattern formations (Schallamach, 1957;Fukahori and Yamazaki, 1994), and nonlinear vibrations (Nakano and Maegawa, 2009;Yamaguchi et al., 2011;Nakano et al., 2019) have attracted the interest of many scientific and engineering researchers. Among various elastomers, rubber is the material that has received the most attention in various practical applications (e.g., tires, seals, and shoes). Since the study done by Grosch (1963) showing the master curve of friction coefficient according to Williams-Landel-Ferry (WLF) theory (Williams et al., 1955), the importance of viscoelasticity has been recognized, and the dependence of the friction coefficient on temperature and velocity has been investigated [e.g., Popov et al. (2018)]. More recently, several swollen polymers showing high elasticity and ultra-low friction (e.g., hydrogels and polymer brushes) have been found (Gong et al., 2001;Nomura et al., 2011), some of which are being developed toward practical applications (Belin et al., 2018;Tadokoro et al., 2020), and their tribological properties have also been discussed in relation to their viscoelastic properties (Mizukami et al., 2019).
To understand the sliding friction of elastomers, various types of modeling methods have been proposed. Among them, those with "viscoelastic foundations" are known to have strong advantages not only of avoiding the difficulties of elastic contact stress theory but also of providing intuitive pictures of how the energy dissipation occurs inside the contact. The first was the extension of the elastic foundation (i.e., the Winkler foundation) from stationary contact problems of thin films to steady sliding contact problems of elastomers (May et al., 1959;Johnson, 1985). Then recently, in the novel framework of the "Method of Dimensionality Reduction" invented by Popov and his colleagues, the modeling method has progressed, enabling us to analyze the true three-dimensional contact with high accuracy (Popov, 2010;Popov and Heß, 2015;Kusche, 2017). For example, Li et al. (2015) numerically and theoretically studied the kinetics of friction coefficient for the sliding contact between a flat elastomer (modeled by a Kelvin-Voigt foundation) and a rough rigid indenter (modeled by a self-affine fractal) under abrupt change in the drive velocity. As a result, they found interesting temporal changes (i.e., jumps and relaxations) of the friction coefficient depending on the drive velocity and the Hurst exponent of the self-affine fractal, under the consideration of quasi-static processes.
Based on the foregoing background, the aim of this study was to find the answer to the following question: What is the minimum requirement for modeling with viscoelastic foundations to describe the transient friction appearing in the onset of sliding between an elastomer and an indenter? In this article, to provide a possible answer in a minimal situation, we consider a simple line contact between a flat elastomer and a cylindrical rigid indenter, focusing on the boundary condition of the rigid indenter in the vertical direction, under the restriction that the elastomer is modeled by the conventional Kelvin-Voigt foundation. Two types of models with different boundary conditions reveal that the vertical dynamics of the rigid indenter strongly affects not only the occurrence of the transient sliding friction but also the characteristics of the steady sliding friction.
Note that the aforementioned idea on the boundary condition in the vertical direction originates from the typical structure of sliding systems under the fluid film lubrication, since when the fluid film lubrication is considered theoretically, it is natural to assume that the "slider" (corresponding to the "indenter" in this study) shifts vertically to find the position to balance with the applied normal load. Also note that in general, it has been known that the friction of elastomers arises by two different mechanisms. One is the "adhesion, " the energy dissipation of which occurs at the contact interface between the elastomer and indenter [i.e., the adhesive friction (e.g., Moore and Geyer (1972))], and the other is the "viscoelasticity, " the energy dissipation of which occurs inside the elastomers [i.e., the hysteresis friction (e.g., Moore and Geyer (1974))]. Although both mechanisms are important and should be interrelated, this study ignores the former and focuses on the latter. Figure 1 shows the two types of models for sliding friction of elastomers considered in this study. The left is termed the "fixed indenter (FI) model, " and the right is termed the "movable indenter (MI) model." They are models describing the two-dimensional sliding contact between a rigid indenter and a viscoelastic foundation in the xz plane, where the x and z axes are taken in the horizontal and vertical directions, respectively.

MODELS Structures
The rigid indenter (mass per unit width: M) has a cylindrical shape (curvature radius: R), the bottom surface of which contacts with the viscoelastic foundation. In the FI model, the rigid indenter is mounted on the rigid walls, not to be able to move in any direction. In the MI model, on the other hand, the rigid indenter is supported by a frictionless linear bearing, to be able to move only in the vertical direction. Note that the boundary condition of the rigid indenter is the single difference between them.
The viscoelastic foundation consists of an infinite number of viscoelastic elements mounted on a rigid base at regular intervals in the horizontal direction. Every viscoelastic element is a one-degree-of-freedom Kelvin-Voigt element consisting of a vertical spring (stiffness per unit width: k) and a vertical damper (damping coefficient per unit width: c) with the same natural length; as a result, its non-disturbed surface is the horizontal plane z = 0. Note that by using the stiffness per unit area (K) and damping coefficient per unit area (C) of the viscoelastic foundation, k and c are given by k = K/N and c = C/N, respectively, where N is the number of viscoelastic elements per unit length. Also note that here we assume that the upper end of every viscoelastic element makes contact with the rigid indenter surface with no adhesion and with no friction.
Owing to the aforementioned structures, the normal load per unit width (W) of the FI model is determined by the indentation depth δ [i.e., W = W(δ), where δ is constant], while that of the MI model is determined by the gravity [i.e., W = Mg, where g is the gravity constant]. In addition, in both models, the rigid base is driven horizontally at the drive velocity V. Therefore, we can say that the FI model represents a "constant-gap sliding system, " while the MI model represents a "dead-weight sliding system." Or, we can also say that the FI and MI models represent the two limiting cases corresponding to a "very stiff apparatus" and a "very soft apparatus, " respectively.

Governing Equations
When the rigid indenter penetrates the viscoelastic foundation by δ > 0, the indenter surfaces in the FI and MI models are given by respectively, where t is the time, and h(x) is the indenter shape function: Let us assume that the ith viscoelastic element at x = x i (t) makes contact with the indenter surface (see Figure 2). The compression and compression rate of the ith viscoelastic element in the FI model are given by  respectively, while in the MI model, respectively, where ( · ) and ( ′ ) are the derivatives with respect to t and x, respectively. Then, as shown in the left of Figure 2, there are four types of local forces acting on the subelement PQR: the normal force f Ni at P (from the rigid indenter), the restoring force ku i at Q (from the spring), the damping force cu i at R (from the damper), and the horizontal constraint force f Ci (from the spring and damper). Considering that the subelement PQR is massless, we obtain the following force balance equations in the vertical and horizontal directions: under -π/2 < θ i < π/2, and from Equation (8), Now, as shown in the right of Figure 2, the normal force acts on the rigid indenter at P from the ith viscoelastic element, which is the reaction force of the normal force f Ni on the left of Figure 2. Therefore, the vertical and horizontal components of the normal force are given by respectively. Note that from Equations (9) and (13), the horizontal constraint force is determined as f Ci = f xi . Consequently, the total vertical force F z and total horizontal force F x acting on the rigid indenter from the viscoelastic foundation are given by respectively. Finally, the vertical position of the indenter bottom A in the FI model is given by while in the MI model, it is determined by the equation of motion: which determines δ andδ in Equations (6) and (7) as follows: Note that when the ith viscoelastic element is not making contact with the indenter surface, it is obvious that f Ni = 0 and therefore f zi = f xi = 0. Then, the compression rate of the non-contacting element is given bẏ where τ is the retardation time of the viscoelastic element, defined as or, by using the macroscopic properties of the viscoelastic foundation, Therefore, when the contacting element violates the following condition: it detaches from the rigid indenter surface. Note that the firstorder differential equation (20) for the non-contacting element has the following solution: where t 0 is the constant.

METHODS
The sets of the governing equations in the previous section were solved numerically for the FI and MI models. To solve the second-order ordinary differential equation (17) for the MI model, the Runge-Kutta method was used, where the time discretization was determined by t = 2 x/V = 2/NV. Standard parameter values for the numerical simulations are listed in Table 1. Note that both models include six independent parameters (R, K, τ , N, δ, and V in the FI model; R, M, K, τ , N, and V in the MI model). However, since N (= 1/ x) is the parameter for the space discretization, the number of essential parameters is five for each model. Also note that if we imagine a "thin" elastic sheet of thickness h = 1 mm as the viscoelastic foundation, it is possible to say that the stiffness per unit area K = 10 GN/m 3 corresponds to the effective elastic modulus E ∼ Kh = 10 MPa, although it depends strongly on the boundary condition of mounting the elastomer on the rigid base [e.g., Popov (2010)]. Figure 3 shows the numerical results of temporal changes in representative variables for the FI model (left column) and MI model (right column) under the standard conditions ( Table 1). Note that as shown by the blue lines in the top row, the drive velocity was abruptly increased from V = 0 to 100 mm/s at t = 0. First, from the results for the FI model, we find that the indenter position provided the initial condition (z A = −0.3 mm, F z ∼ 10 N/mm, and F x = 0). Then, F z and F x were immediately increased at t = 0 and maintained the increased values (z A = −0.3 mm, F z ∼ 35 N/mm, and F x ∼ 5 N/mm).

RESULTS
Then, from the results for the MI model, we find that the gravity force provided the same initial condition as that of the previous model. However, the responses were completely different, showing distinct "transient sliding" toward "steady sliding." The rigid indenter started to move "upward" at t = 0 and then gradually approached z A ∼ −0.1 mm with a time constant in the order of "10 ms" (second row). The response of F z was spiky: it immediately increased to F z ∼ 30 N/mm at t = 0, then rapidly decreased to F z ∼ 9 N/mm with a time constant in the order of "1 ms, " and gradually returned to F z ∼ 10 N/mm (third row). The response of F x was also spiky: it immediately increased to F x ∼ 5 N/mm at t = 0, then rapidly decreased to F x ∼ 2 N/mm with the time constant "1 ms, " and gradually approached F x ∼ 1 N/mm with the time constant "10 ms" (bottom row). It should be noted that the response of F x resembles the typical transition from the static friction to the kinetic friction, although no static friction was considered. This transient behavior will be discussed in the next section. Figure 4 shows the numerical results of spatiotemporal changes in the contact pressure p(x,t) = Nf N (x,t) under the same conditions as those of the previous figure, where the abscissa  Table 1 for parameter values).  Table 1 for parameter values). LE, leading edge of contact; TE, trailing edge of contact.
Frontiers in Mechanical Engineering | www.frontiersin.org is the horizontal position x, the ordinate is the time t, and the magnitude of p is represented by the shade of the red color. The blue lines show temporal changes in the horizontal positions of the leading and trailing edges of the contact area, denoted by x LE (t) and x TE (t), respectively.
First, from the results for the FI model, we again find that the response was not accompanied by any transient processes.
Regarding the edges, they were located at x LE = −2.4 mm and x TE = 2.4 mm as the initial state. After t = 0, the leading edge maintained the position, while the trailing edge immediately moved to the left to x TE = 0.3 mm at t = 0 and maintained the position for t > 0 (which is observed as the peeling of the contact area in the outlet side). Regarding the contact pressure, it was symmetric for t < 0, showing the maximum p max ∼ 3 MPa at the contact center and the minimum p min = 0 at both edges. After t = 0, it showed the maximum p max > 20 MPa at the leading edge and the minimum p min = 0 at the trailing edge.
Then, from the results for the MI model, we again find that the response was accompanied by distinct transient processes.
Regarding the edges, they were located at the same positions as those of the previous model. However, after t = 0, not only the trailing edge but also the leading edge moved (which is observed as the simultaneous peeling of the contact area in both sides). The leading edge started to move to the right at t = 0 and then gradually approached x LE = −1.3 mm with the time constant "10 ms." Meanwhile, the trailing edge immediately moved to the left to x TE ∼ 0 at t = 0, then continued to move to the left to x TE ∼ −1 mm with the time constant "1 ms, " and then gradually approached a limiting value of x TE = 0.1 mm with the time constant "10 ms." Regarding the contact pressure, starting from the same initial state as that of the previous model, after t = 0, it showed the maximum at the leading edge and the minimum at the trailing edge. For example, in the steady sliding, p max ∼ 14 MPa at the leading edge and p min = 0 at the trailing edge.

Transient Sliding Friction
Through direct comparison of the numerical results for the two types of models (Figure 1), we have found that the boundary condition of the rigid indenter is critical for the occurrence of transient sliding. Letting the rigid indenter be vertically movable causes its upward motion (termed the "vertical lift" of the rigid indenter) when the drive velocity is applied (Figure 3), leading to a drastic spatiotemporal change in the area and pressure of the contact (Figure 4).
First, in the FI model, the horizontal positions of the leading and trailing edges are given by x TE = a 2 + l r 2 − l r (FI model) respectively, where a is the half length of the stationary contact, and l r is the retardation length, defined as a = √ 2Rδ (28) respectively. Note that x LE is constant, while x TE is a function of l r . From Equations (26) and (27), we find that there are two limiting cases: where the former means the "non-peeling" of the entire contact area, and the latter means the "complete peeling" of the contact area in the outlet side. Therefore, for example, as the drive velocity is increased, the contact area changes from "symmetric" [Equation (30)] to "asymmetric" [Equation (31)]. Note that x TE given by Equation (27) is immedeately determined when the drive velocity is applied, which is the reason why the response of the FI model is not accompanied by any transient processes.
Then, in the MI model, the horizontal positions of the leading and trailing edges are given by x TE = a 2 1 + τδ δ + l r 2 − l r (MI model) respectively, where δ = δ(t) and therefore a = a(δ) = a(t), which is the reason why the MI model creates transient processes. In Equations (32) and (33), we find two sources creating temporal changes. One is a ∼ δ 1/2 , which shrinks the contact area symmetrically by the vertical lift of the rigid indenter. From the temporal change in x LE in the right of Figure 4, we can say that under the standard condition, the time constant "10 ms" was caused by this effect. (In addition, the temporal change in z A on the right of Figure 3 supports this conclusion.) The other isδ/δ located under the square-root sign in Equation (33), working only for the trailing edge, which deforms the contact area asymmetrically. From the temporal change in x TE on the right of Figure 4, we can say that under the standard condition, the time constant "1 ms" was caused by this effect. Based on the foregoing, let us consider the spiky response of F x on the right of Figure 3. As we saw in the previous section, the temporal change in F x in the MI model is quite similar to the typical transition from the static friction to the kinetic friction, in spite of no static friction considered in the simulation. If we saw this type of response in experiments, we probably believed that this was caused by the typical adhesive friction consisting of two types of friction. In fact, the spiky response consists of the following three parts. The first is the immediate increase responding to the abrupt increase in the drive velocity at t = 0. This is obviously caused by the damping C of the Kelvin-Voigt foundation, which is essentially identical to the response observed in the FI model at t = 0. The second is the rapid decrease with the time constant "1 ms." Considering the discussion in the previous paragraph, we can conclude that the rapid decrease is caused by the rapid motion of the trailing edge to reduce the contact area. The third is the gradual decrease with the time constant "10 ms." Considering the discussion in the previous paragraph, it is natural to say that the gradual decrease in F x is caused by the gradual motion of both edges. However, also considering that the contact pressure takes the maximum p max at the leading edge and the minimum p min = 0 at the trailing edge (Figure 4), we can conclude that the gradual decrease in F x is mainly caused by the gradual motion of the leading edge. Again note that the second and third parts of F x in the MI model never appear in the FI model, which tells us that the vertical lift of the rigid indenter is essential to the spiky response.

Steady Sliding Friction
Through the numerical simulations to examine the response to the abrupt increase in the drive velocity, we have found that the response of the FI model is not accompanied by any transient processes (which means that the steady sliding friction appears immediately), while the response of the MI model is accompanied by distinct transient processes followed by the steady sliding friction. In this section, we focus on the steady sliding friction.
The situation of steady sliding is given bẏ which makes the governing equations for the MI model reduce to those for the FI model. For example, However, as seen in the numerical results for t > 100 ms in Figure 3, the normal load W (= F z in steady sliding) and friction force F (= F x in steady sliding) for the MI model are considerably different from those for the FI model (i.e., W ∼ 35 N/mm and F ∼ 5 N/mm for the FI model, while W ∼ 10 N/mm and F ∼ 1 N/mm for the MI model). Note that the difference is caused by the boundary condition of the rigid indenter: the FI model represents a "constant-gap sliding system" in which the indentation depth δ is controlled, while the MI model represents a "dead-weight sliding system" in which the normal load W (= Mg) is controlled. Therefore, when one tries to measure the sliding friction of elastomers, it is mandatory to pay much attention to the boundary condition of the counter surface: otherwise, measured values could lose their meaning. The foregoing discussion is supported by Figure 5, which shows numerical results in the steady sliding for the FI model (left column) and the MI model (right column). The dependences of δ (top row), W (second row), F (third row), and µ (bottom row) on V in the steady sliding are summarized, where µ is the friction coefficient in steady sliding, defined as The τ -values are 10 −2 s (black), 10 −1 s (blue), and 10 0 s (red), and the other parameter values are the same as those of the standard condition (Table 1). Owing to the difference of boundary conditions, the velocity dependences of µ for the FI and MI models are different from each other: under high-V conditions, the FI model shows the limiting value µ ∼ 0.2, while the MI model shows a negative dependence of µ on V with a slope of −0.5. In addition, the three curves in every graph in Figure 5 are found to be located at regular intervals, which means that the product of V and τ (i.e., l r = Vτ ) is an essential parameter for both models. Through a series of numerical simulations under various sets of parameters, master curves on µ in the steady sliding for the two types of models were obtained (see Figure 6). The red curve in the left graph is the master curve for the FI model, the ordinate and abscissa of which are µ(R/δ) 1/2 and Vτ (Rδ) −1/2 , respectively, while the red curve in the right graph is the master curve for the MI model, the ordinate and abscissa of which are µ(KR 2 /W) 1/3 and Vτ (K/RW) 1/3 , respectively. It should be stressed that every quantity assigned to the axes in Figure 6 is dimensionless.
To examine the asymptotes of the master curves, we again consider the two limiting cases: where again, the former means the "non-peeling" of the entire contact area, and the latter means the "complete peeling" of the contact area in the outlet side. First, when l r ≪ a, W and F can be estimated by respectively, where α = 2 √ 2/3. Therefore, µ for the first limiting case [Equation (38) Then, when l r ≫ a, W and F can be estimated by respectively. Therefore, µ for the second limiting case [Equation (39)] is given by or, by using Equation (43), The black broken lines in Figure 6 are the above asymptotes. Now we find that excellent agreement of the asymptotes with the master curves. Note that Equations (42) and (46) are the same as those shown by Popov (2010).  Table 1 for other parameter values).

Velocity Dependences of Friction Coefficient
In this section, based on Figure 6, let us consider the velocity dependences of the friction coefficient. In Figure 6, the drive velocity V is included only on the abscissa. Therefore, from the intersection of the two asymptotes, we introduce the critical velocity V * defined as By using V * , the formulas for estimating friction coefficient are summarized in Table 2. First, we consider the case of V ≪ V * . From the third row of Table 2, we find that the formulas for both models are the same, which means that when V ≪ V * , the vertical lift effect is negligible, although the rigid indenter could move vertically in the MI model. In addition, considering the assumption given by Equation (38), we find that when V ≪ V * , the peeling of the contact area in the outlet side is also negligible in both models. The velocity dependence of µ is µ ∼ V 1 , which means that it is caused by the damping of the Kelvin-Voigt foundation.
Then, we consider the case of V ≫V * . From the bottom row of Table 2, we find that the formulas for the two types of models are different, which means that when V ≫ V * , the vertical lift effect strongly appears in the MI model. From Equation (43), we find that the vertical lift effect appears according to which is confirmed by the numerical results shown in the upper right graph of Figure 5, where a decrease in δ means an increase in z A (= -δ): that is, the vertical lift of the rigid indenter. In addition, considering the assumption given by Equation (39), we find that when V ≫ V * , the complete peeling of the contact area in the outlet side occurs in both models. Regarding the velocity dependence of µ in the FI model, it is µ ∼ V 0 (i.e., µ is constant). This is because when V ≫ V * , the damping becomes dominant, and therefore the restoring becomes negligible, which leads to the situation that both of W and F are proportional to V, as shown in Equations (43) and (44). Regarding the velocity dependence of µ in the MI model, on the other hand, it shows the negative dependence µ ∼ V −1/2 , which is caused by the vertical lift of the rigid indenter.
As a result, we find that in the MI model, the function µ = µ(V) has a local maximum at V ∼ V * . Note that for many decades, this type of velocity dependence has been observed and discussed by a number of researchers for the sliding friction of elastomers, which seems to be basically understood as the frequency dependence of viscoelasticity (Persson, 2001;Momozono et al., 2010;Carbone and Putignano, 2013). However, the MI model produces a qualitatively similar dependence, although the frequency dependence of the Kelvin-Voigt foundation has no local maximum, where it is caused by the vertical lift of the rigid indenter.
Also note that the negative dependence µ ∼ V −1/2 means which gives us an idea of the ultra-low friction. If we try to embody the situation of Equation (50) in real systems, not only increasing V but also decreasing V * is promising, the method of which is shown by Equation (48). Probably, the most effective parameter is the damping C (or the retardation time τ = C/K).
According to the equation, increasing C leads to decreasing V * , which leads to approaching the situation of Equation (50). It should be noted that the concept can never be embodied in the FI model, because increasing C under a constant δ just linearly increases F, as shown in Equation (44). Therefore, the key is the vertical lift effect arising with the movable boundary condition of the rigid indenter.
Recently, several swollen polymers showing ultra-low friction (e.g., hydrogels and polymer brushes) have attracted attention, with inspirations from natural tribosystems in human bodies (e.g., eyes and joints) (Klein et al., 1994;Lee and Spencer, 2008). An important aspect is obviously the low-adhesive properties of surfaces caused by their "microscopic" structures. However, as another important aspect, proper "macroscopic" structures are also needed to utilize them properly under various conditions. At least, recalling the results on the transient sliding (Figure 3), we can say that in the lubricated sliding of elastomers, macroscopic structures utilizing the vertical lift effect by viscoelasticity seem to have strong advantages for smooth transition to the fluid film lubrication regime, especially at the onset of sliding.

CONCLUSIONS
In this study, the sliding friction of elastomers was investigated numerically and theoretically for the line contact between a cylindrical rigid indenter and a "frictionless" Kelvin-Voigt foundation. The onset of sliding under an abrupt increase in the drive velocity was simulated with different boundary conditions of the rigid indenter. The main conclusions are as follows: 1. When the rigid indenter is not allowed to move in any direction, just an abrupt change in the friction force appears, which is not accompanied by any transient processes. However, when the rigid indenter is able to move in the vertical direction, the transient sliding friction including three different time constants appears, resembling the typical transition from the static friction to the kinetic friction, in spite of no static friction considered in the simulation. The aforementioned drastic difference is caused by the "vertical lift" of the rigid indenter induced by the damping of the Kelvin-Voigt foundation. 2. When the drive velocity is sufficiently low, the vertical lift effect is negligible, where the restoring is more dominant than the damping, which leads to little peeling of the entire contact area and the friction coefficient proportional to the drive velocity. On the other hand, when the drive velocity is sufficiently high, the vertical lift effect becomes strong, where the damping is more dominant than the restoring, leading to the complete peeling of the contact area in the outlet side. The vertical lift of the rigid indenter strongly affects the characteristics of the steady sliding friction, which is explained well by using the critical velocity determined from the asymptotes in the master curve of friction coefficient.
At the end, it is noted again that the foregoing conclusions are obtained for the Kelvin-Voigt foundation against a cylindrical rigid indenter. In general, since behaviors are greatly affected by the rheology of the viscoelastic foundation or the shape of the rigid indenter (Popov and Heß, 2015), further investigations are needed.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/supplementary material.