It is shown in this paper that for bilinear systems with a stable drift term the natural space of input-output models is that of Bochner-integrable causal Fock functionals, the latter being a causal version of the states of non-self-interacting quanta fields. Moreover, the functional spline interpolators or smoothers thereby obtained are found to have time-varying but asymptotically stable and decoupled bilinear realizations, with dimension equal to the number of test signals used for model matching. Motivation is given by the need for order reduction and input-output analysis of Carleman bilinearization, but the results apply to general bilinear systems.