Symmetries and Dynamics of Generalized Biquaternionic Julia Sets Defined by Various Polynomials

Abstract

Higher-dimensional hypercomplex fractal sets are getting more and more attention because of the discovery of more and more interesting properties and visual aesthetics. In this study, the attention was focused on generalized biquaternionic Julia sets and a generalization of classical Julia sets, defined by power and monic higher-order polynomials. Despite complex and quaternionic Julia sets, their biquaternionic analogues are still not well investigated. The performed morphological analysis of 3D projections of these sets allowed for definition of symmetries, limit shapes, and similarities with other fractal sets of this class. Visual observations were confirmed by stability analysis for initial cycles, which confirm similarities with the complex, bicomplex, and quaternionic Julia sets, as well as manifested differences between the considered formulations of representing polynomials.

1. Introduction

The studies of Mandelbrot in the 1970s on the iterative equation:

$$z\leftarrow z^2+c, z,c\in ℂ,$$

(1)

defined previously by Julia and Fatou in their investigation of complex dynamic systems, resulted in the discovery of one of the most geometrically complicated fractal sets, which is known nowadays as the Mandelbrot set, or simply M-set. By variation of $c$, the M-set consists of an infinite number of Julia (usually marked as J sets) and Fatou sets, varied by their connectivity. These studies of Mandelbrot were summarized in his famous book [1]. More information on the Julia and Mandelbrot sets, including the most important theorems, can be also found in [2]. However, in addition to the classical complex iteration mapping given by (1), a generalization of (1) can be found in quaternions $ℍ$, the 4-dimensional hypercomplex algebra, i.e., in this case in (1) $z,c\in ℍ$. This is owing to the cooperation of Mandelbrot with Norton, who defined the quaternionic Mandelbrot-Julia (M-J) sets in terms of quaternions [3,4]. In parallel, Holbrook [5] performed further investigations in terms of their symmetries and connectivity.

In the following decades, numerous generalizations of M-J sets in terms of hypercomplex number spaces were developed. They include the next logical step in the hypercomplex generalization of the M-J sets, that is, their generalization in terms of octonions $\mathbb{O}$ (the 8-dimensional hypercomplex algebra), proposed by Griffin and Joshi [6] with further investigation of their properties [7,8], namely, the structural transitions characteristics for a large class of octonionic maps. Further generalizations of M-J sets were not possible, since the hypercomplex algebra of 16-dimensional sedenions $\mathbb{S}$, being the next generalization after octonions, does not meet any of fundamental properties necessary for addition and multiplication in an iteration process (see (1)), namely, sedenions are not commutative, associative, and alternative. The same is applicable to all subsequent generalizations. A detailed analysis of these properties can be found in [9].

Later, several other interesting generalizations were presented. In particular, Dixon et al. [10] investigated higher-dimensional generalizations of quadratic polynomial M-J sets defined in hypercomplex Clifford algebras. In this study, the authors demonstrated the fractal nature and symmetry properties of higher-dimensional hypercomplex M-J sets. Wang and Sun [11] proposed a generalization in terms of the power of the polynomial (1) in quaternions, namely:

$$z\leftarrow z^p+c, z,c\in ℍ, p\in \mathbb{N}.$$

(2)

In this study, the authors demonstrated the symmetry properties of quaternionic M-J sets, including those with negative power p and the stability analysis of these sets. In later studies, Wang’s team also investigated hypercomplex generalizations defined in terms of power p as in (2) [12], which was an extension of previous studies by Dixon et al. [10].

The new family of generalizations of M-J sets was established by Rochon and his team. Their studies focused on generalizations in terms of tensor product algebras: he introduced bicomplex ($ℂ\otimes ℂ$ or ${ℂ}_2$) M-J sets in [13], later studied by Wang and Song [14] as well, and further generalizations up to multicomplex ${ℂ}_n$ M-J sets [15,16,17]. This short overview shows that the hypercomplex generalizations of M-J sets are currently being intensively studied by numerous research teams.

Proceeding to the generalization of M-J sets in terms of tensor product algebras, it was found that there exist M-J sets defined in biquaternions $ℂ\otimes ℍ$. The first introduction of biquaternionic M-J sets was given by Gintz [18], where he described the fundamentals of algebra of biquaternions and showed several renderings of 3D projections of these fractal sets. Later, Bogush et al. [19] presented their studies on the symmetry of biquaternionic J sets. Studies by the first author on biquaternionic fractal sets were established in [20], where their graphical representation and preliminary analysis based on comparison with other complex and hypercomplex J sets were presented. Furthermore, a more detailed analysis of these fractal sets was performed and described in [9,21], where some unique properties for these sets were demonstrated.

This paper aimed to perform a systematic analysis of the symmetry properties of biquaternionic J sets of type (2) and their additional polynomial variation, as well as to investigate their stability to confirm the observations made during morphological analysis of these sets. Morphological analysis and evaluation of stability for biquaternionic J sets are presented for the first time, which is the main novelty of the study. Moreover, biquaternionic J sets of type (2) were also defined and investigated for the first time.

The rest of the paper is organized as follows. Section 2 consists of fundamental properties of biquaternions, as well as basic operations necessary for constructing biquaternionic J sets, which was also presented in this section. Further, in Section 3, the results of a morphological analysis are presented, including the appearing symmetries in both types of the investigated J sets. The results of stability analysis of initial cycles are presented and discussed in Section 4. The main conclusions are derived in Section 5.

