Show that $M_X(t) = 2 E \left[ e^{tX} \Phi( aX-t) \right], \forall t \in \mathbb{R}$ iff $X$ is Gaussian Let $M_X(t)$ denote the moment generating function of a random variable $X$.  Now suppose that the following expression holds: for a given $a>0$
\begin{align}
M_X(t) = 2 E \left[ e^{tX} \Phi( aX-t) \right], \forall t \in \mathbb{R}
\end{align}
where $\Phi$ is cdf of the standard normal.
Can we show that the only random variable that satisfies the above is Gaussian? The forward direction (i.e., evaluating with Gaussian) is not very difficult to show.  However, the backward direction has been challenging.
Proof of the direct part:  Choose $X$ to be zero mean Gaussian where variance $\sigma^2$ will be selected a bit later. Then, let $Z$ be standard normal, and note
\begin{align}
2 E \left[ e^{tX} \Phi( aX-t) \right]
&= 2 E \left[ e^{tX}  1_{\{Z \le aX-t\}} \right]\\
&= 2 E \left[ E[ e^{tX}  1_{\{Z \le aX-t\}}|Z] \right]\\
&= 2 E \left[ E \left[ e^{tX}  |Z, \{Z \le aX-t\}\right] P[Z \le aX-t|Z] \right]\\
\end{align}
Now, the following is just the moment generating function of a truncated Gaussian
\begin{align}
E \left[ e^{tX}  |Z, \{Z \le aX-t\}\right]= e^{\sigma^2\frac{t^2}{2}} \frac{1-\Phi(\frac{Z+t}{a\sigma}-\sigma t)}{1-\Phi(\frac{Z+t}{a\sigma})}
\end{align}
Also note that  $P[Z \le aX-t|Z]=1-\Phi(\frac{Z+t}{a\sigma})$. Therefore, we arrive at
\begin{align}
2 E \left[ e^{tX} \Phi( aX-t) \right]&=2 E \left[ e^{\sigma^2\frac{t^2}{2}} \left(1-\Phi(\frac{Z+t}{a\sigma}-\sigma t) \right)\right]\\
&=2   e^{\sigma^2\frac{t^2}{2}} E \left[ 1-\Phi \left(\frac{Z+t}{a\sigma}-\sigma t\right)\right]\\
&=2   e^{\frac{t^2}{2a}} E \left[1-\Phi ( \frac{Z}{\sqrt{a}}) \right] \text{  choose } \sigma=\frac{1}{\sqrt{a}}\\
&=2   e^{\frac{t^2}{2a}} \frac{1}{2} \text{ by symmetry}\\
&=   e^{\frac{t^2}{2a}}
\end{align}
Finally, note that the $X$ with $\sigma=\frac{1}{\sqrt{a}}$ has $M_X(t)= e^{\frac{t^2}{2a}}$
 A: Not a complete answer, but too long for a comment.

Let me reformulate$^\ast$ the question in the OP as follows: Consider an integrable, non-negative function $p(x)$, which satisfies the identity
$$\int_{-\infty}^\infty e^{tx}\,\text{erf}\,\bigl((x-t)/\sqrt 2\bigr)\,p(x)\,dx\equiv 0,\;\;\forall t\in\mathbb{R}.\qquad\qquad(\ast)$$
Can we then conclude that $p(x)\propto e^{-x^2/2}$ is a Gaussian? In other words, is the measure such that the convolution of exponential and error function vanishes unique?
If $p_0(x)$ is a solution of $(\ast)$, then also $p_b(x)=e^{bx}p_0(x+b)$ is a solution, for any real $b$. The key point to note is that the Gaussian solution has the unique property that it is invariant under this transformation. If there exists a non-Gaussian solution, then we have a class of independent solutions formed by linear combinations of $p_b(x)$, $b\in\mathbb{R}$.
The argument would be complete if this class of solutions can serve as a basis for any $p(x)$, since the identity $(\ast)$ is definitely false in general.


$^\ast$ To reformulate the identity I used that $\Phi(x)=\tfrac{1}{2}+\tfrac{1}{2}\,\text{erf}(x/\sqrt{2})$. The coefficient $a$ can be set to unity without loss of generality, by rescaling $x$ and $t$. Since the identity is linear in $p(x)$, we need not require that this function is normalized to unity.
