Complex dynamical behaviors in a Bertrand game with service factor and differentiated products

In this paper, taking the factor of product service provided by the manufacturers into consideration, a static Bertrand duopoly game with service factor is studied first, in which these two oligarchs produce differentiated products. A dynamic Bertrand duopoly game with bounded rationality is established by using the gradient mechanism. Keeping the adjustment speeds in a relatively small range may help the long-term stable operation of the market. It is found that there is another 2-cycle, different from the one flip bifurcated from the fixed point, which may appear through a saddle-node bifurcation. The unstable set of the saddle cycle, connecting the saddle to the node, gives a closed invariant curve. In addition, the emergence of intermittent chaos implies that the established economic system has the capability of self-regulating, where PM-I intermittency and crisis-induced intermittency have been studied. With the help of the critical curves, the qualitative changes on the basin of attraction are investigated.


Introduction
When the consumers enter the shopping malls to buy goods, they will feel an increasingly common phenomenon, that is, the difference in function and utility of same kind of products on the shelf is getting smaller and smaller. The products themselves tend to be homogeneous. One of the main reasons for the gradual homogenization between the products is the rapid improvement in enterprise production technology. A traditional view on the degree of product differentiation, proposed by Shy [1], is that the reduction in product differentiation may lead to the increase of market competition intensity and then lead to price war among enterprises. Although the price competition is the key method for oligopoly market competition to enlarge market shares, no one can deny that conducting a long-term price war will lead to vicious price competition among enterprises, which will further damage the profits of enterprises [2]. This low-level price competition is not conductive to stable operation in the long run. Caplin and Nalebuff [3] investigated a price competition model with differentiated products and proposed a new method to analyze the imperfect competition mode. Liu et al. [4] studied a mixed duopoly game with endogenous product differentiation and found that the difference between products arises only when the production cost is very low. The differentiation between products can attract the consumers who love a certain feature of the product, and make the consumers distinguish it from the enterprises selling the same kind of products, so that the enterprises producing differentiated products can gain competitive advantages in market competition.
In many literature on oligopoly market competition, a large number of scholars focus on the effect of product differentiation on the stability and dynamics of the dynamic system. Zhang et al. [5] exploited the form of linear inverse demand function [6] including the product differentiation parameter, and then established a dynamic duopoly model with differentiated products and bounded rationality. They also investigated the dynamic characteristics of this dynamic model. Other scholars, as for example, Zhang et al. [5], Ahmed et al. [7] and Anufriev et al. [8], further studied the price competition game with differentiated products and also analyzed the dynamic characteristics of the built system. Zhang et al. [9] studied complex dynamics in a quantum Bertrand duopoly with differentiated products. It was found that the heterogeneity of products does affect the stability of the model. In addition, Fanti and Gori [10] investigated a Cournot game with two heterogeneous oligarchs and analyzed the influence of changes in consumer demand for differentiated products on the stability of the model. They observed that the increase in the degree of differentiation between products can reduce the intensity of competition. Same conclusion could be found in the research of Agliari et al. [11], that is, the higher the differentiation level between products, the smaller the intensity of market competition.
In the process of consumption, consumers pay more and more attention to some non-price factors of products [12], especially the factor like product service. An increasing number of enterprises and manufacturers begin to attach importance to the aspect of product service. For instance, the clothing manufacturers will present small items, such as buttons, when selling clothes, online shopping on the e-commerce platform can be delivered to your home directly (amazon.com can deliver free of charge when the consumption limit is reached), and electronic product manufacturers, such as Apple and Huawei, set up product experience stores. Both enterprises and manufacturers expect to win the favor of consumers by providing product services, and then they can obtain more market shares. A lot of researchers have made many related researches on this topic. Giri and Sarker [12] considered a situation that two independent retailers compete in both price level and service level and analyzed the equilibrium behaviors of a Stackelberg model. Shin [13] pointed that the retailer who provide product service has the advantage of after-sale service. Wu [14] built a competitive game model for sellers to provide information services for horizontal differentiated products. The study found that the motivation to provide information services not only affected the profits of sellers, but also affected social welfare.
However, a feature of providing product service that cannot be ignored is service intangibility. After providing product services, it is inevitable that consumers will experience product services from one manufacturer and then purchase products from other manufacturers. This is a notable phenomenon of ''free riding'' between manufacturers proposed by Telser [15] as early in 1960s. In the reference [13], the scholar Shin put forward that the phenomenon of ''free riding'' can make free riders reduce the motivation of price competition with service providers. Bernstein [16] also believed that ''free riding'' among manufacturers is widespread and common, but it does not necessarily have a negative impact on competition. Furthermore, the phenomenon of ''free riding'' could lead service spillover effect between different manufacturers. The scholars Tu et al. [17] incorporated service factors into a dynamic hybrid supply chain model and studied the impacts of the price adjustment speed and the service level on the model. In the references [18][19][20][21], the factors of service level or service spillover effect were considered in the establishment of supply chain model. Unfortunately, only few scholars investigate the dynamic competition with service factor by using the knowledge of nonlinear dynamics.
In recent years, a lot of researchers have studied the complexity of dynamic competition between the enterprises and have made great achievements. Fanti et al. [22][23][24] established the nonlinear duopoly game with producing differentiated products and carrying on management entrustment agent. It could be found that the topological structures of basin of attraction have great changes after some global bifurcations, leading to fractal structure of the basins, also in the presence of coexisting attractors. Bischi et al. [25][26][27][28] found the phenomenon of coexistence in the phase space, which may imply that the dynamic evolutions have the path dependence (that is, the historical dependence). Bischi et al. also found that the global bifurcation occurring in the system may be due to contacts of the critical curves with the borders of basin of attraction, which could cause the dramatic changes in size, structure and shape of the basin of attraction. Zhang et al. [29] studied the complex dynamic behaviors appearing in a semicollusion Cournot duopoly game model, in which the enterprises cooperate in production but compete in process of research and development. Zhou et al. [30] observed intermittency phenomenon in the built twostage dynamic game, in which both two enterprises carried on research and development activities. Zhou et al. [31] not only found the phenomenon of coexisting attractors, but also used critical curves to analyze the reason for the structural changes of the basin of attraction. More related investigations can be found in the references [32][33][34].
The rest parts of this research are as follows: a static Bertrand duopoly game is studied in Sect. 2. In Sect. 3, a dynamic Bertrand duopoly game with bounded rationality and the gradient mechanism is established. The analyses on local stability of four equilibrium points are given. The basin of attraction and the critical curves are characterized in Sect. 4. Some local bifurcations appeared in the system are studied in Sect. 5. In Sect. 6, the intermittency phenomena are discussed. The global bifurcations, which cause qualitative change in the topological structure of basin of attraction, are analyzed by the critical curves in Sect. 7. At last, some valuable conclusions are primarily given in Sect. 8.

