Questions about equi-angular sampling

Question

Equiangular sampling is a sampling method for point (sphere) emitter in scattering medium. While it claims to be very efficient when dealing with point (sphere) emitter in scattering medium (see this shadertoy example), I find most of the reference and implementations to have only done about two bounces and for higher order scattering, equiangular sampling fails to distinguish itself from the standard exponential distance sampling:

Equiangular 2 bounces	Exponential 2 bounces	Equiangular 400 bounces	Exponential 400 bounces

As can be seen above, when we have only 2 bounces the equiangular sampling is clearly better but if we do not clip the bounces to just 2 (the number of bounces up-to 400) the result is just ordinary. So why would this happen? Does this mean that equiangular sampling can't actually handle high-order scattering well? Or is there anything wrong with my code? My implementation is embedded in PBRT-v3, simply add a function in media/homogeneous.cpp and pass the position of the point source in integrator/volpath.cpp.

The scene is cornell box scene with point source being LightSource "point" "rgb I" [10. 10. 10.] "point from" [2.779 4.8 2.4] and the parameters for scattering medium is: "rgb sigma_s" [0.4 0.4 0.4] "rgb sigma_a" [0.005 0.005 0.005]. There is a homogeneous medium box enclosing the whole scene. Camera params: LookAt 2.78 2.73 -8.00 2.78 2.73 0 0 1 0 Camera "perspective" "float fov" [39.3077]. The equiangular implementation is simple:

Spectrum HomogeneousMedium::equiangular_sample(
    const Ray &ray, Sampler &sampler, MemoryArena &arena, MediumInteraction *mi, const Point3f& target_p
) const {
    ProfilePhase _(Prof::MediumSample);
    // Sample a channel and distance along the ray
    int channel = std::min((int)(sampler.Get1D() * Spectrum::nSamples),
                           Spectrum::nSamples - 1);
    // Sample a scattering event: surface or medium?
    Float dist = -std::log(1 - sampler.Get1D()) / sigma_t[channel];
    Float d_length = ray.d.Length(), t = std::min(dist / d_length, ray.tMax);
    bool sampledMedium = t < ray.tMax;

    Float theta_a = 0, theta_b = 0, D = 0, sampled_t = 0, sample_pdf = 1.0;
    if (sampledMedium) {
        // equiangular sampling: overriding the sampled `t`
        Vector3f diff_vec = target_p - ray.o;
        Float dot_val = Dot(diff_vec, ray.d);
        Float proj_len = dot_val / d_length;
        D = sqrtf(diff_vec.LengthSquared() - proj_len * proj_len);
        
        theta_a = atan2f(1e-4 - proj_len, D);
        theta_b = atan2f(ray.tMax - proj_len, D);
        Float sample = sampler.Get1D();
        sampled_t = D * tanf(theta_a * (1 - sample) + sample * theta_b);
        t = (sampled_t + proj_len) / d_length;
        
        sample_pdf = D / fabsf(theta_b - theta_a) / (D * D + sampled_t * sampled_t);
        *mi = allocate_it(ray, arena, t);
    }
    // Compute the transmittance and sampling density
    // Note that we can not modify the sigma_t here, since it's calculating contribution
    Float actual_d = std::min(t, MaxFloat) * d_length;
    Float wall_pdf = expf(-sigma_t[channel] * std::min(ray.tMax, MaxFloat) * d_length);
    Spectrum Tr = Exp(-sigma_t * actual_d);
    
    Float pdf = (1.f - wall_pdf) * sample_pdf;
    return sampledMedium ? (Tr * sigma_s / pdf) : (Tr / wall_pdf);
}