Sistema de partículas interactuantes: un modelo resoluble para el estudio en física estadística

Interacting particle system: a solvable model for study in statistical physics

Resumen

Se propone un sistema unidimensional simple de dos partículas interactuantes cuyo comportamiento puede calcularse analíticamente. Se muestra en detalle el cálculo de la función de partición así como también el de cantidades fundamentales como la energía interna y la capacidad calorífica. La estructura espacial del sistema de dos partículas es analizada en términos de las distribuciones que describen la distancia de separaración entre las partículas y la posición de su centro de masa. La evolución del sistema hacia el equilibrio es descrita usando ecuaciones de Langevin en conjunto con las respectivas ecuaciones de Fokker-Planck. Los resultados analíticos son contrastados con los que se obtienen por medio de simulaciones computacionales sencillas basadas en dinámica browniana. En todos los casos el comportamiento del sistema se explica de manera simple y detallada de tal forma que la discusión y el modelo presentado pueden ser usados como un recurso pedagógico en cursos intermedios de mecánica estadística.

Palabras-clave:
Física estadística; ecuaciones de Langevin y Fokker Planck; partículas interactuantes

Abstract

A simple one-dimensional system of two interacting particles whose behavior can be calculated analytically is proposed. The calculation of the partition function is shown in detail as well as that of fundamental quantities such as internal energy and heat capacity. The spatial structure of the two-particle system is analyzed in terms of the distributions that describe the separation distance between the particles and the position of their center of mass. The evolution of the system towards equilibrium is described using Langevin equations in conjunction with the respective Fokker-Planck equations. The analytical results are compared with those obtained using simple computer simulations based on Brownian dynamics. In all cases, the behavior of the system is explained in a simple and detailed manner such that the discussion and the presented model can be used as a pedagogical resource in intermediate courses on statistical mechanics.

Keywords:
Statistical physics; Langevin and Fokker Planck equations; interacting particles

1. Introduction

Uno de los desafíos de la enseñanza de la física estadística reside en encontrar modelos de sistemas físicos que puedan ser resueltos exactamente y que permitan ejemplificar los diversos conceptos teóricos que se quieren transmitir. Un ejemplo de encuentra en el caso de los sistemas de partículas interactuantes en los que, como consecuencia de dichas interacciones, no es posible en general contar con una solución analítica conocida, o en caso de tenerla, frecuentemente está dada por medio de técnicas avanzadas o expresiones matemáticas complejas que no son fáciles de analizar a un nivel introductorio. Si bien las propiedades de un modelo general pueden ser evaluadas a través de simulaciones numéricas usando, por ejemplo, los métodos de dinámica molecular o de Monte Carlo, es muy común que los estudiantes de física requieran enfrentarse a un buen número de problemas concretos para para afianzar y dominar los conceptos y las técnicas de solución. De esta manera, el aprendizaje basado en problemas es considerado un método eficiente de enseñanza en este contexto [¹,²,³] y ha sido el enfoque utilizado en varios libros de texto de física [⁴,⁵,⁶,⁷,⁸].

En este trabajo proponemos un sistema modelo sencillo pero no trivial con solución exacta, que permite transmitir conceptos fundamentales relacionados con sistemas de partículas interactuantes tales como la aproximación al equilibrio, la interacción de las partículas con un reservorio, las correlaciones espaciales de las partículas debidas a las interacciones, entre otros. El modelo propuesto está pensado para introducir a los estudiantes en el estudio teórico, tanto analítico como computacional, de sistemas de partículas interactuantes clásicas dentro y fuera del equilibrio.

El artículo está organizado de la siguiente manera. En la Sec. 2 se presenta el modelo propuesto, en la Sec. 3 se describe el estado de equilibrio del sistema y en la Sec. 4 se plantean las ecuaciones de Langevin y las de Fokker Planck, las cuales se solucionan en el régimen de fricción fuerte. En particular se muestra como el sistema evoluciona hacia el equilibrio. En la Sec. 5 se describe un procedimiento numérico para el estudio de este tipo de sistemas. Finalmente en la Sec. 6 se comparan los resultados numéricos con los analíticos y se proporcionan algunos comentarios finales.

2. Descripción del sistema

Consideremos un sistema unidimensional formado por dos partículas puntuales de masa M restringidas a moverse a lo largo de una línea recta, con posiciones x₁ y x₂. El sistema se encuentra en contacto con un reservorio a temperatura absoluta T, con el cual puede intercambiar energía en forma de calor. Las partículas interactuan entre sí por medio de un potencial logarítmico y con un campo externo que crece cuadráticamente con la distancia al origen. De esta forma, la energía potencial total del sistema de partículas está dado por

(1)$$\mathcal{V}(x_1,x_2)=-B\ln \left|\frac{x_2-x_1}{\xi }\right|+A\left(x_1^2+x_2^2\right),$$

con B y A constantes positivas que determinan la intensidad de la interacción entre las partículas y con el campo externo, respectivamente, mientras que ξ es una constante que determina la escala de longitud de la interacción entre las partículas. Es importante notar que B tiene unidades de energía mientras que A se mide en unidades de energía sobre longitud al cuadrado. De ahora en adelante y sin pérdida de generalidad, tomamos x₂ ≥ x₁, esto es, la partícula 2 siempre estará a la derecha de la partícula 1. Notemos que hay una competencia entre el primer término del potencial que es repulsivo y el último que tiende a confinar a las partículas cerca del origen de coordenadas.

El potencial dado por la Ec. (1) no es arbitrario. Por ejemplo, el potencial armónico puede representar el confinamiento producido por una trampa óptica [⁹,¹⁰,¹¹,¹²], mientras que el término logarítmico está asociado con la solución de la ecuación de Poisson para una partícula cargada en dos dimensiones. Además, éste puede representar la interacción electrostática de dos “varillas” cargadas [⁵, ¹³], de tal manera que esta interacción puede usarse para modelar polielectrolitos y segmentos de ADN, en cuyo caso ξ está relacionado con la longitud de Bjerrum [¹⁴]. Un potencial logarítmico ha sido asociado también con la interacción efectiva entre dos polímeros de estrella, siendo ξ una medida relacionada con el radio de giro de los polímeros [¹⁵]. De esta manera, ξ tiene como rol el establecer la escala de longitud característica del sistema y en los sucesivo se tomará como unidad de longitud (ξ = 1).

Finalmente, debemos recalcar que si bien existen sistemas físicos cuya interacción se puede modelar a través de una interacción logarítmica como la mostrada en la Ec. (1), la principal motivación para la elección de esta forma funcional son sus propiedades matemáticas, las cuales facilitan el cálculo de las magnitudes de interés.

3. Estado de equilibrio

3.1. Función de partición canónica

En ausencia de interacción con el reservorio, el punto de equilibrio del sistema $(x_1^{\text{eq}},x_2^{\text{eq}})$ está dado por el mínimo del potencial $\mathcal{V}$(x₁, x₂) [¹⁶, ¹⁷]. Por lo tanto, derivando el potencial e igualando a cero se obtiene que el punto de equilibrio está dado por las posiciones $x_2^{\text{eq}}=-x_1^{\text{eq}}=\sqrt{B/4A}$. Como es de esperar, las posiciones de equilibrio de las partículas son simétricas con respecto al origen; además, debido al balance entre las intensidades de las contribuciones repulsiva y atractiva del potencial, su valor depende del cociente entre B y A (Ec. (1))

Sin embargo, como consecuencia de la interacción con el reservorio, las variables x₁ y x₂ dejan de ser deterministas y pasan a ser estocásticas. En este caso, el comportamiento del sistema en equilibrio puede obtenerse a través de la función de partición canónica [¹⁸,¹⁹,²⁰,²¹]. La densidad de probabilidad de que el sistema esté en una micro-configuración μ = {p₁, p₂, x₁, x₂} está dada por

