# Trying to implement an algorithm to solve a self-consistent Bloch-Schrodinger Equation

I am trying to understand and to implement an algorithm described by Kunimi-Kato in their supplementary material to solve the self-consistent equation below

$$\left[\frac{\hbar^{2}}{2m}|\boldsymbol{G}|^{2}+n_{0}\tilde{V}\left(0\right)\right]C_{\boldsymbol{G}}+\sum_{\Delta G\neq0}S_{\Delta\boldsymbol{G}}C_{\boldsymbol{G}+\Delta G}=\mu C_{\boldsymbol{G}}\qquad(1)$$

in Mathematica, the equation above is the Gross-Pitaevskii Equation defined in Bloch waves, where $$\boldsymbol G$$ is the reciprocal lattice, $$n_0=N/s=\sum_{\boldsymbol{G}}|C_{\boldsymbol{G}}|^2$$, with $$s=\sqrt{3}/2 \lambda ^2$$, the $$S_{\Delta G}$$ is defined as

$$S_{\Delta G} = \tilde{V}(\Delta G)\sum_{G^\prime}C^{*}_{\boldsymbol{G^\prime}+\Delta G}C_{\boldsymbol{G^\prime}},\qquad (2)$$ with $$\tilde{V}(\boldsymbol{k})=J_{1}(k)/k$$. To run the algorithm we will also need the equation that describe the total energy per particle

$$\frac{E}{N}=\frac{1}{n_{0}}\sum_{\boldsymbol{G}}\left|\boldsymbol{G}\right|^{2}\left|C_{\boldsymbol{G}}\right|^{2}+\frac{1}{2n_{0}}\sum_{\boldsymbol{G}_{1},\boldsymbol{G}_{2},\boldsymbol{G}_{3}}\tilde{V}\left(\boldsymbol{G}_{1}-\boldsymbol{G}_{3}\right)C_{\boldsymbol{G}_{1}+\boldsymbol{G}_{2}-\boldsymbol{G}_{3}}^{*}C_{\boldsymbol{G}_{3}}^{*}C_{\boldsymbol{G}_{2}}C_{\boldsymbol{G}_{1}}. \qquad (3)$$

• What did I understand of the algorithm? As described in the reference, we have the following steps:
1. Choose appropriate value of $$\lambda$$.
2. Choose appropriate {$$C_\boldsymbol{G}$$} as an initial condition.

In this step, I know that for the Eq. (1), when $$\Delta G=1$$, we have a coupled equation between $$C_0$$, which I known, and $$C_1$$ that would can be calculated, since we also fix $$S_{\Delta G=1}$$ with appropriate initial condition.

1. Substitute {$$C_\boldsymbol{G}$$} into $$S_{\Delta G}$$ in the Eq. (1) and diagonalize (1) numerically.

Here, for what I understand, by knowing $$C_0$$, $$C_1$$ and $$C_2$$, we can calculate $$S_{\Delta G=2}$$, for instance. And so this could be done recursively.

1. Calculate the total energy per particle for each eigenstate.
2. Choose {$$C_\boldsymbol{G}$$} for the lowest energy state.
3. Iterate (3)-(5) until convergence.
4. Choose different value of $$\lambda$$ to minimize the total energy per particle.
5. Iterate (2)-(7) until convergence.

So what is the point?

Well, I do not have much experience in mathematica and also I am trying to learn it, so I took this problem to guide me, but after spend days searching for references, books or anything else that could help me in understand and learn how to implement this algorithm in mathematica, now all what left to me was writing here hopefully in someone to give me a light.