Gas de BE en una trampa armónica de potencial.

Qué es lo que hace este notebook.

En una trampa armónica de potencial (simétrica y en 3 dimensiones), el número  de partículas en el ensamble gran canónico puede escribirse como 

n = (∑_(m=0)^∞)(((1/2)(m+1)(m+2))/((z^(-1))Exp[ℏω m β]-1)) 

La energía del oscilador, E = ℏ ω(m+(1/2)), ha sido redefinida para eliminar el término constante, quedando solamente E = ℏ ω m . La suma anterior tiene en cuenta que para cada nivel de energía m hay una multiplicidad (1/2)(m+1)(m+2), tal como ustedes pueden demostrar.

Hemos visto que en el límite termodinámico la suma anterior puede reemplazarse por una integral, que termina siendo una de las funciones (Subscript[g, ν]), separando además el término asociado al nivel fundamental:

n = (1/((z^(.1))-1))+ (((((k T)/(ℏ ω))))^3) (Subscript[g, 3])(z).

En el límite termodinámico ocurre una transición de fase en z = 1. Varias funciones termodinámicas con discontinuas cuando la densidad y la temperatura satisfacen la siguiente condición crítica: 

((((n/(ζ(3)))))^(1/3))= ((k T)/(ℏ ω)).

Recordar que (Subscript[g, 3])(1) = ζ(3).

Ahora bien: si el número de partículas se mantiene finito, no hay una transición abrupta, sino que el cambio cualitativo de las funciones termodinámicas se produce dentro de un intervalo finito de los parámetros. 

En este notebook se estudia cómo se acerca uno a la transición de fase propiamente dicha al aumentar el número de partículas. En especial, se ve como las funciones termodinámicas pasan de un comportamiento suave a uno marcado por cambios abruptos. 

La idea es fijar n, la temperatura y la frecuencia ω y encontrar la fugacidad. A partir de ahí, calcular la energía y el calor específico. Toda la gracia del asunto es que para números  pequeños de partículas uno debe hacer la suma explícitamente, sin pasar a la integral. (En realidad se trunca la suma en un determinado número de términos.)

El ejemplo que figura abajo es con n = 100. Una vez que tengan claro cómo funciona el programa, ustedes pueden aumentar el número de partículas. De todas maneras, abajo de todo incluyo los gráficos que obtuve para n = 100, 1000, ... ,(10^7).

Preliminares

Algunas definiciones: 

- la distribución
- la suma para determinar n
- la energía

Por comodidad defino c = (z^(-1)). [Las celdas como la de abajo se ejecutan con shift+enter o con el enter del teclado numérico.]

(*Esta celda debe ejecutarse*)
dist[c_, m_, kT_, hw_]:= (1/(c (E^(hw m/kT))-1)); (*hw es ℏω*)

normsum[c_, sumlim_, kT_, hw_]:= Sum[dist[c, m, kT,hw]  (((m+2)(m+1))/2), {m, 0, sumlim}]; (*Esta suma es n cuando sumlim=∞. Abajo se fija sumlim para tener precisión suficiente*)

energia[ocupakT_, hw_]:= Sum[hw (j-1) ocupakT[[j]], {j, 1, Length[ocupakT]}] (*ocupakT es un vector con la ocupación de los niveles de energía m*)

Una comparación entre la suma con un número finito de términos  (en rojo), la integral (en verde),  y la integral + la contribución del fundamental (en azul). Lo que se grafica es el número total de partículas en función de c. Puede variarse la temperatura; hw está fijo en 1. El eje vertical  da el número de partículas. Tomo hw = 1.

Para el Mathematica las funciones (Subscript[g, ν])(x) se llaman Polylog[ν, x]. Como el gráfico es logarítmico, corro la variable c para no tener problemas en el cero.

(*Esta celda debe ejecutarse*)
Manipulate[LogLogPlot[{normsum[c+1, sumlim, kT, hw] /. hw-> 1,  ((kT^3)/(hw^3))PolyLog[3, 1/(c+1)]/. hw-> 1,(1/c)+((kT^3)/(hw^3))PolyLog[3, 1/(c+1)]/. hw-> 1}, {c, (10^(-6)), 10}, PlotRange-> All, Frame-> True, ImageSize->600, PlotStyle-> {{Thickness[0.005], RGBColor[1,0,0]}, {Thickness[0.005], RGBColor[0,1,0]}, {Thickness[0.005], RGBColor[0,0,1]}}], {{kT, 5}, 1/100, 100, 1/10, Appearance->{"Open"}}, {{sumlim, 100}, 10, 1000, 10, Appearance->{"Open"}}]