(2)$$P_{\mu }=\frac{1}{{\sout{Z}}_2}\exp \left[-\beta \left(\frac{p_1^2+p_2^2}{2M}+\mathcal{V}(x_1,x_2)\right)\right],$$

con ${\sout{Z}}_2\equiv {\sout{Z}}_2\left(\beta ,A,B\right)$ la función de partición del sistema de las dos partículas

(3)$${\sout{\mathit{\text{Z}}}}_2=\int \frac{d\mu }{h^2}\exp \left[-\beta \left(\frac{p_1^2+p_2^2}{2M}+\mathcal{V}(x_1,x_2)\right)\right],$$

en donde p₁ y p₂ representan la cantidad de movimiento de cada una de las partículas, β = 1/(k_BT), h es la constante de Planck, k_B es la constante de Boltzmann. Por simplicidad, en la Ec. (3) se ha usado la notación dμ = dp₁dp₂dx₁dx₂.

Las integrales sobre p₁ y p₂ son gaussianas y pueden realizarse sin dificultad extendiendo los límites de integración sobre todo el eje real. Esto lleva a

(4)$${\sout{\mathit{\text{Z}}}}_2=\mathrm{}\left({\displaystyle \frac{2\pi M}{\beta h^2}}\right){\displaystyle {\int }_{x_2\ge x_1}}\mathit{d}{x}_2\mathit{d}{x}_1{\left({\displaystyle \frac{x_2-x_1}{\xi }}\right)}^{\tilde{B}}\times \exp \left[-\tilde{A}(x_1^2+x_2^2)\right],$$

en donde se han definido $\tilde{A}=\beta A$ y $\tilde{B}=\beta B$. Por otro lado, las integrales sobre las coordenadas se pueden evaluar usando el cambio de variables s = x₁ + x₂ y r = x₂ − x₁. Note que s es proporcional a la posición del centro de masa de las dos partículas, mientras que r es la distancia de separación entre ellas. Con estas nuevas variables, el dominio de integración está dado por r ∈ [0, ∞) y s ∈ (−∞, ∞), de manera que

(5)$${\sout{\mathit{\text{Z}}}}_2=\mathrm{}\left({\displaystyle \frac{\pi M}{\beta h^2}}\right){\displaystyle {\int }_0^{\mathrm{\infty }}}\mathit{d}r{\left({\displaystyle \frac{r}{\xi }}\right)}^{\tilde{B}}\exp \left(-{\displaystyle \frac{1}{2}}\tilde{A}{r}^2\right)\times {\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}ds\exp (-{\displaystyle \frac{1}{2}}\tilde{A}{s}^2),$$

en donde se ha considerado que el Jacobiano del cambio de variable es igual a 1/2.

Finalmente, considerando la definición de la función gamma, $\mathrm{\Gamma }(z)={\int }_0^{\mathrm{\infty }}\mathit{d}t{t}^{z-1}\exp (-t)$, la función de partición se puede escribir como:

(6)$${𝒵}_2=\frac{{\pi }^{3/2}M}{\beta h^2}\frac{2^{\tilde{B}/2}}{{\xi }^{\tilde{B}}}\frac{1}{{\tilde{A}}^{1+\tilde{B}/2}}\mathrm{\Gamma }\left[\frac{\tilde{B}+1}{2}\right].$$

Debido a la naturaleza del potencial logarítmico las partículas no pueden cruzarse para B ≠ 0, por lo que x₂ > x₁. Esta condición hace que las partículas mantengan el ordenamiento inicial y, por lo tanto, el factor 1/N! con N el número de partículas no sea necesario. Por otro lado, dada la naturaleza indistinguible de las partículas, la configuración x₁ > x₂ es equivalente a aquella en la que x₂ > x₁ de tal forma que corresponden al mismo estado. En el caso B = 0, las partículas no interactuan entre sí y la restricción x₂ > x₁ no es válida. En este caso el dominio de integración de las variables x₁ y x₂ comprende todo el eje real lo que impondría un factor 2 adicional en la función de partición; sin embargo, ya que las partículas no se encuentran ordenadas, es necesario incluir el factor 1/2! de tal manera que ambos factores se cancelan entr sí y, por lo tanto, la función de partición calculada es válida para todo B.

Es importante recordar que en principio la función de partición encapsula la información necesaria para evaluar las propiedades termodinámicas de un sistema en equilibrio. Así que contar con una expresión explícita para ésta permite realizar predicciones cuantitativas del comportamiento termodinámico del sistema en diferentes escenarios físicos, así como comparar esas predicciones con resultados experimentales o de simulaciones numéricas. Asimismo, los cálculos presentados en este artículo corresponden a un sistema de dos partículas puntuales pero pueden generalizarse, por ejemplo, para un gas de N moléculas no interactuantes en el que cada molécula diatómica está formada por un par de partículas como las discutidas en este trabajo. En este caso, la función de partición del gas, ${\sout{Z}}_N\left(\beta \right)$, es proporcional a ${\sout{Z}}_2{\left(\beta \right)}^N$, con ${\sout{Z}}_2\left(\beta \right)$ siendo la función de partición de una sola molécula.

3.1.1. Energía interna y capacidad calorífica

Como se mencionó anteriormente, usando la función de partición se puede encontrar el comportamiento termodinámico del sistema. Por ejemplo, la energía interna está dada por ⟨E⟩ = −∂ ln Z/∂β. Para el sistema en cuestión se encuentra

(7)$$⟨E⟩=\frac{2}{\beta }+\frac{\tilde{B}}{2\beta }\left[1+\ln \left(\frac{\tilde{A}{\xi }^2}{2}\right)-{\psi }^{(0)}\left(\frac{\tilde{B}+1}{2}\right)\right],$$

donde ψ⁽⁰⁾ es la primera derivada¹ del término ln Г con respecto a β. Debido a la forma logaritmica del potencial de interacción, la escala de longitud ξ afecta a la energía potencial y a la energía interna del sistema únicamente a través de una constante aditiva. Para el caso más simple, en el que B = 0, la Ec. (7) se reduce a ⟨E⟩ = 2k_B T como es de esperar debido al teorema de equipartición de la energía [¹⁸,¹⁹,²⁰,²¹]. El exceso de energía interna, es decir, la diferencia entre la energía del sistema en comparación con la de caso no interactuante (B = 0) se muestra en la Figura 1. El comportamiento de la energía de exceso no es trivial y depende fuertemente de los parámetros A y B así como también de la temperatura. En el recuadro se puede ver que para valores suficientemente grandes de A, la energía interna en T = 0 no tiene un comportamiento monótono con el parámetro B.

Figura 1
Exceso de energía interna vs temperatura para diferentes valores de B y A = 1. El recuadro muestra el comportamiento de ⟨E(T = 0)⟩ en función de B para diferentes valores de A.

Por su parte, la capacidad calorífica está dada por C_V = ∂⁡⟨E⟩/∂⁡T, lo que conduce a

(8)$$C_V=\frac{k_B}{4}\left[8-2\tilde{B}+{\tilde{B}}^2{\psi }^{(1)}\left(\frac{\tilde{B}+1}{2}\right)\right],$$

en donde ψ⁽¹⁾ es la segunda derivada de ln Г con respecto a β. En la Figura 2 se muestra el comportamiento de la capacidad calorífica para diferentes valores de B; es importante notar que C_V no depende del parámetro de interacción con el campo externo A y en el caso B = 0, es independiente de la temperatura e igual a 2k_B. Además, el comportamiento de C_V para valores grandes de T está dado por:

