Bifurcations analysis in implicit maps through the dynamics of cumulated surface errors in milling

In this contribution, we examine the evolution of surface errors during consecutive milling operations. Its description is based on a nonlinear implicit map, which is suitable to investigate the surface quality. It describes the series of consecutive Surface Location Errors (SLE) in roughing operations. As one of the principal results of the paper, bifurcations related to the fixed point of the implicit map are analyzed via normal form theorem. We determined a formula for the criticality of the bifurcation, which allows the approximate computation of the arising period-two solution. The method is demonstrated for the surface error model of milling, and the results are verified by numerical computations. Although the amplitude of the SLE would be negligible, its derivatives has a great influence in the model, which can cause stability problems.

of the material removal processes, the cutting edges enter and exit into the material periodically. This phenomenon leads to two kinds of vibration, which may limit the achievable productivity. The mitigation of these vibrations is still an open research field [1]. One of them is referred to as chatter, which is a self-excited oscillation relating to the stability loss of the cutting process [2] while the other one is related to large amplitude forced vibration near to resonant spindle speeds [3].
The most accepted phenomenon to describe chatter is the so-called surface regeneration effect or memory effect [4,5], that is the instantaneous chip thickness depends on the present and the delayed relative positions between the cutting tool and the workpiece. It can be modeled by delay-differential equations (DDE) [6], for which the stability limit draws the line between chatter-free (stable) and chatter (unstable) vibrations in the space of the technological parameters. That is referred to as stability chart or stability lobe diagram illustrating chatter-free regions [7]. These charts can be used to optimize roughing operations, where the workpiece oversize is removed with several consecutive immersions. Usually, the most productive chatter-free parameter regions, where the highest material removal can be achieved, are located near the resonant spindle speeds. However, large amplitude vibrations occur there due to the resonant excitation of the periodic milling force [8,9].
This vibration is copied to the surface and creates the so-called surface location error (SLE), which is an off- Fig. 1 Effects of the surface error on the radial depth of cut set error defined by the largest deviation between the machined and the desired surface [10][11][12][13]. This SLE directly reduces the radial depth of cut which modifies the machining process (see Fig. 1). Due to the fact that SLE is relevant at finishing operations [10], during the prediction of this surface error, the widely used methods do not consider that the SLE can influence the radial depth of cut. However, in the case study in [14], the measured vibrations amplitude resulting surface errors were in the range of the radial depth of cut. Therefore, the SLE can still have a significant impact on the surface quality in case of consecutive immersions during roughing. At every immersion, the machined surface differs from the desired one due to the SLE. The offset errors modify the radial depth of cut and generate different cutting force, which leads to a modified subsequent surface error. During this process, the SLE can accumulate and the phenomenon leads to a new surface quality parameter denoted by the cumulative surface location error (CSLE).
The general concept of the CSLE was firstly introduced in [15]. In this model, it is considered that the initial surface position is shifted relative to the desired one by the surface location error generated in the previous cut. This leads to an explicit nonlinear scalar map which determines how one SLE develops into another immersion-by-immersions. It is published first in [16], then further investigation relating to measurementbased model was carried out in [17], while experimental validation of this explicit model can be found in [18].
This explicit model is not enough to describe this complex phenomenon, and it has a significant draw-back, that it considers only the previous SLE error. However, the actual SLE influences the results similarly that leads to completely different results and a much more complex mathematical model. This improved model, where the surface error arises not only from the error in the previous cutting but also from the current one, was first introduced in [19] for turning operation. Here, static deformation error and its evolution were investigated by means of a numerical method, where the corresponding model is based on an implicit nonlinear map. This implicit model describes a similar memory effect that can be found in the general surface regeneration phenomenon [5], and it better describes the process of SLE evolution, since it takes into account not only the events of the past but also the effects of the current error.
In the existing literature for explicit maps, the linear stability analysis and nonlinear investigation [20] such as bifurcation analysis [21] and normal form calculations [22] are extensively investigated while engineering applications are diverse [23][24][25]. A comprehensive literature review about investigation of explicit maps can be found in [26,27] and all the references therein.
Implicit maps are more complex mathematically. On the one hand, they usually used during discretization in numerical simulation of differential equations, such as the implicit Euler, the implicit Runge-Kutta and the Adams-Moulton methods [28], where the main aspect of the investigations is about the performance and convergence of the numerical solver. On the other hand, implicit maps also occur in engineering applications [29,30]. Linear stability analysis of fixed point can be found in [29,31], while the readers can find more information about implicit systems in [31]. Yet, these works have not addressed normal form calculations in implicit systems and the nonlinear analysis in implicit maps is not discussed. To the best knowledge of the authors, so far the normal form analysis has not been extended to implicit maps.
This paper intends to fill this gap and tackle the challenges arising from the combination of implicit scalar systems and bifurcation analysis. Our contributions are twofold. One of the main goals of this research is to extend the implicit model of the CSLE for milling operation, while the other goal is to understand how these errors are evolving. Therefore, we derive the computational steps of a complete bifurcation analysis with linear and nonlinear investigation of the implicit map.
The rest of the paper is organized as follows. In Sect. 2, we determine the stability of the system by analyzing the fixed point of the corresponding implicit map and calculating its critical normal form coefficients. Via these coefficients, we use analytical formulae to obtain the approximate amplitude of the bifurcating solutions as a function of the bifurcation parameter. Sections 2.1 and 2.2 present the general form of the governing implicit scalar map and the corresponding linear stability analysis of the fixed point. Section 2.3 provides the bifurcation investigation, where the coefficients of the normal forms are given. Then, in Sect. 3, the computation of the SLE and the corresponding mechanical model are briefly discussed. In Sect. 4, the implicit mapping of the SLE values is introduced for consecutive roughing operations and the CSLE is defined. In Sects. 4.1 and 4.2, the linear and nonlinear stability investigation are performed parametrically for SLE mapping based on the results presented in Sect. 2. In Sect. 5, a case study for a selected milling operation is conducted providing the resulting stability charts and bifurcation diagrams. Finally, in Sect. 6, we conclude our findings and discuss future directions.