El eje horizontal es 1/z - 1, de manera que el origen corresponde a z = 1. En el eje vertical está graficado el número de partículas, calculado según las distintas aproximaciones.

Cosas para notar: 

1) La curva verde, que es la integral sin el término fundamental, siempre falla cerca cuando la fugacidad tiende a 1. 

2) Suficientemente cerca de z = 1   las otras dos curvas siempre tienden a acercarse. 

3) A bajas temperaturas esas dos curvas difieren apreciablemente lejos de z = 1, ya que el espaciamiento entre niveles de energía es comparable a kT, y entonces no es una buena aproximación transformar la suma en una integral. Por más que suban el número de términos en la suma, las curvas no se aproximan mejor.

4) Si aumentan la temperatura, las curvas se aproximan, pero si la aumentan demasiado, para conservar el acuerdo entre las dos fuciones también necesitan aumentar el número de términos que aparecen en la sumatoria (sumlim), puesto que más y más niveles se van ocupando.

Solución de c en función del número de partículas y de la temperatura

Fijada la temperatura y el número de partículas, es necesario encontrar  c = (z^(-1)). La función definida aquí abajo encuentra esa solución. Las variables que aparecen son: el número n de partículas, el número sumlim de términos con que se calcula la suma, la temperatura kT y el producto hw.

En lugar de resolver para c directamente, resuelvo para x, definido por  c = 1 + (10^x). Esto da mayor velocidad y precisión cuando la fugacidad se acerca a 1.

(*Esta celda debe ejecutarse*)
q[n_, sumlim_, kT_, hw_]:= 1+ Power[10, FindRoot[normsum[1+(10^x), sumlim, kT, hw] == n, {x,-10, 1}, WorkingPrecision-> 12][[1]][[2]]];

Abajo vuelvo a graficar el número de partículas, pero también agrego la solución numérica para la fugacidad, q[n, sumlim, kT, hw], para que vean que funciona. El número de partículas está escrito como (10^M).

(*Esta celda debe ejecutarse*)
Manipulate[LogLogPlot[{normsum[c+1, sumlim, kT, hw] /. hw-> 1,  ((kT^3)/(hw^3))PolyLog[3, 1/(c+1)]/. hw-> 1,(1/c)+((kT^3)/(hw^3))PolyLog[3, 1/(c+1)]/. hw-> 1}, {c, (10^(-6)), 10}, PlotRange-> All, Frame-> True, ImageSize->600, PlotStyle-> {{Thickness[0.005], RGBColor[1,0,0]}, {Thickness[0.005], RGBColor[0,1,0]}, {Thickness[0.005], RGBColor[0,0,1]}}, Axes-> True, AxesOrigin->{q[(10^M), sumlim, kT, 1]-1, (10^M)}], {{kT, 5}, 1/100, 100, 1/10, Appearance->{"Open"}}, {{sumlim, 100}, 10, 1000, 10, Appearance->{"Open"}}, {{M, 3}, 1, 7, 1, Appearance->{"Open"}}]

Un cálculo de prueba.

La temperatura de transición teórica para este sistema en el límite termodinámico es 

(Subscript[T, c])=(( (n^(1/3)) hw)/((Zeta[3])^(1/3)))

El límite termodinámico significa  hacer tender n a infinito y hw a cero, de modo que el producto  n (hw^3) sea finito. Eso da una medida de la densidad dentro de la trampa de potencial. Yo encontré cómodo trabajar con una temperatura de transición con el valor de ((n^(1/3))× hw) fijo en (10^(4/3)), lo que da

(Subscript[T, c]) = 20.2624....

Pares de valores de n y hw que dan esta misma temperatura crítica teórica son {n = 100, hw = (10^(2/3))}, {n = 1000, hw = (10^(1/3))}, ..., {n = (10^M), hw = (10^((4-M)/3))}.

(*Esta celda debe ejecutarse*)
{(((100^(1/3))(10^(2/3)))/((Zeta[3.])^(1/3))),(((1000^(1/3))(10^(1/3)))/((Zeta[3.])^(1/3))),(((10000^(1/3))1)/((Zeta[3.])^(1/3)))}

