Consider the following observations:
If $A$ and $B$ are bounded operators on a Banach space over $\mathbb{C}$, then the Lie product formula implies that $$ e^{it(A+B)} = \lim_{n\to\infty} (e^{itA/n}e^{itB/n})^n. \tag{1} $$ This then allows to extract $A+B$ from the limit in the right-hand side: $$ A + B = \frac{1}{i}\left( \frac{\mathrm{d}}{\mathrm{d}t}\biggr|_{t=0} \right) \lim_{n\to\infty} (e^{itA/n}e^{itB/n})^n. \tag{2} $$
Now, let $A$ and $B$ be self-adjoint operators defined on dense subsets $\mathrm{Dom}(A)$ and $\mathrm{Dom}(B)$ of a Hilbert space $\mathcal{H}$, respectively. Then a version of Trotter product formula tells that, if $\mathrm{Dom}(A) \cap \mathrm{Dom}(B)$ is dense and $A+B$ is essentially self-adjoint on $\mathrm{Dom}(A) \cap \mathrm{Dom}(B)$, then the equality $\text{(1)}$ holds in the strong operator topology.
These led me to the following question:
Question. Is it possible that $A+B$, as an operator defined on $\mathrm{Dom}(A) \cap \mathrm{Dom}(B)$, is not essentially self-adjoint, but $A+B$ can still be defined as a self-adjoint operator via the formula $\text{(2)}$?
If the answer is negative, does it mean that the RHS of $\text{(1)}$ diverges (or does not define a strongly continuous one-parameter unitary group) whenever $A+B$ is not essentially self-adjoint?
The formula $\text{(1)}$ is less affected from the "densely-definedness" issue of self-adjoint operators. So I am curious whether this small advantage actually makes any differences and possibly opens up more rooms for defining $A+B$ or not.
Unfortunately, I have only tangential knowledge on functional analysis and operator theory, which means that I don't have any good intuition regarding this question. I am fairly certain that this kind of question has been already examined by many people, so the fact that I can't seem to find any reference is likely because there are some serious obstructions in this approach (or I am too clumsy at googling ).
Some Attempts. I tried this question with $\mathcal{H} = L^2[0, 1]$, $A = -\frac{\mathrm{d}^2}{\mathrm{d}x^2}$, and $B = -\frac{\mathrm{d}^2}{\mathrm{d}x^2}$ with the respective domains $$\mathrm{Dom}(A) = \{f \in H^2[0, 1] : f(0) = f(1) = 0\}$$ and $$\mathrm{Dom}(B) = \{f \in H^2[0, 1] : f'(0) = f'(1) = 0\}.$$ Then $A + B$ defined on the intersection of their respective domains is not essentially self-adjoint, since both $2A$ and $2B$ are self-adjoint extensions of $A+B$. I tried simulating the RHS of $\text{(1)}$ numerically, and the results seem suggesting that the limit is convergent. However, I am not sure whether this is an artifact of the numerical scheme I adopted or not.