Mitigate floating-point numerical errors for very-low corner frequency low-pass filter with DSP
I think your issue might be coefficient quantization and filter topology. A direct form biquad has poor quantization effects about 0 and π radians. It’s easier to analyze such effects in fixed point, but even though floating point has a much larger range, it still has shortcomings. In particular, if you add a very small number to a very large number, the small number disappears because it can’t be aligned for the operation in the available number of mantissa bits. This can cause the order of operations to affect the result. For instance, where S is a small number and L is large, L - L + S = S, but L + S - L = 0.
Udo Zolzer covers the differences between several filter structures in his book Digital Audio Signal Processing. I borrow direct form quantization effects on pole locations from the book:

See how precision is lost near 0 and π. Other filter topologies might have higher precision near 0, while being much worse near π, which may be a good tradeoff for uses such as yours. The Gold and Rader form has a very even distribution, it looks like a perfect grid.
Another simple and popular filter that has good quantization characteristics at low frequencies is the &quot;Chamberlin&quot; state variable filter. There are improved versions of this filter, as it has problems at higher frequencies (from about one-sixth of the sample rate and up), but the plain Chamberlin is very good at low frequencies, where you need it.
See my article on the Chamberlin state variable filter here:
The digital state variable filter
Zolzer presents modified Chamberlin structures here:
The Modified Chamberlin and Zölzer Filter Structures
In particular, see the graph of quantization effect near zero for the Chamberlin structure—very dense near zero, at the expense of poorer performance at high frequencies, compared to the direct form graph: