Possible Bug or just my usecase? Memory growth out of bounds

Hello!
I found this rapid memory growth in my simulation but i am not quite sure if its a bug of missing garbage collection or just my usecase that calculates unused results.
Recreation and diagnostic in small a example was done with Fable 5, as by no means i would consider myself a software engineer.

I am calling ElectrodeSOHSolver.solve() repeatedly in my fine interval (dt=1.0s), long-running (Months) simulation. Memory grows about 220KB per call and is not reclaimed from garbage collection.

Memory growth seems to caused by interacting ElectrodeSOHSolver and parameter_values.evaluate(). As ElectrodeSOHSolver unconditionally computes theoretical_energy_integral()and builds brand-new symbolic expression on every call: np.linspace(..., num=1000) stoichiometry arrays become fresh pybamm.Vector symbols with new ids each time.

parameter_values.evaluate() routes through ParameterSubstitutor.process_symbol(),
which stores every processed symbol in self._cache.The dict is unbounded and
keyed by symbol id, so these fresh-every-call entries are never hit again and never
evicted. Both the original and the processed expression trees (each holding
1000-element arrays) are pinned forever.

In my specific usecase where i step a DFN-Model cell at 1s resolution and eSOH at every step it produces around 95MB per 300 steps. This would consume my 32GBs of RAM fast if i simulate my actual wanted time period.

As stated im no expert and therefore not sure if its acually a bug, would be thankful for any opinions on this matter.
Thank you!

Reproduction in small example:

import tracemalloc
import pybamm

parameter_values = pybamm.ParameterValues("OKane2022")
param = pybamm.LithiumIonParameters()

Q_n = parameter_values.evaluate(param.n.Q_init)
Q_p = parameter_values.evaluate(param.p.Q_init)
Q_Li = parameter_values.evaluate(param.Q_Li_particles_init)

esoh_solver = pybamm.lithium_ion.ElectrodeSOHSolver(parameter_values)

tracemalloc.start()
for i in range(1, 101):
    # slight decrease of Q_Li
    inputs = {"Q_n": Q_n, "Q_p": Q_p, "Q_Li": Q_Li * (1 - 1e-6 * i),
              "V_min": 2.5, "V_max": 4.2}
    esoh_solver.solve(inputs)
    if i % 20 == 0:
        cache_len = len(parameter_values._processor.cache)
        heap = sum(s.size for s in tracemalloc.take_snapshot().statistics("filename"))
        print(f"solve #{i:3d}: cache entries = {cache_len:5d}, "
              f"python heap = {heap/1024/1024:6.1f} MB")

# demonstrating the cache is the sole retainer:
parameter_values._processor.clear_cache()
heap = sum(s.size for s in tracemalloc.take_snapshot().statistics("filename"))
print(f"after clear_cache(): heap = {heap/1024/1024:.1f} MB")

Workaround right now is: parameter_values._processor.clear_cache() after each ElectrodeSOHSolver.solve()