Let $\mathfrak{g}$ be a Lie algebra, and let $G$ be the (unique up to isomorphism) simply-connected group with Lie algebra $\mathfrak{g}$. Here, $\mathfrak{g}$ and $G$ should be thought as "abstract", with the exponential map $\exp : \mathfrak{g} \rightarrow G$ defined in terms of the 1-parameter sub-groups in $G$. The BCH formula is fulfilled for this "abstract" exponential map, although the completly general proof is a bit lengthy. It involves deriving a formula for the derivative of the exponential mapping, then proving $\exp \big( [X,\cdot] \big) Y = \exp (X) Y \exp (-X)$ and finally deducing the BCH formula. A concise presentation can be found eg. in chapter 3 of Brian Hall's "Lie Groups, Lie Algebras, and Representations". This reference aims at being fairly accessible, at the price of restricting itself to those Lie goups $G$ which can be constructed as matrix subgroups: this is the case of most groups used in practice (in particular the Heisenberg group generated by the position/momentum operators of quantum mechanics, together with the identity operator, can be constructed as a subgroup of 3x3 real matrices).
Next, let $X \mapsto i\hat{X}$ be a representation of this Lie algebra by unbounded symmetric operators on a complex Hilbert space $\mathcal{H}$ (note that we do not even need to assume that $\hat{X}$ be essentially self-adjoint, just symmetric).
Since we are dealing with unbounded operators, we have to be a bit careful as to what we mean by a representation, because a priori $\hat{X}\hat{Y}$, $\hat{Y}\hat{X}$, and therefore $[\hat{X}, \hat{Y}]$ may fail to be well-defined. Specifically, we demand that there exists a common invariant dense domain $\mathcal{D} \subseteq \mathcal{H}$, such that, for all $X \in \mathfrak{g}$, $\hat{X}$ is defined on $\mathcal{D}$ and stabilizes $\mathcal{D}$ (ie. $\hat{X} \left\langle \mathcal{D} \right\rangle \subseteq \mathcal{D}$). Then, we can define arbitrary products of the $\hat{X}$'s on $\mathcal{D}$, and in particular commutators.
The question is then whether this representation can be exponentiated, ie:
for all $X \in \mathfrak{g}$, $\hat{X}$ is actually essentially self-adjoint;
and there exists a unitary representation $U \mapsto \hat{U}$ of the group $G$ on $\mathcal{H}$ such that, for any $X \in \mathfrak{g}$, $\widehat{\exp X} = \widehat{\exp} (i\hat{X})$. Here, I denote by $\widehat{\exp}$ the "operator" exponential map, which is defined by spectral resolution of the essentially self-adjoint operator $\hat{X}$: as you observed, in the case of unbounded operators, the exponential cannot be defined in terms of the exponential series (note that the notation $\widehat{\exp}$ is non-standard, I use it here only to prevent confusion between the various notions of exponential maps).
If this holds, the BCH formula satisfied by the "abstract" exponential map $\exp$ will be inherited by the "operator" exponential map $\widehat{\exp}$.
So, with all these preliminaries in place: when can a Lie algebra representation by symmetric unbounded operators be exponentiated? A sufficient condition (useful in practice, albeit not a necessary condition) is the Nelson criterion (lemma 9.1 of Edward Nelson, "Analytic Vectors", Annals of Mathematics, Second Series, Vol. 70, No. 3 (Nov., 1959), pp. 572-615): it demands that there exists a basis $X_1,\dots,X_n$ of $\mathfrak{g}$, a dense subdomain $\mathcal{D}_o \subseteq \mathcal{D}$, and a real $s > 0$, such that:
$$
\forall \psi \in \mathcal{D}_o,\; \sum_{m=0}^{\infty} \frac{s^m}{m!} \sum_{k_1,\dots,k_m} \left\| X_{k_1} \dots X_{k_m} \psi \right\| < \infty.
$$
This is a fairly technical result, but the intuition behind it is that this condition is precisely what you need to define $\widehat{\exp} (i\hat{X})$ over $\mathcal{D}_o$ via the exponential series, and, from there, deduce that $\hat{X}$ is indeed essentially self-adjoint (using the Stone's theorem discussed below), that this definition of $\widehat{\exp} (i\hat{X})$ coincides over $\mathcal{D}_o$ with the spectral one, and finally, that this indeed gives you a unitary representation of $G$ (matching the "abstract" BCH-formula with the one that can be proven directly using the exponential series).
In the case of the position/momentum operators of quantum mechanics, one can for example take $X_1 = \text{id}$, $X_2 = q$, $X_3 = p$ and take $\mathcal{D} = \mathcal{D}_o$ to be spanned by finite linear combination of the harmonic osciallator energy eigenstates. Then, the condition can be proven using the expression of the position/momentum operators in terms of ladder operators.
However, coming back to your specific motivation, I do not honestly think that you need all this machinery. Instead, you can use the explicit expression of the Weyl operators to prove that $t \mapsto W(tz)$ is a strongly continuous one-parameter unitary group: ie. $W(sz) W(tz) = W\big((s+t)z\big)$ and, for any $\psi \in \mathcal{H}$, $t \mapsto W(tz) \psi$ is continuous (with respect to the norm of $\mathcal{H}$; in this case, the $L_2$-norm). Then, Stone's theorem (theorem VIII.8 of Reed and Simon, "Methods of Modern Mathematical Physics", volume 1) tells you that $W(tz) = \widehat{\exp} (it\hat{X})$ with $\hat{X}$ the self-adjoint operator defined by:
$$i \hat{X} \psi = \left. \frac{d}{dt} W(tz) \psi \right|_{t=0}$$
for any $\psi$ in the dense domain $\mathcal{D}$ where this derivative exists. Using again the explicit expression of $W(tz)$ you can then check that $\hat{X} = \sqrt{2} (y\hat{q} - x\hat{p})$.
I think you are almost certainly misreading the paper to reverse the flow of logic involved.
You simply evaluate the lowest order CBH expansion mindful of the omitted orders,
$$e^{H_1t/n}e^{H_2t/n}= e^{(H_1+H_2)(t/n) +\frac{1}{2} [H_1,H_2](t/n)^2 + O((t/n)^3)} \\
e^{-H_1t/n}e^{-H_2t/n}= e^{-(H_1+H_2)(t/n) +\frac{1}{2} [H_1,H_2](t/n)^2 + O((t/n)^3)} \\
f(t)\equiv (e^{H_1t/n}e^{H_2t/n}e^{-H_1t/n}e^{-H_2t/n})=
e^{[H_1,H_2](t/n)^2 + O((t/n)^3)}.$$
(Observe the $O(t/n)$ terms mutually commute, being opposite, so they do not contribute to $O((t/n)^2)$.)
It is then evident that
$$
f(t) f(-t)= e^{2[H_1,H_2]t^2/n^2 + O((t/n)^4)},
$$
since the $O((t/n)^3$ odd terms cancel each other in the sum, leaving $O((t/n)^4$.
I'm not sure how you concluded your mystery author needed the first equation to derive this.
(Note, instead, that
$
f(t)^m= e^{m[H_1,H_2](t/n)^2 + O((t/n)^3)}
$, but glibly taking $m=n^2/t$ will net an ill-defined error. I suspect there is essential information in your text you omitted.)
Best Answer
I think this should at least be true on any connected neighborhood of $0$ on which the series converges by analyticity: Lie groups are automatically anaytic manifolds, and both multiplication and the exponential map are analytic. Hence both sides of the formula represent analytic functions, with the left hand side defined on $\mathfrak g\times\mathfrak g$ and an the right hand side defines whehere the power series converges. Since these functions agree locally around the origin, uniqueness of analytic continuation shows that they agree on any connected neighborhood of the identity, on which the series converges. I don't know whether there may be other domains in which the series converges and what happens there.