Como se trata de una función definida en términos de FindRoot, el cálculo de c no es muy rápido. Conviene construir tablas con los valores de c en función de la temperatura o del número de partículas, y usar ListPlot en lugar de Plot.

Abajo figura un ejemplo de cómo construir una tabla para ver la dependencia de c con la temperatura, para un número fijo de partículas. Lejos del límite termodinámico, para hacer más rápidos los cálculos tomo un número de partículas bajo,  n = 100. 

Un poco por prueba y error, lo que he visto es que para tener resultados con varias cifras de precisión, el número de términos en la sumatoria debe estar dado apróximadamente por 

sumlim = 12. Log[10] (kT/(10^((4-M)/3)))


Así, el último término de la sumatoria contribuye a n en una cantidad que es del orden de (10^(-12))n. Para hacer los cálculos más rápidos, pero menos precisos, ustedes pueden bajar un poco el número de términos en la sumatoria.

Entonces, fijado n = 100 y hw = (10^(2/3)), abajo defino una tabla con las temperaturas para las que voy a calcular c. Sabemos que lo interesante ocurre en la vecindad de (Subscript[T, c])= 20.26... así que a la tabla la construyo con T entre Tmin = 1/10, y Tmax = 30, a pasos de 1/10. También construyo una tabla con los valores de sumlim.  La idea es tener calculadas la mayor cantidad de cosas antes de resolver la ecuación para c.

(*Esta celda debe ejecutarse*)
nprueba = 100;
hwprueba = ((10^(4/3))/(nprueba^(1/3)));
Tmin = 1/10;
Tmax =30;
Tpaso = 1/10;

tablekT = Table[kT, {kT, Tmin, Tmax, Tpaso}];

tablesumlim =Table[Max[30,IntegerPart[ 12. Log[10] (kT/hwprueba)]],  {kT, Tmin, Tmax, Tpaso}];
 (*En verdad,un poco arbritrariamente yo tomo como sumlim el máximo entre 30 y12. Log[10] (kT/(10^((4-M)/3))),para evitar considerar un número excesivamente pequeño de términos.Es un lujo que uno puede darse.*)

Ahora calculo los valores de c correspondientes. Los guardo en una tabla a la que llamo tableq. Al construir la tabla voy monitoreando el avance en kT y el número de términos que incluye la suma a cada temperatura. A mí me lleva unos  10 segundos construir toda la tabla.

(*Esta celda debe ejecutarse*)
tableq = Monitor[Table[q[nprueba,tablesumlim[[i]],tablekT[[i]], hwprueba], {i, 1, Length[tablekT]}],{1. tablekT[[i]],  tablesumlim[[i]]}];

El gráfico de la fugacitad en función de la temperatura. La línea vertical señala la temperatura crítica teórica. Para este sistema de pocas partículas no hay una transición definida.

(*Esta celda debe ejecutarse*)
ListPlot[Table[{tablekT[[i]], 1/tableq[[i]]}, {i, 1, Length[tablekT]}], Joined-> True, PlotRange->All, Frame-> True, ImageSize-> 600, PlotStyle-> Thickness[0.005],AxesOrigin-> {20.26, 0}, Axes-> {False, True}]

En el eje horizontal figura la temperatura. En el eje vertical la fugacidad.

Ahora construyo una tabla con los números de ocupación relativos de cada nivel de energía m [recordar que hay una degeneración de (1/2)(m+1)(m+2)]. Tengo que usar los resultados obtenidos antes para las fugacidades.  Para cada temperatura hay que calcular sumlim números de ocupación.

(*Esta celda debe ejecutarse*)
ocupa= Monitor[Table[ ((dist[tableq[[i]], m, tablekT[[i]],hwprueba] (((m+2)(m+1))/2))/nprueba),  {i, 1, Length[tablekT]}, {m, 0,tablesumlim[[i]]}], {Length[tablekT] - i}];

Se trata de una tabla que tiene dos variables: la temperatura y el nivel de energía. Abajo grafico en función del nivel de energía y permito variar la temperatura, que se imprime a la izquierda de la figura.

