In this paper we obtain optimal multipolar Rellich inequality in \(\mathbb {R}^{N}\) for biharmonic Schrödinger operator with positive multisingular potentials of the form \(\begin{aligned} H_{n}:=\Delta ^{2}-\frac{N^{2}(N-4)^{2}}{n^{4}}\sum _{1\le i<j\le n}\frac{\left| a_{i}-a_{j}\right| ^{2}}{|x-a_{i}|^{2}|x-a_{j}|^{2}} \left( \sum _{1\le k<l\le n}\frac{\nu _{k,i,j}\nu _{l,i,j}\left| a_{k}-a_{l}\right| ^{2}}{|x-a_{k}|^{2}|x-a_{l}|^{2}}\right) , \end{aligned}\) where \(a_{1},\ldots ,a_{n}\) are n different singular poles, \(\nu _{k,i,j}=\frac{N+2n-4}{N}\) when \(k=i\) or \(k=j\) , otherwise \(\nu _{k,i,j}=\frac{N-4}{N}\) . We also prove that the best constant \(\frac{N^{2}(N-4)^{2}}{n^{4}}\) is attained in \(D^{2,2}(\mathbb {R}^{N})\) for \(n\ge 3\) , but not attained for \(n=2.\) Moreover, we prove the criticality of the biharmonic Schrödinger operator \(H_{n}\) for \(n\ge 2\) . Finally, we get a class of higher-order bipolar Rellich inequality.