2. Preliminaries

2.1. Biquaternions and Their Properties

The algebra of biquaternions, proposed by Hamilton in 1844 as one of the alternatives to quaternions [9], is a tensor product $ℂ\otimes ℍ$ of the field of complex numbers $ℂ$ and quaternions $ℍ$; therefore, it is sometimes called a complexified quaternion. It is a 4-algebra, which is isomorphic to ${ℂ}_2$ and some Clifford algebras [22,23]. In symbolic form, a biquaternion $\tilde{q}$ can be presented in a complexified form:

$$ℂ\otimes ℍ:=\left\{\tilde{q}=a_1+a_2{i}_1+a_3{i}_2+a_4{i}_3|a_n\in ℂ\right\},$$

(3)

or extended alternative form

$$ℂ\otimes ℍ:=\left\{\tilde{q}=\left(g_1+h_{1}j\right)+\left(g_2+h_{2}j\right)i_1+\left(g_3+h_{3}j\right)i_2+\left(g_4+h_{4}j\right)i_3|g_n,h_n\in ℝ\right\},$$

(4)

where $i_n$ and $j$ are imaginary units, i.e., $i_1^2=i_2^2=i_3^2=-1$, $j^2=-1$.

The fundamental operations on biquaternions, such as addition and multiplication, are the same as for quaternions, i.e., the algebra of biquaternions is alternative. This means that multiplication is performed according to special rules according to their non-commutativity. To simplify these fundamental operations, one can apply the idempotent representation of biquaternions proposed in [24]:

$${\tilde{q}}_{1,2}=\frac{1}{2}\pm \frac{1}{2}\xi j,$$

(5)

where ${\tilde{q}}_{1,2}={\omega }_1+\xi {\omega }_2$, ${\omega }_1,{\omega }_2\in ℂ$, $\xi \in ℂ\otimes ℍ$, ${\xi }^2=-1$, which makes it possible to add and multiply biquaternions element-wise, so that two biquaternions ${\tilde{q}}_1=a_1+a_2{i}_1+a_3{i}_2+a_4{i}_3$ and ${\tilde{q}}_2=b_1+b_2{i}_1+b_3{i}_2+b_4{i}_3$, $a_n,b_n\in ℂ$ can be added and multiplied as follows:

$${\tilde{q}}_1+{\tilde{q}}_2:=\left(a_1+b_1\right)+\left(a_2+b_2\right)i_1+\left(a_3+b_3\right)i_2+\left(a_4+b_4\right)i_3,$$

(6)

$$\begin{gathered} {\tilde{q}}_1:={\tilde{q}}_2:=\left(a_1{b}_1-a_2{b}_2-a_3{b}_3-a_4{b}_4\right)+\left(a_1{b}_2+a_2{b}_1+a_3{b}_4-a_4{b}_3\right)i_1+ \\ \left(a_1{b}_3-a_2{b}_4+a_3{b}_1+a_4{b}_2\right)i_2+\left(a_1{b}_4+a_2{b}_3-a_3{b}_2-a_4{b}_1\right)i_3. \end{gathered}$$

(7)

The full list of properties and basic operations on biquaternions can be found in [25].

2.2. Julia Sets within Biquaternions

2.2.1. Power Polynomials

The generalized J sets defined within biquaternions in the form

$$z\leftarrow z^p+c, z,c\in ℂ\otimes ℍ, p\ge 2$$

(8)

were defined and analyzed in one of the previous author’s publications [21]. Moreover, the condition of $p\ge 2$ was discussed there. The main definitions and theorems for biquaternionic J sets are followed below for consistency.

Definition 1.

The generalized “filled” J set, obtained while iterating recursive Equation (8) is mapped into the biquaternionic vector space with a limited trajectory of c, and

$$J_{ℂ\otimes ℍ}^p=\left\{c\in ℂ\otimes ℍ|f_c^{\left(p\right)}\left(z\right)\nrightarrow \infty \mathrm{if} p\to \infty \right\},$$

(9)

And, thus,$J_{ℂ\otimes ℍ}^p$is bounded [21].

Theorem 1.

$J_{c,ℂ\otimes ℍ}^p=J_{\left(c_1-c_{2}j\right){\tilde{q}}_1+\left({\tilde{c}}_1+{\tilde{c}}_{2}j\right){\tilde{q}}_2,ℂ\otimes ℍ}^p=J_{c_1-c_{2}j,ℂ}^p{\times }_e{J}_{{\tilde{c}}_1+{\tilde{c}}_{2}j,ℍ}^p$ [21].

Proof of Theorem 1.

Let $c\in ℂ\otimes ℍ$. For the recursive equation defined by (8) and $f_c^{\left(p\right)}\left(z\right):=(f_c^{\left(p-1\right)}◦f_c)\left(z\right)$, according to Definition 1, $f_c^{\left(p\right)}\left(0\right)$ has a bounded orbit $\forall s\in \mathbb{N}$. Additionally, when $p\ge 2,$ one obtains

