Obtener una densidad de unión estrecha de los estados con mayor precisión

Lo siento por la respuesta tardía, pero espero que esto pueda ser útil para alguien más.

Puede reducir el ruido utilizando una integral elíptica.

$$D(\varepsilon)=\frac{1}{4 \pi ^3 t}\int_{-\pi }^{\pi }d\phi K\left( \sqrt{1-\left(\frac{\varepsilon +2 t \cos\phi }{4 t}\right)^2}\right)$$

Donde K es la integral elíptica completa de primer tipo: http://mathworld.wolfram.com/CompleteEllipticIntegraloftheFirstKind.html .

No es baladí llegar hasta aquí. E incluso a partir de esta expresión, la integral debe hacerse numéricamente con cuidado (tiene singularidades para muchos valores de$\varepsilon$) pero debería dar mejores resultados. Correr durante seis segundos en Mathematica me da (con$t=\frac{1}{2\sqrt{3}}$):

TL;DR

Con Matlabo Mathematikano puedo ayudar, pero Pythonhay una implementación disponible: sc_dos Aquí$D$es la mitad del ancho de banda$D = 6t$.

import numpy as np
import gftool as gt

eps = np.linspace(-1.2, 1.2, num=6001)
dos = gt.sc_dos(eps, half_bandwidth=1)

La evaluación del DOS me toma ~100 ms.

---

Ya diste las expresiones correctas. Tenemos la función de Green

$$G(z) = \frac{1}{N} \sum_{\boldsymbol{k}} \frac{1}{z-\epsilon_{\boldsymbol{k}}}$$ y la densidad (normalizada) de estados (DOS) $$D(\epsilon) = \frac{1}{N} \sum_{\boldsymbol{k}} \delta(\epsilon - \epsilon_{\boldsymbol{k}}) = -\frac{1}{\pi} \Im G(\epsilon+i0^+),$$ dónde$\epsilon$es una variable de energía real. La segunda igualdad es Sokhotski–Plemelj .

Una suma ingenua sobre los puntos es extremadamente exigente, ya que una gran cantidad de$\boldsymbol{k}$puntos es necesario en 3 Dimensiones.

Para suavizar la función, podemos evaluar la función de Green en un contorno paralelo al eje real desplazado por un finito$\eta>0$en semiplano complejo superior:

$$D(\epsilon) \approx -\frac{1}{\pi} \Im G(\epsilon+i\eta).$$

Cuanto más grande elegimos$\eta$más suave se vuelve la función, pero por otro lado perdemos características.

Como solo nos interesa el límite termodinámico$N\rightarrow \infty$, un enfoque más inteligente que el simple muestreo$\boldsymbol{k}$, es reemplazar la suma por la integral. Para integrales, tenemos algoritmos más o menos eficientes.

Así que calculemos

$$G(z) = \int \frac{d^3 k}{{2\pi}^3} \frac{1}{z-\epsilon_{\boldsymbol{k}}}$$ en cambio. El$\epsilon_\boldsymbol{k}$, es simétrico para todo$k_{x_i}$:$\epsilon(k_{x_i}) = \epsilon(-k_{x_i})$, por lo que es suficiente para integrar más de una octava parte de la zona de Brillouin.

Y finalmente podemos usar resultados analíticos para las integrales. Notamos que podemos expresar la función de Green 3D en términos de resultados conocidos de la función de Green 1D y 2D como tenemos

$$\epsilon_{\boldsymbol{k}} = \epsilon_{k_x, k_y, k_z} = \epsilon^{2D}_{k_x, k_y} - 2t \cos(k_z) = \epsilon^{1D}_{k_x} - 2t[\cos(k_y) + \cos(k_z)]$$ y por lo tanto $$G(z) = \int \frac{dk_z}{2\pi} G^{2D} (z - 2t\cos(k_z)) = \int \frac{dk_y}{2\pi} \int \frac{dk_z}{2\pi} G^{1D}(z - 2t[\cos(k_y) + \cos(k_z)]).$$

La función de Green unidimensional$G^{1D}(z)$se puede evaluar fácilmente, la función de Green bidimensional$G^{2D}(z)$se puede expresar en términos de la integral elíptica completa de primer tipo (que se puede encontrar en los libros de texto estándar). Usando$G^{2D}(z)$es básicamente el resultado dado por bRost03 .

Un tipo muy inteligente llamado Joyce incluso encontró una expresión para$G^{3D}(z)$en 1973. Las ecuaciones son un poco largas y complicadas, así que evitaré copiarlas aquí. Pero los implementamos en un Pythonmódulo gftool>=0.8.0, vea sc_dos . También encontrará las referencias relevantes allí.

Tuve el mismo problema. usé la fórmula$\rho(\epsilon)=\sum_{\vec{k}}\delta(\epsilon - \epsilon(\vec{k}))$para calcular numéricamente la densidad de estados. hice un resumen$100$ $k-$valores para cada componente, y usó una distribución gaussiana con$\sigma =0.1$para que la función delta obtenga el siguiente diagrama. Usando más grande$\sigma$termina suavizando la singularidad en la derivada alrededor de$\epsilon=\pm 2$.

El código se escribió en C++ y se ejecutó durante aproximadamente$60$segundos.$y$el eje es$\rho(\epsilon)$,$x$el eje es$\epsilon$y$t=1$.

PD: Hice la suma de la mitad de la zona de Brillouin, que era lo que necesitaba para mi aplicación.