Governing equation
Consider the following governing nonlinear implicit scalar discrete map depending on a parameter in the form where x ∈ R and α ∈ R represent state variable and bifurcation parameter, respectively, and f : R × R × R → R is supposed to be sufficiently smooth (differentiable). Note that this system is implicit, since the state variable appears at the left and right hand side in (1).
Here we use the following implicit equation notion for (1) where F has the same properties as f . According to the Implicit Function Theorem, if ∂ F ∂ x i is nonsingular, then exists a smooth locally defined function where g : R → R and it satisfies that F(x i , g(x i ); α) = 0. Here the degree of smoothness of the function g is the same as that of F and generally g cannot be expressed explicitly in closed form.

Fixed point and its linear stability
Iterating the map (1), the solution may converge to a fixed point, where the state maps to itself, that is x i = x i+1 =: x * , which indeed satisfies Note that as a consequence of the Implicit Function Theorem, the fixed point also satisfies x * = g(x * ).
We are interested in to analyze the local dynamics in the vicinity of the fixed point. First, expanding (3) into Taylor series and eliminating the higher-order terms and by introducing the perturbation u i = x i − x * , the implicit variational map can be defined as Note that u ≡ 0 is the trivial solution of (5), and the bifurcation analysis of x is equivalent to the analysis of the trivial solution u ≡ 0. In (5), g (x) is derived by the total differentiation of (2) after substituting (3) and dropping the notation of the parameter dependence: It is referred to as implicit function derivation theorem in the Appendix of [29]. This form can also be found in [31] after Eq. (4.11). By using (3), one can find identity , and using (2) the partial derivatives The coefficient of the linear term in (5) can be rearranged as For convenient description, we will use a short notation for the partial derivatives with respect to the state variables in the first and second argument of f . Then the characteristic equation of (5) becomes where the characteristic multiplier is If |μ| < 1, then the fixed point is linearly asymptotically stable, otherwise it is unstable [26,27]. By investigating (8), one can show that the system loses stability via flip (period doubling) bifurcation at a critical parameter α cr where μ = −1 which leads to period-2 solutions around the fixed point and via fold (saddle-node) bifurcation at a critical parameter α cr where μ = 1. Note that only flip and fold bifurcations are possible since (1) is a scalar map. Note that for the so-called Neimark-Sacker (torus) bifurcation, one needs at least two-dimensional system in (1) [26]; however, it is out of the scope of the current work since, in our application, the SLE is scalar-valued.