$$f_c\left(z\right)=z^p+c=\left[{\left(z_1-z_{2}j\right)}^p+\left(c_1-c_{2}j\right)\right]{\tilde{q}}_1+\left[{\left({\tilde{z}}_1+{\tilde{z}}_{2}j\right)}^p+\left({\tilde{c}}_1+{\tilde{c}}_{2}j\right)\right]{\tilde{q}}_2,$$

(10)

where $z=\left(z_1-z_{2}j\right){\tilde{q}}_1+\left({\tilde{z}}_1+{\tilde{z}}_{2}j\right){\tilde{q}}_2$, and $c=\left(c_1-c_{2}j\right){\tilde{q}}_1+\left({\tilde{c}}_1+{\tilde{c}}_{2}j\right){\tilde{q}}_2$, thus

$$f_c^{\left(p\right)}\left(z\right)=f_{c_1-c_{2}j}^{\left(p\right)}\left(z_1-z_{2}j\right){\tilde{q}}_1+f_{{\tilde{c}}_1-{\tilde{c}}_{2}j}^{\left(p\right)}\left({\tilde{z}}_1-{\tilde{z}}_{2}j\right){\tilde{q}}_2.$$

(11)

Considering that $f_c^{\left(p\right)}\left(0\right)=f_{c_1-c_{2}j}^{\left(p\right)}\left(0\right){\tilde{q}}_1+f_{{\tilde{c}}_1-{\tilde{c}}_{2}j}^{\left(p\right)}\left(0\right){\tilde{q}}_2$ is bounded when $p\to \infty$, $f_{c_1-c_{2}j}^{\left(p\right)}\left(0\right){\tilde{q}}_1$ and $f_{{\tilde{c}}_1-{\tilde{c}}_{2}j}^{\left(p\right)}\left(0\right){\tilde{q}}_2$ are also bounded when $p\to \infty$. Then, $c_1-c_{2}j\in J_{ℂ}^p$, ${\tilde{c}}_1-{\tilde{c}}_{2}j\in J_{ℍ}^p$, and $c=\left(c_1-c_{2}j\right){\tilde{q}}_1+\left({\tilde{c}}_1+{\tilde{c}}_{2}j\right){\tilde{q}}_2\in J_{ℂ}^p{\times }_{\tilde{q}}{J}_{ℍ}^p$, and $J_{ℂ\otimes ℍ}^p\subset J_{ℂ}^p{\times }_{\tilde{q}}{J}_{ℍ}^p$ [21]. □

2.2.2. Monic Higher-Degree Polynomials

To investigate the homeomorphisms of M-J sets in a complex plane, Branner and Fagella [26] introduced the family of monic higher-degree polynomials in the following form:

$$z\leftarrow z{\left(z+c\right)}^p.$$

(12)

This allowed us to analyze the connectedness and complex dynamics of limbs and external rays in M-J sets. Later, Zireh investigated the properties of a family of polynomials of type (12) defined in bicomplex numbers in terms of connectedness [27]. From the defined theorems and their proofs presented in [27], it can easily be seen that the same applies for (12) defined in biquaternions; therefore, these theorems are omitted in this paper.

The fundamentals on biquaternionic operations and definitions of J sets defined within power and monic higher-degree polynomials are used for their implementation for visualization and morphology analysis purposes, as well as for analysis of stability of these fractal sets presented in Section 3 and Section 4, respectively.

3. Symmetry of Biquaternionic Julia Sets

The first results on the construction of biquaternionic J sets in [18] were an inspiration for further investigation of these fractal sets, and their preliminary analysis was performed in [20,21], which focused on investigation of the symmetry of well-known analogues of J sets on a complex plane as well as the connectedness of these sets. In the following study, a systematic analysis in terms of symmetry of biquaternionic J sets of both types presented in Section 2.2 is presented.

For analysis of the symmetry of biquaternionic J sets, their 3D projections onto ${ℝ}^3$ are used. Following the previous visualization techniques commonly used for the representation of 4D fractal sets, see, e.g., [11,13,14,20,21], where the 3D projection is obtained by setting $a_4$ in the quadruplet (3) to 0. As shown, for example, in [11], this procedure does not cause a loss of generality.

All 3D projections of the biquaternionic J sets presented in this paper were prepared in ChaosPro 4.0 freeware software. To implement addition and multiplication operations on biquaternions according to the rules presented in Section 2.1, which are necessary for the construction of biquaternionic J sets according to (8) and (12), a new library was prepared in the native programming language in ChaosPro, which is based on C/C++.

3.1. Symmetry of Biquaternionic Julia Sets Defined by Power Polynomials

As shown in [21], biquaternionic J sets do not resemble their quaternionic or bicomplex analogues in terms of symmetry. Wang and Sun in [11] demonstrated the rotational symmetry of quaternionic J sets, whereas Rochon showed the quadrilateral symmetry of bicomplex J sets, which was confirmed in [21]. The conjecture of uniqueness of biquaternionic J sets in terms of their symmetry made in [18] was confirmed in further studies presented in [21].

