Studying discrete-time dynamical systems (III): generalization of dynamical systems able to produce spirals and clusters of points

In a previous post I wondered if there were any other dynamical systems able to generate spirals and clusters of points as the found expression:

$$S_{D}=\{(x,y): (x_{n+1},y_{n+1}) = (Im(\psi(x_{n} + y_n  i)\cdot\sin{\frac{D}{f}}, Re(\psi(x_{n} + y_n i))\cdot\cos{\frac{D}{f}})\}$$

 To answer this, first I had a look to the image of the complete plot of the family of spirals: 

We can see that there is a quasi-horizontal symmetry there. Both the "upper body" ($z=a+bi , b \ge 0$) and "lower body" ($z=a+bi , b \lt 0$) are quite symmetrical, but not exactly equal. The structure also has three vertical bodies / structures, a right side arc-like structure, a central body and a left-side body. These are the regions where the density of points is bigger.

The density of points make this structures to look like they are "glowing" (e.g. they remind a bioluminiscent jellyfish). The areas where there is more density of points seems to "glow" compared with other areas of lower density.

First question: who is responsible for making the spirals and who is responsible for deciding the location and distortion of the spirals?

To answer this, first I need to amend the expression:
$$S_{D}=\{(x,y): (x_{n+1},y_{n+1}) = (Im(\psi(x_{n} + y_n  i)\cdot\sin{\frac{D}{f}}, Re(\psi(x_{n} + y_n i))\cdot\cos{\frac{D}{f}})\}$$ And make it more general, it is: $$S_{D}=\{(x,y, f(z)): (x_{n+1},y_{n+1}) = (Im(f(x_{n} + y_n  i)\cdot\sin{\frac{D}{f}}, Re(f(x_{n} + y_n i))\cdot\cos{\frac{D}{f}})\}$$ So indeed: the family of dynamical systems is even bigger. Depending on the complex function that is used we will have a different structure. All of them have similar transformations but the results depend obviously on the selected function. So basically my original question was about the above more general family of dynamical systems when $f(z)=\psi(z)$ is specifically the digamma function.

Now, going back to the digamma expression:
$$S_{D}=\{(x,y): (x_{n+1},y_{n+1}) = (Im(\psi(x_{n} + y_n  i)\cdot\sin{\frac{D}{f}}, Re(\psi(x_{n} + y_n i))\cdot\cos{\frac{D}{f}})\}$$ Let us split the expression into two versions:

Version 1: No sine or cosine are applied:
$$S_{D}=\{(x,y): (x_{n+1},y_{n+1}) = (Im(\psi(x_{n} + y_n  i)), Re(\psi(x_{n} + y_n i)))\}$$and

Version 2: No function is applied, only the crossing of $Re$ and $Im$ (the meaning of "crossing" is as explained in the question):
$$S_{D}=\{(x,y): (x_{n+1},y_{n+1}) = (Im(x_{n} + y_n  i)\cdot\sin{\frac{D}{f}}, Re(x_{n} + y_n i)\cdot\cos{\frac{D}{f}})\}$$ Version 1 will provide the following result:



Version 2 will provide the following result (for the same $n$ iterations and $f$ and $D$ parameters):



So the conclusion is: the sine and cosine expressions provide the observed horizontal quasi-symmetry of the dynamical systems, and the digamma function was providing the spiral-like structures. The combination of both expressions provide a final structure in which the spirals of the digamma function are transformed, displaced and distorted by the sine and cosine transformations.

Second question: what do we obtain if we replace the digamma function by other functions?

My observations are that we obtain different "glowing" systems, but they share similar characteristics: horizontal quasi-symmetry, and three vertical structures (left, center, right), in other words, regions in which the density of points is bigger.

These are some of those beauties:

Scaled complementary error function, $f(z)=exp(z^2) \cdot erfc(z)$. (see Python scipy.special functions list, Steven G. Johnson, Faddeeva W function implementation). $D \in [0,10^4], n \in [0,10^4]$ zoom into region $z=[+/-5]+[+/-5]i$:



$f(z) = $ Exponential integral $E_1$ of complex argument $z$. (see Python scipy.special functions list). $D \in [0,10^3], n \in [0,10^4]$ zoom into region $z=[+/-1]+[+/-4]i$:



$f(z)=$ Exponential integral $E_i$. (see Python scipy.special functions list). $D \in [0,10^3], n \in [0,10^4]$ zoom into region $z=[+/-5]+[+/-5]i$:



Finally, my favorite families so far, $f(z)=\frac{1}{(z+1)^t}, t \in \Bbb R$. $D \in [0,8 \cdot 10^3], n \in [0,8 \cdot 10^3]$ zoom into region $z=[+/-5]+[+/-5]i$::

$f(z)=\frac{1}{(z+1)}$:



$f(z)=\frac{1}{(z+1)^{1.1}}$:



$f(z)=\frac{1}{(z+1)^{1.2}}$:



$\cdots$

$f(z)=\frac{1}{(z+1)^{2}}$:



$\cdots$

$f(z)=\frac{1}{(z+1)^{3}}$:



I just can say that these structures are visually beautiful. I will keep my question regarding these families of dynamical systems open at MSE in hope that somebody will be able to provide more insights about why this transformation / mapping works in this way.