De la ecuación de Kramers a Smoluchowski: límite sobreamortiguado del movimiento browniano

FACULTAD DE FÍSICA DEPARTAMENTO DE FÍSICA ATÓMICA, MOLECULAR Y NUCLEAR De la ecuación de Kramers a Smoluchowski: límite sobreamortiguado del movimiento browniano UNIVERSIDAD DE SEVILLA Carlos Ríos Monje Curso 2021-2022 Doble Grado en Física y Matemáticas Tutorizado por Antonio Prados Montaño y Carlos Alberto Plata Ramos Resumen El movimiento browniano consiste en el movimiento estocástico de partículas suspendidas en un fluido debido a las colisiones de estas con las partículas del mismo. Su estudio es de gran importancia en áreas de conocimiento donde las fluctuaciones térmicas son relevantes, tales como la materia condensada blanda, la biofísica o los nanodispositivos. En este trabajo realizamos un estudio del movimiento browniano a través de dos enfoques distintos, la ecuación de Fokker-Planck y el enfoque de Langevin. Comenzamos por un breve estudio del régimen subamortiguado descrito por la ecuación de Kramers, en el que demostramos una propiedad interesante cuando el sistema está sometido a una fuerza armónica. A continuación, como punto fundamental de este trabajo, probamos que la ecuación de Smoluchowski, que describe el régimen sobreamortiguado, aparece como un límite bien definido de la ecuación de Kramers, usando para ello un desarrollo perturbativo singular. Seguidamente, nos centramos en el estudio de la partícula browniana en el régimen sobreamortiguado, donde analizamos dos situaciones de confinamiento interesantes: armónico y en una caja. Simulamos estos sistemas y comparamos los resultados con las expresiones analíticas para comprobar la equivalencia entre ambos enfoques. Por último, estudiamos el caso de una partícula browniana en un recinto en expansión. Abstract Brownian motion is the stochastic motion of particles suspended in a fluid due to collisions of these particles with those of the fluid. It is of particular interest in research fields where thermal fluctuations are relevant, such as soft condensed matter, biophysics and nanodevices. In this work, we study the Brownian motion by following two approaches, the Fokker- Planck equation and the Langevin description. We start with a brief study of the underdamped regime, governed by the Kramers equation, where we demonstrate an interesting property of its solutions when the Brownian particle is submitted to a harmonic force. Then, we prove, as the main objective of this work, that the Smoluchowski equation can be obtained as a limit case of the Kramers equation, for which we use a singular perturbation method. Subsequently, we focus on the study of the Brownian particle in the overdamped regime, where two different types of confinement are analyzed: harmonic and in a box. Simulations have been carried out, comparing numerical results with analytic expressions to verify the equivalence between both approaches. Finally, we look into a Brownian particle in a growing domain. 1 Índice general 1. Introducción al movimiento browniano 3 1.1. Historia del movimiento browniano . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2. Enfoque moderno: procesos estocásticos markovianos . . . . . . . . . . . . . 5 1.3. La ecuación de Fokker-Planck . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.4. EnfoquedeLangevin ................................ 9 2. Ecuación de Kramers 12 2.1. Resolución para el caso armónico con condiciones iniciales gaussianas . . . 12 3. Límite sobreamortiguado de la ecuación de Kramers: la ecuación de Smoluchowski 17 3.1. Obtención del límite sobreamortiguado por el método de escalas múltiples 17 4. Partícula browniana confinada en el régimen sobreamortiguado 23 4.1. Ecuación de Smoluchowski para el caso armónico . . . . . . . . . . . . . . . 23 4.1.1. Simulaciones ................................. 25 4.2. Ecuación de Smoluchowski en un dominio finito . . . . . . . . . . . . . . . . 30 4.2.1. Simulaciones ................................. 31 5. Movimiento browniano en un recinto en expansión 35 5.1. Simulación....................................... 38 6. Conclusiones 43 A. Simulaciones: Integración numérica de la ecuación de Langevin 46 B. Hermiticidad del operador ∂v∗[v∗+∂v∗]48 Anexos: Códigos de simulación 50 Partícula browniana sobreamortiguada bajo un potencial armónico . . . . . . . . 50 Partícula browniana sobreamortiguada en un dominio finito . . . . . . . . . . . . 53 Bibliografía 55 2 Capítulo 1 Introducción al movimiento browniano 1.1. Historia del movimiento browniano Las primeras observaciones registradas del movimiento estocástico de partículas suspendidas en un fluido se remontan a poco después de la invención del microscopio [1]. No fue hasta el año 1828, cuando el botánico Robert Brown, quien da nombre a este interesante fenómeno, el movimiento browniano, publicó un estudio detallado de dicha dinámica [2]. En su análisis, Brown observó el movimiento irregular de pequeñas partículas de polen suspendidas en agua. En un principio, Brown intentó razonar este comportamiento argumentando que estas partículas estaban vivas. Sin embargo, tras realizar observaciones similares en diversas situaciones, incluso en el interior de cristales de gran antigüedad, donde la hipótesis carecía de sentido, descartó esta idea. La explicación de este movimiento debía tener, por tanto, un origen puramente físico. Hubo numerosos intentos de explicar el origen de este suceso. Algunos trataron de explicarlo centrando su origen en factores totalmente externos al fluido, como diferencias de temperatura con el entorno. No obstante, fue el físico francés Leon Gouy el que demostró que esta hipótesis debía de ser errónea [3]. El movimiento debía provenir de propiedades intrínsecas del sistema. En 1905, Einstein publicó su famoso artículo “Investigaciones en la teoría del movimiento browniano” [4], en el que al fin se ofrecía una explicación físicamente satisfactoria a este fenómeno. En su desarrollo, obtuvo la ecuación que describía la posición de partículas suspendidas en un líquido debido a los choques con las partículas de este último. Para ello, se basó en dos hipótesis: 1. El movimiento de una partícula es independiente del movimiento de las demás. 3 CAPÍTULO 1. INTRODUCCIÓN AL MOVIMIENTO BROWNIANO 4 2. Escogiendo una escala de tiempo adecuada, el estado de una partícula en un instante solo depende del instante inmediatamente anterior. Bajo estas hipótesis, Einstein obtuvo la ecuación de difusión ∂f ∂t =D∂2f ∂x2,(1.1) donde f(x, t)es la densidad de partículas, y Duna constante característica de la dupla partícula-fluido conocida como coeficiente de difusión. Einstein también obtuvo la expresión del coeficiente de difusión, que relaciona las fluctuaciones de las partículas con γ, el coeficiente de amortiguamiento de la velocidad de las partículas en el fluido, D=kBT γ,(1.2) donde kBes la constante de Boltzmann y Tla temperatura absoluta del sistema. A esta expresión se la conoce actualmente como el teorema de fluctuación-disipación [5]. A través de (1.1), Einstein dedujo además cuál era la evolución del desplazamiento cuadrático medio de las partículas ⟨∆x2⟩=2Dt, (1.3) obteniendo una dependencia lineal en el tiempo, sello distintivo del régimen difusivo. La confirmación de los resultados de Einstein llegó de la mano de Jean Perrin en una serie de experimentos realizados en 1908 [6, 7]. Esto supuso un argumento a favor de la descripción discontinua de la materia, la teoría atómica, que había sido objeto de debate durante mucho tiempo. De hecho, este trabajo, junto a otros relacionados con la teoría discontinua, hicieron a Perrin merecedor del Premio Nobel en 1926. Paralelamente a Einstein, Smoluchowski, y posteriormente Langevin, llegaron a la misma solución utilizando enfoques completamente distintos. La forma que utilizó Langevin es particularmente simple, a la par que intuitiva desde el punto de vista físico [8]. Partiendo de la ecuación que describía el movimiento macroscópico de la partícula, Langevin modeló el efecto de los choques de las partículas con las moléculas del líquido mediante la adición de una fuerza aleatoria, L(t): dx dt =v, (1.4a) Mdv dt =−γv +L(t).(1.4b) donde Mes la masa de la partícula browniana. Con esto, Langevin dedujo la misma expresión que Einstein para el desplazamiento cuadrático medio, salvo por un término exponencial en el tiempo que decaía a 0. Este término podía eliminarse escogiendo una CAPÍTULO 1. INTRODUCCIÓN AL MOVIMIENTO BROWNIANO 5 escala de tiempo adecuada, tal como hacía Einstein. De esta manera, el proceso se podía considerar completamente difusivo y, así, ambos enfoques coincidían. El trabajo de Langevin fue relevante más allá del punto de vista puramente físico, fue el origen de una nueva rama de las matemáticas, las ecuaciones diferenciales estocásticas [9]. Una revisión más detallada acerca del origen del movimiento browniano puede encontrarse en la referencia [10]. En la actualidad, el estudio del movimiento browniano es de especial interés en áreas de investigación donde el efecto de las fluctuaciones debidas a la temperatura no es despreciable. Entre sus muchas aplicaciones, encontramos el diseño de nanodispositivos, como los motores microscópicos [11]. 1.2. Enfoque moderno: procesos estocásticos markovianos Actualmente, el movimiento browniano se enmarca dentro de la denominada teoría de procesos estocásticos [12]. En particular, se modela como un proceso de Markov, en el que la probabilidad de observar cualquier estado solo depende del estado inmediatamente anterior. Es decir, para cualquier conjunto de instantes de tiempo t1<t2<... <tn, se tiene que P(yn, tn∣y1, t1;y2, t2;...yn−1, tn−1)=P(yn, tn∣yn−1, tn−1),(1.5) donde P(yn, tn∣y1, t1;y2, t2;...yn−1, tn−1)denota a la probabilidad condicionada de que en el instante tn, el estado del sistema sea yn, habiendo sido, respectivamente, y1, y2,...,yn−1 los estados en los instantes previos t1, t2,...,tn−1. Si recordamos las hipótesis que había realizado Einstein en su artículo, vemos que realmente lo que había hecho, visto desde un punto de vista actual, era considerar la posición de la partícula browniana como un proceso de Markov. La pregunta natural que surge es: ¿es realmente el movimiento browniano un proceso de Markov? Presentamos aquí una breve discusión acerca de la hipótesis markoviana en la línea de la reflexión más detallada que puede encontrarse en la referencia [12]. Si tenemos una partícula inmersa en un fluido, las colisiones con las partículas de este provocarán su movimiento aleatorio. Si en un instante la partícula tiene una velocidad determinada, será mucho más probable que se den colisiones en la dirección y sentido de esta velocidad que en cualquier otra. De esta forma, la velocidad de la partícula en un instante dependerá únicamente de la velocidad en el instante previo, y no de valores previos. Este hecho conlleva a que la velocidad en el movimiento browniano pueda ser tratada como un proceso de Markov. Sin embargo, en su artículo Einstein trataba como proceso markoviano la posición de CAPÍTULO 1. INTRODUCCIÓN AL MOVIMIENTO BROWNIANO 6 t y y1 t1t2 y3 t3 Figura 1.1: Imagen intuitiva de la ecuación de Chapman-Kolmogorov. Para un proceso markoviano, la probabilidad de transición desde un estado y1en el tiempo t1hasta un estado y3en el tiempo t3puede obtenerse sumando sobre todas los estados posibles y2en un tiempo intermedio t2, ya que la probabilidad de transitar desde (y2, t2)hasta (y3, t3)no depende del instante anterior, esto es, de (y1, t1). la partícula. El motivo por el cual esto es posible es que, experimentalmente, la medida de la posición de la partícula se hacía en intervalos de tiempos muy superiores al tiempo de relajación de la velocidad. De esta forma, la velocidad había alcanzado su distribución de equilibrio y, por tanto, la posición de la partícula en un instante solo dependía de su propio estado en el instante anterior. Hay que destacar que el carácter markoviano de la posición en el movimiento browniano solo se manifiesta en escalas de tiempo que son mucho mayores que el tiempo que tarda la velocidad en alcanzar el equilibrio. Si esto no es así, tal como veremos a continuación, el proceso de Markov sería para (x, v). De manera general, para que un proceso físico pueda ser considerado como un proceso de Markov no solo hay que definir bien las variables que estamos considerando, sino también identificar la escala de tiempo relevante en cada caso. A partir de la definición (1.5) se puede llegar a una ecuación que debe satisfacer la probabilidad de transición de cualquier proceso de Markov, conocida como ecuación de Chapman-Kolmogorov:1 P(y3, t3∣y1, t1)=∫P(y3, t3∣y2, t2)P(y2, t2∣y1, t1)dy2,(1.6) donde t1<t2<t3. La ecuación (1.6) tiene una interpretación física sencilla. Para llegar del estado y1al estado y3, puede hacerse a través de cualquier estado intermedio y2. La probabilidad de dicho suceso será la suma de la probabilidad de todos los caminos intermedios posibles, véase la Figura 1.1. No obstante, la ecuación de Chapman-Kolmogorov no es especialmente manejable para realizar predicciones físicas. Es una ecuación integral y su manipulación operativa es complicada, ya que implica dos intervalos temporales: 1La derivación de esta ecuación se puede encontrar en [12]. CAPÍTULO 1. INTRODUCCIÓN AL MOVIMIENTO BROWNIANO 7 [t1, t2]y[t2, t3]. Es por ello que, considerando uno de los pasos temporales infinitesimalmente pequeño, se reformula de manera equivalente en la forma conocida como ecuación maestra:2 ∂P(y, t) ∂t =∫[W(y∣y′)P(y′, t)−W(y′∣y)P(y, t)]dy′,(1.7) donde W(y2∣y1)es la probabilidad de transición por unidad de tiempo de pasar del estado y1ay2. En principio, la ecuación maestra se deduce para la evolución de la probabilidad condicionada P(y2, t2∣y1, t1). Sin embargo, a partir de esta puede demostrarse que la probabilidad P(y, t)en un tiempo obedece la misma ecuación maestra. La diferencia entre las ecuaciones maestras para P(y2, t2∣y1, t1)yP(y, t)radica únicamente en la condición inicial, P(y2, t2=t1∣y1, t1)=δ(y2−y1)mientras que P(y, t =0)=P0(y), siendo P0(y)la distribución de probabilidades inicial para y. Matemáticamente, esto significa que P(y2, t2=t1∣y1, t1)=δ(y2−y1)es la función de Green del problema. Nótese que el carácter markoviano del proceso hace que conociendo únicamente la probabilidad P(y1, t1) yP(y2, t2∣y1t1)es posible conocer P(y1, t1;y2, t2;...;yn, tn). La ecuación (1.7) también tiene una interpretación física clara, no es más que una ecuación de ganancia-pérdida. La probabilidad de estar en el estado yen el instante taumenta por la probabilidad de pasar del estado y′al estado yy disminuye por la probabilidad de abandonar el estado yhacia otro y′. 1.3. La ecuación de Fokker-Planck La ecuación maestra es una ecuación integrodiferencial para P(y, t). Cuando solo son posibles transiciones entre estados cercanos, esto es, W(y∣y′)es no nulo únicamente para ∣y−y′∣“pequeño” (sobre la escala en la que varía P(y, t)respecto de y), la ecuación maestra puede aproximarse por una ecuación diferencial en derivadas parciales que se conoce como la ecuación de Fokker-Planck. El movimiento browniano se modela mediante la ecuación de Fokker-Planck. La forma general de esta ecuación, en el caso de una variable en una dimensión, es ∂P(y, t) ∂t =−∂ ∂y[A(y)P(y, t)]+1 2 ∂2 ∂y2[B(y)P(y, t)],(1.8) donde A(y)yB(y)>0son funciones reales diferenciables. Al primer término del segundo miembro de (1.8) se le suele llamar término de arrastre, mientras que al segundo se le denomina término de difusión. Es posible interpretar la ecuación de Fokker-Planck como 2De nuevo, la obtención de esta ecuación puede encontrarse en [12]. CAPÍTULO 1. INTRODUCCIÓN AL MOVIMIENTO BROWNIANO 8 una ecuación de continuidad para la probabilidad ∂P(y, t) ∂t =−∂J(y, t) ∂y ,(1.9) donde J(y, t)es la corriente de probabilidad, definida como J(y, t)=A(y)P(y, t)−1 2 ∂ ∂y[B(y)P(y, t)].(1.10) Por una parte, la ecuación de Fokker-Planck (diferencial en derivadas parciales) es más manejable desde un punto de vista práctico que la ecuación maestra (integrodiferencial). Además, en muchas situaciones físicas relevantes las funciones coeficiente A(y)yB(y) puedan relacionarse directamente con la física del problema, a partir de la ecuación que gobierna la evolución de ydesde un punto de vista macroscópico y el hecho de que el sistema debe tender al equilibrio termodinámico en el límite de tiempos grandes. Esto es [12], d dt⟨y⟩=⟨A(y)⟩,l´ım t→∞P(y, t)=Pe(y),(1.11) donde Pe(y)es la distribución de yen el equilibrio, proporcionada por la colectividad canónica si el sistema está en equilibrio a la temperatura T. Como se discute en la referencia [12], este enfoque tiene limitaciones cuando A(y)es una función no lineal de y—esta dificultad no es relevante para el trabajo que presentamos aquí. A la hora de particularizar (1.8) para el caso del movimiento browniano, tenemos dos alternativas, ya que dependiendo de la escala de tiempo y las variables escogidas, el proceso será o no de Markov. Lo haremos, por simplicidad, en el caso de una partícula en una dimensión. Si escogemos una escala de tiempo fina, donde la posición por sí misma no sea un proceso de Markov, tenemos que describir el problema considerando tanto la posición como la velocidad. En este caso, usando la versión bivariante de (1.8), la evolución de P(x, v, t)viene dada por la ecuación conocida como ecuación de Kramers ∂P ∂t +v∂P ∂x +F(x) M ∂P ∂v =γ M ∂ ∂v [vP +kBT M ∂P ∂v ],(1.12) donde F(x)es una fuerza determinista aplicada sobre la partícula. Durante el desarrollo de este trabajo, resolveremos esta ecuación para el caso en el que tengamos un potencial armónico y condiciones iniciales gaussianas tanto para la posición como la velocidad, y comprobaremos que mantienen esta distribución a lo largo del tiempo. Si la velocidad relaja sobre una escala de tiempo más pequeña que la posición y consideramos una escala de tiempo más grande, entonces la posición sí será un proceso CAPÍTULO 2. ECUACIÓN DE KRAMERS 15 zamos las variables. Para las unidades de longitud y velocidad, empleamos los valores de la dispersión de xyven el equilibrio. De esta forma, x=x0x∗, x0=¿ Á Á ÀkBT λ;v=v0v∗, v0=¿ Á Á ÀkBT M;t=t0t∗, t0=x0 v0=¿ Á Á ÀM λ.(2.14) Si definimos el parámetro α=γ √λM ,(2.15) el sistema (2.9) adimensionalizado queda d dt∗⟪X∗⟫=⟪V∗⟫,⟪X∗⟫(0)=µx∗,0,(2.16a) d dt∗⟪V∗⟫=−⟪X∗⟫−α⟪V∗⟫,⟪V∗⟫(0)=µv∗,0,(2.16b) cuya solución es ⟪X∗⟫(t∗)=e−α 2t∗[αµx∗,0+2µv∗,0 √α2−4sinh(√α2−4 2t∗)+µx∗,0cosh(√α2−4 2t∗)],(2.17a) ⟪V∗⟫(t∗)=e−α 2t∗[µv∗,0cosh(√α2−4 2t∗)−2µx∗,0+αµv∗,0 √α2−4sinh(√α2−4 2t∗)],(2.17b) Como es de esperar, los valores medios tanto de la posición como de la velocidad tienden al valor medio predicho por la distribución de equilibrio, en ambos casos 0. Por otro lado adimensionalizando el sistema (2.10) obtenemos d dt∗⟪X∗2⟫=2⟪X∗V∗⟫,⟪X∗2⟫(0)=σ2 x∗,0,(2.18a) d dt∗⟪V∗2⟫=−2⟪X∗V∗⟫−2α⟪V∗2⟫−2α, ⟪V∗2⟫(0)=σ2 v∗,0,(2.18b) d dt∗⟪X∗V∗⟫=⟪V∗2⟫−⟪X∗2⟫−α⟪X∗V∗⟫,⟪X∗V∗⟫(0)=0,(2.18c) Al igual que (2.16), la adimensionalización simplifica la resolución del sistema (2.18). Sin embargo, dados su extensión y su poco valor adicional, se ha decidido no incluir las soluciones de forma explícita en este trabajo. Analizando el comportamiento de las soluciones de (2.18) para tiempos grandes, observamos que ⟪X∗V∗⟫tiende a 0 y que tanto ⟪X∗2⟫como ⟪V∗2⟫tienden a 1, justamente los valores adimensionales de la varianza de la distribución de equilibrio.1 Para finalizar este análisis nos gustaría enfatizar que, en variables adimensionales, los sistemas dependen únicamente de un parámetro, también adimensional, α. Este parámetro 1Nótese que hemos adimensionalizado la posición y la velocidad con sus respectivos valores de desviación estándar en el equilibrio. CAPÍTULO 2. ECUACIÓN DE KRAMERS 16 será de gran importancia en el Capítulo 3, pues el desarrollo perturbativo singular que realizaremos en él se realizará respecto a ε=α−1. Capítulo 3 Límite sobreamortiguado de la ecuación de Kramers: la ecuación de Smoluchowski En sistemas donde el papel de la inercia es pequeño, la escala de tiempo en la que la velocidad relaja al equilibrio es mucho más rápida que la de la posición. De esta manera, el estado de una partícula vendrá descrito únicamente por la función de probabilidad de la posición, que se obtiene marginalizando sobre las velocidades. Estas últimas se suponen, de forma ideal, instantáneamente en equilibrio. Este límite, denominado sobreamortiguado, describe el comportamiento de un gran abanico de sistemas realistas, como por ejemplo las partículas coloidales [14]. El objetivo de este capítulo será obtener la ecuación del régimen sobreamortiguado, la ecuación de Smoluchowski, partiendo de la ecuación de Kramers. Para ello, realizaremos un método perturbativo de escalas múltiples respecto a un parámetro, que está relacionado con el cociente de las escalas de tiempo características de la relajación de la velocidad y la posición. Al orden más bajo, obtendremos la ecuación de Smoluchowski. 3.1. Obtención del límite sobreamortiguado por el método de escalas múltiples Nuestro punto de partida es la ecuación de Kramers (1.12). Por simplicidad, lo hacemos en el caso armónico, es decir, F(x)=−λx. En primer lugar, adimensionalizamos las variables, utilizando para ello la adimensionalización (2.14) introducida en el capítulo anterior. La probabilidad en variables adimensionalizadas está relacionada con la 17 CAPÍTULO 3. DE KRAMERS A SMOLUCHOWSKI 18 probabilidad original con la expresión P∗(x∗, v∗, t∗)dx∗dv∗=P(x, v, t)dxdv. (3.1) Como tanto x0como v0son constantes, usando la linealidad de (1.12), podemos intercambiar P∗(x∗, v∗, t∗)yP(x, v, t)sin problemas. Además, ∂ ∂x =1 x0 ∂ ∂x∗,∂ ∂v =1 v0 ∂ ∂v∗,∂ ∂t =1 t0 ∂ ∂x∗,(3.2) de modo que en variables adimensionales (1.12) pasa a ser 1 t0 ∂t∗P+v0v∗ x0 ∂x∗P−λx∗ M x0 v0 ∂v∗P=γ Mv0 ∂v∗[v0v∗P+kBT Mv0 ∂v∗P].(3.3) donde hemos simplificado la notación de las derivadas parciales y omitido ∗de P. Sustituyendo los valores de x0, v0yt0y simplificando, obtenemos ε[∂t∗P+v∗∂x∗P−x∗∂v∗P]=∂v∗P[v∗P+∂v∗P], ε ≡√λM γ,(3.4) es decir, la ecuación de Kramers adimensional solo depende de un parámetro, ε, en analogía con los sistemas de cumulantes adimensionalizados en el capítulo anterior, que dependían de α=ε−1. Llegados a este punto, estamos en condiciones de aplicar el método de las escalas múltiples para obtener la ecuación en el caso sobreamortiguado, usando εcomo parámetro perturbativo. Nótese que el significado físico de εes claro: es el cociente de dos escalas de tiempo, t1=M γ, que es la escala sobre la que decae al equilibrio la velocidad, yt2=√M λ, que es la escala de las oscilaciones debida a la fuerza armónica. Suponer que εes muy pequeño implicará suponer que t1≪t2, o lo que es lo mismo, que el tiempo de relajación de la velocidad es mucho más pequeño que el tiempo característico de las oscilaciones. Esta será la condición física que nos determinará la validez del desarrollo que realizamos a continuación. Admitimos entonces que ε≪1y definimos dos escalas de tiempo t1=t∗, τ =εt∗,(3.5) que tratamos como independientes. Esto es, P(x∗, v∗, t∗)Ð→P(x∗, v∗, t1, τ),∂ ∂t∗Ð→ ∂ ∂t1+ε∂ ∂τ (3.6) CAPÍTULO 3. DE KRAMERS A SMOLUCHOWSKI 19 y la ecuación (3.4) nos queda ε[∂t1+v∗∂v∗+ε∂τ]P(x∗, v∗, t1, τ)=∂v∗[v∗+∂v∗]P(x∗, v∗, t1, τ).(3.7) La idea del método de escalas múltiples es que tratar t1yτcomo variables independientes nos introduce un “grado de libertad” adicional en la teoría perturbativa, en comparación con un desarrollo perturbativo estándar o regular. Este grado de libertad adicional lo usaremos para que la solución perturbativa no tenga divergencias.1 Ahora, procedemos a desarrollar P(x∗, v∗, t1, τ)en serie de potencias de ε P(x∗, v∗, t1, τ)=P(0)(x∗, v∗, t1, τ)+εP(1)(x∗, v∗, t1, τ)+ε2P(2)(x∗, v∗, t1, τ)+... (3.8) Si sustituimos en (3.7), e igualamos los coeficientes para cada potencia de ε, obtenemos las ecuaciones para cada P(i)(x∗, v∗, t1, τ) ε0∶∂v∗[v∗+∂v∗]P(0)=0,(3.9a) ε1∶∂v∗[v∗+∂v∗]P(1)=[∂t1P(0)+v∗∂x∗P(0)−x∗∂v∗P(0)],(3.9b) ε2∶∂v∗[v∗+∂v∗]P(2)=[∂t1P(1)+v∗∂x∗P(1)−x∗∂v∗P(1)+∂τP(0)],(3.9c) ⋮ Para el desarrollo a orden 0, que es el que queremos, nos basta retener las ecuaciones anteriores hasta orden 2. Procedamos a resolver la ecuación para el orden 0 ∂v∗[v∗+∂v∗]P(0)=0Ô⇒[v∗+∂v∗]P(0)=C(x∗, t∗),(3.10) con Cuna función arbitraria. Multiplicando lo anterior por e1 2v∗2 √2π, integrando respecto a v∗y despejando, obtenemos. P(0)(x∗, v∗, t1, τ)=ϕ(x∗, t1, τ)e−1 2v∗2 √2π,(3.11) donde ϕ(x∗, t1, τ)es otra función arbitraria y donde hemos escogido C(x∗, t∗)=0para que P(0)(x∗, v∗, t1, τ)sea normalizable. Reescribiendo (3.11), la probabilidad para el orden 0 1Un ejemplo elemental de la implementación del método de escalas múltiples puede encontrarse en la referencia [13], donde se aplica a un oscilador no lineal (ecuación de Duffing) y se muestra cómo el efecto esencial de la no linealidad es modificar la frecuencia natural del oscilador [15]. CAPÍTULO 3. DE KRAMERS A SMOLUCHOWSKI 20 queda P(0)(x∗, v∗, t1, τ)=ϕ(x∗, t1, τ)PMB(v∗), PMB(v∗)≡e−1 2v∗2 √2π,(3.12) donde PMB(v∗)es la distribución de equilibrio de velocidades (Maxwell-Boltzmann) en variables adimensionales. A priori, no obtenemos ninguna condición sobre ϕ(x∗, t1, τ), más allá de que tiene que estar normalizada al integrar sobre x(nótese que PMB(v∗)está normalizada al integrar sobre las velocidades). Tenemos que recurrir al siguiente orden para obtenerla. Antes de continuar, definimos el producto escalar ⟨φ∣ψ⟩=∫+∞ −∞ dv∗φ(v∗)ψ(v∗) PMB(v∗).(3.13) Este producto escalar, con “función peso” 1/PMB es habitual en la teoría de procesos estocásticos markovianos porque hace que el operador que caracteriza la dinámica sea hermítico. Esto es lo que justamente sucede en nuestro problema: si nos fijamos en las expresiones de orden superior en (3.9), vemos que estas siguen una forma común L∣P(i)⟩=∣g(i)⟩,L≡∂v∗[v∗+∂v∗],(3.14) donde ∣g(i)⟩depende de P(j<i). El operador Les hermítico2bajo el producto escalar que acabamos de definir. De la ecuación de orden 0, obtenemos que L∣P(0)⟩=0. Si multiplicamos por la izquierda por ⟨P(0)∣la ecuación de orden 1, utilizando la hermiticidad de L, 0=⟨P(0)∣L∣P(1)⟩=⟨P(0)∣g(1)⟩,(3.15) obtenemos una condición para que la ecuación sea resoluble: la inhomogeneidad ∣g(1)⟩ tiene que ser perpendicular a la solución de la ecuación homogénea ∣P(0)⟩.3Usando la expresión (3.12), 0=⟨P(0)∣g(1)⟩=ϕ∂t1ϕ 1 ⎧⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎩ ∫+∞ −∞ dv∗PMB(v∗)+ϕ∂x∗ϕ : 0 ∫+∞ −∞ dv∗v∗PMB(v∗) +x∗ϕ2  : 0 ∫+∞ −∞ dv∗v∗PMB(v∗),(3.16) 2Véase Apéndice B. 3Esta condición se conoce a veces como alternativa de Fredholm [16]. CAPÍTULO 3. DE KRAMERS A SMOLUCHOWSKI 21 lo que implica ∂t1ϕ=0⇐⇒ϕ(x∗, t1, τ)=ϕ(x∗, τ).(3.17) Por tanto, para que la ecuación de P(1)tenga solución, ϕdebe ser independiente de t1. Con esta información, podemos resolver la ecuación para P(1)(x∗, v∗, t1, τ) ∂v∗[v∗+∂v∗]P(1)= f(x∗, τ) ⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩ [∂x∗ϕ+x∗ϕ]v∗PMB(v∗).(3.18) Aplicando el método de variación de las constantes (el caso homogéneo es la ecuación de P(0)(x∗, v∗, t1, τ)), P(1)(x∗, v∗, t1, τ)=F(x∗, v∗, t1, τ)PMB(v∗).(3.19) Si sustituimos (3.19) en (3.18) e integramos en [−∞, v∗], se obtiene ∂v∗F(x∗, v∗, t1, τ)=−f(x∗, τ)⇒F(x∗, v∗, t1, τ)=−v∗f(x∗, τ)+ψ(x∗, t1, τ),(3.20) donde ψes una función arbitraria. De esta manera, P(1)(x∗, v∗, t1, τ)=[(−x∗ϕ−∂x∗ϕ)v∗+ψ]PMB(v∗).(3.21) De nuevo, no obtenemos condiciones sobre ψ, por lo que hay que recurrir al orden 2. Conociendo (3.12) y (3.21), podemos sustituir en la ecuación de orden 2 LP(2)=[∂t1ψ+v∗2∂x∗(−x∗ϕ−∂x∗ϕ)+v∗∂x∗ψ−x∗v∗ϕ+x∗(x∗ϕ+∂x∗)(1−v∗2)+∂τϕ]PMB(v∗).(3.22) Procedemos como antes para obtener la condición de solubilidad de ψ 0=⟨P(0)∣L∣P(2)⟩=⟨P(0)∣g(2)⟩.(3.23) que conduce a ∂t1ψ=−∂x∗[−x∗ϕ−∂x∗ϕ]−∂τϕ. (3.24) Como el segundo miembro de (3.24) no depende de t1, esta ecuación implica ∂x∗[−x∗ϕ−∂x∗ϕ]+∂τϕ=0,(3.25) Si no fuera cero, integrando (3.24) obtendríamos que ψes lineal en t1y divergiría en el límite t1→∞, lo cual no es físicamente aceptable. Notemos que la ecuación anterior puede reescribirse como ∂τϕ=∂x∗[x∗ϕ+∂x∗ϕ].(3.26) CAPÍTULO 3. DE KRAMERS A SMOLUCHOWSKI 22 Integrando P(0), dado por la ecuación (3.12), sobre las velocidades, es claro que ϕ(x, τ)es la probabilidad marginal de la posición al orden más bajo. La probabilidad marginal ϕsolo evoluciona sobre la escala lenta τ, como expresa la ecuación (3.17). Por tanto, la ecuación (3.26) no es otra cosa que la ecuación de Smoluchowski en variables adimensionales. De esta manera, hemos probado que la ecuación de Smoluchowski no es más que el límite sobreamortiguado de la ecuación de Kramers. Nótese además que hemos obtenido de manera natural la escala de tiempo más “lenta” en la que evoluciona la ecuación de Smoluchowski, que se traduce en la escala de tiempo a considerar para que la posición pueda considerarse un proceso de Markov por sí sola. Capítulo 4 Partícula browniana confinada en el régimen sobreamortiguado En el capítulo anterior, obtuvimos la ecuación de Smoluchowski como el límite sobreamortiguado de la ecuación de Kramers. En este capítulo vamos a resolver la ecuación de Smoluchowski para una partícula confinada. En concreto, consideraremos dos situaciones físicas de interés: (i) confinamiento armónico y (ii) confinamiento en una caja. Para el caso de un confinamiento armónico, ya vimos que la ecuación de Kramers poseía la interesante propiedad de conservar el carácter gaussiano de sus soluciones. Es de esperar que las soluciones de la ecuación de Smoluchowski tengan la misma propiedad, y probaremos que efectivamente esto es así en este capítulo. El caso de una partícula browniana confinada en una caja también admite solución analítica, pero en este caso involucra una suma infinita sobre los modos de relajación del sistema. Para ambos tipos de confinamiento, presentamos simulaciones de la dinámica del sistema para comprobar que nuestras soluciones analíticas coinciden con las obtenidas en la simulación, utilizando para esto último el enfoque de Langevin. 4.1. Ecuación de Smoluchowski para el caso armónico El problema que nos planteamos esta vez es la resolución de la ecuación (1.13) con una fuerza armónica F(x)=−λx y condiciones iniciales gaussianas, ∂P(x, t) ∂t =∂ ∂x ⎡⎢⎢⎢⎢⎢⎣ λx γP(x, t)+D∂P(x, t) ∂x ⎤⎥⎥⎥⎥⎥⎦, P (x, 0)=1 σ0√2πe−(x−µ0)2 2σ2 0.(4.1) De manera análoga a lo que hicimos con la ecuación de Kramers, recurrimos al espacio de momentos, lo que simplifica la resolución. Podemos pasar al problema equivalente para la 23 CAPÍTULO 4. PARTÍCULA BROWNIANA SOBREAMORTIGUADA 24 función característica G(k, t)=∫+∞ −∞ eikxP(x, t)dx (4.2) tomando la transformada de Fourier en (4.1), ∂G(k, t) ∂t =−λk γ ∂G(k, t) ∂k −k2DG(k, t), G(k, 0)=eiµk−1 2σ2k2.(4.3) En consecuencia, la ecuación que obedece Z(k, t)=lnG(k, t), la función generatriz de cumulantes, es ∂Z(k, t) ∂t =−λk γ ∂Z(k, t) ∂k −k2D, Z(k, 0)=iµ0k−1 2σ2 0k2.(4.4) Para obtener las ecuaciones para los cumulantes, realizamos un desarrollo en serie de de potencias de kcon coeficientes (cumulantes) dependientes del tiempo Z(k, t)=∞ ∑ n=1 κn(t)(ik)n n!.(4.5) Igualando los coeficientes para las distintas nobtenemos las ecuaciones de evolución para los cumulantes. Para n=1, se tiene dκ1(t) dt =λ γκ1(t), κ1(0)=µ0.(4.6) Para mayor generalidad, podemos suponer que λ=λ(t), esto es, la “constante elástica” puede depender del tiempo. La solución de (4.6) es µ(t)=κ1(t)=µ0e−1 γ∫t 0λ(t′)dt′ .(4.7) Particularizando para el caso λconstante, (4.7) queda µ(t)=µ0e−λ γt,(4.8) es decir, el valor medio de la posición decae de forma exponencial, propio del régimen sobreamortiguado, al mínimo de potencial. Para n=2podemos escribir dκ2(t) dt +2λ γκ2(t)=2D, κ2(0)=σ2 0.(4.9) De igual forma, podemos considerar que D=D(t), por ejemplo, si la temperatura es CAPÍTULO 4. PARTÍCULA BROWNIANA SOBREAMORTIGUADA 31 atravesando los bordes del dominio. Sabiendo esto, el problema a resolver es ∂P(x, t) ∂t =D∂2P(x, t) ∂x2,(4.22a) J(0, t)=J(L, t)=0⇒∂P(0, t) ∂x =∂P (L, t) ∂x =0,(4.22b) P(x, 0)=P0(x).(4.22c) La condición inicial P0(x)es una distribución de probabilidad arbitraria. La ecuación (4.22) se resuelve con el método de separación de variables. La solución es P(x, t)=1 L+∞ ∑ n=1 Ane−(nπ L)2Dt cos(nπ Lx),(4.23) donde An=2 L∫L 0 P0(x)cos(nπ Lx)dx. (4.24) La forma de la solución consta de una suma infinita de modos normales que decaen exponencialmente a 0, más una constante, 1/L, que es precisamente el valor que tiene la densidad de probabilidad uniforme en el intervalo [0, L]. Esto es razonable desde un punto de vista físico, ya que en el límite de tiempos grandes esperamos que ese sea precisamente el estado de equilibrio en un proceso difusivo en un dominio finito. En particular, consideramos la condición inicial correspondiente a que la partícula se encuentra muy próxima a un extremo de la caja, por ejemplo x=0+, para t=0. Esto es, consideramos P0(x∗)=δ(x∗+).1Para esta condición inicial, la solución general (4.23) se reduce a P(x, t)=1 L+2 L ∞ ∑ n=1 e−(nπ L)2Dt cos(nπ Lx).(4.25) 4.2.1. Simulaciones A continuación, comprobamos nuestras predicciones teóricas mediante la simulación del sistema. Como hicimos anteriormente, comenzamos por adimensionalizar las variables del problema. En este caso, al no tener un potencial actuando sobre el sistema, no podemos adimensionalizar la posición como hacíamos antes. Ahora escogemos x∗=x x0 , x0=L;t∗=t t0 , t0=L2 D.(4.26) 1La notación δ(x∗+)es una abreviatura para δ(x∗−ϵ), con 0<ϵ≪1. Esto evita problemas con la definición de la delta de Dirac en un extremo del intervalo. CAPÍTULO 4. PARTÍCULA BROWNIANA SOBREAMORTIGUADA 32 De esta forma, la ecuación a simular es (4.17), la misma que en el caso anterior, salvo por la ausencia de la fuerza determinista y las condiciones de contorno. Tal como hemos adimensionalizado las variables, la partícula solo tiene permitido su movimiento en el intervalo [0,1]. La implementación de las condiciones de contorno en la simulación es un problema delicado, véanse las referencias [17, 18] para un mayor detalle. Diferentes implementaciones microscópicas pueden dar lugar al mismo comportamiento macroscópico. Por simplicidad, en este trabajo hemos utilizado el siguiente criterio: Si xi>1⇒xi→1−(xi−1). Si xi<0⇒xi→−xi. Realizamos una simulación con 5000 trayectorias en [0, t∗ máx], con t∗ máx =1y un paso de tiempo ∆t∗=10−3. La condición inicial corresponde a la que hemos comentado antes, la partícula se encuentra perfectamente localizada en un punto arbitrariamente cerca de la pared izquierda, esto es, en x=0+, de modo que P0(x∗)=δ(x∗+). De igual forma que en el caso de un confinamiento armónico, obtenemos el valor medio y la varianza (adimensionales) de la posición como función del tiempo en la simulación, para compararlos con sus respectivos valores analíticos. Para calcular estos últimos, partimos de (4.25) adimensionalizada, P∗(x∗, t∗)=1+2∞ ∑ n=1 e−(nπ)2t∗cos(nπx∗).(4.27) La media y varianza analíticas resultan ser µ∗(t∗)=1 2+2∞ ∑ n=1(−1)n−1 n2π2e−(nπ)2t∗,(4.28a) σ2∗(t∗)=1 3+2∞ ∑ n=1(−1)n2 n2π2e−(nπ)2t∗−µ∗2(t∗).(4.28b) Estas soluciones analíticas involucran sumatorios infinitos sobre todos los modos n= 1,2,.... En la práctica, esto significa que para evaluar las soluciones analíticas hay que introducir un “corte” en n,nm´ax, de manera que el sumatorio se hace realmente entre n=1 ynm´ax. Esperamos que el error introducido por este corte disminuya al aumentar nm´ax. En el caso que nos ocupa, tomar nm´ax =40 es más que suficiente para obtener buenos resultados. En la Figura 4.4 podemos observar el comportamiento del valor medio y varianza frente al tiempo adimensionalizado. Al igual que en el caso armónico, podemos hacer el mismo desarrollo para obtener información sobre la distribución µsim y así determinar que en el área sombreada correspondiente a µ±3σµsim deben caer los puntos en un 99.73% de los casos. Notemos que µsim no sigue una distribución gaussiana en el caso que nos ocupa, salvo en el límite de un gran número de trayectorias, en el que podemos invocar CAPÍTULO 4. PARTÍCULA BROWNIANA SOBREAMORTIGUADA 33 0 0.2 0.4 0.6 0.8 1 0 0.02 0.04 0.06 0.08 0.1 Teoría Simulación Figura 4.4: Media (arriba) y varianza (abajo) de la posición adimensional de una partícula browniana en un dominio finito. Hemos promediado sobre NTray =5000 trayectorias de la ecuación de Langevin, integrada con un paso de tiempo ∆t∗=10−3. El sombreado corresponde a un rango de µ±3σµsim , donde σµsim =σ∗/√NTray. CAPÍTULO 4. PARTÍCULA BROWNIANA SOBREAMORTIGUADA 34 Figura 4.5: Distribución de probabilidad de la posición adimensionalizada de una partícula browniana en un dominio finito para (a) t∗=0.01, (b) t∗=0.05, (c) t∗=0.15, (d) t∗=0.80. El número de trayectorias y el paso de tiempo son los mismos que en la Figura 4.4. el Teorema del Límite Central. Tal como cabría esperar teniendo en cuenta la forma de (4.28), el sistema evoluciona hacia valor medio 1/2y varianza 1/12, que son los propios de la distribución uniforme en [0,1]. Asimismo, comprobamos que, en efecto, el sistema evoluciona hacia la distribución uniforme atendiendo a los histogramas de la Figura 4.5. El desarrollo que hemos realizado en esta sección nos será de gran utilidad para el capítulo siguiente, donde trataremos con un dominio finito en expansión, ya que seremos capaces de reducir ese problema al analizado en este capítulo. Capítulo 5 Movimiento browniano en un recinto en expansión En el capítulo anterior, hemos estudiado el movimiento browniano en un dominio finito [0, L]estático, esto es, con Lindependiente del tiempo. Para ello, hemos resuelto analíticamente la ecuación de Smoluchowski en ese caso y comparado estos resultados con la integración numérica de la ecuación de Langevin. En este capítulo, vamos a considerar el movimiento browniano en un recinto en expansión, esto es, consideramos de nuevo un dominio finito [0, L]pero ahora L=L(t), donde L(t)es una función arbitraria del tiempo. El estudio de este tipo de dominios tiene un gran interés en diversas áreas como la biología, donde algunos procesos involucran difusión junto con expansión del medio, o la cosmología, donde el estudio de sistemas difusivos en la expansión del universo puede tener aplicaciones a los rayos cósmicos [19]. Dependiendo del comportamiento del fluido en el que está inmerso la partícula browniana bajo la expansión, es decir, si el fluido no se ve afectado en absoluto por la expansión o si este se ve “arrastrado” por ella, tendremos dos casos bien diferenciados. En general, el problema que debemos resolver puede escribirse como ∂P(x, t) ∂t =−∂J(x, t) ∂x , J(x, t)=F(x) γP(x, t)−D∂P(x, t) ∂x ;(5.1a) J(0, t)=0, J(L(t), t)=˙ L(t)P(L(t), t),(5.1b) P(x, 0)=P0(x).(5.1c) Ahora, el extremo superior L(t)es una función dada del tiempo y el término F(x)incluye tanto las posibles fuerzas externas actuando sobre el sistema como el posible efecto del arrastre debido a la expansión de la caja. Nótese que la condición de contorno en el extremo móvil garantiza la conservación de la probabilidad en el interior de la caja en 35 CAPÍTULO 5. RECINTOS EN EXPANSIÓN 36 expansión, ∫L(t) 0P(x, t)dx =1,∀t.1 Para afrontar este problema de manera más sencilla, realizamos un cambio de variable y=x L(t)L0,(5.2) donde L0=L(0). A diferencia de casos anteriores donde el cambio de variable era independiente del tiempo, aquí debemos actuar con cautela a la hora de obtener la ecuación en la nueva variable. Primero, de la relación ˜ P(y, t)dy =P(x, t)dx obtenemos P(x, t)=˜ P(y, t) L(t)L0.(5.3) Para obtener la ecuación para la nueva variable debemos tener presente qué variables dejamos constante cuando tomamos derivadas parciales. Las relaciones entre estas son (∂ ∂t)x=(∂ ∂t)y+(∂y ∂t )x(∂ ∂y)t ,(5.4) (∂ ∂x)t=(∂y ∂x)t(∂ ∂y)t (5.5) Sustituyendo las relaciones anteriores en (5.1a) tenemos que el problema (5.1) para la nueva variable y∈[0,1]puede escribirse como ∂˜ P(y, t) ∂t =−∂˜ J(y, t) ∂y ,˜ J(0, t)=˜ J(L0, t)=0,˜ P(y, 0)=˜ P0(y),(5.6) donde ˜ J(y, t)=(˜ F(y) γ−˙ L(t) L0 y)L0 L(t)˜ P(y, t)−DL2 0 L2(t)∂˜ P(y, t) ∂y .(5.7) donde ˜ F(y)=F(x). Nótese que el efecto del cambio de variable sobre la corriente es hacerla nula en los extremos de la caja, ya que ahora estos están fijos. Analicemos ahora los dos posibles casos que tenemos, esto es, (i) que el fluido no se vea afectado por la expansión y por tanto no arrastre a la partícula y (ii) que el fluido se expanda “síncronamente” con la pared y arrastre en consecuencia a la browniana inmersa en el mismo. En el caso (i), si asumimos que no existen fuerzas externas, esto es, F(x)=0, la ecuación para ˜ P(y, t)es ∂˜ P(y, t) ∂t =∂ ∂y [˙ L(t) L(t)y˜ P(y, t)+DL2 0 L2(t)∂˜ P(y, t) ∂y ].(5.8) 1Por concreción en el lenguaje, hablamos de un recinto en expansión, pero nuestros resultados son válidos tanto para ˙ L>0(expansión) como ˙ L<0(compresión), mientras que L(t)permanezca positiva. CAPÍTULO 5. RECINTOS EN EXPANSIÓN 37 Si atendemos más en detalle a los términos de (5.8), vemos esta ecuación es equivalente a la de una partícula browniana en una caja de longitud L0sometida a un potencial armónico cuya “constante elástica”, en este caso γ˙ L(t)/L(t), y coeficiente de difusión, DL2 0/L2(t) dependen del tiempo. Este caso, aunque para un dominio infinito, lo resolvimos en el Capítulo 4 con una condición inicial de tipo gaussiana. En el caso (ii), cuando el movimiento de la pared “arrastra” al fluido, tratamos independientemente el movimiento de la partícula debido a su dinámica natural incluyendo las colisiones con el fluido, es decir, el movimiento browniano, y el debido exclusivamente al efecto de arrastre de la expansión. Así, la ecuación de movimiento de la partícula es dx dt =(dx dt )browniano +(dx dt )arrastre .(5.9) El término debido al movimiento browniano ya lo conocemos, es la ecuación de Langevin (1.19) donde asumiremos de nuevo F(x)=0, es decir, difusión pura. Para calcular el de arrastre recurrimos al cambio de variable (5.2). En la variable yno hay arrastre, por lo que deducimos (dx dt )arrastre =˙ L(t) L0 y=˙ L(t) L(t)x, (5.10) es decir, que el efecto del arrastre sobre la partícula es la aparición de una fuerza proporcional a la posición con el mismo sentido que la expansión de la caja. El efecto de esta fuerza en la ecuación de la nueva variable yes el de contrarrestar el primer término del segundo miembro de (5.8). De esta manera, la ecuación para ˜ P(y, t)queda ∂˜ P(y, t) ∂t =˜ D(t)∂2˜ P(y, t) ∂y2,˜ D(t)=DL2 0 L2(t).(5.11) En la variable y, el problema pasa a ser un proceso puramente difusivo donde el coeficiente de difusión depende del tiempo. Tal como comentamos en el párrafo final del capítulo anterior, podemos reducir este problema al caso donde la difusión no depende del tiempo si definimos una nueva escala temporal τcomo ∂t=˜ D(t)∂τ⇒τ=∫t 0 ˜ D(t′)dt′=∫t 0 DL2 0 L2(t′)dt′.(5.12) De esta forma, la ecuación (5.11) queda ∂˜ P(y, t) ∂τ =∂2˜ P(y, t) ∂y2,(5.13) que coincide con (4.22) tomando D=1. La solución a esta ecuación la tenemos en (4.23). CAPÍTULO 5. RECINTOS EN EXPANSIÓN 38 Si deshacemos el cambio de variable (5.12), obtenemos la solución para (5.11). 5.1. Simulación En esta sección simulamos el sistema correspondiente al caso (ii), es decir, donde la expansión de la caja afecta a la partícula “arrastrándola”. En particular, simulamos, partiendo de la distribución de equilibrio, la expansión de una caja unidimensional hasta que esta tenga una longitud que duplique su valor inicial, dejando un tiempo para que el sistema alcance el equilibrio tras esto. Para la expansión, escogemos L(t)de manera que su variación sea suave, utilizando para ello un polinomio de grado 3 en t texp , donde texp es el tiempo que va a estar la caja en expansión. Definir así L(t)nos permitirá controlar de manera directa la expansión a través del tiempo total simulado sin tener que cambiar el polinomio utilizado. Para t=0,L(0)=L0y para t=texp,L(texp)=2L0, valor que permanecerá constante para t>texp. Además, exigimos que tanto t=0como t=texp sean extremos locales.2Con todo esto, L(t)=⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩ L0[1+3(t texp )2−2(t texp )3],0<t≤texp. 2L0, t >texp. (5.14) Tenemos dos formas de simular este problema. Podemos optar por simular directamente la ecuación (5.9) o partir de la difusión pura en una caja y deshacer los cambios de variable tanto en el tiempo (5.12) como en la posición (5.2). Ambas formas deben ofrecer las mismas soluciones, por lo que para mayor simplicidad empleamos la primera. Sin embargo, a la hora de calcular la solución analítica del problema, usamos la segunda, pues ya tenemos la solución de la ecuación de difusión en una caja estática (4.23). Antes de comenzar, conviene adimensionalizar nuestras variables. Escogiendo la adimensionalización x∗=x x0 , x0=L0;t∗=t t0 , t0=L2 0 D,(5.15) la ecuación de Langevin a simular es dx∗ dt∗=˙ L(t∗) L(t∗)x∗+√2η(t∗)(5.16) Realizamos 2500 trayectorias en el intervalo [0, t∗ máx]con t∗ máx =5,∆t∗=10−3yt∗ exp = t∗ máx −2=3. De nuevo, comparamos el valor medio y varianza de la posición adimensional en distintos instantes con los valores analíticos obtenidos de la solución de (5.1). Para obtenerlos, 2Esto es, que ˙ Lse anule en los tiempos inicial y final de la expansión. CAPÍTULO 5. RECINTOS EN EXPANSIÓN 39 partiendo de (4.23) para D=1y˜ P0(y)=1 L0(recordemos que partimos de la distribución de equilibrio para la longitud inicial) y deshaciendo el cambio de variable τobtenemos la solución de (5.11), ˜ P(y, t)=1 L0+2 L2 0 ∞ ∑ n=1 : 0 [∫L0 0cos(nπ L0 y′)dy′]exp[−(nπ)2D∫t 0 dt′ L2(t′)]cos(nπ L0 y)=1 L0 . (5.17) Deshaciendo el cambio de variable y, mediante (5.3), la probabilidad para xes P(x, t)=1 L(t), x ∈[0, L(t)],(5.18) es decir, esperamos que el sistema, conforme se expanda, no abandone su distribución de equilibrio instantánea. Si adimensionalizamos (5.18), P∗(x∗, t∗)=1 L∗(t∗), x∗∈[0, L∗(t∗)],(5.19) donde L∗(t∗)=1+3(t∗ t∗ exp )2−2(t∗ t∗ exp )3 .(5.20) Por tanto, el valor medio y varianza de la posición adimensionales son los correspondientes a la distribución uniforme para cada instante de tiempo, es decir, µ∗(t∗)=L∗(t∗) 2, σ∗2(t∗)=L∗2(t∗) 12 .(5.21) En la Figura 5.1, observamos que la evolución de la posición media y varianza adimensional de las trayectorias simuladas y las expresiones analíticas (5.21) coinciden salvo las fluctuaciones habituales. Esto nos confirma que, en efecto, ambas maneras de resolver el problema ofrecen los mismos resultados. También observamos lo que predice (5.19), que el sistema siempre se encuentra en su distribución de equilibrio instantáneo. Esto se observa claramente en la Figura 5.2 Para finalizar el estudio de este sistema, por analogía con el capítulo anterior, donde simulamos una partícula browniana en un caja de extremos fijos partiendo desde el extremo izquierdo, simulamos este sistema bajo esa misma condición inicial, es decir, ˜ P0(y)=δ(y∗+). Procediendo como antes, usando (4.25) y deshaciendo el cambio de variable (5.12), la solución analítica de este sistema en la variable yes ˜ P(y, t)=1 L0+2 L0 ∞ ∑ n=1 exp[−(nπ)2D∫t 0 dt′ L2(t′)]cos(nπ L0 y),(5.22) CAPÍTULO 5. RECINTOS EN EXPANSIÓN 40 012345 0 0.1 0.2 0.3 0.4 0.5 Teoría Simulación Figura 5.1: Media (arriba) y varianza (abajo) de la posición adimensional de la partícula browniana en un dominio finito en expansión. Hemos promediado sobre NTray =2500 trayectorias con un paso de tiempo ∆t∗=10−3y un tiempo de expansión del recinto t∗ exp =3. El sombreado corresponde a un rango de µ±3σµsim , donde σµsim =σ∗/√NTray. APÉNDICE A. SIMULACIONES 47 Conociendo esto, expresamos (A.3) como x(t+∆t)=x(t)+F(x, t) γ∆t+√2D(t)∆t N(t),(A.5) donde ahora N(t)es una variable aleatoria gaussiana de media nula y varianza unidad.1 Esta última ecuación nos proporciona un método iterativo para poder hallar la solución del problema de manera numérica. Para construir una trayectoria en el intervalo de tiempo [0, tmáx], hallamos la posición en distintos instantes de tiempo equiespaciados por una cantidad ∆t, que denominamos paso de tiempo. Si llamamos xj ial valor de x(ti)en la trayectoria j-ésima, donde ti=i∆t, podemos hallar xj i+1como xj i+1=xj i+F(xj i, ti) γ∆t+√2D(t)∆tNj i,(A.6) donde Nj ies un número aleatorio que extraemos de una distribución gaussiana de media nula y varianza unidad. Este método de resolución numérica para la ecuación de Langevin se conoce como método de Euler. El valor medio de cualquier función de xen cada instante tise obtiene como un promedio sobre trayectorias. Esto es, en la simulación realizamos un número de trayectorias NTray suficientemente grande e identificamos ⟨f(x, ti)⟩sim =1 NTray NTray ∑ j=1 f(xj i, ti).(A.7) A través de (A.7), calculamos la media y varianza de las trayectorias simuladas en función del tiempo. También hemos calculado la distribución de probabilidad a partir de los datos obtenidos en la simulación. Para ello, utilizamos histogramas. El primer paso para construir los histogramas es definir un intervalo de estudio lo suficientemente grande para asegurarnos que la totalidad de los puntos de las trayectorias caen en su interior. A continuación, realizamos una equipartición de este intervalo, de manera que tendremos muchos subintervalos con la misma anchura ∆hist. En cada instante ti, contamos el número de puntos xj ide las distintas trayectorias que caen en cada subintervalo, y tomamos como la altura del histograma en dicho subintervalo este número dividido por NTray∆hist. De esta manera, el área del histograma (sumatorio de la altura de cada subintervalo por la anchura del mismo) está normalizada a la unidad, y su representación nos ofrece la estimación numérica de la distribución de probabilidad buscada. 1El error cometido al escribir (A.5) es de orden (∆t)3/2. Apéndice B Hermiticidad del operador ∂v∗[v∗+∂v∗] En este apéndice probamos que, bajo el producto escalar definido en (3.13), el operador L=∂v∗[v∗+∂v∗](B.1) es hermítico. Recordemos que un operador es hermítico si cumple ⟨φ∣L∣ψ⟩=⟨ψ∣L∣ϕ⟩,∀ψ, φ ∈X. (B.2) donde, en el caso que nos ocupa, Xes el espacio de funciones de cuadrado integrable con peso 1/PMB(v∗), con PMB(v∗)=(2π)−1/2exp(−v2/2), lo suficientemente regulares. Probemos que Les un operador hermítico. Sean φ, ψ funciones arbitrarias de X: tenemos que ⟨φ∣L∣ψ⟩=∫+∞ −∞ φ∂v∗[v∗+∂v∗]ψ PMB dv∗=∫+∞ −∞ φψ +φv∗∂v∗ψ+φ∂2 v∗ψ PMB dv∗.(B.3) Ahora integramos por partes una vez, ⟨φ∣L∣ψ⟩=∫+∞ −∞ ψφ PMB dv∗+∫+∞ −∞ φv∗∂v∗ψ PMB dv∗+∫+∞ −∞ φ∂2 v∗ψ PMB dv∗ =∫+∞ −∞ ψφ PMB dv∗+∫+∞ −∞ φv∗∂v∗ψ PMB dv∗+ * 0 ∂v∗φ∂v∗ψ PMB ∣+∞ −∞ −∫+∞ −∞ ∂v∗(φ PMB )∂v∗ψdv∗. (B.4) A continuación, utilizando la forma explícita de la distribución de Maxwell-Boltzmann, ∂v∗PMB(v∗)=−v∗PMB(v∗),(B.5) 48 APÉNDICE B. HERMITICIDAD DEL OPERADOR ∂V∗[V∗+∂V∗]49 y entonces ⟨φ∣L∣ψ⟩=∫+∞ −∞ ψφ PMB dv∗+  ∫+∞ −∞ φv∗∂v∗ψ PMB dv∗−  ∫+∞ −∞ φv∗∂v∗ψ PMB dv∗−∫+∞ −∞ ∂v∗φ∂v∗ψ PMB dv∗ =∫+∞ −∞ ψφ PMB dv∗− >0 ∂v∗φψ PMB ∣+∞ −∞ +∫+∞ −∞ ∂v∗(∂v∗φ PMB )ψdv∗.(B.6) Finalmente, usando de nuevo (B.5), tenemos que ⟨φ∣L∣ψ⟩=∫+∞ −∞ ψφ +ψv∗∂v∗φ+ψ∂2 v∗φ PMB dv∗=⟨ψ∣L∣φ⟩,(B.7) esto es, Les hermítico. Anexos: Códigos de simulación En este anexo incluimos algunos de los códigos en el programa MATLAB para la simulación de los sistemas estudiados. Partícula browniana sobreamortiguada bajo un potencial armónico 1%% Simulacion numerica de movimiento browniano sobreamortiguado de una particula sometida a un potencial armonico en una dimension. 2% A(x)=-k/( gamma )x, B(x)=2D, donde D= k_BT/ gamma . 3% La ecuacion que describe la posicion de la particula segun el enfoque de 4% Langevin es dx /dt=A(x)+B(x) ^(1/2) * eta (t) , donde eta (t) es ruido blanco 5% gaussiano . 6% Utilizando el metodo de Euler podemos calcular la 7% trayectoria de la particula : 8% x(t+ Deltat )=x(t)+A(x) Deltat +( B(x)* Deltat ) ^(1/2) N(t) , donde ahora N(t) es 9% un numero aleatorio que sigue una distribucion normal de medio 0 y 10 % varianza 1. 11 12 %% Adimensionalizacion 13 % Antes de simular , es conveniente adimensionalizar la ecuacion con la que 14 % trabajaremos: 15 % x*=x/x_o , x_o =raiz(k/( k_BT)), t*=tk /( gamma ) 16 % dx*/ dt*=-x*(t*)+raiz (2) eta(t*) 50 CÓDIGOS DE SIMULACIÓN 51 17 % con x*(0) siguendo una distribucion gaussiana de media mu *(0) y varianza var(x*) (0). 18 % 19 %% Simulacion 20 % Escogemos Delta_t =1e -3. 21 % Simularemos cada trayectoria con 5000 puntos cada 22 % una , es decir , llegando a t*=5. Realizaremos un numero elevado de trayectorias , y para 23 % cada t* hallaremos la media de x*, y compararemos con la expresion 24 % teorica <x*(t*) >=mu *(0) exp(-t*). 25 % Tambien la compararemos con la expresion teorica para la varianza 26 % var (x*) =( -1+ var(x*) (0) exp (-2t*) +1. 27 Delta_t =1e -3; 28 Tmax =5; 29 Trayectorias=5000; 30 NPasos = Tmax / Delta_t ; 31 medias = zeros (1 , NPasos /100+1) ; 32 momento2 = zeros (1 , NPasos /100+1) ; 33 % Histogramas equiespaciados : 34 NHistogramas=100; % Sin contar el inicial . 35 xmin = -15; 36 xmax =15; 37 Anchura =0.25; 38 Bordes = xmin : Anchura : xmax ; % Elegimos estos como " cajas" globales . 39 HistData = zeros ( NHistogramas +1, length ( Bordes ) -1); 40 % Media y desviacion estandar de x* en t=0 41 mu_o =3; 42 sigma =3; 43 for i=1:Trayectorias 44 x_o = normrnd (mu_o , sigma ); 45 medias (1) = medias (1) +x_o ; 46 momento2 (1) =momento2 (1)+x_o ^2; 47 N0= histcounts (x_o , Bordes ); 48 HistData (1 ,:)= HistData (1 ,:) +N0; CÓDIGOS DE SIMULACIÓN 52 49 for j=1: NPasos 50 x_1 =x_o *(1 - Delta_t )+sqrt (2* Delta_t )* normrnd (0 ,1); 51 x_o = x_1; 52 if mod(j ,100) ==0 53 medias ((j) /100+1) = medias (( j) /100+1) +x_1 ; 54 momento2 ((j) /100+1) = momento2 (( j) /100+1) + x_1 ^2; 55 end 56 if mod (j , NPasos / NHistogramas ) ==0 57 if x_o >xmax 58 HistData (j/( NPasos / NHistogramas )+1 , end )= HistData (j /( NPasos / NHistogramas )+1,end)+1; 59 elseif x_o <xmin 60 z=x_o - xmin; 61 HistData (j/( NPasos / NHistogramas )+1 ,1)= HistData (j/( NPasos / NHistogramas ) +1 ,1) +1; 62 else 63 z=x_o - xmin; 64 N= ceil (z/ Anchura ); 65 HistData (j/( NPasos / NHistogramas )+1 ,N)= HistData (j /( NPasos / NHistogramas )+1,N) +1; 66 end 67 end 68 end 69 end 70 medias = medias / Trayectorias ; 71 varianza =( momento2 / Trayectorias - medias .^2) ; CÓDIGOS DE SIMULACIÓN 53 Partícula browniana sobreamortiguada en un dominio finito 1%% Simulacion del movimiento monodimensional de una particula browniana sobreamortiguada en una caja de longitud L debido unicamente a la difusion . 2% La ecuacion que rige este comportamiento es: 3% dx/dt=raiz (2* k_BT /( gamma ))*eta(t) 4% Si adimensionalizamos x*=x/L, t*=t/(L^2/D) 5% dx */ dt*= raiz (2)*eta(t*) 6% Supondremos que la condicion inicial es x*(t*=0) =0. Por encontrarse en 7% una caja , tendremos por condiciones de contorno 0<x<L => 0< x*<L*, donde 8% L*=1). 9 10 11 %% Simulacion 12 % Vamos a simular la ecuacion anterior de la siguiente manera : 13 % x(t*+ Delta_t *)=x(t*)+raiz (2* Delta_t *)N(t) donde N(t)= normrnd (0 ,1) . 14 Delta_t =1e -3; 15 Tmax =1; 16 Trayectorias=5000; 17 NPasos = Tmax / Delta_t ; 18 NPuntos =100; % Sin contar el inicial 19 medias = zeros (1 , NPuntos +1) ; 20 momento2 = zeros (1, NPuntos +1) ; 21 % Media y varianza de x* en t=0 22 x_0=0; 23 NHistogramas=100; 24 Anchura =5e -2; 25 Bordes =0: Anchura :1; 26 HistData = zeros ( NHistogramas +1, length ( Bordes ) -1); % Sin contar el inicial CÓDIGOS DE SIMULACIÓN 54 27 for i=1:Trayectorias 28 x_o = x_0; 29 medias (1) = medias (1) +x_o ; 30 momento2 (1) =momento2 (1)+x_o ^2; 31 N0= histcounts (x_o , Bordes ); 32 HistData (1 ,:)= HistData (1 ,:) +N0; 33 for j=1: NPasos 34 x_1 =x_o + sqrt (2* Delta_t )* normrnd (0 ,1) ; 35 if x_1 >1 36 x_1 =1 -(x_1 -1) ; 37 end 38 if x_1 <0 39 x_1 =-x_1; 40 end 41 x_o = x_1; 42 if mod (j , NPasos / NPuntos ) ==0 43 medias (j /( NPasos / NPuntos ) +1) = medias (j /( NPasos / NPuntos ) +1)+ x_1; 44 momento2 (j/( NPasos / NPuntos )+1) = momento2 (j/( NPasos / NPuntos ) +1) +x_1 ^2; 45 end 46 if mod (j , NPasos / NHistogramas ) ==0 47 N= histcounts (x_o , Bordes ); 48 HistData (j/( NPasos / NHistogramas )+1 ,:)= HistData (j/( NPasos / NHistogramas )+1 ,:)+N; 49 end 50 51 end 52 end 53 medias = medias / Trayectorias ; 54 varianza =( momento2 / Trayectorias - medias .^2) ; Bibliografía [1] S. Gray. «III. Several microscopical observations and experiments, made by Mr. Stephen Gray». En: Phil. Trans. R. Soc. 19 (1696), págs. 280-287. [2] R. Brown. «XXVII. A brief account of microscopical observations made in the months of June, July and August 1827, on the particles contained in the pollen of plants; and on the general existence of active molecules in organic and inorganic bodies». En: The philosophical magazine 4(1828), págs. 161-173. [3] L. Gouy. «Note sur le movement brownien». En: Journal de Physique 7(1888), págs. 561-564. [4] A. Einstein. Investigations on the Theory of the Brownian Movement. Dover Corporation, 1956. [5] H. B. Callen y T. A. Welton. «Irreversibility and generalized noise». En: Physical Review 83 (1951), pág. 34. [6] J. Perrin. Brownian movement and molecular reality. Dover Corporation, 2013. [7] J. Perrin (17 de abril de 2013). Mouvement brownien et molécules. Laboratoire de chimie physique de la Sorbonne. url:https://images.cnrs.fr/video/2080. [8] D. S. Lemons y A. Gythiel. «Paul Langevin’s 1908 paper “On the theory of Brownian motion”[“Sur la théorie du mouvement brownien,” cr Acad. Sci.(Paris) 146, 530–533 (1908)]». En: American Journal of Physics 65 (1997), págs. 1079-1081. [9] C. W. Gardiner. Stochastic Methods: A Handbook for the Natural and Social Sciences. Springer, 2009. [10] M. D. Haw. «Colloidal suspensions, Brownian motion, molecular reality: a short history». En: Journal of Physics: Condensed Matter 14 (2002), pág. 7769. [11] C. A. Plata, D. Guéry-Odelin, E. Trizac y A. Prados. «Building an irreversible Carnot-like heat engine with an overdamped harmonic oscillator». En: Journal of Statistical Mechanics: Theory and Experiment 2020 (2020), pág. 093207. [12] N. G. Van Kampen. Stochastic Processes in Physics and Chemistry. Elsevier, 1992. [13] C. M. Bender y S. A. Orszag. Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory. Springer Science & Business Media, 1999. [14] S. Ciliberto. «Experiments in stochastic thermodynamics: short history and perspectives». En: Phys. Rev. X 7(2017), pág. 021051. [15] L. D. Landau y E. M. Lifshitz. Mechanics: Volume 1. Elsevier Science, 1982. [16] A. G. Ramm. «A simple proof of the Fredholm alternative and a characterization of the Fredholm operators». En: The American Mathematical Monthly 108 (2001), págs. 855-860. 55 BIBLIOGRAFÍA 56 [17] A. V. Skorokhod. «Stochastic equations for diffusion processes in a bounded region». En: Theory of Probability & Its Applications 6(1961), págs. 264-274. [18] A. V. Skorokhod. «Stochastic equations for diffusion processes in a bounded region. II». En: Theory of Probability & Its Applications 7(1962), págs. 3-23. [19] F. Le Vot, S. B. Yuste y E. Abad. «Standard and fractional Ornstein-Uhlenbeck process on a growing domain». En: Physical Review E 100 (2019), pág. 012142.