A simple way to show this uniqueness is to set $c=0$ in (8). The resulting recursion equation leads to the appearance of the analog of the so-called “unit disc” in the case of the $ℂ$-plane [4], i.e., $\left\{z\in ℂ:\left|z\right|<1\right\}$. Previous studies on quaternionic and bicomplex J sets for such conditions [4,28,29] show that in the first case when $p\to \infty ,$ one obtains a hypersphere (which was observed also in [11]), whereas in the second case, the resulting shape is a Steinmetz hypersolid. This can be generalized for their higher-dimensional analogues as well. To analyze the biquaternionic J sets for $c=0$ and increasing $p$, their representative 3D projections are shown in Figure 1.

From the above projections, one can make the following observations on the properties of biquaternionic J sets.

Remark 1.

For $c=0,$ J sets have nine planes of symmetry regardless of the value of $p$.

From the visualizations shown in Figure 1, the symmetry with respect to the plane ${\Im }_1{\Im }_2$ in the coordinates defined as $\Re {\Im }_1{\Im }_2$ for a 3D projection (${\Im }_3$ is not taken into consideration) is well visible. To show the other symmetry planes, it is essential to visualize biquaternionic sets in other projections. Let us analyze the limit cases presented in Figure 1a,f. The orthogonal 3D projections are presented in Figure 2 and Figure 3, respectively, for $p=2$ and $p=1000$. Despite the increase in $p$, the symmetry of biquaternionic J sets is preserved. One can observe that the symmetry of biquaternionic J sets is similar to that of bicomplex J sets (cf. the results presented in [29]).

Remark 2.

The shape of the biquaternionic J set tends to be a unique nontrivial set when $p\to \infty$.

As discussed earlier, for $p\to \infty ,$ every M-J set of type (8) defined in a given number space tends to a set with a unique geometry characteristic for a given number space. Naturally, with an increase of $p$, the value of $c$ has less influence and for $p\to \infty ,$ it loses this influence completely. In the case of biquaternionic J sets, one can observe a unique shape of the J sets tend when $p\to \infty$, as seen in Figure 3. Because of the high value of $p,$ one can consider this shape as very similar to the limit shape. The top 3D projection presented in Figure 3a depicts that this set contains characteristic horn-like structures, which are formed at the very beginning of the increase of $p$ (see Figure 1). Moreover, one can observe that the rotational symmetry in the $\Re {\Im }_1$ plane is broken, i.e., the shape presented in Figure 3a is not the same if it is rotated by 90°. This is a direct consequence of the definition of these fractal sets in a hypercomplex number space composed from a tensor product of two dimensionally unequal algebras, $ℂ\otimes ℍ$.

Remark 3.

By reducing the number of non-zero elements in $z$ and $c$, it is possible to obtain the analogues of the J sets defined in the $ℂ$-plane.

To visualize the relationship between biquaternionic and classical J sets defined in the $ℂ$-plane, it is essential to analyze sets with a specific geometry and values of $c$, which demonstrate a similarity between the mentioned groups of fractal sets. The biquaternionic analogs of famous Dendrite and San Marco J sets are presented in Figure 4.

The Dendrite presented in Figure 4a shows a geometric structure characteristic for all types of number spaces used for construction of fractal sets (e.g., it looks like similar to J sets constructed in the $ℂ$-plane and ${ℂ}_n$-spaces as well as in the $ℍ$-space). The San Marco J set presented in Figure 4b–d in various projections reveals a much more complicated shape and it is difficult to observe its $ℂ$-plane analogue. The characteristic shape of the San Marco J set is quite visible in the top view (Figure 4c), which can be directly obtained by cutting its biquaternionic analogue in the $\Re {\Im }_1$ plane. However, it is interesting to observe that in the $\Re {\Im }_2$ plane (Figure 4d), this set has the form of a square. The lack of rotational symmetry typical for quaternionic or quadrilateral symmetry of bicomplex J sets as well as the above-discussed differences in various orthogonal projections is the result of a tensor product of dimensionally unequal algebras in the number space, where these sets were defined. However, the examples presented in Figure 4 proved that the biquaternionic J sets remain fractals.

3.2. Symmetry of Biquaternionic Julia Sets Defined by Monic Higher-Degree Polynomials

To investigate the homeomorphisms of J sets defined in biquaternions, the recurrence in Equation (12) is applied to generate visualizations with the same assumptions made in the preamble of Section 3 and Section 3.1. This allowed for observing the variability of the shapes of particular biquaternionic J sets defined by (12) with increasing $p$ as well as determining the limit shape when $p\to \infty$. The 3D projections of the J sets constructed using (12) for selected values of $p$ are presented in Figure 5.

From the above-presented visualizations, it can be observed that the resulting J sets for various values of $p$ reveal some differences, especially for low values of $p$; however, for higher values of $p$, these differences become less and less recognizable. Finally, when $p\to \infty$, the resulting limit figure is the same as for the classical J set defined by (8), which is evident when substituting to $p=\infty$ to (8) and (12). Due to this, all symmetries characteristic for biquaternionic J sets defined by (8) are preserved.

Nevertheless, the dynamics of the biquaternionic J sets defined by (12) are completely different, which can be observed based on previously analyzed 3D projections of the analogues of Dendrite and San Marco fractals; see Figure 6. One can observe that the resulting shapes do not correspond to the biquaternionic J sets defined by (8), nor to their analogues defined in the $ℂ$-plane.

The J sets of both types, investigated from the point of view of their morphology and appearance symmetries, were further investigated in light of their stability, in particular, the determination of the boundaries of fixed points for variable $p$, which is the subject of the next section.

