In the statistical physics research, I encountered the following inequality that I need to prove: \begin{aligned} & \int_{[0,1]^n} \frac{\partial f(a; x_1, \cdots, x_n)} {\partial a}\mathrm{d}x_1\cdots \mathrm{d}x_n \int_{[0,1]^n} f^2(a; x_1, \cdots, x_n) \mathrm{d}x_1\cdots \mathrm{d}x_n \geq \\ & \int_{[0,1]^n} \frac{\partial f(a; x_1, \cdots, x_n)} {\partial a} f(a; x_1, \cdots, x_n) \mathrm{d}x_1\cdots \mathrm{d}x_n \int_{[0,1]^n} f(a; x_1, \cdots, x_n) \mathrm{d}x_1\cdots \mathrm{d}x_n, \end{aligned} where $$ f(a; x_1, \cdots, x_n) = \frac{\sum_{i=1}^n x_i e^{a x_i}}{\sum_{i=1}^n e^{a x_i}}. $$ It is clear that $f(a; x_1, \cdots, x_n)$ represents the mean with respect to the Gibbs distribution $\frac{e^{a x_i}}{\sum_{i=1}^n e^{a x_i}}$, and that $\frac{\partial f(a; x_1, \cdots, x_n)}{\partial a}$ corresponds to the variance under this distribution. I have performed Monte Carlo simulations with $n$ ranging from $2$ to $10$ over the interval $a\in[-20, 20]$, simulation results suggest that the inequality holds. For example, for $n=2$, the figure of $\text{LHS} - \text{RHS}$ is as follows: Simulation results
I would appreciate any insights or suggestions on how to approach the proof of this inequality.