(9)$$C_V=2k_B-\frac{1}{2}{k}_B\tilde{B}+O\left({\beta }^2\right),$$

de donde se concluye que C_V tiende a 2k_B cuando T → ∞. En este límite, la interacción logarítmica entre las partículas puede despreciarse y el sistema se comporta como un gas ideal confinado a través de un potencial cuadrático, y por lo tanto, la capacidad calorífica es constante.

Figura 2
Comportamiento de la capacidad calorífica en función de la temperatura para diferentes valores de B. En todos los casos se ha usado A = 1.

3.1.2. Entropía

De manera análoga a la energía interna, el comportamiento de la entropía se puede evaluar usando, por ejemplo, la expresión S = −k_Bβ²∂ [ln Z₂(β)/β]/∂β. De esta forma se encuentra

(10)$$S(\beta ;A,B)=S(\beta ;A,0)+k_B\ln \left(\frac{1}{\sqrt{\pi }}\mathrm{\Gamma }\left(\frac{\beta B+1}{2}\right)\right)+\frac{k_B}{2}\beta B\left(1-{\psi }^{(0)}\left(\frac{\beta B+1}{2}\right)\right),$$

donde se ha definido

(11)$$S(\beta ;A,0)=2k_B\left(1-\ln \left(\beta A{\xi }^2\right)\right)+k_B\ln \left({\displaystyle \frac{{\pi }^{2}MA{\xi }^4}{h^2}}\right).$$

Es importante resaltar que, debido a la aproximación clásica que se ha empleado para resolver el problema, al calcular la entropía se encuentra que ésta diverge para T → 0. Esto es evidente en el término S(β; A, 0) cuya contribución logarítmica diverge en el límite β → ∞. Este resultado es inconsistente con lo esperado de acuerdo con las leyes de la termodinámica, puesto que en este límite no hay fluctuaciones térmicas y en ausencia de degeneramiento, solo hay un estado posible para el sistema, lo que implica entropía nula. La Figura 3 muestra el exceso de entropía en función de la temperatura que, como se puede observar, es una función bien comportada.

Figura 3
Exceso de entropía vs temperatura para diferentes valores de B.

3.2. Distribuciones espaciales de las partículas

En ausencia total de interacciones, las partículas pueden ocupar cualquier sitio con la misma probabilidad. Sin embargo, las interacción entre las partículas así como la interacción con el campo externo restringen las configuraciones espaciales posibles para las partículas.

3.2.1. Función de correlación de las partículas

La densidad de probabilidad de encontrar las partículas en las posiciones x₁ y x₂ está dada por

(12)$$𝒫(x_1,x_2)=\int \mathit{d}{p}_1\mathit{d}{p}_2{P}_{\mu }=\frac{2^{1-\frac{\tilde{B}}{2}}{\tilde{A}}^{\frac{\tilde{B}}{2}+1}}{\sqrt{\pi }\mathrm{\Gamma }\left(\frac{\tilde{B}+1}{2}\right)}{\left(x_2-x_1\right)}^{\tilde{B}}\times \exp \left[-\tilde{A}\left(x_1^2+x_2^2\right)\right]$$

Claramente la distribución no es uniforme y depende de una forma no trivial de los parámetros A y B, así como de la temperatura T. Por otro lado, se debe notar que, debido a la interacción entre las partículas, la distribución $𝒫$(x₁, x₂) no se puede factorizar de la forma $\widehat{𝒫}(x_1)\widehat{𝒫}(x_2)$. En otras palabras, las variables x₁ y x₂ no son independientes debido a que están correlacionadas por la interacción logarítmica. Sin embargo, para analizar la estructura formada por el sistema de dos partículas es más conveniente usar distribuciones reducidas. La probabilidad de que la distancia entre las dos partículas sea r se obtiene a partir de

(13)$$p(r)=\frac{1}{{𝒵}_2}\int \mathit{d}\mu P_{\mu }\delta \left(r-(x_2-x_1)\right)=\frac{2^{\frac{1-\tilde{B}}{2}}{\tilde{A}}^{\frac{\tilde{B}+1}{2}}}{\mathrm{\Gamma }\left(\frac{\tilde{B}+1}{2}\right)}{r}^{\tilde{B}}\exp \left(-\frac{1}{2}\tilde{A}{r}^2\right),$$

en donde δ(x) es la función delta de Dirac. El comportamiento de esta distribución se muestra en la Figura 4 para distintos valores de temperatura. La presencia de una interacción repulsiva entre las partículas se refleja en el hecho de que p(r) → 0 cuando r → 0; que este decaimiento sea de la forma de potencia $r^{\tilde{B}}$ es debido al potencial de tipo logarítmico. Además, que p(r) → 0 para valores grandes de r es consecuencia directa del potencial externo que confina a las partículas. En este caso, la forma cuadrática del potencial de confinamiento da lugar a un decaimiento de tipo gaussiano, $\exp \left(-\tilde{A}{r}^2/2\right)$, para distancias grandes. Además, es importante notar que en el caso de B → 0, la distribución p(r) es una gaussiana cuyo valor medio y varianza están determinadas por $\tilde{A}$. Recordemos que este límite corresponde al de un gas ideal confinado por un potencial cuadrático.

Figura 4
Comportamiento de p(r) para diferentes valores de T. En todos los casos se tomó A = 1, B = 1, ξ = 1 y k_B = 1.

La distancia de separación promedio entre las partículas está dada por

(14)$$⟨r⟩={\int }_0^{\mathrm{\infty }}\mathit{d}rrp(r)={\left(\frac{2}{\tilde{A}}\right)}^{\frac{1}{2}}\frac{\mathrm{\Gamma }\left(\frac{\tilde{B}+2}{2}\right)}{\mathrm{\Gamma }\left(\frac{\tilde{B}+1}{2}\right)}.$$

Considerando que $\mathrm{\Gamma }\left(x+s\right)\sim x^s\mathrm{\Gamma }\left(x\right)$, para valores pequeños de la temperatura T, o equivalentemente β ≫ 1, la distancia promedio toma la forma

(15)$$⟨r⟩=⟨x_2-x_1⟩=\sqrt{\frac{B}{A}}\left[1+\frac{1}{4\tilde{B}}+\frac{1}{2{(4\tilde{B})}^2}\right]+𝒪\left({\beta }^{-3}\right).$$

Note que en el límite T → 0 (β → ∞) la expresión anterior se reduce a $〈r〉=x_2^{eq}-x_1^{eq}=\sqrt{B/A}$. En este límite las fluctuaciones térmicas desaparecen y se recupera el resultado determinista encontrado minimizando el potencial de interacción. Por su parte, en el límite T → 0 la distribución p(r) tiende a $\delta (r-\sqrt{B/A})$. Como se observa en la Figura 4, al disminuir T, la distribución p(r) tiende a una delta de Dirac marcando la desaparición de las fluctuaciones térmicas. Por el contrario, a temperaturas altas, p(r) tiende a ser uniforme puesto que la energía térmica es mucho más grande que la asociada al potencial de interacción, de tal forma que las partículas tiene a distribuirse homogéneamente.

El “tamaño” típico del sistema está dado por el valor medio ⟨r⟩. El trabajo isotérmico de compresión, W, realizado sobre el sistema al variar el parámetro A está dado por

(16)$$W=k_{B}T{\int }_{A_1}^{A_2}\frac{\partial \ln {𝒵}_2}{\partial A}\mathit{d}A=\frac{k_{B}T}{2}(2+\tilde{B})\ln \left(\frac{A_2}{A_1}\right)=k_{B}T(2+\tilde{B})\ln \left(\frac{{⟨r⟩}_1}{{⟨r⟩}_2}\right).$$

