From Osculating to Mean Elements with Brouwer-Lyddane Propagator, negative eccentricity

Hi team as discussed in this thread https://forum.orekit.org/t/from-osculating-to-mean-elements-with-brouwer-lyddane-propagator/1867 There is a potential issue of negative eccentricity on the BL propagator. Here the code to reproduce the error:

long start = 1613870730000L;

/* TLEs of the satellite */
TLE tleOrbit = new TLE
            "1 43196U 18015E   21055.59816856  .00000894  00000-0  38966-4 0  9996",
            "2 43196  97.4662 188.8169 0016935 299.6845  60.2706 15.24746686170319");
Propagator propagator = TLEPropagator.selectExtrapolator(tleOrbit);
    
//Get state at initial date and 3 days before
SpacecraftState tleState = propagator.getInitialState();
SpacecraftState tleStateAtDate = propagator.propagate(Time.getAbsoluteDateFromTimestamp(start));
    
//BL mean orbit
double epsilon = 1.0e-12;
int maxIter = 500;
UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(7, 7);
Orbit meanOrbitOK = BrouwerLyddanePropagator.computeMeanOrbit(tleState.getOrbit(), provider,
provider.onDate(tleState.getDate()), BrouwerLyddanePropagator.M2, epsilon, maxIter);
Orbit meanOrbitKO = BrouwerLyddanePropagator.computeMeanOrbit(tleStateAtDate.getOrbit(), provider,
provider.onDate(tleStateAtDate.getDate()), BrouwerLyddanePropagator.M2, epsilon, maxIter);

Thanks for your support! Alberto