Bifurcation analysis
The linear stability analysis only provides information about the local behavior around the fixed point, but does not give information about how the system behaves farther from it. Furthermore, it does not provide information on how close to the fixed point the local stability behavior is valid.
In this section, we perform the bifurcation analysis of the implicit system in (1) which can approximate the domain of attraction providing insight into the global behavior of the system. Based on the theory of nonlinear dynamical systems, the dominant dynamics can be investigated on an invariant manifold. In case of scalar system, this manifold is described equivalently by the map in (3). Setting the system parameter α close to the critical bifurcation parameter α cr , one can approximate the essential dynamics by the normal form equations on the manifold. For the normal form calculation and bifurcation analysis, the third-order approximation of the map (3) is needed [26] In this way, arising solutions and their stability can be characterized where the fixed point undergoes a bifurcation.
Since the governing equation is implicit in the state variable, the derivation of its third-order approximation is based on the total differentiation (implicit differentiation) of (2). The derivation steps are similar as presented in Sect. 2.2, but higher-order total derivatives are needed. The detailed derivation of the third-order approximation is discussed in Appendix A.

Flip bifurcation
Using the nonlinear near-identity transform in [26,32,33] i for flip-type bifurcation μ = −1, the second-order term can be eliminated yielding the normal form where the leading coefficient reads as (see the formulae after Eq. (5.69) in [26]) where the complete formula for the derivatives of g can be found in Appendix A.
In case of flip bifurcation, in the vicinity of the bifurcation point α cr the steady-state oscillations can be approximated as a function of the bifurcation parameter α [26] aŝ where γ cr and δ cr are the so-called root tendency and leading coefficient, respectively. By taking the derivative of the characteristic equation (8) at the critical point α cr , one can obtain the root tendency where f a;α : ∂α denote the mixed partial derivatives with respect to the state variables and the bifurcation parameter. This root-tendency represents the speed by which the critical characteristic multiplier crosses the −1 during the flip bifurcation leading to stability loss.
Based on the above derivation, we can define the closed-form formula for the leading coefficient in case of implicit scalar system: where the second-and third-order partial derivatives of f with respect to the state variables are defined based on the indices, e.q.: The general conditions for the existence and uniqueness of the arising oscillations are provided by means of the transversality condition γ cr = 0 and the nondegeneracy condition δ cr = 0, [26]. Here, we assume that these conditions are fulfilled.
The criticality of the flip bifurcation is determined by the sign of the leading coefficient in (15). The bifurcation is subcritical and the arising oscillations in (13) are unstable for δ cr > 0 while the bifurcation is supercritical and the oscillations are stable for δ cr < 0. Note that the arising oscillation of the state variable can be approximated as x i = x * + u i .

Fold bifurcation
In case of fold bifurcation, the critical characteristic multiplier is μ cr = 1 and the corresponding normal form reads as which can be obtained by the second-order approximation of map (3). Here, the genericity conditions for the existence and uniqueness are provided by the nondegeneracy condition g (x * ) = 0 and the transversality condition ∂g ∂α s * = 0 [26]. In the followings, the evaluation of the series of the surface location error (SLE) values in milling operations and the corresponding implicit map is investigated to demonstrate the above derivations.

