Problema 7: Defectos de Frenkel
-------------------------------

Una red cristalina perfecta está formada por :math:`N` átomos de la
misma especie.  Si se extraen :math:`n` átomos de sus lugares en la
red (con :math:`1\ll n\ll N` ) y se los coloca en posiciones
intersticiales, se obtienen :math:`n` defectos de tipo Frenkel. El
número :math:`N'` 0 de posiciones intersticiales en la red es del
orden de magnitud de :math:`N`. Sea :math:`W` la energía necesaria
para producir un defecto. Halle el valor de :math:`\langle E\rangle =
W \langle i\rangle` y de allí muestre que 

.. math::
   \langle n\rangle \approx \sqrt{N\,N'}\text{e}^{-\beta W/2}

Grafique cualitativamente :math:`\Omega(n)\text{e}^{-\beta nW}` en
función de :math:`n`. Resuelva este problema tanto en el ensamble
microcanónico como en el canónico.


Ensamble microcanónico
^^^^^^^^^^^^^^^^^^^^^^

En el ensamble microcanónico tenemos que calcular :math:`\Omega(N,
E)`, y a partir de esta superficie obtenemos la entropía:

.. math::
   S(\Omega) = k_B\,\log\Omega

Queremos hallar los estados compatibles con una energía :math:`E =
W\,n`. Si quitamos :math:`n` átomos de sus posiciones de la red,
tenemos :math:`{N \choose n}` formas de quitar los átomos y :math:`{N'
\choose n}` formas de escoger qué posiciones intersticiales
ocupamos. Así, la superficie es:

.. math::
   \Omega(n) &= {N \choose n}{N' \choose n} \\
   &= \frac{N!}{(N-n)!n!}\frac{N'!}{(N'-n)!n!}

y calculamos la entropía:   

.. math:: 
   S(n) = k_B\,\log{\left(\frac{N!}{(N-n)!n!}
   \frac{N'!}{(N'-n)!n!}\right)}

Con la aproximación en la que :math:`1\ll n\ll N`, todos los términos
son muy grandes, y podemos aplicar la aproximación de Stirling.

.. math::
   S(n) \approx S_S(n) = k_B\,[&N \log N - N - n \log n + n - (N - n)
   \log (N - n) + (N - n) + \\ 
    +&N' \log N' - N' - n \log n + n - (N' - n)
   \log (N' - n) + (N' - n)]

con las simplificaciones típicas de números combinatorios (ver
:ref:`apendice`), definiendo :math:`x = n/N`, :math:`\alpha = N'/N`, resulta

.. math::
   S_S(x)  = - k_B[&N(x \log{x} + (1 - x) \log{(1-x)}) \\
   +&N'\left(\frac{x}{\alpha} \log{\frac{x}{\alpha}} + \left(1 -
   \frac{x}{\alpha}\right) \log\left(1-\frac{x}{\alpha}\right)\right)]
   :label: s_approx


Si queremos encontrar la temperatura, tenemos que calcular la
derivada de la entropía respecto de la energía. Como hicimos en
:ref:`dosniveles`, 

.. math::
   \frac{1}{T} &= \frac{dS}{dE}\\
   &= \frac{dS}{dn} \frac{dn}{dE} \\
   &= k_B\log{\frac{\Omega(n+1)}{\Omega(n)}}\, 
   \frac{dn}{dE}
   :label: domega

En nuestro caso, :math:`dn/dE = 1/W` y 

