## Appendix A  
We present the steps from the Schrodinger system (1), (2) for the variables \(u\), \(\psi\) to the infinite system (30) for their Wannier coefficients, see also [1],  
We first consider the second nematicon Equation (2). We multiply (2) by \(\phi^{*}_{n^{\prime},k^{\prime}}(y)\), integrate over \(y\in\mathbb{R}\), and use the Schrodinger Equation (12) to obtain  
\[\int_{\mathbb{R}}\psi(y)(E^{*}_{n^{\prime},k^{\prime}}+g^{2})\phi^{*}_{n^{ \prime},k^{\prime}}(y)\,dy=\alpha\int_{\mathbb{R}}|u(y)|^{2}\phi^{*}_{n^{\prime },k^{\prime}}(y)\,dy,\] (A1)  
hence  
\[\int_{\mathbb{R}}\psi(y)\phi^{*}_{n^{\prime},k^{\prime}}(y)\,dy=\alpha(E_{n^{ \prime},k^{\prime}}+g^{2})^{-1}\int_{\mathbb{R}}|u(y)|^{2}\phi^{*}_{n^{\prime },k^{\prime}}(y)\,dy.\] (A2)We then multiply (A2) by \(e^{imk^{\prime}b}\) and integrate both sides over \(k\in[-\pi/b,\pi/b]\). Interchanging the order of integration in \(k,y\), and using the definition of Wannier functions we obtain that \(\psi\) is given by (19) with  
\[c_{n,m}=\alpha\sqrt{\frac{b}{2\pi}}\sum_{n_{1},n_{2}\in\mathbb{N}}\sum_{m_{1},m _{2}\in\mathbb{Z}}K_{n_{1},n_{2},m}^{m_{1},m_{2},m}u_{n_{1},m_{1}}u_{n_{2},m_{2 }}^{*},\] (A3)  
where  
\[K_{n_{1},n_{2},m}^{m_{1},m_{2},m}=\int_{\mathbb{R}}w_{n_{1}}^{m_{1}}(y)w_{n_{2 }}^{m_{2}*}(y)\left\{\int_{-\pi/b}^{\pi/b}\frac{\phi_{n,k}^{*}(y)e^{imbk}}{E_{n, k}+g^{2}}\,dk\right\}dy.\] (A4)  
Expanding \(\phi_{n,k}\) in the Wannier basis as in (17) we have  
\[K_{n_{1},n_{2},m}^{m_{1},m_{2},m} = \sqrt{\frac{b}{2\pi}}\int_{\mathbb{R}}w_{n_{1}}^{m_{1}}(y)w_{n_{2 }}^{m_{2}*}(y)\left\{\sum_{m^{\prime}\in\mathbb{Z}}\left(\int_{-\pi/b}^{\pi/b }\frac{e^{i(m-m^{\prime})bk}}{E_{n,k}+g^{2}}\,dk\right)w_{n}^{m^{\prime}*}(y) \right\}dy\] (A5) \[= \sqrt{\frac{b}{2\pi}}\int_{\mathbb{R}}w_{n_{1}}^{m_{1}}(y)w_{n_{ 2}}^{m_{2}*}(y)\left\{\sum_{m^{\prime}\in\mathbb{N}}\sum_{m^{\prime}\in \mathbb{Z}}\delta_{n,n^{\prime}}\left(\int_{-\pi/b}^{\pi/b}\frac{e^{i(m-m^{ \prime})bk}}{E_{n^{\prime},k}+g^{2}}\,dk\right)w_{n^{\prime}}^{m^{\prime}*}(y) \right\}dy\] \[= \sum_{n^{\prime}\in\mathbb{N}}\sum_{m^{\prime}\in\mathbb{Z}} \sqrt{\frac{b}{2\pi}}\delta_{n,n^{\prime}}\left(\int_{-\pi/b}^{\pi/b}\frac{e^{ i(m-m^{\prime})bk}}{E_{n^{\prime},k}+g^{2}}\,dk\right)\int_{\mathbb{R}}w_{n_{1}}^{m _{1}}(y)w_{n_{2}}^{m_{2}*}(y)w_{n^{\prime}}^{m^{\prime}*}(y)\,dy.\]  
By (A2), (A5) we then have  
\[c_{n,m}=\alpha\sum_{n_{1},n_{2},m^{\prime}\in\mathbb{N}}\sum_{m_{1},m_{2},m^ {\prime}\in\mathbb{Z}}G_{n,n^{\prime}}^{m,m^{\prime}}l_{n_{1},\mathbb{H}_{2}, \mathbb{H}^{\prime}}^{m_{1},\mathbb{H}_{2},\mathbb{H}^{\prime}}u_{n_{1},m_{1 }}u_{n_{2},m_{2}}^{*},\] (A6)  
with  
\[G_{n,n^{\prime}}^{m,m^{\prime}}=\frac{b}{2\pi}\delta_{n,n^{\prime}}\left( \int_{-\pi/b}^{\pi/b}\frac{e^{i(m-m^{\prime})bk}}{E_{n^{\prime},k}+g^{2}}\,dk \right),\quad l_{n_{1},\mathbb{H}_{2},\mathbb{H}^{\prime}}^{m_{1},\mathbb{H} _{2},\mathbb{H}^{\prime}}=\int_{\mathbb{R}}w_{n_{1}}^{m_{1}}(y)w_{n_{2}}^{m_ {2}*}(y)w_{n^{\prime}}^{m^{\prime}*}(y)\,dy.\] (A7)  
To expand the first nematicon Equation (1) in coefficients of Wannier functions, we substitute the series expansions (19), (20) into (1), multiply (1) by \(w_{n^{\prime}}^{m^{\prime}*}(y)\), and integrate over \(y\in\mathbb{R}\).