Computation of surface location error
In this section, the mechanical model of the milling process and the main steps of the SLE computation are summarized based on [14]. Here, the mechanical model is kept as simple as it is possible but detailed enough to contain all relevant dynamics for the SLE computation. Therefore, the quality of the machined surface property is determined numerically in case of straight fluted tool. The workpiece is considered as a rigid body, and the milling tool is described as a flexible part (presented in Fig. 2). In case of chatter-free machining process, the forced stationary periodic motion can be calculated as the solution of the following ordinary differential where F is the periodic directional force acting on the workpiece, d is the number of degrees of freedom, and the modal parameters such as ζ k , ω n,k and U are the natural angular frequency, the corresponding relative damping ratio and the modal transformation matrix, respectively [34]. In this study, linear force model is considered [35] in which the resultant cutting force is linearly proportional to the chip width w and the chip thickness h. The chip thickness consists of two parts, stationary and dynamic chip thickness [5]. The stationary one is resulted from the projection of the feed motion in the direction of the cutting edge. The dynamic one corresponds to the surface regenerative effect of the machining process [4]. Now only the periodic motion is analyzed which creates the surface error for which the dynamical components cancel out; hence, only the stationer component of the chip thickness is used. For the SLE computation, only the forced periodic vibration of the tool-tip is needed in the direction of the machined surface normal (perpendicular to the tool path). Therefore in our model, only the dynamics along the y direction are analyzed, and only the y component of the resultant cutting force acting on the tool-tip is included. The radial and the tangential cutting force component projected to the y direction is given by where Z is the number of the cutting edges, a p is the axial depth of cut, f Z is the feed per tooth, K r and K t are the radial and tangential cutting force coefficients, respectively [14]. These parameters are usually identified by several cutting experiments for given cutting parameters and tool geometry [36]. The current angular position of the jth cutting edge is denoted by ϕ j (t) = Ωt + 2π( j − 1)/Z , where Ω is the angular spindle speed in rad/s. Note that the spindle speed in round-per-minute (rpm) is usually denoted with n in the machining society. The so-called screen function g indicates that the jth edge is in contact with the material between the enter ϕ enter and exit angle ϕ exit only These where a e is the radial depth of cut and D is the diameter of the tool. In order to avoid uncertainties resulting from the modal fitting, frequency response function H (ω) can be directly used to compute the forced stationary vibration where φ y (ω) denotes the Fourier transformed the cutting force (18) and F −1 denotes the inverse Fourier transformation. The motion of the j th cutting edge is described by means of the superposition of the rotary motion of the edge and the forced stationary vibration of the tool centre point y(t), that is Then the SL E is determined as the extremum of the cutting edge motion: The SLE depends on various cutting parameters. It is linearly proportional to the feed per tooth f Z , and in case of straight-edged-tool, it is linearly proportional to the axial depth of cut a p in case of linear force model. The spindle speed n has a strong nonlinear influence since the natural frequencies of the system could be excited by the Fourier harmonics of the cutting force, which can lead to large resonant vibrations and large SLE values at certain resonant spindle speeds [37]. The SLE is also a nonlinear function of the radial depth of cut, too, as seen in (20). This influence has a key role in the following analysis; therefore, this function is denoted by s = f SLE (a), where s = SL E/D denotes the dimensionless surface location error. This nonlinear function demonstrated in Fig. 3 near the first resonant spindle speed for up-milling (the parameters detailed in Sect. 5). Note that since the function f SLE is interpreted in the range of a ∈ [0, 1], where there is no milling for a = 0 and full immersion for a = 1. We extend the validity range of the function as follows This mathematical extension is valid in an engineering point of view, if a is larger than 1, it represents a situation when the tool is completely in the material and if a is less than 0, it represents a case, when the tool does not touch the surface and there is no machining at all.
Note that in this calculation method, first, we define the radial immersion, then we compute the resulting cutting force and the corresponding vibration, which finally leads to an SLE value [10][11][12]. In this way, we do not take into account that the SLE can affect its radial depth of cut itself. To describe this phenomenon, we introduce the implicit model, which closes this computational loop in the next section.