Así, el trabajo de compresión depende del cociente de las intensidades del campo externo al final (A₂) y al inicio del proceso (A₁). Es importante notar que el trabajo de compresión puede escribirse en términos del valor medio de r puesto que una compresión implíca una disminucion en ⟨r⟩.

Una forma de cuantificar las desviaciones de una variable de su valor medio es a través de las fluctuaciones. Para el caso partícular de la variable r se tiene que

(17)$$⟨\mathrm{\Delta }{r}^2⟩=⟨{(r-⟨r⟩)}^2⟩=⟨r^2⟩-{⟨r⟩}^2=\frac{2}{\tilde{A}}\left(\frac{1+\tilde{B}}{2}-\frac{{\mathrm{\Gamma }}^2\left(\frac{\tilde{B}+2}{2}\right)}{{\mathrm{\Gamma }}^2\left(\frac{\tilde{B}+1}{2}\right)}\right).$$

Como es de esperar, las fluctuaciones dependen de los parámetros A, B y T, como se muestra en la Figura 5. En el límite T → ∞ la Ec. (17) se reduce a

(18)$$⟨\mathrm{\Delta }{r}^2⟩=\frac{B(\pi -2\ln (4))}{\pi A}+\frac{\pi -2}{\pi A\beta }+O\left(\beta B^2\right),$$

lo que implica que ⟨Δr²⟩ crece proporcionalmente con la temperatura. En el límite T → 0, la Ec. (17) se reduce a ⟨Δr²⟩ ≈ k_BT/A, por lo que las fluctuaciones en la variable r no dependen de la interacción logaritmica a bajas temperaturas.

Figura 5
Comportamiento de ⟨Δr²⟩ para diferentes valores de T. En todos los casos se tomó A = 1, ξ = 1 y k_B = 1.

A través de un procedimiento análogo al empleado anteriormente para determinar p(r), se encuentra que la distribución de la variable s, q(s), es

(19)$$q(s)=\sqrt{\frac{\tilde{A}}{2\pi }}\exp \left(-\frac{1}{2}\tilde{A}{s}^2\right),$$

con lo que se puede concluir que la ubicación del centro de masa de las dos partículas sigue una distribución gaussiana con media cero debido a la simetría del potencial externo (ver Figura 6). Es importente notar que la distribución del centro de masa no depende de la interacción entre las partículas, esto es, no depende de B. Lo anterior es de esperar puesto que las fuerzas internas de un sistema no influyen en el comportamiento de su centro de masa. Por el contrario, p(r) depende de A y de B y por lo tanto su forma funcional depende explícitamente del balance entre las fuerzas de repulsión y attracción. A diferencia de lo que ocurre con las variables x₁ y x₂, las variables r y s son independientes entre sí, y la distribución conjunta puede escribirse como el producto de las distribuciones de r y s.

Figura 6
Distribución de la variable s = x₂ + x₁ para diferentes valores de T. En todos los caso se tomó A = 1, B = 1, ξ = 1 y k_B = 1. La línea punteada representa el potencial de confinamiento que ha sido incluida como referencia. El eje izquierdo representa los valores de la distribución q(s) mientras que el derecho el del potencial cuadrático medido en unidades de k_BT.

3.2.2. Distribución de densidad de una partícula

La densidad local del sistema, ρ(x), está relacionada con la probabilidad de encontrar cualquiera de las dos partículas en la posición x. Si bien está cantidad solo depende de una partícula y no de dos como p(r) y q(s), para este problema en particular, la evaluación de ρ(x) tiene un grado mayor de dificultad que las mencionadas anteriormente. Por definición, ρ(x) está dada por

(20)$$\rho (x)=\frac{1}{{𝒵}_2}\int \mathit{d}\mu P_{\mu }\left[\delta (x-x_2)+\delta (x-x_1)\right]=2^{1-\frac{\tilde{B}}{2}}\sqrt{\frac{\tilde{A}}{\pi }}\exp \left(-2\tilde{A}{x}^2\right)ℐ(\tilde{A},\tilde{B},x),$$

en donde

(21)$$ℐ(\tilde{A},\tilde{B},x)=\gamma (\tilde{B}){\int }_0^1\mathit{d}t\frac{t^{\frac{\tilde{B}-1}{2}}}{{(1-t)}^{1+\frac{\tilde{B}}{2}}}\exp \left(\tilde{A}{x}^{2}t\right),$$

y

(22)$$\gamma (\tilde{B})=\frac{\mathrm{\Gamma }\left(\frac{1}{2}\right)}{\mathrm{\Gamma }\left(-\frac{\tilde{B}}{2}\right)\mathrm{\Gamma }\left(\frac{\tilde{B}+1}{2}\right)}.$$

Para el desarrollo de la Ec. (20) se expresó la densidad de probabilidad $𝒫$(x₁, x₂) en términos de las variables r y s y posterioremente, se evaluó una de las integrales. Es importante recalcar que la densidad está normalizada de tal forma que cuenta el número de partículas del sistema, esto es, ∫ρ(x)dx = N_part = 2.

En el caso B = 0, la densidad se reduce a

(23)$$\rho (x)=2\sqrt{\tilde{A}/\pi }\exp (-\tilde{A}{x}^2),$$

mientras que la expresión general dada por la Ec. (20) puede evaluarse numéricamente para B > 0 a fin de determinar el comportamiento de la densidad para diferentes valores de A y T. Los resultados encontrados se muestran en la Figura 7. Independientemente de los valores de A y B, en el límite de temperaturas altas (β → 0), la densidad tiende a ser uniforme alrededor del origen debido a que la energía térmica es mucho mayor que la energía de interacción. De hecho, cerca del origen $\tilde{A}(\tilde{B}-1)x^2\ll 1$ y se tiene que

(24)$$\rho (x)\approx \sqrt{\frac{\tilde{A}}{\pi }}{2}^{1-\frac{\tilde{B}}{2}}\left(\tilde{A}(\tilde{B}-1)x^2+1\right)+O(x^3).$$

Figura 7
Densidad ρ(x) en función de la temperatura T para diferentes valores de los parámetros de interacción A y B.

En el caso de fluctuaciones térmicas pequeñas, para un valor fijo de A, las partículas tienden a ubicarse cada vez mas cerca una de la otra a medida que B disminuye, hasta que eventualmente los dos máximos de la densidad desaparecen dando lugar a un único máximo cuando B = 0. Esto se debe a que B controla la interacción repulsiva que separa las partículas; así, a medida que B disminuye las partículas pueden encontrarse cada vez más cerca la una de la otra. De forma análoga, si T es pequeña, para un valor fijo de B > 0 las partículas estarán ubicadas alrededor de sus posiciones de equilibrio pero cada vez más lejos la una de la otra a medida que A disminuye. En este caso, la fuerza de repulsión domina el comportamiento del sistema manteniendo las partículas separadas una mayor distancia.

La condición para la cual ρ(x) pasa de tener dos a un solo máximo en el origen de coordenadas a medida que la temperatura aumenta, se puede evaluar explícitamente para valores generales de A y B > 0 analizando el cambio de concavidad de ρ(x), es decir, cuando ∂² ρ(x)/∂x² cambia de signo para x = 0. Con la ayuda de la Ec. (24) se puede verificar que para B > k_BT la densidad ρ(x) tiene dos máximos. De lo contrario, la energía asociada a la interacción logarítmica es menor que la asociada a las fluctuaciones térmicas y por lo tanto ρ(x) solo tiene un máximo. En el caso marginal B = k_BT la densidad es localmente plana alrededor de x = 0.

