ISJ bandwidths consistently >~ 2 times smaller than ArviZ's

I was confused why KDEpy's ISJ bandwidth produced such noisy results, especially compared to ArviZ. Turns out its calculation is consistently larger than KDEpy's (reproducing code below):

Here are comparisons for unimodal and bimodal distributions (including a comparison with KDEpy's bandwidth fixed to that calculated by ArviZ):

I can't weigh in on the correctness of one implementation over the other, but in these two examples, ArviZ's ISJ result seems unbiased compared to KDEpy's, just without the wiggliness.

import numpy as np
from scipy.stats import norm
from KDEpy import FFTKDE
from arviz_stats.base import array_stats
import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
ns = 2**np.arange(9, 21)
kdepy_bws = []
az_bws = []
for n in ns:
    data = norm.rvs(size=n)
    kdepy_bws.append(FFTKDE(bw="ISJ").fit(data).bw)
    az_bws.append(array_stats.kde(data, bw="isj")[-1])

axes[0].loglog(ns, kdepy_bws, ".-", label="KDEpy")
axes[0].loglog(ns, az_bws, ".-", label="ArviZ")
axes[0].set_ylabel("bandwidth")
axes[1].semilogx(ns, np.array(az_bws) / np.array(kdepy_bws), ".-")
axes[1].set_ylabel("ArviZ/KDEpy")
for ax in axes.flat:
    ax.set_xlabel("sample size")
axes[0].legend()

def plot_comparisons(ax, data):
    _fit = FFTKDE(bw='ISJ').fit(data)
    ax.plot(*_fit(), label=f"FFTKDE, ISJ bw={_fit.bw:.4g}")
    
    x, y, bw = array_stats.kde(data, bw="isj")
    ax.plot(x, y, "--", label=f"arviz, ISJ bw={bw:.4g}")
    
    _fit = FFTKDE(bw=float(bw)).fit(data)
    ax.plot(*_fit(), ":", label=f"FFTKDE, bw={_fit.bw:.4g}")

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
dist = norm()
n = 10**4
x = np.linspace(-4, 4, 1000)
axes[0].plot(x, dist.pdf(x), color="k", label="truth")
data = dist.rvs(n)
plot_comparisons(axes[0], data)
axes[0].legend(loc="lower center", bbox_to_anchor=(1/2, 1), ncol=2)
axes[0].set_xlim(-3, 3)

dist1 = norm(loc=-2)
dist2 = norm(loc=2)
x = np.linspace(-6, 6, 1000)
axes[1].plot(x, (dist1.pdf(x) + dist2.pdf(x))/2, color="k", label="truth")
n = 10**4
data = np.hstack([dist1.rvs(n), dist2.rvs(n)])
plot_comparisons(axes[1], data)
axes[1].legend(loc="lower center", bbox_to_anchor=(1/2, 1), ncol=2)
axes[1].set_xlim(-5, 5)