Is there an algorithm for second order filtering of data where the filter frequency is not known a priori? It is well known that a moving average algorithm done in the time domain is equivalent to a filter with frequency response $\mathrm{sinc}(\omega\tau)$ where $\tau$ is the averaging time. (see this related answer)

This has the following beneficial property: you are streaming a time series of data ${x_n}$, and the average at any point ($a_n$) is just:
$$ a_n = a_{n-1} \frac{n-1}{n} + \frac{x_n}{n} .$$

Thus you may apply the above recursive algorithm for an arbitrary amount of time ($\tau$), and when you stop, the value you have is filtered by $\mathrm{sinc}(\omega\tau)$, and has a correspondingly reduced variance. Now the $\mathrm{sinc}$ function is a first order low pass, modulated by a $\sin$ envelope. So in effect you have done a first order low pass where the characteristic low pass time constant $\tau$ is equal to the length of the data stream, and $\tau$ was not necessarily known before you started.

My question is: is there some analogous procedure which allows for an (approximate) second order low pass where the time constant is not known a priori? 

A possibility is to "average the averages" but that requires keeping all the averages in memory. Is there some law preventing such a procedure with small memory requirements?
The formula $$ a_n = a_{n-1} \frac{n-1}{n} + \frac{x_n}{n} .$$
is a way to compute the average of $n$ samples, and update this average whenever a new sample is obtained, without storing all of the samples individually.

Ignoring the part about online computation, I think your question is equivalent to "is there a way to estimate the mean of N identically distributed random variables that has lower standard deviation than the arithmetic average?"  If so, the answer is: no, there isn't.

I would also worry about the numerical behavior of this formula for large $n$.