4. Aproximación al Equilibrio

Para estudiar la evolución del sistema hacia el equilibrio es necesario describir más detalladamente la interacción de las dos partículas con el reservorio. Una aproximación posible es modelar el reservorio como un solvente a temperatura T, el cual está constituido por partículas puntuales de masa m ≪ M. Las partículas de la suspensión colisionan continuamente con las partículas del sistema (masa M), las cuales pueden ser asociadas a partículas coloidales, dando lugar a pequeños saltos aleatorios en sus trayectorias. Este comportamiento da lugar al movimiento browniano de las partículas del sistema [⁷, ¹⁸, ²²]. La primera evidencia de este fenómeno viene de la observación de partículas de polvo en la superficie del alcohol; sin embargo, el nombre es debido a Robert Brown quien lo estudió en granos de polen sobre el agua [⁷].

4.1. Ecuación de Langevin

Uno de los enfoques más utilizados para estudiar la dinámica de partículas inmersas en un reservorio como el descrito anteriormente está dado por las ecuaciones de Langevin. En este enfoque, se supone que el movimiento de las partículas microscópicas (masa m) es mucho más rápido que el de las partículas mesoscópicas (masa M), de tal forma que se puede separar la dinámica de ambos sistemas. Además, la interacción entre la suspensión y las partículas se modela por medio de una fuerza estocástica y una fuerza de fricción. De esta forma, el solvente se modela de forma implícita usando una aproximación de fondo contínuo que simplifica los cálculos numéricos.

De acuerdo con las consideraciones anteriores, la evolución temporal de las partículas bajo la dinámica de Langevin está dada por

(25)$$M{\displaystyle \frac{d^2{x}_1}{d{t}^2}}=-{\displaystyle \frac{B}{x_2-x_1}}-2A{x}_1-\gamma {\displaystyle \frac{d{x}_1}{dt}}+{\eta }_1(t)$$$\mathrm{}M{\displaystyle \frac{d^2{x}_2}{d{t}^2}}$$=+{\displaystyle \frac{B}{x_2-x_1}}-2A{x}_2-\gamma {\displaystyle \frac{d{x}_2}{dt}}+{\eta }_2(t),$

donde γ es el coeficiente de fricción y η(t) es la fuerzas estocásticas aplicada por el solvente sobre cada partícula. Esta fuerza satisface las condiciones [¹⁸, ²²]

(26)$$⟨\eta (t)⟩=0\text{y}⟨\eta (t)\eta (t^{\prime })⟩=g\delta (t^{\prime }-t).$$

La primera condición garantiza que la fuerza sea isotrópica mientras que la segunda considera que las fuerzas aplicadas en cada una de las diferentes colisiones sean independientes entre sí. Por su parte, la constante g, que indica la intensidad de la fuerza estocástica, está relacionada con la intensidad de la fuerza viscosa a través de g = 2γk_BT [⁷, ²⁰, ²²]. Los dos primeros términos de las Ecs. (25) se calculan por medio del gradiente del potencial dado por la Ec. (1) y corresponden a la interacción entre las partículas y con el campo externo, respectivamente.

Para avanzar analíticamente es conveniente considerar el régimen de fricción fuerte. En este régimen, las partículas alcanzan rápidamente su velocidad terminal y por tanto, sus aceleraciones pueden despreciarse en comparación con los otros términos de las ecuaciones de Langevin. Matemáticamente, esto permite describir el movimiento de las partículas usando dos ecuaciones diferenciales de primer orden en lugar de dos ecuaciones de segundo orden. Para desacoplar el sistema de ecuaciones (25) se introducen nuevamente las variables r = x₂−x₁ y s = x₂ + x₁. De esta manera, en el régimen de fricción fuerte las nuevas variables r y s evolucionan de acuerdo con

(27)$$\gamma \frac{dr}{dt}=\frac{2B}{r}-2Ar+{\eta }_1(t)-{\eta }_2(t),$$

y

(28)$$\gamma \frac{ds}{dt}=-2As+{\eta }_1(t)+{\eta }_2(t),$$

respectivamente.

La Ec. (28) es la misma ecuación de un oscilador armónico estocástico sobreamortiguado y la función s(t) puede ser encontrada facilmente multiplicando por el factor integrante exp⁡(2At/γ). Esto da como resultado

(29)$$s(t)=\mathrm{}s(0)\exp \left(-{\displaystyle \frac{2At}{\gamma }}\right)+\exp \left(-{\displaystyle \frac{2At}{\gamma }}\right){\displaystyle {\int }_0^t}{\displaystyle \frac{\eta (t^{\prime })}{\gamma }}\exp \left({\displaystyle \frac{2A{t}^{\prime }}{\gamma }}\right)\mathit{d}{t}^{\prime },$$

en donde η(t) = η₁(t) + η₂(t) satisface las condiciones (26) con una una intensidad equivalente al doble de la original (g). Tomando el promedio sobre el ensamble en la Ec. (29) y tomando en cuenta la primera condición de la Ec. (26), se obtiene que

(30)$$⟨s(t)⟩=s(0)\exp \left(-\frac{2At}{\gamma }\right).$$

Por lo tanto, el sistema tiende exponencialmente al equilibrio con un tiempo típico τ_s = γ/(2A), también llamado tiempo de relajación. El estado de equilibrio se recupera tomando el límite t → ∞ lo que implica que ⟨s⟩→0. Este resultado coincide con lo encontrado en la sección anterior para la variable s puesto que implica que, en promedio, las dos partículas están situadas simétricamente alrededor del origen.

La solución para la variable r presenta un grado mayor de dificultad debido al término proporcional a 1/r, asociado a la interacción entre las partículas. Tomando el promedio en la Ec. (27) y considerando el límite t → ∞,

se encuentra que la relación

(31)$$⟨\frac{B}{r(t\to \mathrm{\infty })}⟩=⟨Ar(t\to \mathrm{\infty })⟩.$$

es satisfecha. La Ec. (31) implica que, en el equilibrio, el valor medio de la fuerza de repulsión entre las partículas es igual al promedio de la componente restauradora asociada al campo externo cuadrático.

Si bien la Ec. (27) no puede solucionarse analíticamente, es posible determinar la evolución temporal de desviaciones pequeñas alrededor de la distancia de equilibrio r^eq. Para este fin se define r = r^eq + δr, con δr una pequeña desviación de la distancia. Usando esta definición en la Ec. (27) y tomando en cuenta que $r^{\text{eq}}=\sqrt{B/A}$, se encuentra que

(32)$$\gamma \frac{d\delta r}{dt}=-4A\delta r+{\eta }_1(t)-{\eta }_2(t),$$

cuya solución está dada por

(33)$$\delta r(t)=\mathrm{}\exp \left(-{\displaystyle \frac{4At}{\gamma }}\right)\delta r(0)+{\displaystyle \frac{1}{\gamma }}\exp \left(-{\displaystyle \frac{4At}{\gamma }}\right)\times {\displaystyle {\int }_0^t}d{t}^{\prime }\exp \left({\displaystyle \frac{4A{t}^{\prime }}{\gamma }}\right)({\eta }_1(t^{\prime })-{\eta }_2(t^{\prime })).$$

Tomando el promedio de la expresión anterior se encuentra que

(34)$$⟨\delta r(t)⟩=\delta r(0)\exp \left(-\frac{4At}{\gamma }\right).$$

Claramente la variable r(t) = r^eq + δr(t) decae exponencialmente hacia el equilibrio con una tasa dos veces mayor que la variable s(t), esto es, τ_r = γ/(4A) = τ_s/2.

