Problem. Be $f\in \mathcal{C}(\mathbf{R})$ and fix $c_0\in\mathbf{R}$ and $c_2<0$. Consider the integral equation $$ f(y) = c_0 + c_2 y^2 + \int dz\, f(z) \exp(f(y-z)) $$ where the integral is over $(-\infty,\infty)$. Under some conditions on $c_0$ and $c_2$, it is possible to find $a_0$ and $a_2$ such that $f(y)=a_0+a_2 y^2$ is a solution. The problem is: Is this solution unique?
Let's be free to choose specific (nontrivial) values for $c_0$ and $c_2$ if that simplifies the analysis; in particular though, let's keep $c_2\neq 0$. Also, let's be free to assume that $f$ is smooth, e.g. $f\in \mathcal{C}^2(\mathbf{R})$.
Roughly speaking, the problem is how to prove that there is only one solution to an equation like $$ f * \exp f(y) = c_2 y^2 + f(y). $$
Attempt 1. From the form of the equation, it is natural to define the operator $$ F_r(f)(y) \triangleq c_0 + c_2 y^2 + \int_{-r}^r dz\, f(z) \exp(f(y-z)) $$ and regard the solution as a fixed point of $F_r$ for any fixed $r>0$. So my first attempt would be to apply Banach fixed point theorem to $F_r$. I don't know if this works at all, even in principle, since $f\in \mathcal{C}(\mathbf{R})$, but for now this is what I have: \begin{align} \| F_r(f)-F_r(g) \|_\infty & \leq \sup_{y\in\mathcal{B}_r(0)} \int_{-r}^r dz \left| f(z)\exp(f(y-z)) - g(z)\exp(g(y-z)) \right| \end{align} where $\mathcal{B}_r(0)$ is the ball of radius $r$ and center the origin. From here, I should try to end up with something like $C\|f-g\|_\infty$ for some $C>0$...
Attempt 2. A different approach is to show that, if $f^*$ and $g^*$ are two solutions, i.e., $F(f^*)=f^*$ and $F(g^*)=g^*$, then $f^*=g^*$. I tried the following: Since by assumption $f^*$ and $g^*$ are two solutions, then $$ f^*(y) - g^*(y) = F(f^*)(y) - F(g^*)(y) = \int dz\, \big[f^*(y-z) \exp(f^*(z))-g^*(y-z) \exp(g^*(z))\big]. $$ I should show that the above quantity is identically equal to zero (for all $y$). I would proceed by contradiction. Suppose that $f^*\neq g^*$, that is, there exists $h=f^*-g^*$ that is not identically zero. Suppose that $h(y)\neq 0$ in a neighborhood of $y=y_0$. Then I try to expand $h$ in Taylor series around $y=y_0$. For any $y\neq y_0$ it results \begin{multline} 0 \neq h(y_0) + h'(y_0)(y-y_0) + \cdots \\ = \int dz\, \bigg[ ( f^*(y_0-z) + (f^*)'(y_0-z)(y-y_0+z) + \cdots) \exp(f^*(z)) \\ - ( g^*(y_0-z) + (g^*)'(y_0-z)(y-y_0+z) + \cdots) \exp(g^*(z)) ) \bigg]. \end{multline} In the end, I should lead to some contradiction, e.g. the integrand is zero...
Attempt 3. I don't know if it is just a coincidence or not, but the solution $f(y)=a_0 + a_2 y^2$ has the same form of the term $c_0+c_2 y^2$. In other words, if we define the operator
$$ G(f) \triangleq f * \exp f - f $$
then the problem is to find $f$ such that $$ G(f)(y) = c_0 + c_2 y^2 .$$
Therefore, an approach may be to find the solution to a “Green equation” like $$G(\varphi)(y)=\delta_y$$ where $\delta_y$ represents a Dirac mass at $y$. Such a $\varphi$ would (?) unlock the possibility of solving the “inhomogeneous” case.
If $f(y) = c_0 + c_2 y^2 + \int f(z) e^{f(y-z)}\mathrm{d}z$ differentiating three times gives $f'''(y)=\int f'''(z)e^{f(y-z)}\mathrm{d}z$. Thus, $f'''=f'''*e^f= f'''*\underbrace{e^f \ast \cdots \ast e^f}_n$, so that $f'''(y)=0$. (Maybe a bit more here.) Higher derivtives are also zero so that $f$ is quadratic. Let $f=a_0+a_2(y-\mu)^2$. If $f(y) = c_0 + c_2 y^2 + \int f(z) e^{f(y-z)}\mathrm{d}z$, then \begin{align*} -a_2\mu& = \frac12 f'(0) = a_2\int(z-\mu) e^{a_0+a_2(z-\mu)^2}\mathrm{d}z %% = 0 \implies \mu=0\\[3mm] %% a_0 & = f(0) = c_0 + %% e^{a_0}\int (a_0+a_2z^2) e^{a_2z^2}\mathrm{d}z\\ %% & = c_0 + a_0e^{a_0}\sqrt{-4\pi a_2}+ %% a_2e^{a_0}\sqrt{-4\pi a_2}/(-2a_2) \\ %% & = c_0 + 2a_0e^{a_0}\sqrt{-\pi a_2} %% -e^{a_0}\sqrt{-\pi a_2 }\\[3mm] %% a_2& = \frac12 f''(0) = c_2 + a_2\int e^{a_0+a_2z^2}\mathrm{d}z %% = c_2 + 2a_2e^{a_0}\sqrt{-\pi a_2}. %% \end{align*} Thus the only solution is $f(y)=a_0+a_2y^2$, where $a_0$ and $a_2$ give the unique (and here) solution to \begin{align*} a_0& = c_0 + 2a_0e^{a_0}\sqrt{-\pi a_2} %% -e^{a_0}\sqrt{-\pi a_2 }\\ %% a_2&=c_2 + 2a_2e^{a_0}\sqrt{-\pi a_2}. %% \end{align*}