ordered exponential of unbounded operators Let $H$ be a Hilbert space,
and let $A_t$ be a family of unbounded positive (self-adjoint) operators on $H$ parametrized by $\mathbb t\in R_{\ge 0}$. Consider the ordinary differential equation
$$
\qquad\qquad\qquad\qquad\qquad\qquad\qquad\frac{d}{dt} E_t = -A_tE_t
\quad\qquad\qquad\qquad\qquad\qquad\qquad(1)
$$
that defines the ordered exponential of the family $A_t$.
If $A_t=A$ is independ of $t$, then the solution of the ODE is the usual exponential $E_t=e^{-tA}$.
Note that the above operators $E_t$ are bounded.
I suspect that, if I put appropriate hypotheses on $A_t$, (such as having a common dense domain, depending continuously on $t$, whatever that might mean, etc.) the solution of (1) will also be bounded. Intuitively, it's kind of clear:
$$
E_t = \lim_{N\to\infty} \Big(e^{-\frac t N  A_t} \cdot e^{-\frac t N A_{t(1-1 /N)}} \cdot e^{-\frac t N A_{t(1-2 /N)}}\cdots  \cdot e^{-\frac t N A_{t(3 /N)}} \cdot e^{-\frac t N A_{t(2 /N)}} \cdot e^{-\frac t N A_{t(1 /N)}}\Big)
$$
and each little exponential in the above product has norm $\le 1$.
Q:
Which properties should I impose on the family $A_t$ in order for the solution of (1) to be well defined and bounded?
 A: There are zillions of articles and books on this topic, for more and more general families of operators, where usually you can choose a specific trade-off between regularity of solutions, variability of domains of the operators, and smoothness of the dependence on time.
In the "variational" setting you consider, it is probably most efficient to introduce a weak formulation based on a family of quadratic forms. 
In this case, there is a classical theory by Lions (see e.g. Chapt. 3 of Equations Differentielles Opérationelles et Problèmes aux Limites, Springer 1961), which essentially says that if the family of forms is equi-continuous and equi-coercive, and the dependence on time is merely measurable, then you have well-posedness (in a certain weak sense) of your equation (1) (even if you add an inhomogeneous term in the equation), as well as boundedness. Also Chapt. 4 of Tanabe's Equations of Evolutions (Pitman 1979) is a good reference.
A: This is usually called the magnus expansion method and has a nice literature in numerical analysis. Kato also used this method to show the existence of solutions in the hyperbolic case. 
I would say that strong resolvent continuity and a sufficiently big common domain is sufficient in your case. See Section 5.3 in Pazy.
I can also give a related self-reference, where also the investigation of this product appears and some of the ideas are explained in a simpler situation.
ADDED: My answer concentrates on the method you propose to converge to the solution. To make the content of the references short: yes. A common dense domain and continuity of hte maps $t\mapsto A_tx$ implies the convergence of the product to the solution of the differential equation.
A: Certainly, one should pay attention to the domain of the operator. However, the following argument should survive a reasonable assumption. We have
$$
\frac{d}{dt}(E^\ast(t)E(t))=-2  E^\ast(t) A(t) E(t)
$$
and thus for a fixed $u$
$$
\frac{d}{dt}(\Vert E(t)u\Vert^2)=-2\langle A(t)E(t)u, E(t)u(t)\rangle\le 0\Longrightarrow
\Vert E(t)u\Vert^2\le \Vert E(0)u\Vert^2,
$$
implying boundedness of $E(t)$, provided $E(0)$ is bounded. Assuming $A(t)\le -c_0Id$, the same inequality along with Gronwall gives
$$
\Vert E(t)u\Vert^2\le e^{-2c_0t}\Vert E(0)u\Vert^2.
$$