La solución para las variables x₁ y x₂ se obtiene directamente de las relaciones anteriores

(35)$$⟨x_1(t)⟩=\frac{1}{2}\left[⟨s(t)⟩-⟨\delta r(t)⟩-r^{\text{eq}}\right],$$

y

(36)$$⟨x_2(t)⟩=\frac{1}{2}\left[⟨s(t)⟩+⟨\delta r(t)⟩+r^{\text{eq}}\right].$$

De esta forma, se satisfacen las condiciones iniciales ⟨x₁(t → 0)⟩ = x_1,0 y ⟨x₂(t → 0)⟩ = x_2,0, mientras que para tiempos grandes se tiene ⟨x₁(t → ∞)⟩ = −r^eq/2 y ⟨x₂(t → ∞)⟩ = r^eq/2.

4.2. Ecuaciones de Fokker Planck

Un enfoque alternativo para estudiar la aproximación del sistema hacia el equilibrio está dado por la ecuación de Fokker-Planck (FP). Como se discute detalladamente en las Refs. [¹⁸, ²²], es posible genéricamente asociar una ecuación de FP a una ecuación de Langevin. A diferencia del enfoque de Langevin, en el que se analiza directamente la evolución temporal de las variables s y r, en el enfoque de FP se determina la evolución de la distribución de probabilidad de dichas variables. De forma general, una variable estocástica x que evoluciona según la dinámica de Langevin satisface la ecuación

(37)$$\frac{dx}{dt}=\mathrm{\Gamma }(x,t)+\mathrm{{\rm Y}}(x,t)\kappa (t),$$

donde Г(x, t) se denomina comúnmente término de arrastre mientras que Υ(x, t) está relacionado con el coeficiente de difusión. Por su parte κ(t) es una variable aleatoria con media cero y varianza uno. La Ec. (37) tiene asociada una distribución dependiente del tiempo, P(x, t), que satisface

(38)$$\frac{\partial P(x,t)}{\partial t}=\left(-\frac{\partial }{\partial x}{D}^{(1)}(x,t)+\frac{{\partial }^2}{\partial x^2}{D}^{(2)}(x,t)\right)P(x,t).$$

Los coeficientes D^(1,2) vienen dados por

(39)$$D^{(1)}(x,t)=\mathrm{\Gamma }(x,t)+{\mathrm{{\rm Y}}}^{\prime }(x,t)\mathrm{{\rm Y}}(x,t),$$

con Υ′(x,t) la derivada de Υ(x,t) con respecto a x, y

(40)$$D^{(2)}=\frac{1}{2}{\mathrm{{\rm Y}}}^2(x,t).$$

Aplicada a la evolución de la variable s descrita por la Ec. (28), la ecuación de Fokker-Planck para la distribución de probabilidad Q(s,t) resulta

(41)$$\frac{\partial Q(s,t)}{\partial t}=\frac{\partial }{\partial s}\left[\frac{2As}{\gamma }Q(s,t)\right]+\frac{g^{\prime }}{2\gamma }\frac{{\partial }^{2}Q(s,t)}{\partial s^2}.$$

Es importante notar que en este caso g′ = 2g/γ puesto que el término estocástico de la ecuación diferencial para s es la suma de η₁(t) y η₁(t). La Ec. (41) se puede resolver mediante separación de variables considerando Q(s, t) = X(s)θ(t), que empleada en la Ec. (41) da lugar a las dos ecuaciones siguientes:

(42)$$\frac{1}{\theta (t)}\frac{d\theta (t)}{dt}=-\lambda ,$$

y

(43)$$\frac{d}{ds}\left[\frac{2As}{\gamma }X(s)\right]+\frac{g^{\prime }}{2\gamma }\frac{d^{2}X(s)}{d{s}^2}=-\lambda X(s),$$

donde el signo negativo que acompaña a los valores propios se elige a conveniencia de tal forma que λ ≥ 0. Los valores propios de la Ec. (43) están determinados por las condiciones de frontera Q(s → ∞, t) → 0 y dQ(s → ∞, t)/ds → 0 y resultan λ_n = 2An/γ, con n ∈ {0, 1, 2, …}; por su parte, las funciones propias correspondientes están dadas por

(44)$$F_n(z=s/\sigma )=\frac{1}{\sigma \sqrt{2\pi n!}}\frac{d^n}{d{z}^n}{e}^{-z^2/2}.$$

Más detalles sobre la solución de este tipo de ecuaciones diferenciales está disponible en [²³, ²⁴]. De esta forma, para la condición inicial Q(s, 0) = δ(s − s₀), la solución de esta ecuación diferencial está dada por [²⁵]

(45)$$Q(s,t)=\mathrm{}\sqrt{2\pi {\sigma }^2}\exp \left[{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{s_0}{\sigma }}\right)}^2\right]\times {\displaystyle \sum _{n=0}^{\mathrm{\infty }}}{F}_n\left({\displaystyle \frac{s}{\sigma }}\right)F_n\left({\displaystyle \frac{s_0}{\sigma }}\right)\exp (-{\displaystyle \frac{2An}{\gamma }}t),$$

en donde se ha definido ${\sigma }^2=g/(2A\gamma )=k_{B}T/A=1/\tilde{A}$, que introduce la dependencia de la temperatura. En el límite t → ∞, se alcanza el estado estacionario cuya distribución de probabilidad, Q_est(s), está determinada por el primer término de la expansión, esto es, para n = 0. De esta manera, la distribución se reduce a

(46)$$Q_{\text{est}}(s)=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left[\frac{1}{2}{\left(\frac{s}{\sigma }\right)}^2\right]=\sqrt{\frac{\tilde{A}}{2\pi }}\exp \left[\frac{1}{2}\tilde{A}{s}^2\right],$$

que es equivalente al resultado obtenido en equilibrio a través del uso de la función de partición (ver Ec. (19)).

Por su parte, la ecuación de Fokker-Planck para la distribución de la variable r, P(r, t), esta dada por

(47)$${\displaystyle \frac{\partial P(r,t)}{\partial t}}=\mathrm{}{\displaystyle \frac{1}{\gamma }}{\displaystyle \frac{\partial }{\partial r}}\left[\left(2Ar-{\displaystyle \frac{2B}{r}}\right)P(r,t)\right]+{\displaystyle \frac{g^{\prime }}{2\gamma }}{\displaystyle \frac{{\partial }^{2}P(r,t)}{\partial r^2}},$$

que al igual que en el caso de Q(s, t), se puede resolver mediante separación de variables pero ahora con las condiciones de frontera P(r → 0, t) = P(r → ∞, t) = 0, es decir, las dos partículas no pueden estar juntas ni infinitamente separadas. Es impotante notar que la primera condición de frontera es válida únicamente si B ≠ 0. La evolución temporal está dada por la Ec. (42) mientras que el comportamiento asociado a la variable r satisface la ecuación

(48)$$\left[\frac{\partial }{\partial r}\left(Ar-\frac{B}{r}\right)+\frac{g^{\prime }}{4}\frac{{\partial }^2}{\partial r^2}\right]X(r)=-\frac{\gamma \lambda }{2}X(r).$$

Al igual que antes, λ_n = 4An/γ y n ∈ {0, 1, 2, …}, en tanto que las funciones propias son proporcionales a

(49)$$z^{\mu +\frac{1}{2}}{e}^{-z}{L}_n^{(2\mu )}(z),$$

en donde z = r²/(2σ²) y $4\mu =\tilde{B}-1$. Las funciones $L_n^{(2\mu )}(z)$ corresponden a los polinomios generalizados de Laguerre, que están definidos a través de la formula de Rodrigues [²⁵]. La solución de la Ec. (47) para la condición inicial P(r, 0) = δ(r − r₀) está dada por

