Struggling to Understand Algorithm for Displaying Polynomial Matings with Julia Sets

From what you are saying I am unfortunately almost certain that you do not have enough mathematical background to really understand the theory of matings (which is something quite advanced, not realistically understandable before grade school). You typically want to have taken an advanced course in complex analysis, topology, and some background in complex dynamics.

It is difficult to explain shortly the notion of matings, but I'll try to give a quick summary of my own understanding of the topic. An equipotential is a line where the Green function is constant and strictly positive (if you don't know what the Green function is, you need background in complex dynamics). If the Julia set is connected then this equipotential is homeomorphic to a circle. In particular, it splits the sphere into two components, one containing the Julia set and the other containing $\infty$, and each of those components are homeomorphic to a disk. If you take two different connected Julia sets and two such equipotentials, you can glue the components containing the Julia sets by identifying the two equipotentials (if you don't understand this sentence you need background in topology). What you get then is something homeomorphic to a sphere (let's call it $S_1$, and you get a continuous map $f$ defined on $S_1$ that coincides with the restrictions of both polynomials outside of the glueing line. However the range of $f$ is not $S_1$ but rather a similar object obtained by glueing two different equipotentials (the ones that are images of the previous ones by the two polynomials). So you get a continuous map $f: S_1 \to S_2$, where $S_1$ and $S_2$ are topological spaces homeomorphic to spheres.

With some work and a really deep theorem that I won't even try to explain here, called the measurable Riemann mapping theorem, you can somehow get a holomorphic map $g: \hat{\mathbb C} \to \hat{ \mathbb C}$ from this whole business. The map $g$ is conjugated to $f$ by homeomorphisms that map the $S_i$ to the Riemann sphere. However, you shouldn't consider that $g$ is a dynamical system since it is conjugated to $f$, which has different domain and range. But if you used equipotentials $G=t$ for $S_1$, then you used equipotentials $G=dt$ for $S_2$, and as $t \to 0$ the difference between these two equipotentials shrinks to zero. So you want to prove that the map $g_t$ that you get with this whole procedure has a limit when $t \to 0$. This is not true in general, but when it is the case that limit is exactly what's called the mating (in one sense) between the two polynomials.

The theory is also beyond my level of education, but a practical implementation to make pictures is relatively straight forward following Chapter 5 of Wolf Jung's paper "The Thurston Algorithm for Quadratic Matings" (preprint linked in question). However, an important thing missing in my code is detecting homotopy violations, so there is no proof that the images are correct.

In my code, the slow mating is computed as per Wolf Jung's chapter 5, pulling back the critical orbits using continuity to choose the sign of the square root. Pulling back an orbit means the next orbit's $z_n$ depends in some way on the previous orbit's $z_{n+1}$. The process has a sequence of orbits, where the orbit at time $t+s+1$ depends on the orbits at time $t + s$ (for choosing roots by continuity) and time $t + 1$ (for the input to the square root function). Increasing the granularity $s$ presumably makes the continuity test more reliable.

To render images, the pull-back in Wolf Jung's paper is inverted to a sequence of functions of the form $z\to\frac{az^2+b}{cz^2+d}$, which are composed in reverse order starting from the desired pixel coordinates. Then choose the hemisphere based on $|z|<1$ or $|z|>1$, find $w=Rz$ or $w=R/z$ and $c=c_1$ or $c=c_2$ depending on hemisphere, and continue to iterate $w→w^2+c$ until escape (or maximum iteration count is reached).

Here is a scrappy diagram of the process that I made, which is how I initially understood how it works. The top left triangle (green) is computed for the critical orbits, with the aim to calculate the coefficients of the inverse of the bottom diagonal. Then the red path is computed per pixel. The subdiagram at the right shows the continuity checking process.

For colouring filaments with distance estimation I use dual complex numbers for automatic differentiation, for colouring interior I use a function of the final $w$ to adjust hue. To keep images stable for animations, the same number of total iterations is needed for the interior pixels each frame.