It is possible to render 3D fractals by ray-marching with distance estimators. The distance estimate tells you, given a point, how far the nearest point in the fractal is. The ray-marching algorithm calculates the distance from a point on the ray (in any direction, not necessarily the ray direction), and steps the ray forward by this amount, thus the ray gets closer and closer to the object (or passes by it and hits the background).
Syntopia has a blog with a series of posts on Distance Estimated 3D Fractals, applicable within the software Fragmentarium and its updated fork FragM. FragM 2.0.0 has a Complex.frag with complex dual numbers for automatic differentiation. Using FragM's DE Raytracer, a "Multibrot Stack" distance estimator can be implemented like this:
uniform float Iterations; // e.g. 100
uniform float EscapeRadius; // e.g. 30
uniform float DistanceFactor; // e.g. 0.5
float DE(vec3 position)
{
vec4 c = vec4(position.xy, 1.0, 0.0);
float d = position.z;
vec4 z = c;
int i = 0;
float r = length(z.xy); // r = |z|
while(r < Bailout && (i < Iterations))
{
z = cExp(cLog(z) * d) + c; // z = z^d + c
r = length(z.xy); // r = |z|
i++;
}
float dr = length(z.zw); // dr = |dz/dc|
// compute scaled distance estimate
return 0.5 * log(r) * r / dr * DistanceFactor;
}
Note the use of the DistanceFactor variable to scale the distance estimate. This is to try to mitigate the inherent problem with this approach: the distance estimate is only valid in the plane where the power $d$ is constant, but the ray direction in general has varying $d$. Small changes in $d$ might lead to large changes in the distance estimate, but hopefully the change is bounded above by a constant factor.
There is another problem: non-integer powers of complex numbers are multivalued, because complex $\log$ is multi-valued (the imaginary part can have arbitrary integer multiples of $2 \pi$ added to it). This leads to the "cut" that appears along the negative X axis. A proper rendering of the fractal (even in 2D with constant $d$) should better take into account the multiple values, but I don't know the best way to do this.
How does it look? This animation took about 45 minutes to render at an acceptable quality (with DistanceFactor = 0.5):
