In Lorenzo Bergomi, Stochastic Volatility Modeling, Chapter 5 Appendix A.1, Equation (5.64), as shown below, seems to assume $\hat\sigma$ to be constant. If that is the case, why do we bother to invoke the Feynman-Kac formula? We can simply use the Black-Scholes formula.
Just want to add the observation that the pricing PDE solution can be formally written as $$ C(\tau) = e^{\tau \mathcal H} C(0) \quad (*) $$ where $\tau$ is time to maturity and $\mathcal H$ is a differential operator. For example, in the BS world with zero interest rate it is $$ \mathcal H = \tfrac12 \sigma^2 S^2 \frac{\partial^2}{\partial S^2} $$ Thus $U(\tau) = e^{\tau \mathcal H}$ is an 'evolution operator'.
In the BS case $\sigma$ is not a variable but a parameter. So you can differentiate both sides of equation (*) to very quickly obtain the vega gamma relation by noting that the operator $U(\tau)$ depends on the parameter $\sigma$. You could reinsert dividends and rates to also obtain sensitivities to $r$ and $q$ in the same manner.
In stochastic volatility models you can similarly show that the sensitivity of the option price to the correlation parameter is the stochastic volatility model vanna.
Even though it is true that the volatility is constant in this setting, the relationship is valid for all terminal condition or pay-off function -- beyond the typical $(\pm(S-K))_+$ -- so long as the pay-off function is independent of the volatility. We can certainly write out the integral expressions of vega and gamma (of arbitrary pay-off functions) and find their relationship. But it seems simpler dealing with the PDE directly. Moreover, this methodology can be used to find other high order partial derivatives.
