What I already know:
Say we have a Gaussian distribution $X \sim N(0,\sigma^2)$. We know that $\text{Pr}[a\le x\le b]$ is a unimodal function of $\sigma$. The reason is as follows. First, we define $f(\sigma) = \text{Pr}[a\le x\le b]$, and it can be written by $$f(\sigma) = \int_a^b \frac{1}{\sigma\sqrt{2 \pi} } e^{-\frac{1}{2}\left(\frac{x}{\sigma}\right)^2} \ dx=\frac{1}{2} \left(\text{erf}\left(\frac{b}{\sqrt{2}\sigma}\right)-\text{erf}\left(\frac{a}{\sqrt{2}\sigma}\right)\right),$$ where $\text{erf}(\cdot)$ is Error function. And the derivative of $f(\sigma)$ is $$f'(\sigma)=\sqrt{\frac{2}{\pi}} \frac{\left(a e^{-\frac{a^2}{2\sigma^2}} - b e^{-\frac{b^2}{2\sigma^2}}\right)}{2\sigma^2}.$$ Clearly, this derivate has only one zero point, so $f(\sigma)$ is unimodal.
What I want to know:
Can we extend the unimodality to bivariate (or multivariate if possible) Gaussian?
More specifically, we have a bivariate Gaussian distribution $(X,Y) \sim N(0,\sigma^2)$. (I get the bivariate Gaussian function of variables $x$ and $y$ from this paper, i.e., $X$ and $Y$ are independent, and $\sigma_1 = \sigma_2 = \sigma$ in our problem.) I would like to know if the double integral over an area as a function of $\sigma$ is also unimodal. Define function $g(\sigma)$ as $$g(\sigma) = \iint_D \frac{1}{2\pi \sigma^2}e^{-\frac{x^2+y^2}{2\sigma^2}}dxdy ,$$ where $D$ is an area for point $(x,y)$.
What I have tried:
If area D is an annular sector whose center is origin, we can know $g(\sigma)$ is unimodal by the conclusion from univariate Gaussian. However, for other shapes, the double integral is not easy to get an expression consisting of elementary functions to analyze.
I have tried some areas (e.g., $D_1:1\le x\le2 \wedge 1\le y\le2$, $D_2:(x-2)^2+y^2 \le 1$) via Mathematica, and the plots show that $g(\sigma)$ is unimodal. I guess that $g(\sigma)$ is unimodal if area $D$ is a connected set (or more strictly: convex set).
Anyone has ideas to prove it?
Update:
By the counter-example from @MathWonk, a connected area $D$ is not enough to make $g(\sigma)$ unimodal. Then what if a convex area $D$?
Consider this counter-example: a ball centered at the origin and a concentric annular ring. The annular subregion is defined by $2<r<3$ and the ball has radius 1/2. The union of these two subregions is the total region of integration. The integral depends on the parameter $t=\sigma$. The graph is shown below:
P.S. Note that you can connect the the two subregions with a narrow neck that has negligible area to create a connected region that has similar properties.
The counter-example is based on the idea that (i) the integral over the ball is initially close to 1 but is monotone decreasing, but (ii) the integral over the annulus is initially close to zero, then rises and falls. The sum of these two effects creates a double spike. The underlying intuition that explains this counter-example is that the probability distribution describes diffusion via a random walk, and (i) as time progresses the walker who is initially almost certainly near the origin will be progressively less likely to be in the ball, and (ii) the walker will initially not be in the annulus but will eventually reach it and then gradually wander farther out beyond it. The second picture shows the graphs of (i), (ii) and (i)+(ii).
I still need to think about the convex case!