\mathbf{k}}}({\mathbf{r}})={\mathrm{e}^{\mathrm{i} \phi_{n{\mathbf{k}}}}}{\mathrm{e}^{\mathrm{i} \mathbf{k} \mathbf{r}}}{{\mathbf{u}}_{n \mathbf{k}}}({\mathbf{r}})
\label{eq:bloch_def}$$ with respect to the wave vector ${\mathbf{k}}$, $${{\mathbf{W}}_{n \mathbf{R}}}({\mathbf{r}})=\frac{1}{\sqrt{N}}\sum_{{\mathbf{k}}\in\text{BZ}}{\mathrm{e}^{-\mathrm{i} \mathbf{k} \mathbf{R}}}{{\mathbf{B}}_{n \mathbf{k}}}({\mathbf{r}}).\label{eq:wan_def}$$ The Bloch mode is a solution of the corresponding wave equation in a periodic medium, where ${{\mathbf{u}}_{n \mathbf{k}}}({\mathbf{r}})$ is a lattice-periodic envelope function, and we have explicitly denoted the arbitrary phase $\phi_{n{\mathbf{k}}}$ of the Bloch function, the Bloch phase. $n$ is the band index. For simplicity, in the present work we do not consider mixing of different bands [@Jiao2005; @john] in the construction of Wannier functions from Bloch modes.

From the definitions (\[eq:bloch\_def\]), (\[eq:wan\_def\]) it is seen that the Wannier functions would be $\delta ({\mathbf{r}})$-localized in space only if the envelopes ${{\mathbf{u}}_{n \mathbf{k}}}({\mathbf{r}})$ were constant. By contrast, for arbitrary Bloch phases $\phi_{n{\mathbf{k}}}$ the Wannier functions have actually a rather large spatial extension due to the oscillatory character of the ${{\mathbf{u}}_{n \mathbf{k}}}({\mathbf{r}})$, a feature which becomes more and more pronounced for higher band indices $n$. However, the gauge freedom in the Bloch phases may be employed to adjust, for each ${\mathbf{k}}$, the value of $\phi_{n{\mathbf{k}}}$ so as to construct maximally localized Wannier functions with respect to some locality criterion. This constitutes a complex, multidimensional optimization problem. The choice of the locality measure is 