Implicit model of cumulative surface location error
The dimensionless SLE for the i th consecutive immersion s i can be calculated by a map as where the current dimensionless radial depth of cut a i is composed by the pre-set radial depth of cut a 0 , the actual error s i and the previously resulted one s i−1 in the form as as shown in Fig. 4a, consequently, the enter and exit angles in (20) can change in each step. Note that there is an analogy between this model and the surface regenerative effect [5], where the actual chip thickness depends on the feed, the present and the delayed positions of the tool. It can be modeled by delay differential equations, while the evaluation of the SLE is based on stationary solution. Thus, it can be described as a difference equation realized by the following implicit map where f SLE ∈ R → R with the state variable s and it is assumed to be smooth. This implicit map can be presented in the same form as the governing equation (1) by considering the following form f (s i , s i+1 ; a 0 ) = f SLE (a 0 + s i − s i+1 ). For the special pattern of derivatives, which will be needed for the linear and nonlinear investigations, see detailed derivation in Appendix B. The implicit map (26) determines how an SLE develops into another SLE over immersion-by-immersion. The series of the SLEs may converge to a fixed point, which is called cumulative surface location error (CSLE = lim i→∞ SLE i ). In other words, fixed point is a state that dynamics maps to itself s i = s i+1 ≡ s * , that is, in the special implicit map in (26), the fixed point can be determined simply by It means that every point of the function f SLE is also a fixed point, but still, the surface errors do not always

Linear stability analysis
In this section, we analyze the local dynamics in the vicinity of the fixed point. The corresponding characteristic multiplier can be calculated using (9), which simplifies to . (28) As it is shown in (28), μ can be expressed as a hyperbolic function of the derivatives of the f SLE (a 0 ) at the fixed point. By investigating (28), one can show that the system loses stability via flip bifurcation μ = −1 (period-2 solutions around the fixed) at a critical pre-set radial depth of cut a 0,cr where f SLE (a 0,cr ) = −0.5.
Note that usual fold (or saddle-node) bifurcation cannot occur in this system since it is not possible that μ reaches 1.

Nonlinear investigation
In the followings, we derive the required formulas for bifurcation analysis. The root tendencies at the critical point a 0,cr can be calculated according to (14) and using that in case of flip bifurcation (29) holds, which simplifies (14) to The leading coefficients can be calculated for the implicit CSLE model by the formula in (15) considering the special patterns in the derivatives presented in Appendix B Finally, the approximated oscillations of the period-2 solution in the vicinity of the fixed point based on (13) can be can be given by Note that the analysis based on normal form calculation is valid only in the vicinity of the bifurcation point and can catch neither the trend of the period-2 solution farther from the critical point nor the presence of fold bifurcation of period-2 solutions.

Case study: cumulative surface error in milling
In this section, we demonstrate numerical bifurcation analysis of the implicit scalar system provided by the evolution of surface errors in case of milling operations. In the conducted case study, all the computations were performed for the parameters presented in Table  1, while the measured tool-tip FRF of a long peripheral milling tool was used, which is shown in Fig. 5. Figure 6 shows the surface location errors with gray colormap in the plane of the spindle speed n and the dimensionless radial depth of cut a. Here, down-milling and up-milling are represented above each other, as a = 0 represents no machining for both cases. It is well-known that large resonant stationary vibrations can occur if one of the natural frequencies is excited by one of the Fourier harmonics of the cutting force.  Table 1 [17] As a result, considerable SLE is generated at resonant spindle speeds (e.g., n ≈ 14.2, 7.1, 3.55 krpm).
According to (27), these errors (SLE) are also the fixed points of the implicit model of the surface evolution process (CSLE= Ds * ). The stability properties can be determined with further computation steps, as it is presented in Sect. 2 to ensure that surface error remains in the vicinity of the predicted one after several consecutive immersions.
It is important to note that large surface errors (around 1 in the dimensionless error) are not realistic. In that case, this model is likely to be inadequately, because the corresponding large amplitude resonant vibration can cause the cutting edges to exit the material on a cut-by-cut basis and can influence the enter and exit angles. Furthermore, the circular approximation of the tool path is not valid anymore and some other linear approximations during the modeling lose their validity and nonlinear investigation would be needed [38,39]. In practice, significant SLE value (s 0.2) should be avoided, this the precise calculation this large amplitude vibration is out of scope of this paper.

