How do I improve my Buddhabrot rendering?

Question

I am trying to create a beautiful render of the Buddhabrot fractal with a mathematical background and I want to make it faster and more importantly, prettier. I have a script (buddhabrot_functions) and a main program. It basically splits the image into chunks, which I then later glue together. I have found the code to be rather slow, is there some way that I am missing to make it faster? Also, when I'm running more chunks there seems to be large yellow squares towards the outside of the image and rather noticable inconsistencies at the borders of the chunks.

I have a script (buddhabrot_functions) and a main program. It basically splits the image into chunks, which I then later glue together. I have found the code to be rather slow, is there some way that I am missing to make it faster? Also, when I'm running more chunks there seems to be large yellow squares towards the outside of the image and rather noticable inconsistencies at the borders of the chunks.

This is the script:

def buddhabrot(c, max_iter):
    z = c
    trajectory = [c]
    for i in range(max_iter):
        if abs(z) > 4.0:
            return trajectory
        z = z*z + c
        trajectory.append(z)
    return []

def buddhabrot_helper(inputs):
    c, max_iter = inputs
    return buddhabrot(c, max_iter)

and this is the main program:

from multiprocessing import Pool
from PIL import Image
from buddhabrot_functions import buddhabrot_helper
import numpy as np
from matplotlib import cm


RESOLUTION =10000
xmin, xmax, ymin, ymax = -2.0, 1.0, -1.5, 1.5
width = xmax - xmin
height = ymax - ymin
CHUNK_SIZE = 5000
max_iter = 100000

def calculate_and_save_chunk(i, j):
    xmin_chunk = xmin + i * width * CHUNK_SIZE / RESOLUTION
    xmax_chunk = xmin + (i + 1) * width * CHUNK_SIZE / RESOLUTION
    ymin_chunk = ymin + j * height * CHUNK_SIZE / RESOLUTION
    ymax_chunk = ymin + (j + 1) * height * CHUNK_SIZE / RESOLUTION

    r1 = np.linspace(xmin_chunk, xmax_chunk, CHUNK_SIZE)
    r2 = np.linspace(ymin_chunk, ymax_chunk, CHUNK_SIZE)
    meshgrid = np.meshgrid(r1, r2)

    pixels = [(complex(r, i), max_iter) for r, i in zip(np.ravel(meshgrid[0]), np.ravel(meshgrid[1]))]
    
    pool = Pool()
    trajectories = pool.map(buddhabrot_helper, pixels)
    pool.close()

    # Initialize density map
    density = np.zeros((CHUNK_SIZE, CHUNK_SIZE))

    for trajectory in trajectories:
        for z in trajectory:
            x, y = (z.real - xmin_chunk) * CHUNK_SIZE / (xmax_chunk - xmin_chunk), (z.imag - ymin_chunk) * CHUNK_SIZE / (ymax_chunk - ymin_chunk)
            if 0 <= x < CHUNK_SIZE and 0 <= y < CHUNK_SIZE:
                density[int(x), int(y)] += 1

    density = np.sqrt(density) 

    density = (density - np.min(density)) / (np.max(density) - np.min(density))
    
    color_data = cm.viridis(density)
    img = Image.fromarray((255 * color_data).astype(np.uint8))

    
    img.save(f"buddhabrot_{i}_{j}.png")


for i in range(RESOLUTION // CHUNK_SIZE):
    for j in range(RESOLUTION // CHUNK_SIZE):
        calculate_and_save_chunk(i, j)

I definitely want to stick to this project and Python, can someone help me out? Also, I can not use GPU for my code.