.. math::
   \Omega(n + 1) &= {N \choose (n + 1)}{N' \choose (n + 1)} \\
   &= \frac{N!}{(N-n-1)!n+1!}\frac{N'!}{(N'-n-1)!n+1!}\\
   &= \frac{N!}{(N-n)!n!}\frac{N-n}{n+1}\,
   \frac{N'!}{(N'-n)!n!}\frac{N'-n}{n+1}\\
   &= \Omega(n)\frac{(N-n)(N'-n)}{(n+1)^2}
   :label: nextomega

Queda la temperatura, entonces
   
.. math::
   \frac{W}{k_BT} = \log{\frac{(N-n)(N'-n)}{(n+1)^2}}
   
De aquí despejamos :math:`n`, definiendo :math:`\omega =
\text{e}^{W/k_BT}`.

.. math::
   \frac{(N-n)(N'-n)}{(n+1)^2} & = \omega \\
   n^2 - (N' + N) n + N'N & = \omega n^2 + 2\,\omega n + \omega\\
   (1-\omega) n^2 - (N' + N - 2\omega) n + N'N -\omega & = 0
   :label: caract


Y obtenemos

.. math:: 
   n = \frac{N' + N - 2\omega \pm \sqrt{(N' + N - 2\omega)^2 - 4
   (1-\omega) (N'N -\omega)}}{2-2\omega}


En el límite termodinámico, :math:`N \gg \omega`, y la expresión para
:math:`n` queda

.. math:: 
   n = \frac{N' + N + \sqrt{(N' + N)^2 - 4 (1-\omega)
   (N'N)}}{2-2\omega}

Finalmente, veamos cómo se comporta :math:`n` el régimen que nos piden
en el enunciado, :math:`1 \ll n \ll N`. Parece intuitivo notar que esta
aproximación es para temperaturas bajas. En este caso, podemos hacer
la aproximación en :eq:`caract`

.. math::
   \omega = \frac{(N-n)(N'-n)}{(n+1)^2} \approx \frac{N'N}{n^2} \\
   :label: caract_approx
   
y despejando, llegamos al resultado

.. math::
   \langle n\rangle \approx \sqrt{N\,N'}\text{e}^{-\beta W/2}

Podemos ver en qué régimen estas dos expresiones :eq:`caract` y
:eq:`caract_approx` son parecidas, con :math:`N = 10000`, y
:math:`\alpha = 0.8`:

.. plot:: guia4/problema7/defectos.py
   :include-source:

Aparece nuevamente un régimen de temperaturas negativas. Podemos ver
además que la aproximación no es muy buena para :math:`x \approx 0.1`
[#frenk]_. ¿Cuándo comienzan las temperaturas negativas? Cuando se
anula el argumento del logaritmo, es decir

.. math::
   n = \frac{N'N}{N'+N}

Ensamble canónico
^^^^^^^^^^^^^^^^^

Resolvamos ahora este problema en el ensamble canónico, suponiendo que
:math:`N < N'`.

.. math::
   Z(N, T) &= \sum_{n=0}^{N} \Omega(n)\,\text{e}^{-\beta E} \\
   &= \sum_{n=0}^{N}\frac{N!}{(N-n)!n!}\frac{N'!}{(N'-n)!n!},\text{e}^{-n\beta W}\\
   &= \sum_{n=0}^{N} \zeta_n
   :label: zeta


Esta suma no parece nada fácil de calcular, así que vamos a ver cómo
la podemos aproximar. Vamos a gráficar :math:`\zeta_n` para :math:`N = 500`

.. plot:: guia4/problema7/particion.py
   :include-source:

Vemos que hay claramente un máximo, y tiene una forma de campana. Por
lo pronto, entonces, vamos a tratar de aproximarla por una
gaussiana. Para eso, analicemos los términos que la componen,
:math:`\zeta = \text{e}^\lambda`

.. math::
   \text{e}^{\lambda} &= \Omega(n)\,\text{e}^{-n\beta W} \\
   \lambda &= S(n) - n\beta W

El máximo de :math:`\lambda` se da para

.. math::
   \frac{d\lambda}{dn} &= 0 \\
   \frac{dS}{dn} - \beta W &= 0

Utilizando :eq:`domega` y :eq:`nextomega`, podemos hallar la derivada
de la entropía. De esta manera, el máximo se da en

.. math::
   \log{\frac{(N-n)(N'-n)}{(n+1)^2}} = \beta W

A partir de aquí vamos a analizar (para que las cuentas sean más
sencillas) el caso :math:`1\ll n\ll N` [#aprox]_. Así, resulta

.. math::
   \log{\frac{N'N}{n^2}} &= \beta W \\
   n_{\text{máx}} &= \sqrt{N'N} \text{e}^{-\beta W/2} \\

Y la derivada segunda en el máximo es

.. math:: 
   \left.\frac{d^2\lambda}{dn^2}\right|_{n = n_{\text{máx}}} &= \frac{d\log{\frac{N'N}{n^2}}}{dn} \\
   & = -2 / n_{\text{máx}} \\
   & = -2 \frac{\text{e}^{\beta W/2}}{\sqrt{N'N}} 

Tenemos entonces, a orden 2 en Taylor, la expresión de los términos de
la suma:

.. math:: 

   \lambda \approx \lambda_{\text{máx}} - \frac{(n - n_{\text{máx}})^2}{\sqrt{N'N}\text{e}^{-\beta W/2}}

Volviendo a la expresión original, :eq:`zeta`, 

.. math:: 
   \zeta \approx \zeta_{\text{máx}} \text{e}^\frac{(n -
   n_{\text{máx}})^2}{2\sigma^2}
   :label: part_approx

con

.. math::

   \sigma^2 = \sqrt{N'N}\text{e}^{-\beta W/2}

Bueno, entonces, ¿para qué dimos tantas vueltas? En la ecuación
:eq:`part_approx` aproximamos la función de partición cerca del máximo
como una distribución gaussiana. Podemos calcular la sumatoria como el
área bajo esta gaussiana:

.. math::
   Z(N, T) &= \sum_{n=0}^{N} \Omega(n)\,\text{e}^{-\beta E} \\
   &\approx \int_{-\infty}^\infty dn \zeta_{\text{máx}} \text{e}^\frac{(n -
   n_{\text{máx}})^2}{2\sigma^2} \\
   &\approx \zeta_{\text{máx}} \,\sqrt{2\,\pi}\sigma

Recordemos que al fin y al cabo, lo único que nos importa en la
termodinámica es el logaritmo de la función de partición:

.. math::
   \log{Z(N, T)}&\approx \lambda_{\text{máx}} + \frac{\sqrt{2\,pi}}{4}(\log{N'N} -
   \beta W/2) \\
   &\approx S(n)/k_B - n\beta W - \gamma\log{N'N}

Como :math:`S(n)` y :math:`n` son variables extensivas, crecen mucho
más rápido que cualquier término :math:`\log{N}`, por lo que 

.. math::
   \log{Z(N, T)}\approx S(n_{\text{máx}})/k_B - n_{\text{máx}}\beta W

Con buen tino pueden pensar que hicimos demasiadas aproximaciones como
para que esto sea verdad. Vamos a ver el caso general. Queremos
estudiar específicamente cuál es el error relativo que obtenemos en la
termodinámica al aproximar :math:`Z(N, V, T)` por
:math:`\zeta_{\text{máx}}(N, V, T)`. Como ya dijimos, en la
termodinámica nos interesan los logaritmos, por lo que queremos
calcular entonces

.. math::
   \Delta = \frac{\log{Z} - \log{\zeta_{\text{máx}}}}{\log{Z}}
   :label: delta

Partimos entonces de la expresión de la función de partición

.. math::
   Z(N, V, T) = \sum_{n=0}^{N} \zeta_n

Como regla general, vale

.. math:: 
   \zeta_{\text{máx}} \leq \sum_{n=0}^{N} \zeta_n \leq
   N\zeta_{\text{máx}}

Tomemos logaritmo de la expresión anterior, y obtenemos

.. math:: 
   
   \log{\zeta_{\text{máx}}} \leq \log{Z} \leq
   \log{N} + \log{\zeta_{\text{máx}}}

Y ahora restamos :math:`\log{\zeta_{\text{máx}}}`, para acercarnos a
la expresión :eq:`delta`:

.. math:: 
   0 \leq \log{Z} - \log{\zeta_{\text{máx}}} \leq \log{N}

Finalmente, nosotros sabemos que estamos tratando sistemas
termodinámicos. De ese modo, sabemos que :math:`\log{Z} = \beta F(N,
V, T)`.  Como :math:`F` es un potencial termodinámico, es homogéneo de
primer orden y vale :math:`F(N, V, T) = N\,F(1, V, T) = N\,f(V,
T)`. Así, dividiendo por :math:`\log{Z}` la expresión anterior tenemos
   
.. math:: 
   0 \leq \frac{\log{Z} - \log{\zeta_{\text{máx}}}}{\log{Z}} leq
   \frac{1}{\beta f(V, T)} \frac{\log{N}}{N}

Usando :eq:`delta`, 

.. math:: 
   0 \leq \Delta leq \frac{1}{\beta f(V, T)} \frac{\log{N}}{N}

y, entonces, en el límite :math:`N\rightarrow\infty`, 

.. math:: 
   0 \leq \Delta &\leq 0\\
   \Delta &= 0

.. [#frenk] De todos modos, para valores típicos de :math:`W =
            1\,\text{eV}`, :math:`T = 1000\,\text{K}` y :math:`\alpha
            \approx 1`, resulta :math:`x \approx 0.001`

.. [#aprox] El razonamiento vale igual para los otros casos, pero las
            cuentas se vuelven engorrosas. Pero pronto vamos a ver que
            no importa exactamente la forma de la función de partición.