(50)$$P(r,t)=\mathrm{}{2}^{\frac{1-\tilde{B}}{2}}{\tilde{A}}^{\frac{\tilde{B}+1}{2}}{e}^{-\frac{\tilde{A}{r}^2}{2}}{r}^{\tilde{B}}\mathrm{\Gamma }\left({\displaystyle \frac{\tilde{B}+1}{2}}\right)\times {\displaystyle \sum _{n=0}^{\mathrm{\infty }}}{\displaystyle \frac{\exp \left(-\frac{4An}{\gamma }t\right)L_n^{\frac{\tilde{B}-1}{2}}\left(\frac{\tilde{A}{r}^2}{2}\right)L_n^{\frac{\tilde{B}-1}{2}}\left(\frac{\tilde{A}{r}_0^2}{2}\right)}{n!\mathrm{\Gamma }\left(\frac{\tilde{B}+1}{2}\right)\mathrm{\Gamma }\left(\frac{\tilde{B}+1}{2}+n\right)}},$$

En el estado estacionario (t → ∞), la distribución de probabilidad está dada por primer término de la suma en la Ec. (50), recuperando la Ec. (13).

5. Simulaciones de Dinámica Browniana

En el método estándar de dinámica molecular el movimiento de las partículas se calcula integrando directamente las Ecs. 25. Sin embargo, en el presente caso se está interesado en analizar el límite de amortiguamiento fuerte, que da lugar a una descripción conocida como dinámica browniana. En este límite, la posición de la i-ésima partícula puede escribirse a primer orden en δt como

(51)$$x_i(t+\delta t)\approx x_i(t)+\frac{\delta t}{\gamma }\left[-\frac{\partial U}{\partial x_i}\right]+\frac{1}{\gamma }{\int }_t^{t+\delta t}{\eta }_i(t^{\prime })\mathit{d}{t}^{\prime }.$$

En este esquema, las posiciones de las partículas en el instante t + δt pueden estimarse a partir de sus posiciones y fuerzas netas aplicadas en el instante t. Este método tiene una precisión muy baja dada la aproximación de orden δt, y el error en la determinación de la posición es de orden δt². Este problema de precisión puede subsanarse usando métodos de integración más sofisticados, pero dada la simplicidad del sistema tratado no es necesario recurrir a métodos de integración más sofisticados. La nueva fuerza estocástica ϑ(t) dada por

(52)$$\vartheta (t)={\int }_t^{t+\delta t}\mathit{d}{t}^{\prime }\eta (t^{\prime }),$$

contiene el efecto neto de la variable η(t) en el intervalo de tiempo δt. Recordemos que las variables x_i(t) y v_i(t) cambian lentamente mientras que η(t) lo hace rápidamente. Para determinar la intensidad de la fuerza ϑ(t) consideremos la cantidad

(53)$$⟨\vartheta {(t)}^2⟩={\int }_t^{t+\delta t}\mathit{d}{t}^{\prime }{\int }_t^{t+\delta t}\mathit{d}{t}^{″}⟨\eta (t^{\prime })\eta (t^{″})⟩={\int }_t^{t+\delta t}\mathit{d}{t}^{\prime }{\int }_t^{t+\delta t}\mathit{d}{t}^{″}2\gamma k_{B}T\delta (t^{″}-t^{\prime })=2\gamma k_{B}T\delta t.$$

De esta forma, la posición de las partículas está dada por:

(54)$$x_i(t+\delta t)\approx x_i(t)+\frac{\delta t}{\gamma }\left[-\frac{\partial U}{\partial x_i}\right]+\sqrt{\frac{2k_{B}T\delta t}{\gamma }}{\vartheta }^{\prime }(t),$$

con ϑ′(t) la nueva variable estocástica que satisface

(55)$$⟨{\vartheta }^{\prime }(t)⟩=0\text{y}⟨{\vartheta }^{\prime }(t){\vartheta }^{\prime }(t^{\prime })⟩=\delta (t-t^{\prime }).$$

Es importante anotar que dadas las propiedades estadísticas de η(t), su integral evaluada para el intervalo δt (Eq. 52) tiene una varianza proporcional a este último. Por lo tanto, el término aleatorio en la solución numérica debe tener una magnitud que escale como $\sqrt{\delta t}$ para mantener las propiedades estadísticas correctas del ruido térmico.

La aproximación descrita anteriormente permite estudiar sistemas en los que no es posible encontrar una solución analítica. Las simulaciones realizadas en este trabajo se ejecutaron tomando k_BT = 1/5, γ = 10, B = 4, A = 1 y δt = 0.0005, y promediando sobre 10⁵ repeticiones. El mínimo número de pasos usados para llegar al equilibrio es de 2×10⁴.

La Figura (8) muestra el promedio de las posiciones de las dos partículas para dos condiciones iniciales diferentes (x₁(0), x₂(0)) = (− 1.5, 1.5) (líneas superior e inferior) y (x₁(0), x₂(0)) = (− 0.7, 0.8) (líneas intermedias). Las líneas corresponden a la solución analítica dada por la Eq. (45), mientras que los puntos a los resultados numéricos. Como se puede apreciar, la concordancia es excelente para los parámetros y condiciones iniciales considerados. Como es de esperar, el estado de equilibrio es independiente de la condicion inicial. Para el conjunto de parámetros y condiciones iniciales, la aproximación dada por la Ec. (34) proporciona una buena estimación de los resultados numéricos obtenidos por medio de simulaciones computacionales, siempre y cuando se satisfaga $\delta r/r^{\text{eq}}=\sqrt{A/B}\delta r\ll 1$. El recuadro de la Figura (8) muestra el error porcentual $E_p=[x_2^{\text{teo}}(t)-x_2^{\text{sim}}(t)]/x_2^{\text{sim}}(t)$ como función de t, el cual es menor al 3% en el periodo de tiempo considerado.

Figura 8
Posición promedio de las partículas como función del tiempo. Los parámetros usados son A = 1, B = 4, γ = 10, k_BT = 1/5. Los puntos corresponden a los resultados de las simulaciones de dinámica browniana mientras que las líneas corresponden a las Ecs. 35 y 36. El recuadro muestra el error porcentual en el cálculo de las trayectorias de la partícula 2 como función del tiempo. La línea negra representa la condición inicial x₂(0) = 1.5 mientras que la roja x₂(0) = 0.8.

En la Figura 9 se muestran los resultados para la distribución Q(s, t) (Ec. (45)) para diferentes valores del tiempo con condiciones iniciales x₁(0) = 0.7 y x₂(0) = 0.8. A medida que el tiempo transcurre, la distribución se desplaza hacia la izquierda hasta que el valor promedio de la variable s es igual a cero. De igual forma, la varianza de la distribución aumenta hasta alcanzar el valor de equilibrio. Por su parte, la Figura 10 muestra el comportamiento de P(r, t) para la misma condición inicial. La distribución evoluciona hasta que el sistema alcanza el equilibrio en donde ⟨r⟩ = 2. Nuevamente, los resultados numéricos concuerdan con los analíticos. Esto resalta no solo la importancia de las ecuaciones de Langevin para describir la aproximación al equilibrio en un sistema de partículas interactuantes acopladas a un baño térmico sino que también pone de manifiesto que el método de simulación propuesto a pesar de ser simple, permite una fácil implementación y una descripción adecuada de la dinámica del sistema.

