Dividing by extremely small numbers in expectation value

Hi, I hope this isn’t a duplicate question.
I’m trying to calculate quantities of the form

\frac{\bra{\psi}e^{-H\tau_2}\hat{n}_xe^{-H\tau_1}\ket{\psi}}{\bra{\psi}e^{-H(\tau_2+\tau_1)}\ket\psi}

by using TDVP to evolve the initial nonequilibrium state \ket{\psi} in imaginary time. My naive approach is to obtain a vector of (normalized) wavefunctions at each time step and evaluate the inner product:

inner(npsis[t2],apply(op("n",s,x),npsis[t1]))
/inner(npsis[t2],npsis[t1])

Which appears to work as expected, but for large time intervals the overlap between states is vanishingly small on the top and bottom and I get NaN results. If instead I don’t normalize the wavefunctions at each step I quickly hit the float64 limit and time evolution fails to converge. Is there a more direct way of calculating these sort of correlators over a partition function, similar to expect but when \tau_1 \neq \tau_2, or failing that some workaround to handle states with very large/small norms?
Thanks!

Just in case it wasn’t clear, npsis is a vector of matrix product states, and I’m working with fermionic site types but the issue isn’t exclusive to them. I’ve also tried loginner everywhere I can but that still returns -Inf in places where it shouldn’t.