If you'd like to avoid geoplot or are experiencing crashing, you should be able to also get the desired effect using seaborn and geopandas (though not with just geopandas to my knowledge).

with data like this:

Polygon GeoDataFrame:
          id                                           geometry
0  polygon_1  POLYGON ((-101.7 21.1, -101.6 21.1, -101.6 21....

Points GeoDataFrame:
   PointID                     geometry
0        0  POINT (-101.61326 21.14453)
1        1  POINT (-101.66483 21.18465)
2        2  POINT (-101.61764 21.11646)
3        3  POINT (-101.65355 21.12132)
4        4  POINT (-101.68183 21.17071)
5        5  POINT (-101.61088 21.14948)
6        6  POINT (-101.66336 21.17007)
7        7   POINT (-101.6774 21.14027)
8        8  POINT (-101.66757 21.13169)
9        9  POINT (-101.66333 21.12997)

in any crs:

fig, ax = plt.subplots(figsize=(3, 3))
assert polygon_gdf.crs == points_gdf.crs
print("crs is", polygon_gdf.crs.coordinate_system)
polygon_gdf.plot(ax=ax, facecolor='none')
points_gdf.plot(ax=ax)
ax.set_axis_off()

crs is ellipsoidal

you can just re-project it to a cartesian coordinate reference system, like Albers Equal Area:

polygon_gdf.crs = 'EPSG:9822'
points_gdf.crs = 'EPSG:9822'
print("crs is", polygon_gdf.crs.coordinate_system)

crs is cartesian

Then you can plot as follows:

import seaborn as sns
fig, ax = plt.subplots(figsize=(3, 3))
sns.kdeplot(
    x=points_gdf.geometry.x,
    y=points_gdf.geometry.y,
    fill=False,
    cmap="viridis",
    bw_adjust=0.5,
    thresh=0.05,
    ax=ax,
)
polygon_gdf.plot(ax=ax, facecolor="none")
ax.set(xlabel="", ylabel="")
fig.tight_layout()