# Numerical problem in solving the Bogoliubov de Gennes equations- methods to solve?

I am trying to solve an assignment on solving the Bogoliubov de Gennes equations self-consistently in Matlab. BdG equations in 1-Dimension are as follows:-

$$\left(\begin{array}{cc} -\frac{\hbar^{2}}{2m}\frac{\delta^{2}}{\delta z^{2}}-\mu+V\left(z\right) & \triangle(z)\\ \triangle(z) & \frac{\hbar^{2}}{2m}\frac{\delta^{2}}{\delta z^{2}}+\mu-V(z) \end{array}\right)\left(\begin{array}{c} u_{n}(z)\\ v_{n}(z) \end{array}\right)= \epsilon_{n}\left(\begin{array}{c} u_{n}(z)\\ v_{n}(z) \end{array}\right)$$ along with the equations for gap function $\triangle(z)$ and number density $n(z)$. $$\triangle(z)=U\sum_{n}\left(1-2f_{n,}\right)u_{n}(z)v_{n}^{\star}(z)$$ and $$n(z)=2\sum_{n}|{u_{n}(z)}|^{2}f_{n}+|{v_{n}(z)}|^{2}\left(1-f_{n}\right).$$

For the case of solving the BdG equations in Fourier space in Matlab for the case of a periodic potential and periodic gap function (assumed), we can take $$u_{n}(z)=\sum_{k}\exp\left[ikz\right]U_{n,k},$$ $$\triangle(z)=\sum_{K}\exp(iKz)T_{K},$$ and $$V(z)=\sum_{K}\exp(iKz)P_{K}$$ where the sum is over the reciprocal lattice vectors $K$ leading to the number equation $$N=2\sum_{n,k}\left[f_{n}|{U_{n,k}}|^{2}+\left(1-f_{n}\right)|{V_{n,k}}|^{2}\right]$$ with $f_{n}$ as the Fermi distribution function. Solving the set of equations self-consistently for a fixed $N$, I am trying to get a value of chemical potential from the number equation each time after solving the eigenvector components $U_{n,k}$ and $V_{n,k}$, but due to the form of the exponentials in the number equation and sum over large number of them, I am unable to get a correct value of chemical potential out of them using Matlab routines as the root of the equation to put it back into the equations for eigenvector components.

In most cases, I get random values of chemical potential since the equation is more or less insoluble. How can I avoid this error ? Is there a better way to numerically solve the BdG equations self-consistently ? I also want to do this assignment in real space avoiding finite size effects but started with the Fourier space case to avoid errors associated with discretizing the differential. Please guide and ask for any details you might need.

Following is my MATLAB code to solve the equations in real space but the code does not work as fsolve does not find the mu value.

http://postimg.org/image/v7amx5vd9/full/

## http://postimg.org/image/do9pyyk7z/full/

If you have solved BdG equations numerically, please tell about the method and steps you used such that the above problems are eliminated.

• Would Computational Science be a better home for this question? – Qmechanic Mar 10 '14 at 8:19
• I have posted this in CompSci but I am only getting downvotes and not getting any comments or answers there. – user38579 Mar 10 '14 at 8:43
• I will not help you since I have no idea on how to resolve your problem. But still, I'm really wondering why your equations are written in dimension-full variables... You may never wonder about the chemical potential if you cancel it consistently in writing everything in dimensionless variables. Also, to go to the real space with boundary conditions, you should NOT use Fourier transform at all I believe. It works only for periodic boundary conditions. In self-consistent problem, you need to know a bit about the $\Delta$. Most of the solutions are for $\Delta=0$, the trivial one. – FraSchelle Mar 17 '14 at 12:24
• (cont.) The interesting solutions are obviously for $\Delta\neq 0$, but I'm not sure you started with such a possibility. You may find interesting the following solution physics.stackexchange.com/questions/54200 as well regarding the self-consistent problem. – FraSchelle Mar 17 '14 at 12:26
• @FraSchelle : I keep the chemical potential because it remains in the number equation even if you manage to remove it from the matrix equation and these equations have to be solved together self-consistently. Secondly, I go to the Fourier space only for the case of periodic potential, as mentioned in my question and I am doing that only to avoid discretization errors encountered if I directly solve the matrix equations in real space. Though I agree FFT will give errors of its own. – user38579 Mar 17 '14 at 19:58