Cartopy has a demo that addresses this issue here: https://scitools.org.uk/cartopy/docs/v0.15/examples/always_circular_stereo.html?highlight=set_extent

Basically, make a clip path around the border of your map. The clip path is defined underneath your call to generate the figure, and there are two set_boundary calls for the maps with the limited extents.

The output (the automatic gridlines are a little funky but you can always make your own):

Here's your modified code:

from cartopy import crs
from math import pi as PI
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import numpy as np
import matplotlib.path as mpath

CEL_SPHERE = crs.Globe(
    ellipse=None,
    semimajor_axis=180/PI,
    semiminor_axis=180/PI,
)
PC_GALACTIC = crs.PlateCarree(globe=CEL_SPHERE)

def render_map(path, width, height):
    fig = plt.figure(layout="constrained", figsize=(width, height))
    theta = np.linspace(0, 2*np.pi, 100)
    center, radius = [0.5, 0.5], 0.5
    verts = np.vstack([np.sin(theta), np.cos(theta)]).T
    circle = mpath.Path(verts * radius + center)
    try:
        gs = GridSpec(2, 2, figure=fig)

        axN1 = fig.add_subplot(
            gs[0, 0],
            projection=crs.AzimuthalEquidistant(
                central_latitude=90,
                globe=CEL_SPHERE,
            )
        )
        axN1.gridlines(draw_labels=True)

        axS2 = fig.add_subplot(
            gs[0, 1],
            projection=crs.SouthPolarStereo(globe=CEL_SPHERE)
        )
        axS2.gridlines(draw_labels=True)

        axN2 = fig.add_subplot(
            gs[1, 0],
            projection=crs.AzimuthalEquidistant(
                central_latitude=90,
                globe=CEL_SPHERE,
            )
        )
        axN2.set_extent((-180, 180, 70, 90), crs=PC_GALACTIC)
        axN2.gridlines(draw_labels=True)
        axN2.set_boundary(circle, transform=axN2.transAxes)

        axS2 = fig.add_subplot(
            gs[1, 1],
            projection=crs.SouthPolarStereo(globe=CEL_SPHERE)
        )
        axS2.set_extent((-180, 180, -90, -70), crs=PC_GALACTIC)
        axS2.gridlines(draw_labels=True)
        axS2.set_boundary(circle, transform=axS2.transAxes)

        fig.savefig(path)
    finally:
        plt.close(fig)

if __name__ == "__main__":
    render_map("map_test.pdf", 12, 12)