(*Esta celda debe ejecutarse*)
Manipulate[{1.tablekT[[i]], ListLogLogPlot[ocupa[[i]],  PlotRange-> All, Joined-> True, PlotMarkers->Automatic,ImageSize-> 500]}, {i, 1, Length[ocupa], 1, Appearance-> "Open"}]

F

El eje horizontal es el nivel m de energía. En el eje vertical figura la fracción del total de partículas que tiene cada nivel.

Abajo, grafico la ocupación de un dado nivel en función de la temperatura, empezando con el fundamental. 

Notar el cambio cualitativo entre el nivel fundamental y los siguientes.

(*Esta celda debe ejecutarse*)
Manipulate[ListPlot[Table[{tablekT[[i]], ocupa[[i]][[m]]}, {i, 1, Length[tablekT]}],  PlotRange-> All, Joined-> True,ImageSize-> 500], {m, 1, Length[ocupa[[1]]], 1, Appearance-> "Open"}]

Población vs temperatura para un nivel de energía fijo.

Ahora calculo la energía como función de la temperatura; construyo una tabla de valores {T, U}.

(*Esta celda debe ejecutarse*)
tableU = Table[{tablekT[[j]], energia[ocupa[[j]], hwprueba]}, {j, 1, Length[ocupa]}];

(*Esta celda debe ejecutarse*)
ListPlot[tableU, ImageSize-> 600, Axes-> {False, True}, AxesOrigin-> {20.26, 0}, Frame-> True, Joined-> True, PlotMarkers->None, PlotStyle-> Thickness[0.003]]

La energía en función de la temperatura. La línea vertical representa la temperatura de transición teórica de un sistema con infinitas partículas.

Para calcular el calor específico es necesario interpolar una función suave por los pares de puntos de tableU.

(*Esta celda debe ejecutarse*)
tableUinterpol[x_]:= Evaluate[Interpolation[tableU][x]];

(*Esta celda debe ejecutarse*)
Plot[tableUinterpol[kT], {kT, 1, 30}, ImageSize-> 600, Frame-> True,  PlotStyle-> Thickness[0.003],  Axes-> {False, True}, AxesOrigin-> {20.26, 0}]

Abajo, el gráfico del calor específico:

(*Esta celda debe ejecutarse*)
Plot[tableUinterpol'[kT], {kT, 1, 30},ImageSize-> 600, Frame-> True,  PlotStyle-> Thickness[0.003],  Axes-> {False, True}, AxesOrigin-> {20.26, 0}]

El calor específico vs la temperatura.

Resultados para otros números de partículas

Yo repetí los cálculos anteriores tomando un mayor número de puntos en el intervalo de temperaturas entre 10 y 24  y con n tomando los valores  (10^2), (10^3), ... (10^7). Aquí abajo sólo incluyo los gráficos. Los archivos con los datos son un poco grandes para incluirlos aquí. Toma un buen rato calcular las cosas con un número alto de partículas, ya que hay que incluir un gran número de términos en la sumatoria para n.

La fugacidad. La curva más suave corresponde a N = 100. La más angulosa, a N = (10^7).

(*Esta celda debe ejecutarse*)
grafq ={,,,,,};

Show[grafq, ImageSize-> 600]

Aquí abajo pueden ir viendo cada curva por separado, entre N = 100 y N = (10^7).

(*Esta celda debe ejecutarse*)
Manipulate[{N[(10^(i+1)),1],Show[grafq[[i]], ImageSize-> 600]}, {i, 1, Length[grafq], 1}]

La energía:

(*Esta celda debe ejecutarse*)
grafU ={,,,,,};
Show[grafU, ImageSize-> 600]

Aquí abajo pueden ir viendo cada curva por separado, entre N = 100 y N = (10^7).

(*Esta celda debe ejecutarse*)
Manipulate[{(10^(i+1)),Show[grafU[[i]], ImageSize-> 600] }, {i, 1, Length[grafU], 1}]

Notar como a medida que aumenta el número de partículas el cambio en la curva se vuelve más pronunciado y se produce más cerca de la temperatura crítica teórica.

Ídem para los calores específicos:

(*Esta celda debe ejecutarse*)
grafCV = {,,,,,};
Show[grafCV, ImageSize-> 600]

(*Esta celda debe ejecutarse*)
Manipulate[Show[grafCV[[i]], ImageSize-> 600], {{i,1}, 1, 6,1}]