Figura 9
Distribución de la variable s para diferentes tiempos. Los parámetros usados son los mismos de la Figura 8 con condiciones iniciales s(0) = 0.1 y r(0) = 1.5. Los puntos corresponden a los resultados de las simulaciones mientras que las líneas corresponden al resultado de la Eq. (45). El máximo error porcentual encontrado en la región central de la distribución s ∈ (−0.5,05) es menor a 5% para todos los valores de t considerados. Los errores en las colas de la distribución son mayores y pueden disminuirse incrementando el número de repeticiones.

Figura 10
Distribución de la variable r para diferentes tiempos, los parámetros usados son los mismos de la Figura 8 con condiciones iniciales s(0) = 0.1 y r(0) = 1.5. Los puntos corresponden a los resultados de las simulaciones mientras que las líneas corresponden al resultado de la Eq. (50). El máximo error porcentual encontrado en la región central de la distribución, r ∈ (1.5,3.75) para t = 1 y r ∈ (1.1,2.1) para t = 10, es menor a 5%. Los errores en las colas de la distribución son mayores y pueden disminuirse incrementando el número de repeticiones.

Finalmente, es importante recalcar que el tiempo que tarda el sistema en alcanzar el equilibrio depende de los parámetros A y γ, esto es, depende de la intensidad del potencial de confinamiento parabólico y del coeficiente de fricción. Sin embargo, el tiempo de relajación al equilibrio de la variable s es dos veces mayor que el de la variable r. Por otro lado, en el equilibrio la distribución Q(s) es simétrica con respecto al origen debido a la simetría del potencial parabólico; por su parte, P(r) no es simétrica en general.

6. Conclusiones

En este trabajo se propone un sistema sencillo de dos partículas clásicas interactuantes el cual es analíticamente resoluble. Es importante resaltar que para la mayoría de sistemas de partículas interactuantes es necesario recurrir al uso de simulaciones computacionales debido a la ausencia de soluciones analíticas. Sin embargo, se debe tener en cuenta que, aunque las simulaciones de dinámica molecular permiten reproducir detalladamente la evolución de los sistemas, también presentan algunas dificultades que a veces hacen inconveniente su implementación. Una de ellas es su costo computacional, lo que implica que no sea actualmente posible explorar sistemas con muchas partículas en escalas de tiempo grandes. En estos casos usualmente es conveniente el uso de otro tipo de simulaciones como aquellas basadas en Monte Carlo cinético. Estas simulaciones permiten explorar sistemas grandes en escalas de tiempo mucho mayores a las que permite el método de dinámica molecular [²⁶, ²⁷].

Por otro lado, del modelo presentado se pueden resaltar varios aspectos. En primer lugar, con este modelo es posible calcular de forma explícita la función de partición (Ec.(3)), a través de la cual se pueden evaluar propiedades termodinámicas en diferentes escenarios, así como sus propiedades estructurales. Además, permite predecir la respuesta del sistema en el equilibrio ante cambios de los parámetros que lo definen, por ejemplo, la temperatura, las intenisdades de las interacciones, el acomplamiento con el baño térmico, etc. El sistema propuesto es sencillo pero no trivial y sirve como un ejemplo concreto para comprender algunos principios fundamentales. Además, permite al estudiante familiarizarse con técnicas analíticas y numéricas que son de uso común en la física estadística.

Agradecimientos

El trabajo de D.L.G fue apoyado por la Vicerrectoría de Investigaciones de la Universidad del Valle por medio del proyecto C.I. 71369. M.C. agradece el apoyo de la Vicerrectoría de Ciencia, Tecnología e Innovación de la Universidad Antonio Nariño a través del proyecto interno No. 2024209.

Referencias

• ^[1] P. Heller, R. Keith and S. Anderson, American Journal of Physics 60, 627 (1992).
• ^[2] P. Heller and M. Hollabaugh, American Journal of Physics 60, 637 (1992).
• ^[3] E. Koerfer and B. Gregorcic, Phys. Rev. Phys. Educ. Res. 20, 010105 (2024).
• ^[4] J.D. Wall, Introductory Physics: A Problem-solving Approach (Heath, Massachusetts, 1977).
• ^[5] M. Zahn, Electromagnetic Field Theory: A Problem Solving Approach (Krieger Publishing Company, Florida, 2003).
• ^[6] V.K. Thankappan, Quantum Mechanics: Through Problems (New Age International (P) Limited, New Delhi, 2003).
• ^[7] P.L. Krapivsky, S. Redner and E. Ben-Naim, A Kinetic View of Statistical Physics (Cambridge University Press, Cambridge, 2010).
• ^[8] M. Cahay and S. Bandyopadhyay, Problem Solving in Quantum Mechanics: From Basics to Real-World Applications for Materials Scientists, Applied Physicists, and Devices Engineers (Wiley, New Jersey, 2017).
• ^[9] D.C. Appleyard, K.Y. Vandermeulen, H. Lee and M.J. Lang, American Journal of Physics 75, 5 (2007).
• ^[10] G. Pesce, P.H.T. Jones, O.M. Maragò and G. Volpe, The European Physical Journal Plus 135, 949 (2020).
• ^[11] C. Zhang, J. Muñetón Díaz, A. Muster, D.R. Abujetas, L.S. Froufe-Pérez and F. Scheffold, Nature Communications 15, 1020 (2024).
• ^[12] J. Muñetón-Díaz, C. Zhang, P.M. Priti and F. Scheffold, arXiv:2410.02641 (2024).
• ^[13] D.J. Griffiths, Introduction to electrodynamics (Pearson, London, 2013).
• ^[14] Y. Levin, Reports on Progress in Physics 65, 1577 (2002).
• ^[15] C.N. Likos, Soft Matter 2, 478 (2006).
• ^[16] D. Kleppner and R.J. Kolenkow, An Introduction to Mechanics (Cambridge University Press, Cambridge, 2010).
• ^[17] H. Goldstein, Classical Mechanics (Addison-Wesley, San Francisco, 1980).
• ^[18] L.E. Reichl, A modern course in statistical physics (John Wiley & Sons, New Jersey, 2016).
• ^[19] M. Kardar, Statistical Physics of Particles (Cambridge University Press, Cambridge, 2007).
• ^[20] D. Chandler, Introduction to Modern Statistical Mechanics (Oxford University Press, Oxford, 1987).
• ^[21] K. Huang, Statistical Mechanics (John Wiley & Sons, New Jersey, 1987), 2 ed.
• ^[22] H. Risken, The Fokker-Planck Equation: Methods of Solution and Applications (Springer-Verlag, Berlin, 1989).
• ^[23] A. Zettl, Sturm-Liouville Theory (American Mathematical Society, Rhode Island, 2005).
• ^[24] M. Al-Gwaiz, Sturm-Liouville Theory and its Applications (Springer Undergraduate Mathematics Series, New York, 2008).
• ^[25] R.L. Stratonovich, Topics in the Theory of Random Noise (Gordon and Breach, New York, 1963), v. 1.
• ^[26] M.E. Tuckerman, Statistical Mechanics: Theory and Molecular Simulation (Oxford University Press, Oxford, 2010).
• ^[27] T. Tomé and M.J. Oliveira, Stochastic dynamics and irreversibility (Springer, New York, 2015).

• 1
  Las funciones ψ⁽ⁿ⁾ son comunmente conocidas como funciones Poligama. Como es usual, los paréntesis angulares representan promedio sobre el ensamble.

Editado por

• Editor-Chefe:
  Marcello Ferreira https://orcid.org/0000-0003-4945-3169

Fechas de Publicación

• Publicación en esta colección
  31 Mar 2025
• Fecha del número
  2025

Histórico

• Recibido
  05 Dic 2024
• Revisado
  03 Feb 2025
• Acepto
  21 Feb 2025