4. Stability of Biquaternionic Julia Sets

In this section, the stability of generalized Julia sets defined in biquaternions is analyzed. The analysis was performed for both investigated recurrence Equations (8) and (12) for 1-, 2-, and 3-cycle stability and various values of $p$, which made it possible to observe variation of stability for these sets when $p\to \infty$.

4.1. Stability of Biquaternionic Julia Sets Defined by Power Polynomials

4.1.1. 1-Cycle Stability

There are numerous similarities between biquaternions, quaternions, and bicomplex numbers; that is, all of these number spaces represent 4-dimenional hypercomplex numbers with numerous similar properties (see Section 2.1 for more details). Moreover, $ℂ\otimes ℍ$ is isomorphic to ${ℂ}_2$ [22]. Thus, the dynamics and stability of Julia sets defined as biquaternions are expected to reveal some similarities as in the case of their morphological analysis performed in Section 3.

Following [11,30], let us represent the n-th generation biquaternion in terms of Pauli matrices:

$$z_n=x_{n}I+i\mathsf{\sigma }{V}_n,$$

(13)

where $\mathsf{\sigma }=\left({\mathsf{\sigma }}_1,{\mathsf{\sigma }}_2,{\mathsf{\sigma }}_3\right)$ is a set of 2 × 2 Pauli matrices:

$${\mathsf{\sigma }}_1=\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right], {\mathsf{\sigma }}_2=\left[\begin{array}{cc}0& -i\\ i& 0\end{array}\right], {\mathsf{\sigma }}_3=\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right],$$

(14)

completed by 2 × 2 invariant unit matrix

$$I=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right],$$

(15)

and $\left\{x,V\right\}\equiv \left\{x,y,v,w\right\}$, i.e., the biquaternion can be represented by its scalar and vector parts:

$$S_n=\left\{x_n\right\}, V_n=\left\{y_n,v_n,w_n\right\}, x_n,y_n,v_n,w_n\in ℂ.$$

(16)

Note that for quaternions $\mathsf{\sigma }=-i\left({\mathsf{\sigma }}_3,{\mathsf{\sigma }}_2,{\mathsf{\sigma }}_1\right)$, i.e., the quaternion needs to be multiplied by $-i$ and indices of ${\mathsf{\sigma }}_n$ need to be replaced (see, e.g., [25]), whereas for biquaternions, Pauli matrices (14) can be used directly without additional operations [31].

The recurrence relation for the biquaternionic quadratic map can be presented as:

$$z_{n+1}=z_n^2+c,$$

(17)

where

$$c=d_{1}I+i\mathsf{\sigma }D,$$

(18)

$$D=\left\{d_2,d_3,d_4\right\}, d_1,d_2,d_3,d_4\in ℂ.$$

(19)

Taking into account (13) and (16), we may write

$$z^2={\left(SI+i\mathsf{\sigma }V\right)}^2=\left(S^2-V^2\right)I+2i\mathsf{\sigma }SV.$$

(20)

The iterations of (17) considering (18) and (20) are represented by:

$$S_{n+1}=S_n^2-V_n^2+d_1,$$

(21)

$$V_{n+1}=2S_n{V}_n+D.$$

(22)

To determine fixed points and boundary regions for the biquaternionic quadratic map, it is necessary to determine the eigenvalues of the matrix with perturbed terms. Introducing the perturbation of $z$, we have:

$$S=S_n+\epsilon S_n,$$

(23)

$$V=V_n+\epsilon V_n,$$

(24)

where $\epsilon$ is the so-called small perturbation parameter. This leads to:

$$\epsilon S_{n+1}=S-S_{n+1}=S^2-V^2-S_n^2+V_n^2,$$

(25)

$$\epsilon V_{n+1}=V-V_{n+1}=2\left(SV-S_n{V}_n\right),$$

(26)

which can be simplified to the following form:

$$\epsilon S_{n+1}=2\left(S_n\epsilon S_n-V_n\epsilon V_n\right)+O\left(\epsilon S_n^2,\epsilon V_n^2\right),$$

(27)

$$\epsilon V_{n+1}=2\left(S_n\epsilon V_n-V_n\epsilon S_n\right)+O\left(\epsilon S_n^2,\epsilon V_n^2\right),$$

(28)

where $O\left(\epsilon S_n^2,\epsilon V_n^2\right)$ is the Landau symbol, which denotes truncation in terms of the order higher than the first one in solutions (27) and (28). Expanding $\epsilon V_n$ in terms of the orthogonal triplet $\left(i_1,i_2,i_3\right),$ we have:

$$\epsilon V_n=\epsilon {\alpha }_n{i}_1+\epsilon {\beta }_n{i}_2+\epsilon {\gamma }_n{i}_3$$

(29)

and omitting the Landau symbol in (27) and (28), we obtain:

$$\epsilon S_{n+1}=2\left(S\epsilon S_n-V\epsilon {\alpha }_n{i}_1\right),$$

(30)

$$\epsilon V_{n+1}=2\left(S_n\epsilon {\alpha }_n{i}_1-V\epsilon S_n\right).$$

(31)

From the matrix

$$\left(\begin{array}{c}\epsilon S_{n+1}\\ \epsilon {\alpha }_{n+1}\end{array}\right)=\Lambda \left(\begin{array}{c}\epsilon S_n\\ \epsilon {\alpha }_n\end{array}\right),$$