Linear investigation
The linear stability limit of the fixed point is calculated by the formula in (29). It can be directly computed in the plane of technological parameters by computing a contour line over a dense grid or by using the socalled multi-dimensional bisection method (MDBM) [40], which is a more efficient algorithm. Note that here we use numeric derivation to calculate f SLE in (29). These curves correspond to flip bifurcation.
Boundaries of the linearly unstable parameter domains are visualized with green/red curves in Fig. 6. Selecting technological parameters from these domains, the final surface error cannot be predicted at the end of the roughing operation due to the fact that the SLE diverges from the fixed point at each immersion. Thus, the oversize of the workpiece will be uncertain before the finishing operation, therefore the corresponding technological parameters may be unfavourable during the process design. Note that chatter vibrations due to classical regenerative effect may occur along with the CSLE stability problem [17]; however, it is out of the scope of the current work.
Note that chatter and CSLE stability are not mutually exclusive; however, typically one can be found where the other is not. Unfortunately, the CSLE stability problem is more typical at favorable pockets of the stability lobes of classical regenerative vibrations (chatter) near the resonant spindle speeds.

Nonlinear investigation
The bifurcation analysis is carried out based on Sect. 4.2. The criticality of the flip bifurcation is investigated based on the sign of the normal form coefficient δ cr in (31). In case of supercritical/subcritical bifurcation (which are marked by green and red dots in Fig. 6, respectively), a branch of stable/unstable periodic-2 solution emerges from the flip points. Note that where green and red curves are met, the bifurcation is nongeneric since the normal form coefficient is zero. To get a clear view of the global dynamics, Fig. 7 present the bifurcation diagrams along horizontal cross sections at a = 0.3 for down-milling and a = 0.54 for up-milling (see dotted lines in Fig. 6). Here, the approximated amplitude of the surface error oscillations around the fixed point are plotted with black dashed curves as a function of the bifurcation parameter n (spindle speed). Note that for the fixed point, this amplitude is zero corresponding to the horizontal axis. For down-milling, as shown in Fig. 7a the fixed point is unstable for spindle speeds n ∈ [13, 14.5] krpm, while for up-milling operation (see Fig. 7b), it is unstable between n ∈ [6.3, 7] krpm. The flip bifurcations, denoted in both panels by F 1 , F 3 and F 2 , F 4 , are supercritical and subcritical, respectively.
From the engineering point of view, the subcritical case is the more dangerous, since it causes an unstable  period-2 solution around the linearly stable fixed point. It also indicates the presence of a stable period-2 solution by folding back this unstable one, which creates a bistable parameter region (see a similar phenomenon in [41]). Here, the final surface location error can converge either to a stable fixed point or to a stable period-2 solution depending on initial conditions. Consequently, this bistable parameter domain increases the unfavorable parameter range even further.

Two-step mapping and verification
The verification of the analytical results of the emerging branches could be performed by means of long numerical iterations for different bifurcation parameter values. This is computationally inefficient and leads to errors close to bifurcation points due to the very slow converging solutions (|μ| ≈ 1). However, a computationally more efficient technique can be used based on the direct computation of all the period-2 solutions along the bifurcation parameter. Two-step mapping of (26) must be considered by substituted successively into itself . (33) In case of period-2 solutions, all the even and odd values are the same s 0 = s 2 = . . . =: s * 0 and s 1 = s 3 = . . . =: s * 1 . Considering these properties in (33), the period-2 solutions s * 0,1 can be calculated by means of solving the following equations One can directly compute these solutions within a selected domain of spindle speed by means of the MDBM, which can find the solution curves for three variables (n, s * 0 , s * 1 ) for two implicit nonlinear equations. The maximum amplitude of u of the calculated period-2 solutions are plotted with continuous curves in Fig. 7. As excepted, the approximated analytically calculated bifurcation curves (black dashed curves) well follow the trend of the numerical solution near the fixed point.
The stability of the period-2 solutions provided by the two-step map can be analyzed by the implicit function derivation theorem presented in [29]. In Fig. 7, stable and unstable period-2 solutions are colored with green and red, respectively. Note that in contrast to the fixed point where only flip bifurcation occurs, here, fold type stability loss can also happen. This fold point (denoted by black stars) takes place where stable and unstable period-2 solutions merge.
Note that periodp solution can be constructed in a similar way which is briefly presented in Appendix C, for which the linear stability behavior can be determined based on [42].