Static Bertrand duopoly game with service factor
With the continuous improvement in production technology, the products produced by different manufacturers tend to be homogeneity-oriented. That is to say, the differences in quality, function and utility between products are very small. On the basis that there is no restriction on the production capacity of manufactures (namely, the manufactures do not compete in outputs), if one manufacturer attempts to cut the price of its own products so as to increase sales, then other manufacturers will follow the price reduction measures one after another, which will lead vicious price competition among manufacturers [1]. The final result of vicious price war can only be that the selling price of products will drop again and again until it is equal to the marginal cost of production and then the manufacturers cannot make a profit anymore. On the contrary, the case that one manufacturer occupies a large market share due to its price reduction measures will not arise, if there are differences between the products produced by different manufacturers. This is because there will always be loyal consumers who have a preference for a certain feature of the product and will not be shaken by price changes [10]. Hence, the manufacturers should consider increasing product differentiation to ease the dilemma caused by vicious price competition [11]. Furthermore, more and more consumers pay attention on nonprice factors during the process of consumption, such as product service. The product service is an effective way to avoid price wars, compete for consumer preferences and gain larger market share, and then enhance the competitive advantage of manufacturers in the market. This service measure also can turn manufacturers from a single-price competition mode to a multi-dimensional competitive mode [12]. Based on the above analyses, it is assumed that there are two oligarchs existing in an industry, labeled by M i (i ¼ 1; 2). These two manufacturers produce differentiated products and provide product services. The selling price p i and the product outputs q i of the manufacturer M i satisfy the following form: where a [ 0 denotes the largest scale of the market, b [ 0 represents price sensitivity coefficient, and c represents the product differentiation coefficient. Moreover, the parameter c belongs to ½0; 1Þ and 0\c\b. The smaller the value of c, the greater the degree of differentiation between products, and c ¼ 0 denotes that the products produced by different manufacturers M i are totally heterogeneous. Since providing product services have the characteristic of service intangibility, it is impossible for manufacturers to avoid the problem that consumers experience product service in their own retail stores but buy products from their competitors. That is to say, there will be spillover effect of service, which will lead to the ''free riding'' phenomenon [15]. The service level that the consumers experience from the manufacturer M i is marked ass i in Eq. (1). It satisfies the equation where s i [ 0 and s j [ 0 represent the service level provided by two manufacturers, respectively. The parameter h 2 ð0; 1 denotes the intensity of service spillover effect. In addition, the cost functions of these two oligarchic manufacturers M i have linear forms, which can be expressed as follows, in which c i [ 0 is the marginal production cost that producing unit product for manufacturer M i . The profit functions with service factor of two manufacturers can be expressed as following, where S i ¼ 1 2 hs 2 i represents the service cost that the manufacturer M i needs to spend when providing a service level as s i , and h [ 0 is the service cost coefficient. Based on above assumptions, a static duopolistic game can be established, in which the manufacturers compete in the price and both of them pursue profit maximization. Under this duopolistic setting, the following proposition is tenable.
Both two manufacturers want to maximize their respective profit, which means that the first-order conditions should be satisfied. In addition, the corresponding marginal profit function of manufacturer M i can be given by To further facilitate calculation, two auxiliary variables can be introduced as A ¼ a þ s 1 þ hs 2 þ bc 1 and B ¼ a þ s 2 þ hs 1 þ bc 2 . Then, the Nash equilibrium can be obtained by calculating the algebraic system as following, From the first equation in Eq. (5), it can be easily got that p 1 ¼ Aþcp 2 2b . By substituting this relation into the second equation of Eq. (5), we can obtain that p Ã 2 ¼ 2bBþcA 4b 2 Àc 2 . Finally, the value of p Ã 1 can be gained by retrieving the value of p Ã 2 , namely p Ã 1 ¼ 2bAþcB 4b 2 Àc 2 . Furthermore, since the second-order conditions are satisfied, that is, 4b 2 Àc 2 is really a point of maximum, and it is also the unique Nash equilibrium point of this static Bertrand game. This proves the following, The unique Nash equilibrium of the static Bertrand game exists and it can be expressed as 3 Dynamic game with the gradient mechanism In real market, gathering market information is costly, so it is difficult for enterprises to grasp all the market information and make competitive decisions completely rationally [35]. In order to accord with actual market situation, we further suppose that these two manufacturers considered in this economic system are boundedly rational and will adopt the gradient mechanism [35][36][37] to adjust their pricing decisions. That is to say, each manufacturer will adjust the selling price of products in period t þ 1 by estimating its marginal profit in period t. Moreover, both these two manufacturers make long-term product pricing decisions through repeated competitive games. Therefore, the specific form of the dynamical gradient mechanism with time evolution can be expressed as where k i [ 0 denotes the adjustment speed of pricing decision of manufacturer M i , which also implies that the difference between the price in period t þ 1 and the price in period t is directly proportional to its marginal profit in period t. Finally, the dynamic adjustment model of two boundedly rational oligarchs who produce differentiated products and provide services could be expressed as T : in which T is a time evolution operator.

Local stability analysis
The fixed points of system T can be defined as the solutions of the following algebraic system, That is to say, these two manufacturers may think the prices of their opponents in the next period as equal to the selling prices in current period, p 1 ðt þ 1Þ ¼ p 1 ðtÞ and p 2 ðt þ 1Þ ¼ p 2 ðtÞ. After some simple calculation steps, it is not hard to find out from Eq. (8) that system T includes four different equilibrium points, Due to the nonnegativity of parameters and the condition that 0\c\b, it is clear that these four equilibrium points are located in the positive quadrant. In addition, the points E 1 , E 2 and E 3 are boundary equilibrium points of system T.
The local stability of these four equilibria can be determined through analyzing the eigenvalues of the Jacobian matrix. The corresponding Jacobian matrix of system T is given by, where From the nonnegativity of the above parameters, the following propositions hold.
Proposition 3 When the adjustment speed k 2 of manufacturer M 2 satisfies the condition that Analogously, some similar conclusions with respect to the boundary equilibrium point E 3 can be given in the following proposition.
Proposition 4 When the adjustment speed k 1 of manufacturer M 1 satisfies the condition that The corresponding Jacobian matrix at the Nash 4b 2 Àc 2 can be written as The characteristic equation of J E Ã N À Á is given by in which TrðE Ã N Þ and DetðE Ã N Þ are the trace and the determinant of the matrix J E Ã N À Á , respectively. The specific expressions of TrðE Ã N Þ and DetðE Ã N Þ are Proposition 5 The Nash equilibrium E Ã N will never lose its stability through a Neimark-Sacker bifurcation.
Proof The discriminant of Eq. (11) can be given by Since the nonnegativity of the parameters and 4b 2 À c 2 [ 0, it is quite clear that the discriminant D is always greater than zero. This means that it is impossible for the characteristic equation PðkÞ to have complex roots. Namely, the Nash equilibrium E Ã N will never lose its stability through a Neimark-Sacker bifurcation. h When k 1 and k 2 satisfy the additional conditions that Þ , a sequence of flip bifurcations may occur.
Proof Since it is difficult to obtain the analytic expressions of the eigenvalues corresponding to J E Ã N À Á , we cannot analyze the asymptotical stability of the Nash equilibrium through judging the positions of eigenvalues on complex plane. With the help of Jury criterion [35] and references therein, the sufficient conditions for local stability of the Nash equilibrium can be given. At the Nash equilibrium 4b 2 Àc 2 , the specific forms corresponding to the Jury criterion can be written as the following three conditions, From the nonnegativity of the parameters, we have that the condition (i) always holds. Namely, there is no bifurcation related to an eigenvalue equal to ? 1 involving the Nash equilibrium. As for the conditions (ii) and (iii), they will be satisfied when the following conditions hold, : That is to say, if the above two conditions are satisfied, the Nash equilibrium E Ã N is asymptotically stable.
In addition, if 1 þ Tr þ Det ¼ 0, the mathematical expression of the flip bifurcation curve can be given by which means that if the values of k 1 and k 2 cross this bifurcation curve, namely, sequence of flip bifurcations may occur. From Proposition 5, it is clear that the eigenvalues of J E Ã N À Á are both real. It is also possible that two eigenvalues cross -1. When k 1 and k 2 satisfy the additional conditions that From the perspective of qualitative analysis, the local stability of four equilibrium points has been studied above and the corresponding stability conditions for these equilibria also have been given. Next, from the perspective of numerical simulation, we will further investigate the dynamic characteristics of this economic system when the parameters take different values in their possible ranges.

Basin of attraction and critical curves
Suppose that A is an attractor of the map T, where A may be a Nash equilibrium, a cycle or other complex attractor. The basin of attraction of A contains all points that generate bounded trajectories converging to A, which can be given as  [38]). We denote the basin of infinity as Bð1Þ, which is the set of points having unbounded trajectories. Its basin boundary is denoted as oB 1 , and we have oB 1 ¼ oBðF Þ, where oBðF Þ is the boundary of the closure of F . In addition, the two coordinate axes are invariant because the points on them will be mapped into the same axes by map T. The boundary of the set of bounded trajectories may be composed of segments belonging to the invariant coordinate axes and their higher rank preimages.
Apparently, the map T is a noninvertible map, since taking a point ðp 1 ; p 2 Þ from R 2 þ arbitrarily, we can get a unique point ðp 0 1 ; p 0 2 Þ through the map T but this point may have more than one preimage. Furthermore, from calculating a fourth degree algebraic system (7) with respect to p 1 and p 2 , we can get four preimage points of the origin E 1 ¼ ð0; 0Þ. That is, the rank-1 preimages of the origin can be obtained by the following equations, It is clear that the origin point E 1 itself is a rank-1 preimage, and there are two other rank-1 preimages located on two invariant coordinate axes. One is 2 and another is E And one more preimage exists in the positive quadrant.
The line segments E 1 E ð2Þ À1 and E 1 E ð1Þ À1 have been remarked as w 1 and w 2 for simplicity. The corresponding rank-1 preimages of these two line segments are denoted as w À1 1 and w À1 2 , respectively. Moreover, it is easy to obtain w À1 1 and w À1 2 . Taking any point p; 0 ð Þ on w 1 arbitrarily, then its rank-1 preimages need to satisfy the following equations, From the second formula of Eq. (15), we can know that the rank-1 preimages of the invariant line w 1 are located on itself or on the straight line 1 þ k 2 ½B À 2bp 2 þ cp 1 ¼ 0. Analogously, the rank-1 preimages of the invariant line w 2 may be located on the straight line 1 þ k 1 ½A À 2bp 1 þ cp 2 ¼ 0 or on w 2 itself.
Furthermore, the fourth rank-1 preimage of the origin point E 1 is at the intersection of two equations that satisfying that is, The basin boundary oB 1 is included in that set of preimages of two invariant line segments w 1 and w 2 , that is, the basin boundary oB 1 is included in that set of preimages, An example is shown in Fig. 1, in that case the basin of attraction of the Nash equilibrium E Ã N is a regular convex polygon with four vertices E 1 , E  Fig. 1, the feasible set F E Ã N À Á of E Ã N is depicted in white, defined as the set of points that can generate bounded trajectories. The basin of infinity Bð1Þ is depicted in gray, and the points in this region will generate unbounded trajectories.
In this situation, the boundary of this basin of attraction can be given as oB 1 2 . Under this set of parameters, the boundary equilibrium E 2 is a saddle point with locally stable manifold along p 2 , because its corresponding eigenvalues are k 1 j j ¼ 2:1237 [ 1 and k 2 j j ¼ 0:1817\1. Similarly, another boundary equilibrium E 3 is also a saddle point with locally stable manifold along p 1 , because its two eigenvalues are k 1 j j ¼ 0:1817\1 and k 2 j j ¼ 2:1237 [ 1, respectively.
From both perspectives of mathematical and economic, it is worth studying and analyzing the structure of BðAÞ, because the long-run evolution of the economic system (7) may have path dependence. A comprehensive knowledge of the properties of feasible region is of vital importance to predict a system's feasibility. In addition, with the help of the critical curves, the global dynamics and bifurcations of the noninvertible map can be better studied.
The critical curves of a two-dimensional noninvertible map is denoted as LC, and the preimage of LC (two of which must be merging, by definition, see Bischi et al. [26] and Mira et al. [41]) can be denoted as LC À1 . Furthermore, the preimages LC À1 of the critical curves LC must belong to the locus of points where the determinant of the Jacobian matrix (9) vanishes. After some simple calculations, we can have that the preimages LC À1 are formed by two branches, marked as LC ðaÞ À1 and LC ðbÞ À1 in Fig. 2a. The specific expressions of LC ðaÞ À1 and LC ðbÞ À1 can be given as following, in which m and n are auxiliary variables, m ¼ k 1 c þ k 1 k 2 cB À 4k 2 b À 4k 1 k 2 Ab and n ¼ k 2 c þ k 1 k 2 cA À 4k 1 b À4k 1 k 2 Bb. The corresponding critical curves LC are also formed by two branches, denoted as LC ðaÞ and LC ðbÞ in Fig. 2b, where TðLC ðaÞ À1 Þ ¼ LC ðaÞ and TðLC ðbÞ À1 Þ ¼ LC ðbÞ . Moreover, as shown in Fig. 2b, two branches of the critical curves, LC ðaÞ and LC ðbÞ , separate the two-dimensional phase plane into three areas, labeled by Z 0 , Z 2 and Z 4 , Z 0 is the set of points having no preimage, the points belonging to Z 2 have two different preimages, while the points in Z 4 have four distinct preimages.

Local bifurcation
In a nonlinear dynamical system, a slight change of a system parameter may lead to a great change in the dynamic behavior of the system. Next, the system parameters, such as the adjustment speed k 1 , the service level s 1 , the service spillover effect h and the price sensitivity coefficient b, are, respectively, selected as bifurcation parameters. As shown in Fig. 3, the parameters are fixed as a ¼ 4:6198, c ¼ 0:3668, c 1 ¼ 1:0524 and c 2 ¼ 1:0079. It is of interest to study the effects of these bifurcation parameters on the stability of system T, especially the influences on manufacturer M 2 .
The bifurcation diagrams with respect to p 2 are depicted in blue and the corresponding largest Lyapunov exponents diagrams are remarked as black. The two 1-D bifurcation diagrams with respect to adjustment speed k 1 are shown together in Fig. 3a. The other parameters are chosen as b ¼ 2:6707, h ¼ 0:3639, s 1 ¼ s 2 ¼ 4:0000, and it is worth noticing that the values of k 2 are chosen as k 2 ¼ k 1 in upper part of Fig. 3a and k 2 ¼ 0:1496 in lower part of Fig. 3a, respectively. Furthermore, the upper part in Fig. 3a, which shows that a 2-cycle that flip bifurcated from the fixed point may undergo a Neimark-Sacker LC ðaÞ À1 : and LC ðbÞ À1 : bifurcation, as the value of k 1 is gradually increasing. It is quite remarkable that the values of black curves marked by pink box fluctuate up and down near zero, which means that system T is at a quasi-periodic state. The lower figure of Fig. 3a shows that a 2-cycle will experience a flip bifurcation with continuous increase of the value of k 1 . Finally, system T will drop into chaos when k 1 % 0:1788. It can be seen that the economic system (7) may keep in a steady state when the value of k 1 is relatively small. Once the value of k 1 crosses the stable threshold, the economic system may uncontrollably drop into a disorderly state, and then it will be quite difficult to predict the trends of this economic system. The two different bifurcation types with respect to s 1 are shown together in Fig. 3b. As shown in the upper figure of Fig. 3b, the fixed point will bifurcate to a 2-cycle, and then this 2-cycle will experience a Neimark-Sacker bifurcation with constant increase of s 1 . The corresponding largest Lyapunov exponents, marked by a pink box, fluctuate up and down near zero, which means that system T is at a quasi-periodic state. The bifurcation diagram corresponding to another bifurcation type of system T is displayed in the lower figure of Fig. 3b, and system T is at a stable period-1 state when s 1 2 ½0; 2:0280Þ. Similarly, with increasing the value of s 1 , system T will experience a sequence of flip bifurcations and eventually enter into chaotic motion when s 1 ¼ 6:8160.
The 1-D bifurcation diagrams with regard to h and b are given in Fig. 3c, d, respectively. We can clearly see from Fig. 3c that system T will keep in a stable state only when the value of the spillover effect h is relatively small. Once the value of h exceeds 0.6970, then system T will become unstable. Furthermore, it is quite worth to notice that the bifurcation shown in Fig. 3d is different from the others in Fig. 3. Hence, the effect of the parameter b on the stability of system T needs to be further investigated and analyzed. When b 2 1:46; 2:7 ð Þ , it is visible that the fixed point is attracting, and then the 2-cycle appears both increasing and decreasing the parameter b. In With the same parameters as in Fig. 3d, but we choose the value of h as 0.3639. The corresponding bifurcation diagram can be seen in Fig. 4a, which only shows a 2-cycle when b 2 0:9775; 4:862 ð Þ . It is necessary to investigate the 2-cycle appeared in the simulations and there may be other 2-cycles existing in the basin of attraction, which are different from the one flip bifurcated from the fixed point. In Fig. 4b, there is only an attracting 2-cycle when b ¼ 4:0000, whose feasible region is depicted in white and the basin of infinity is denoted as gray. As shown in Fig. 4c, a pair of attracting 2-cycles will coexist in the basin of attraction when b % 4:226055, and the second one may appear through a saddle-node bifurcation. In addition, under this set of parameters, the Nash equilibrium is marked as a green point. It is clear that the second attracting 2-cycle is located on the boundary of its basin of attraction, which is depicted in red. With gradually increasing b to 4.5, the basin of attraction corresponding to the 2-cycle appeared via a saddle-node bifurcation becomes bigger, see in Fig. 4d. When b ¼ 5, the first 2-cycle will bifurcate into two closed invariant curves via a Neimark-Sacker bifurcation, as shown in Fig. 4e. As the value of b further increases to 5.1656, the size of these two closed invariant curves will become bigger, see in Fig. 4f.
When b ¼ 5:1657, it is worth noting that there exists an attracting 4-cycle a 1 ; a 2 ; a 3 ; a 4 f g (which is also four fixed points for the map T 4 , the fourth iteration of the map T) coexisting with an attracting 2cycle, as shown in Fig. 4g. In fact, the two closed attracting curves (due to the Neimark-Sacker bifurcation of the two-cycle focus b 1 ; b 2 f g) shown in Fig. 4f, as the parameter b increases, undergo a phase locking, leading to a pair of 4-cycles, an attracting node (with periodic points a 1 ; a 2 ; a 3 ; a 4 f g ) and a saddle 4-cycle (with periodic points c 1 ; c 2 ; c 3 ; c 4 f g ) and the two closed invariant curves still exist as saddle node connections (for map T 2 ) of the pair of 2-cycle leading to two attracting closed curves. This is confirmed by the eigenvalues that we have computed at the pair of 2cycle visible in Fig. 4g, which are, one very close to ? 1, but smaller than 1 and the second one positive and smaller than 1. Figure 4g and the enlargement in Fig. 4 show the basins for the map T 4 . For the second iterate T 2 the unstable set of the saddle 2-cycle forms the attracting closed curve on which the attracting set is the 2-cycle a 1 ; a 2 ; a 3 ; a 4 f g , while the internal branch of the stable set of the saddle 2-cycle is issuing from the repelling fixed point (for map T 2 ) b 1 ; b 2 f g. The eigenvalues of Jacobian matrix for the pair of 2-cycles of map T 2 , a 1 ; a 2 ; a 3 ; a 4 f g , are also calculated, leading to j¼ 0:21083962 7868446\1. Moreover, when b % 5:16566457, the values of k 1 a 1 ð Þk 1 a 3 ð Þ j j and k 1 a 2 ð Þk 1 a 4 ð Þ j j will cross ? 1, that is, a saddle-node bifurcation occurs. Figure 4h is a local enlargement of Fig. 4g, showing the structure of the attracting closed invariant curve, saddle-node connections, as explained above.
The system T may have bistability, which can be observed from the 1-D bifurcation diagrams with different initial conditions in When k 1 ¼ 0.3278665 (case shown in Fig. 5c), there exists an attracting 4-cycle a 1 ; a 2 ; a 3 ; a 4 f g which, as in the case shown in Fig. 4, coupled with a saddle-4 cycle c 1 ; c 2 ; c 3 ; c 4 f gforms a saddle-node connection for map T 2 , leading to two attracting closed curves. The unstable set of the saddle goes to the attracting node forming a closed curve around the repelling focus 2-cycle (which had a Neimark-Sacker bifurcation), while the internal branch of the 4-cycle saddle is issuing from the 2-cycle repelling focus. Similarly in Fig. 5d and its enlargement in Fig. 5e, where, calculating the eigenvalues of the Jacobian matrix in this 4-cycle a 1 ; a 2 ; a 3 ; a 4 f g , we obtain and k 2 ða 1 Þk 2 ða 3 Þ j j¼ k 2 ða 2 Þk 2 ða 4 Þ j j¼ 0:999999747 372979\1. It is clear that a 1 ; a 3 f g and a 2 ; a 4 f g are two 2-cycles for the second iteration of the map T and they are still attracting, but close to a saddle-node bifurcation with the saddle 4-cycle, leading to a reverse phase locking, as shown in the enlarged part in Fig. 5e. The basins are shown for map T 4 .
When we slightly increase the value of k 1 from Fig. 5d, due to reverse phase locking, the two 4-cycles merge in a saddle-node bifurcation and disappear, leading to a pair of attracting closed invariant curves, as shown in Fig. 5f, on which we observe quasiperiodic trajectories.
When k 1 ¼ 0:3400, there have a 2-cycle that coexists with two closed invariant curves in the set of bounded trajectories, see in Fig. 5f. The 2-cycle shown in Fig. 5f is not the one flip bifurcated from the fixed point, it appears, as in the case shown in Fig. 4, by a saddle-node bifurcation at a value of k 1 2 0:3355; 0:3356 ð Þ . This 2-cycle coexists with another attracting set associated with the sequence of bifurcations of the fixed point, that is, two closed invariant curves which have experienced a Neimark-Sacker bifurcation. With gradual increase of k 1 , the 2cycle will disappear via another reverse saddle-node bifurcation at a value of k 1 2 0:3509; 0:3510 ð Þ , as shown in Fig. 5g, h. In Fig. 5g, it is clear that each point of the 2-cycle is located on the boundary of its absorbing area, when k 1 ¼ 0:3509. In Fig. 5h, when k 1 ¼ 0:3510, there is a unique attracting set, a twopiece chaotic attractor for the map T 2 , exists in the basin of attraction, which is preexisting far from the 2-cycle (see Fig. 5g).

Intermittency
Another notable dynamic phenomenon taking place in system T is the appearance of intermittency. The main characteristic of intermittent chaos in the time series plot is that the amplitude of state variable changes alternately between regular oscillation and irregular oscillation. Corresponding to the real economic system, if the phenomenon of intermittent chaos occurs in the system, it means that the system has self-regulation ability. That is to say, without the help of external force, the economic system can adjust and reverse the market state from disorderly competition to orderly competition by itself. In addition, back in the early 1980s, two French scholars, Pomeau and Manneville, had subdivided the phenomena of intermittent chaos into three types, see in the reference [39]. Nowadays, the phenomena of intermittent chaos include five types in general, which are specifically divided into PM-type (PM-I intermittency, PM-II intermittency, PM-III intermittency), On-off and In-out intermittency, and crisis-induced intermittency [40].
The parameters are keeping unchanged as in Fig. 5, we further fix the value of k 1 as 0.3747 and choose the service level s 1 provided by M 1 as bifurcation parameter. The corresponding 1-D bifurcation diagram with respect to s 1 is shown in Fig. 6a, in which the bifurcation diagram about p 1 and the bifurcation diagram about p 2 are depicted in red and blue, respectively. When s 1 ¼ 1:0060, there is an attracting 4-cycle existing in the set of bounded trajectories, which can be seen in Fig. 6b, where the basin is given in 4 colors for the map T 4 . The fractal structure of those basins shows the existence of a repelling chaotic set inside the set of bounded trajectories that will appear decreasing the parameter (as can be seen in Fig. 6a, when the attraction 4-cycle and a saddle one merge in a saddle-node bifurcation and disappear. This occurs at s 1 ¼ 1:0050, leading to a one-piece chaotic attractor, as shown in Fig. 6c, where the place of the pre-existing 4-cycle is marked by red points, and it is the place where a trajectory spends a lot of time (i.e., iterations). Therefore, it is a PM-I intermittency, due to ''the ghost'' of the 4-cycle. The corresponding time series plot is shown in Fig. 6d. Figure 6e shows a local enlargement of Fig. 6a in the interval s 1 2 0:75; 0:8 ð Þ , and only the bifurcation diagram about p 2 is shown in Fig. 6e, which is depicted in blue. From Fig. 6e, it can be clearly seen that there is a ''periodic window'' of period six, appearing inside a 2-cyclic chaotic attracting set. Increasing the parameter from the two-piece chaotic set an attracting 6-cycle appears via saddle-node bifurcation, and the attracting 6-cycle leads to a 6-piece chaotic attractor, shown in Fig. 6f when it is already close to a contact bifurcation of the chaotic set with its immediate basin of attraction, leading to an expansion bifurcation which reveals again the 2-cyclic chaotic attractor, as shown in Fig. 6g. This transition is due to the periodic window occurring inside an absorbing region where a 6-piece chaotic attractor becomes a chaotic repeller due to a contact bifurcation with its immediate basin boundary. Therefore, this is a crisis-induced intermittency and the corresponding time series plot is shown in Fig. 6h. b Fig. 4 The parameters are same as in Fig. 3d, but the value of h is chosen as 0.3639. a Another 1-D bifurcation diagram with respect to b. It is clear to show a 2-cycle when b 2 0:9775; 4:862 ð Þ . b There is only an attracting 2-cycle when b ¼ 4:0000, whose feasible region is depicted in white and the basin of infinity is depicted in gray. c A pair of attracting 2cycles will coexist in the set of bounded trajectories, and the second attracting 2-cycle will appear from a saddle-node bifurcation when b % 4:226055. In addition, the basin of attraction of the second 2-cycle is depicted in red and the Nash equilibrium is denoted as green. d It is clear that the basin of attraction of the second 2-cycle will become bigger when b ¼ 4:5. e The first 2-cycle will bifurcate into two closed invariant curves via a Neimark-Sacker bifurcation when b ¼ 5. f The size of the two closed invariant curves will become bigger when b ¼ 5:1656. g An attracting 4-cycle a 1 ; a 2 ; a 3 ; a 4 f g and a saddle 4-cycle (here too close to the attracting one to be distinguished) give a saddle-node connection for map T 2 on the two attracting closed curves at b In addition, it is clear that the frontier of these two pieces of chaotic attractor almost contacts with the boundary of their immediate basin of attraction for map T 2 . As the value of s 1 increases to 1.4040, we can find from Fig. 7c, d that a contact bifurcation has occurred. After the contact between the two pieces chaotic attractor and the boundary of its immediate basin of attraction, the two pieces chaotic attractor will expand in a unique chaotic set. It is also shown in Fig. 7c, d that most time is spent in the two previous intervals in the unique chaotic set, due to the ghost of the 2-piece chaotic set.

Global bifurcation
In the previous sections, the basins of attraction are used to comment the local bifurcation appeared in the system. Slight changes of the values of parameters will not only change the number and structure of the attracting sets of system T, but also cause the changes in topological structure of the basin of attraction, and even the ''gaps'' (i.e., changes in the boundary of the set of divergent trajectories) will appear in the basin boundary. Then, the analyses with respect to the global bifurcations that may lead qualitative changes in the topological structure of the basin of attraction will be carried out with the help of critical curves in this section.
Under the case that the parameters are chosen as a ¼ 1:1673, b ¼ 0:7101, c ¼ 0:5302, h ¼ 0:9303, k 2 ¼ 0:4139, s 1 ¼ 0:7033, s 2 ¼ 0:6205, c 1 ¼ 0:1121 and c 2 ¼ 1:9901, we will be interest in investigating the ''gaps'' appearing on the basin boundary. As shown in Fig. 8a, when k 1 ¼ 0:2916, a 12-cycle coexists with a 8-piece chaotic attractor, and their basins of attraction are depicted in white and red, respectively. The basin of infinity B 1 ð Þ is depicted in gray. Furthermore, the two-dimensional phase plane has been divided into three areas (Z 0 , Z 2 and Z 4 ) by two branches of the critical curves, LC ðaÞ and LC ðbÞ . It can be clearly seen from Fig. 8a that the upper end of LC ðbÞ has crossed the line segment w À1 1 and the intersection is denoted as the point C. Another branch of the critical curves, LC ðaÞ , is completely contained in the basin of attraction. In this situation, the rank-1 preimage w À1 1 of w 1 can be seen as the combination of the line segments E ð1Þ À1 C and CE ð3Þ À1 , that is, À1 . It is worth noting that the line segment CE ð3Þ À1 locates in Z 0 area, while the line segment E ð1Þ À1 C locates in Z 2 area. The rank-1 preimage w À1 2 of w 2 can be given as w À1 À1 . Since one branch of the critical curves, LC ðbÞ , has crossed the basin boundary oB 1 , there is a ''gap'' that appears on the right boundary of the basin of attraction. The boundary of this ''gap'' can be represented by a curve DE. Hence, when this occurs, the basin boundary oB 1 can be obtained as where the curve DE is the rank-1 preimage of the line segment E ð1Þ À1 C and it is also the rank-2 preimage of the line segment w 1 . Moreover, the specific expression of the curve DE can be given as When the value of k 1 increases to 0.3066, the 8piece chaotic attractor will become the 18 pieces of b Fig. 5 The parameters are fixed at a ¼ 2:9147, b ¼ 2:8236, c ¼ 0:3552, h ¼ 0:8802, k 2 ¼ 0:3744, s 1 ¼ 0:8872, s 2 ¼ 0:2240, c 1 ¼ 0:8811 and c 2 ¼ 0:7042. a The 1-D bifurcation diagram when the initial conditions are fixed as p 1 ¼ 0:7694 and p 2 ¼ 0:0802. b Another 1-D bifurcation diagram when the initial conditions are fixed as p 1 ¼ 1:5787 and p 2 ¼ 0:7353. c There exists a 4-cycle a 1 ; a 2 ; a 3 ; a 4 f gin the basin of attraction of the attracting 4-cycle for map T 4 , belonging to two attracting closed invariant curves when k 1 ¼ 0:3278665. d The two 4-cycles node and saddle lead to two closed invariant curves when k 1 ¼ 0:3278675, as shown in the enlargement in e. The basins are shown for map T 4 . f An attracting 2-cycle coexists with two closed invariant curves when k 1 ¼ 0:3410 and a saddle 2-cycle belongs to the boundary of its basin of attraction. g When k 1 ¼ 0:3509, the pair of 2cycles (the attracting node and the saddle) are now very close to each other and close to the reverse saddle-node bifurcation leading to their disappearance. The basins are shown for map T 4 . h There only attractors are now two pieces of chaotic attractor exist in the set of bounded trajectories, when k 1 ¼ 0:3510. The basins are shown for map T 2 attractors, as shown in Fig. 8b. A change which is worth to note is that another end of LC ðbÞ has crossed the line segment E ð2Þ À1 D, which is denoted as the line segment HI marked by a red box. That is to say, now, the rank-1 preimage w À1 2 of w 2 can be expressed as À1 , and the line segment HI enters into Z 2 area from Z 0 area. This contact between the critical curve LC ðbÞ and the basin boundary will lead a ''gap'' appearing on the upper basin boundary, and the boundary of this ''gap'' is remarked as the curve FG. It can be known that the curve FG is the rank-1 preimage of the line segment HI and it is also the rank-2 preimage of the line segment w 2 . The rank-1 preimage w À1 1 of w 1 can be expressed as b Fig. 6 The parameters are chosen as a À1 . In addition, the specific expression of the curve FG can be given as Therefore, the basin boundary oB 1 under this situation can be obtained by From Fig. 8c, we can clearly see that the line segment HI becomes longer and the ''gap'' on the upper basin boundary also becomes bigger when k 1 ¼ 0:3594. It is clear that the 18 pieces of chaotic attractor disappears and a new chaotic attractor appears in the basin of attraction. With gradual increasing the value of k 1 to 0.3874, the size of this new chaotic attractor apparently becomes bigger and its boundary has contacted with the segment CF of the basin boundary, as shown in Fig. 8d. As a result, there is a new ''gap'' that appears on the lower basin boundary, whose frontier is denoted as the curve JK. In addition, the curve JK is the rank-1 preimage of the curve CF and it is also the rank-3 preimage of w 2 . Therefore, the basin boundary oB 1 can be seen as a combination of crossed the line segment w À1 2 when k 1 ¼ 0:3066. c Two ''gaps'' appear on the basin boundary and a new chaotic attractor appears when k 1 ¼ 0:3594. d The size of chaotic attractor has become bigger when k 1 ¼ 0:3874 À1 . With further increasing of k 1 , the final bifurcation [25,26] will appear in system T, and then the chaotic attractor will disappear. However, the skeleton of its attracting domain will still exist in the basin of attraction, which is consisted by dense scattered points. This is the so-called the ''ghost.''

Conclusions
In this research, we first analyze a static Bertrand duopoly game taking differentiated products and service factor into consideration. It is concluded that a unique Nash equilibrium exists. By using the gradient mechanism, a dynamic Bertrand duopoly game is established. The local stability analyses on the four equilibria can be further obtained with the help of the knowledge of nonlinear dynamics and the Jury criterion. It can be shown that all the boundary equilibrium points are unstable. About the Nash equilibrium point, we show that it cannot lose stability via a Neimark-Sacker bifurcation.
In the part of local bifurcation, we can see that the system parameters, such as the adjustment speeds, service levels, service spillover effect and price sensitivity coefficient, do have significant impacts on steady properties of the system. With constant increase of the parameters, the system will become disorderly. In addition, the two manufacturers should keep the values of adjustment speeds in a relatively small range, which is conducive to steady operation of two manufacturers in the long run. The 2-cycle, flip bifurcated from the fixed point, may coexist with other 2-cycles, appearing by saddle-node bifurcations. When this occurs as a phase locking on closed invariant curves, the unstable set of the saddle cycle, connecting the saddle to the node, gives an attracting closed invariant curve. In addition, with different initial conditions, the dynamic properties may be quite different.
The phenomena of intermittent chaos appearing in the system are special dynamic behaviors that are worthy to further study and analyze. The emergence of intermittency implies that this established economic system has good capability of self-regulating. That is to say, without the help of external force, the economic system can adjust and reverse the market state form disorderly competition to orderly competition by itself. Two different kinds of intermittency phenomenon have been discussed, namely PM-I intermittency and crisis-induced intermittency.
Moreover, it is found that changing the values of the parameters will not only change the number of attracting sets, but also cause the changes in structure of the basins of attraction. There may be many ''gaps'' appearing in the basin boundary. The global bifurcations which may cause qualitative changes in the topological structure of the basin of attraction of divergent trajectories are analyzed by the critical curves. It can be found that the basin boundary is composed of the invariant axes and their higher rank preimages. with eigenvector r ð2Þ ¼ 0; 1 ð Þ. From the nonnegativity of the adjustment speeds and auxiliary variables, we can know that the boundary equilibrium point E 1 is a repelling node with eigen directions along the coordinate axes. h Proof (proof for Proposition 3) The corresponding Jacobian matrix J E 2 ð Þat the boundary equilibrium point E 2 ¼ 0; B 2b À Á can be written as The matrix J E 2 ð Þ is a lower triangular matrix, and the corresponding two eigenvalues can be obtained directly as k 1 ¼ and k 2 ¼ 1 À k 2 B with eigenvector r ð2Þ ¼ 0; 1 ð Þ along p 2 -axis, respectively. From the nonnegativity of the parameters, it is clear that one eigenvalue k 1 is located out of the unit circle. When k 2 \ 2 B , the boundary equilibrium point E 2 is a saddle point. Otherwise, the boundary equilibrium point E 2 is a unstable node. The two eigenvalue can be given by the diagonal entries, k 1 ¼ 1 À k 1 A with eigenvector r ð1Þ ¼ 1; 0 ð Þ along p 1 -axis and k 2 ¼ It is clear that when 0\k 1 \ 2 A , the boundary equilibrium point E 3 is a saddle point. Otherwise, E 3 is a unstable node. h