(32)

we can determine $\Lambda$ considering (30) and (31):

$$\Lambda =\left[\begin{array}{cc}2S& -2V{i}_1\\ 2V{i}_1& 2S\end{array}\right],$$

(33)

which makes it possible to obtain the eigenvalues:

$${\lambda }_{1,2}=2\left(S\mp V{i}_1\right).$$

(34)

Considering that $\left|{\lambda }_n\right|=2\left(S+V\right)$ and noting that $\left|\lambda \right|<1$, the boundary region of the biquaternionic quadratic map is determined by

$$S+V=\frac{1}{2},$$

(35)

which leads to:

$$\Vert z\Vert =\frac{1}{2}.$$

(36)

According to the trigonometric representation proposed in [11]:

$$z=\Vert z\Vert \left(\cos \theta +m\sin \theta \right),$$

(37)

where $m=K/\left|K\right|$, we have:

$$S=\frac{1}{2}\cos \theta , V=\frac{1}{2}m\sin \theta , 0\le \theta <2\pi .$$

(38)

The boundary of the region of fixed points for the biquaternionic quadratic map is thus defined by:

$$d_1=\frac{1}{4}\left(2\cos \theta -\cos 2\theta \right), D=\frac{m}{4}\left(2\sin \theta -\sin 2\theta \right), 0\le \theta <2\pi .$$

(39)

In the same way, we can obtain solutions for higher-degree biquaternionic maps. Let us analyze the next cases discussed previously for $p=3$ and $p=5$. The boundary region of fixed points for these cases are, respectively, as follows:

$$p=3: \Vert z\Vert =\frac{1}{\sqrt{3}}, d_1=\frac{1}{\sqrt{3}}\left(\cos \theta -\frac{1}{3}\cos 3\theta \right), D=\frac{m}{\sqrt{3}}\left(\sin \theta -\frac{1}{3}\sin 3\theta \right),$$

(40)

$$p=5: \Vert z\Vert =\frac{1}{\sqrt[4]{5}}, d_1=\frac{1}{\sqrt[4]{5}}\left(\cos \theta -\frac{1}{5}\cos 5\theta \right), D=\frac{m}{\sqrt[4]{5}}\left(\sin \theta -\frac{1}{5}\sin 5\theta \right).$$

(41)

For detailed calculations for these cases, see Appendix A.

Taking into account (37) and (39)–(41), one can observe that these formulas lead to their general form proposed in [11] for quaternions:

$$\Vert z\Vert ={\left|p\right|}^{-\frac{1}{p-1}},$$

(42)