Numeric iteration
In order to illustrate different kind of behavior during the evolution process of the surface error, four points were selected in Fig. 7a and successive implicit iteration of (26) based on [43,44] were performed. These iterations are not straightforward because multiple solution can exist for implicit equations (1) or (26), thus the uniqueness of the solution is not guaranteed. We bypassed this problem in the implementation of the numerical computation: we calculated the new solutions by a Newton-Raphson method starting from the previous solution which provides only a single solution. This issue could be solved with dynamic time-domain (A) (B) Fig. 8 The evolution of the dimensionless surface location errors immersion-by-immersion in case of linearly stable and unstable fixed point for parameter points A (n = 12.8 krpm) and B (n = 13.6 krpm) in Fig. 7a simulation which takes into account complex interaction between tool-workpiece engagement during forced vibration [38,39], but it would require excessive computational resources and loses the advantages of the relatively simple solution offered by the implicit model. The surface errors are plotted immersion-byimmersion in Figs. 8 and 9. A globally stable case is illustrated in Fig. 8a, where after an initial oscillation, the surface location error converges to the stable fixed point (green dash-dotted line). Panel b represents the case for a linearly unstable fixed point, where alternating surface errors emerge due to the small initial perturbation of the fixed point, which tends to the period-2 solution (green dotted lines).
The behavior in case of bistable parameter region is illustrated in Fig. 9. Here, we use the same technological parameters and start the iteration from two different initial conditions, which highlights the role of the unstable period-2 solution (red dotted lines). In our case study, if we select an initial value between the unstable period-2 solution (s * U 0 , s * U 1 ), the surface error converges to the stable fixed point (panel a), otherwise, it goes to the stable period-2 solution (s * S 0 , s * S 1 ) (panel b). In other words, the unstable period-2 solution determines a threshold that separates the attraction zone of the stable solutions.
From the viewpoint of the CSLE stability problem, during the parameter design of roughing operation, the globally stable fixed point (point A) can be acceptable because in all other cases, the oversize of the workpiece may not be precisely predicted before the finishing process.

Conclusion
In this study, we introduce an approach to analyze the bifurcations of implicit scalar maps. The presented (A) (B) Fig. 9 The evolution of the dimensionless surface location errors immersion-by-immersion in case of bistable parameter point at n = 14.52 krpm and a = 0.3 during down-milling for different initial conditions. The initial conditions for point a is s 0 = −0.2783 while for point b is s 0 = −0.2784. See the location C and D in Fig. 7a method uses the Implicit Function Theorem and its differentiation form. Analytical formulae for normal form coefficients are derived, from which approximated amplitude and stability of the arising solutions and criticality of the bifurcation point can be determined. The derived theory for implicit maps was applied in an engineering problem. The evolution of surface errors during roughing operations was investigated. The essence of the proposed method is the introduction of the novel implicit map, which takes into account not only the effects of the past but also the effects of the current state, while it describes how one error develops into another. The stability analysis of this problem shows that the final surface location error cannot be predicted at the end of the roughing process, which can affect the finishing operation. In addition, we have shown that stability loss can not only occur in resonant spindle speeds (where a large surface location error occurs), but it can also happen where the surface error itself is small, but has a large derivative along with the radial depth of cut. Also, we proved that only flip-type bifurcation could occur for the SLE model, which may lead to alternating final surface errors. With the proposed method, the traditional stability chart relating to the regenerative effect can be extended, from which chatter-free and CSLE stable technological parameters can be selected in order to support the technological design.
During the nonlinear investigations of the special milling case, we utilized further simplifications to determine the coefficients of the normal form. It has shown that bistable parameter ranges can also occur due to the presence of subcritical flip bifurcation. From the engineering point of view, it must also be avoided since these domains further increase the range of the uncertain SLE prediction. These analyses may open ways for efficient cutting parameter optimization strategies crossing the border toward reliable and optimal cutting. Our further goals are to perform measurements for the surface error problems to validate and support the theoretical results.
Although, the proposed bifurcation analysis of nonlinear implicit maps is applied only on milling example, this could be used generally in any application field. To support general usage, it is still left for future work to extend the normal form computation for systems with higher dimensions.