Given a complex-valued function $f: \mathbb{R}^N\rightarrow\mathbb{C}$ with $f\in L^2$, which vanishes at the domain boundaries, is twice differentiable, and antisymmetric upon two-coordinate exchange—i.e. $f(x_i,x_j)=-f(x_j,x_i)$—, can we prove that the following term is real?
$$\frac{\nabla^2 f}{f}\stackrel{?}{\in}\mathbb{R}$$
This is trivially true, if $f$ is real or an eigenfunction of $\nabla^2$, but is this true in general?
I can only show, that its imaginary parts have to cancel out upon integration: we know, that $\nabla^2$ is a Hermitian operator. Thus, its expectation value $\langle\nabla^2\rangle$ is real. $$ \langle\nabla^2\rangle=\int\frac{\nabla^2 f}{f}f^\ast f\,\mathrm{d}\tau\in\mathbb{R} $$ The probability density in this expected value is $f^\ast f$, which is real. Thus, if $\nabla^2f/f$ is not real somewhere, the imaginary parts of the whole integrand have to cancel out exactly, since the result is real.
Another approach is to simpliy insert $f=a+ib$ with $a,b\in\mathbb{R}$ and getting the denominator real: $$\frac{\nabla^2f}{f}=\frac{\nabla^2a+i\nabla^2b}{a+ib}=\frac{a\nabla^2a+b\nabla^2b+i(a\nabla^2b-b\nabla^2a)}{a^2+b^2}$$ Thus, the question really is: $$a\nabla^2b-b\nabla^2a\stackrel{?}{=}0$$ Is there any way to prove or disprove this with the properties of a function in $L^2$? A counter-example is also sufficient.
It is false. Let us first construct a counterexample without the symmetry assumption, then we will modify it to make it symmetric, using an idea suggested in comments by Hussein.
Consider the following function, defined for $x\in\mathbb R^d$; $$ f(x)= e^{2i\xi\cdot x}+ e^{i\xi\cdot x}, $$ where $\xi\in \mathbb R^d$ is a unit vector. Then $$ -\frac{\nabla^2 f}{f}=\frac{2e^{2i\xi\cdot x} + e^{i\xi\cdot x}}{e^{2i\xi\cdot x}+e^{i\xi\cdot x}}=\frac{2 e^{i\xi\cdot x} +1}{ e^{i\xi\cdot x} +1}, $$ whose imaginary part equals $\frac12\sin(\xi\cdot x)\ne 0$.
Remark. There is nothing special about $\xi$ and $2\xi$. The idea is to sum two different complex eigenfunctions of the Laplacian $\nabla^2$.
Now we make this example symmetric. Unfortunately the argument looks a bit messy, but the idea is really simple: we are going to localize the previous $f$ in a small ball and then extend it by symmetry.
To make this formal, let $\sigma_{ij}\colon \mathbb R^d\to \mathbb R^d$ be the permutation map $$ \sigma_{ij}(x_1, \ldots, x_i, \ldots, x_j, \ldots, x_d)=(x_1, \ldots, x_j, \ldots, x_i, \ldots, x_d).$$ There are $N=\binom{d}{2}$ such permutations. Let $B$ denote a ball in $\mathbb R^d$, with center at the point $(1, 2, \ldots, d)$ and radius so small that $\sigma_{ij}(2B)$ is disjoint from $2B$ for all $i$ and $j$. (The choice of the center is mostly arbitrary. Any point that does not fix any $\sigma_{ij}$ will do). Consider a smooth real function $\phi$ that equals $1$ on $B$ and is compactly supported inside $2B$, which is often written $$ B\prec\phi\prec 2B.$$ And now iteratively construct $$ \begin{cases} f_0:=f\phi, \\ f_1:=f_0-f_0\circ \sigma_{12}, \\ f_2:=f_1-f_1\circ \sigma_{13}, \\ \vdots\\ f_N:=f_{N-1}-f_{N-1}\circ \sigma_{(d-1)\, d}. \end{cases} $$ The function $f_N$ has the desired symmetry $f_N\circ \sigma_{ij}=-f_N$ for all $i, j\in\{1, \ldots, d\}$ with $i\ne j$, and it equals $f$ on $B$. Thus, in particular, $\nabla^2 f_N/f_N$ is not purely real. $\Box$
Remark. The original question contains an argument based on $\langle \nabla^2\rangle$ being real. That shows that $\nabla^2 f/f$ is real on average, that is, that its imaginary part must have a vanishing integral. However, as this counterexample shows, it is not possible to conclude from that argument that such imaginary part is zero at every point.