$$\{\begin{array}{c}{d}_1={\left|p\right|}^{-\frac{1}{p-1}}\left(\cos \theta -\cos p\theta \right),\\ D={\left|p\right|}^{-\frac{1}{p-1}}\left(\sin \theta -\sin p\theta \right).\end{array}$$

(43)

Using (42), one can get the 1-cycle stability for the limit shape considered in the previous section. Substituting $p=\infty$ to (42), the direct result will have an indeterminate form; however, when representing it in terms of limits, we have

$$\underset{p\to \infty }{\lim }{\left|p\right|}^{-\frac{1}{p-1}}=1.$$

(44)

This confirms the convergence of the 1-cycle stability region of the generalized biquaternionic J set to the unit disc.

4.1.2. 2-Cycle Stability

The 2-cycle for the biquaternionic quadratic map is defined by the equation:

$$z_{n+2}={\left(z_n^2+c\right)}^2+c=z_n^4+2z_n^{2}c+c^2+c,$$

(45)

where $c$ is given by (18) and (19). Using the notation given by (13) and (16), the powers of $z$ can be given by (20) and the following expression:

$$z^4={\left(SI+i\mathsf{\sigma }V\right)}^4=\left(S^4-6S^2{V}^2+V^4\right)I+4i\mathsf{\sigma }\left(S^{3}V-S{V}^2\right).$$

(46)

The iterations of (45) considering (18), (20) and (46) are represented by:

$$S_{n+2}=S_n^4-6S_n^2{V}_n^2+V_n^4+2a\left(S_n^2-V_n^2\right)+d_1^2+d_1,$$

(47)

$$V_{n+2}=4S_n{V}_n\left(S_n^2-V_n+D\right)-D^2+D.$$

(48)

When introducing the perturbation of $z$ by (23) and (24) as in Section 4.1.1, one obtains:

$$\epsilon S_{n+2}=S-S_{n+2}=S^4+V^4-6\left(S^2{V}^2-S_n^2{V}_n^2\right)+2\left(S^2-V^2-S_n^2-V_n^2\right)-S_n^4-V_n^4,$$

(49)

$$\epsilon V_{n+2}=V-V_{n+1}=4\left(S^{3}V-S{V}^2+SV\right)-4\left(S_n^3{V}_n-S_n{V}_n^2+S_n{V}_n\right).$$

(50)

After simplification of (49) and (50) and omitting the Landau symbol as in Section 4.1.1, we have:

$$\epsilon S_{n+2}=4S\epsilon S_n\left(S^2-3V^2+1\right)+4V\epsilon {\alpha }_n{i}_1\left(V-3S^2-1\right),$$

(51)

$$\epsilon V_{n+2}=4S\epsilon {\alpha }_n{i}_1\left(S^2-2V+1\right)-4V\epsilon V_n\left(V+3S^2+1\right).$$

(52)

From matrix (32) we can determine $\Lambda$:

$$\Lambda =\left[\begin{array}{cc}4S\left(S^2-3V^2+1\right)& 4V{i}_1\left(V-3S^2-1\right)\\ -4V{i}_1\left(V+3S^2+1\right)& 4S\left(S^2-2V+1\right)\end{array}\right],$$

(53)

and the corresponding eigenvalues:

$${\lambda }_{1,2}=4\left(S^3-1.5S{V}^2-SV+S\right)\pm 2Vi\mathsf{\sigma }\sqrt{36S^4+9S^2{V}^2-12S^{2}V+28S^2-4V^2+4}.$$

(54)

Taking the module of (54), we have:

$$\Vert z\Vert =\frac{1}{\sqrt[3]{4}},$$

(55)

which leads to the determination of the boundary of the region of fixed points for the 2-cycle using (37) and (38).

In a similar way, we can determine the boundary regions of fixed points for the 2-cycle for $p=3$ and $p=5$, where:

$$p=3: \Vert z\Vert =\frac{1}{\sqrt[8]{9}},$$

(56)

$$p=5: \Vert z\Vert =\frac{1}{\sqrt[24]{25}}.$$

(57)

The presentation of the detailed calculations was omitted in this paper because of its severe extension.

From the obtained values, one can derive the general formula for $\Vert z\Vert$ for the 2-cycle:

$$\Vert z\Vert ={\left|p\right|}^{-\frac{2}{p^2-1}}.$$

(58)

Assuming $p=\infty$ in (58), one leads to the conclusion that the 2-cycle stability region of the generalized biquaternionic J set also leads to the unit disc.

4.1.3. 3-Cycle Stability

Starting again from the biquaternionic quadratic map, the 3-cycle analogue is defined by the equation:

$$z_{n+3}={({(z_n^2+c)}^2+c)}^2+c=z_n^8+4z_n^{6}c+6z_n^4{c}^2+2z_n^{4}c+4z_n^2{c}^3+4z_n^2{c}^2+c^4+2c^3+c^2+c,$$

(59)

where $c$ is given by (18) and (19). Using the notation given by (13) and (16), the powers of $z$ can be given by (20), (46), and the following expressions:

$$z^6=\left(S^6-15S^4{V}^2+15S^2{V}^4-V^6\right)I+2i\mathsf{\sigma }\left(3S^{5}V-10S^3{V}^3+3S{V}^5\right),$$

(60)

$$z^8=\left(S^8-28S^6{V}^2+70S^4{V}^4-28S^2{V}^6+V^8\right)I+8i\mathsf{\sigma }\left(S^{7}V-7S^5{V}^3+7S^3{V}^5-S{V}^7\right).$$

(61)

Introducing perturbation of $z$ and simplifications, which were omitted here due to complex notation, one obtains the matrix from which the eigenvalues are calculated. Taking the module of the obtained eigenvalues, we have:

$$\Vert z\Vert =\frac{1}{\sqrt[7]{4}}.$$

(62)

Using a similar procedure, the value of $\Vert z\Vert$ was calculated for $p=3$:

$$\Vert z\Vert =\frac{1}{\sqrt[26]{27}}.$$

(63)

Calculations for the higher values of $p$ were not possible because of the hardware limitations. Because of this, the determination of the general formulas for the 3-cycle and higher cycles is not possible through the computational way. Nevertheless, it is expected that the dynamics are different for higher cycles, which finally leads to the limit shape observed in Section 3.1.

4.2. Stability of Biquaternionic Julia Sets Defined by Monic Higher-Degree Polynomials

4.2.1. 1-Cycle Stability

The recurrence relation for the biquaternionic quadratic map defined in terms of monic higher-degree polynomials (12) can be presented as:

$$z_{n+1}=z_n{\left(z_n+c\right)}^2=z_n^3+2z_n^{2}c+z_n{c}^2.$$

(64)

Using the notation given by (13) and (16), the powers of $z$ can be given by (20) and (A1). Applying the same procedure as in Section 4.1.1, one can calculate the iterations of (64) considering (18), (20) and (46), which can be represented by:

$$S_{n+1}=S_n\left(S_n^2+2S_n{d}_1-3V_n^2-4D{V}_n-D^2+d_1^2\right)-2V_n{d}_1\left(V_n+D\right),$$

(65)

$$V_{n+1}=V_n\left(-V_n^2+4S_n{d}_1+S_n^2-2D{V}_n-D^2+d_1^2\right)+2S_{n}D\left(S_n+d_1\right).$$

(66)

Introducing the perturbation of $z$ and simplifications, we can obtain the eigenvalues in the following form:

$${\lambda }_{1,2}=\frac{3}{2}{S}^2-3SV-4V-\frac{3}{2}{V}^2\pm \frac{1}{2}i\mathsf{\sigma }\sqrt{\begin{array}{c}9S^4+108S^{3}V+96S^3+18S^2{V}^2+240S^{2}V+164S^2-108S{V}^3 \\ -192S{V}^2+24SV+64S+9V^4-48V^3-100V^2-32V+4 \end{array}}$$

(67)

Taking the module of (67), we have:

$$\Vert z\Vert =\frac{\sqrt{6}}{3},$$

(68)

which leads to the determination of the boundary of the region of fixed points for a 1-cycle using (36) and (37).

Similarly, we can determine the boundary regions of fixed points for a 1-cycle for $p=3$ and $p=5$, where:

$$p=3: \Vert z\Vert =\frac{1}{\sqrt[3]{2}},$$

(69)

$$p=5: \Vert z\Vert =\frac{1}{\sqrt[5]{3}}.$$

(70)

4.2.2. 2-Cycle Stability

The 2-cycle for the biquaternionic quadratic map defined in terms of monic higher-degree polynomials (12) is given by equation:

$$z_{n+2}=z_n{(z_n{(z_n+c)}^2+c)}^2=z_n^8+4z_n^{7}c+6z_n^6{c}^2+4z_n^5{c}^3+z_n^4{c}^4+2z_n^{4}c+4z_n^3{c}^2+2z_n^2{c}^3+c^2.$$

(71)

The powers of $z$ are given by (20), (46), (60), (61), (A1), (A3), and the following expression:

$$z^7=\left(S^7-21S^5{V}^2+35S^3{V}^4-7S{V}^6\right)I+i\mathsf{\sigma }\left(7S^{6}V-35S^4{V}^3+21S^2{V}^5-V^7\right).$$

(72)

Omitting the steps presented above, one leads to the determination of modules of eigenvalues for 2-cycle stability, which in this case have the following values for the considered $p$:

$$p=2: \Vert z\Vert =\frac{1}{\sqrt[7]{4}},$$

(73)

$$p=3: \Vert z\Vert =\sqrt[14]{\frac{2}{15}},$$

(74)

4.2.3. 3-Cycle Stability

Finally, the 3-cycle for the biquaternionic quadratic map defined in terms of monic higher-degree polynomials (12) is given by the equation:

$$\begin{array}{l}{z}_{n+3}=z_n{(z_n{(z_n{(z_n+c)}^2+c)}^2+c)}^2\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=z_n^{10}+4z_n^{9}c+6z_n^8{c}^2+4z_n^7{c}^3+z_n^6{c}^4+2z_n^{6}c+4z_n^5{c}^2+2z_n^{5}c\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+2z_n^4{c}^3+4z_n^4{c}^2+2z_n^3{c}^3+z_n^2{c}^2+2z_n{c}^2+c^2.\end{array}$$

(75)

In this case, the powers of $z$ are given by all expressions mentioned in Section 4.2.2 and the following ones:

$$\begin{array}{l}{z}^9=\left(S^9-36S^7{V}^2+126S^5{V}^4-84S^3{V}^6+9S{V}^8\right)I\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+i\mathsf{\sigma }\left(9S^{8}V-84S^6{V}^3+126S^4{V}^5-36S^2{V}^7+V^9\right).\end{array}$$

(76)

$$\begin{array}{l}{z}^{10}=\left(S^{10}-45S^8{V}^2+210S^6{V}^4-210S^4{V}^6+45S^2{V}^8-V^{10}\right)I\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+2i\mathsf{\sigma }\left(5S^{9}V-60S^7{V}^3+126S^5{V}^5-60S^3{V}^7+5S{V}^9\right).\end{array}$$

(77)

Using the same computational procedure as previously, the modulus of the eigenvalue for 3-cycle stability was determined for $p=2$:

$$\Vert z\Vert =\frac{1}{\sqrt[17]{9}}.$$

(78)

Because of the increasing computational complexity (for $p=3$ the maximal power within the polynomial is 48, while for $p=5$ this value equals 180), the determination of modules of eigenvalues for higher values of $p$ were not possible because of hardware limitations.

5. Conclusions

This paper discusses the generalized biquaternionic Julia sets defined by power and monic higher-degree polynomials in terms of their symmetry and stability. The 3D projections of both types of fractal sets were constructed and discussed. The symmetry planes for these fractal sets were defined together with the definition and analysis of the limit shapes (when $p\to \infty$). The similarities between fractal sets were constructed based on power polynomials and its analogues in other hypercomplex number spaces with the same dimension, such as fractal sets defined in quaternions and bicomplex numbers. Moreover, it was shown that all of the mentioned fractal sets have the same origin as complex Julia sets. The biquaternionic Julia sets constructed based on monic higher-degree polynomials were presented for the first time, and the morphological analysis of these sets shows that their shapes and dynamics are completely different with respect to the above-discussed fractal sets defined in power polynomials, again underlining the multimodality and increasing complexity of Julia sets defined in hypercomplex number spaces.

To support the observations made during the morphological analysis, stability calculations were performed for both investigated types of fractal sets. It was shown that the dynamics of biquaternionic Julia sets defined by power polynomials reveal similarities to the dynamics of fractal sets defined in quaternions, at least in the initial cycles of stability. The severe differences between these two types of fractal sets are expected to appear at high cycles of stability; however, because of hardware limitations, these differences have not been identified yet. In contrast, the biquaternionic fractal sets defined by monic high-order polynomials reveal some similarities but different and more complex behavior, which confirms the observations made for the morphology of these sets.