The Wannier center was chosen as ’on-site’ for the 1st band and as ’between’ for the 2nd, 4th and 5th band. []{data-label="fig:MLWF_airRods_te"}](mlwf_airRods_te){width="0.95\columnwidth"}

Sometimes, it is beneficial to use the time-reversal symmetry of wave operators $\mathcal{L}_{E}$ and $\mathcal{L}_{H}$ to construct real Wannier functions. This is possible, since the envelope functions of the Bloch modes transform as $${u_{n \mathbf{-k}}}({\mathbf{r}})={u_{n \mathbf{k}}}^{*}({\mathbf{r}})$$ under time reversal (inversion of the reciprocal space), and the Bloch phases can be constrained by the condition $$\phi_{n{\mathbf{k}}}=-\phi_{n{\mathbf{-k}}}\ .\label{eq:phaseInversion}$$ By restricting the parameter space of the optimization problem in such a way, we found that both SM and IM optimization criteria possess multiple extrema, making the use of local optimization method not efficient without special choice of the initial Bloch phases. However, using real-valued Wannier functions is of advantage for the considerations of the next section.

Choice of initial conditions\[sec:BC\]
======================================

Here we propose an analytical expression for a generic set of Bloch phases to be used as a starting point for the optimization procedure. It is based on the following theorem.

[*Theorem:*]{} Suppose that (i) the Wannier functions are real-valued, i.e., $\phi_{n{\mathbf{k}}}=-\phi_{n{\mathbf{-k}}}$, and the Bloch functions ${{\mathbf{B}}_{n \mathbf{k}}}({\mathbf{r}})={\mathrm{e}^{\mathrm{i} \phi_{n{\mathbf{k}}}}}
{{\mathbf{\tilde B}}_{n \mathbf{k}}}({\mathbf{r}})$ conform to the following conditions: $\text{Re}({{\mathbf{B}}_{n \mathbf{k}}}({\mathbf{r}}))$ has 