Introducción al control de sistemas estocásticos

INTRODUCCIÓN AL CONTROL DE SISTEMAS ESTOCÁSTICOS FACULTAD DE FÍSICA FÍSICA ATÓMICA, MOLECULAR Y NUCLEAR David Molina Barrera Curso 2023-2024 Tutor: Antonio Prados Montaño Índice general 1. Introducción 2 1.1. Motivaciones .......................................... 2 1.2. Objetivos ............................................ 3 1.3. Movimiento Browniano: Ecuaciones de Langevin y de Fokker-Planck . . . . . . . . . . . . 3 1.3.1. EcuacióndeLangevin ................................. 4 1.3.2. Ecuación de Fokker-Planck (FP) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4. Termodinámicaestocástica................................... 8 1.5. Transformaciones estado-estado rápidas . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.5.1. Ingenieríainversa.................................... 12 1.5.2. Controlóptimo..................................... 12 2. Movimiento sobreamortiguado de una partícula browniana en una trampa armónica 14 2.1. Solución de la ecuación de Fokker-Planck para una trampa armónica . . . . . . . . . . . . . 15 3. Ingeniería inversa para una partícula browniana en una trampa armónica 17 3.1. Ingenieríainversa........................................ 17 3.2. Conexión ESE polinómica de tercer grado . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.2.1. Simulaciones numéricas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.3. Conexión ESE polinómica de cuarto grado . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.3.1. Simulaciones numéricas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 4. Conexiones óptimas sin restricciones 28 4.1. Caso general: Transporte óptimo de masa . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.2. Casoarmónico ......................................... 34 4.2.1. Conexión óptima sin restricciones . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.2.2. Simulaciones numéricas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.2.2.1. Protocolo óptimo en expansión . . . . . . . . . . . . . . . . . . . . . . . 39 5. Conclusiones 43 A. Simulaciones numéricas con Langevin y determinación de la función de distribución 45 Nomenclatura ESE Engineered Swift Equilibration FB Fast-Backward FDP Función de Densidad de Probabilidad FF Fast-Forward FP Fokker-Planck SST Swift State-to-state Transformation STEP Protocolo de salto abrupto Resumen Controlar los procesos termodinámicos que conectan estados en sistemas mesoscópicos es posible, dentro del marco teórico de la Termodinámica Estocástica. Estos sistemas son aquellos en los que el contacto con un baño térmico hace aparecer un ruido o fuerza aleatoria en su dinámica; un ejemplo paradigmático es una partícula browniana sumergida en un fluido. En este trabajo, tras una introducción teórica sobre el movimiento browniano y los conceptos termodinámicos en la mesoescala, controlamos la evolución temporal de la constante elástica de un potencial armónico para conectar dos estados de equilibrio, correspondientes a valores distintos de la constante elástica, de una partícula browniana confinada en este potencial. En primer lugar, haciendo uso de la técnica de ingeniería inversa, mostramos que la conexión acelerada entre estados es factible, alcanzándose un estado de equilibrio con una constante elástica que cuadriplica la inicial en un tiempo 50 veces inferior al de relajación natural. Esto se consigue a cambio de efectuar un trabajo medio sobre el sistema unas 30 veces superior a la diferencia de energía libre entre los estados. También tratamos el problema de minimización del trabajo efectuado en una conexión isoterma, siguiendo para ello dos enfoques bien diferenciados: uno numérico, que no lleva a ninguna mejora significativa; y uno formal, por medio de un método de mapeado a problema de transporte óptimo de masa—que describimos para un potencial arbitrario. El protocolo que minimiza el trabajo isotermo también puede obtenerse, para el potencial armónico, usando directamente un principio variacional; recurrimos aquí al método de mapeado simplemente para ilustrar una de las herramientas de optimización disponibles dentro del campo del control. El protocolo óptimo de variación de la rigidez conlleva una reducción de, aproximadamente, un 20% en el trabajo irreversible—entendiéndolo como la diferencia entre trabajo medio y variación de energía libre entre los estados—que los protocolos subóptimos anteriormente obtenidos. Todos los protocolos anteriores se implementan en simulaciones numéricas para comprobar que se han determinado correctamente, consiguiendo un buen acuerdo con las predicciones teóricas. 1 Capítulo 1 Introducción 1.1. Motivaciones Los avances a lo largo de los últimos 20 años en el campo de la física estadística de no equilibrio permiten el control de los procesos termodinámicos que tienen lugar en los sistemas mesoscópicos. El calificativo mesoscópico se refiere a sistemas con pocos grados de libertad1en contacto con un baño térmico. La interacción con el baño introduce fluctuaciones en las variables dinámicas del sistema [1]. Esta interacción puede modelarse como una fuerza aleatoria o estocástica, que se añade a las ecuaciones de movimiento (imagen de Langevin). Cada trayectoria de un mismo experimento es diferente de las demás, dependiendo de cada realización de la fuerza aleatoria. Una descripción equivalente propone tratar al sistema como un colectivo de configuraciones, con una función de densidad de probabilidad que obedece una cierta ecuación de evolución (ecuación de Fokker-Planck) [2]. Como ejemplo paradigmático de sistema mesoscópico, consideremos una partícula browniana inmersa en un fluido en equilibrio a la temperatura T. La partícula es pesada comparada con las demás partículas del fluido, y sigue un movimiento aleatorio al estar sujeta a las colisiones con estas mismas partículas [3]. Este llamado movimiento browniano juega un papel destacado en pequeños sistemas biológicos, y en micro y nanodispositivos. Esto es así porque—en estos sistemas—las energías involucradas en sus procesos termodinámicos son del orden de las fluctuaciones energéticas2[4]. Una partícula browniana puede ser confinada por "pinzas ópticas", que usan la presión de radiación y las fuerzas de gradiente3de un haz de luz fuertemente focalizado para crear un poten- 1Es significativo únicamente el estudio de sistemas con pocos grados de libertad, ya que el ruido térmico que introduce el baño depende en un factor 1/√Ndel número de grados de libertad, N.En el límite macroscópico N→∞ las fluctuaciones desaparecen. 2Típicamente estas fluctuaciones son del orden de kBT,con kBla constante de Boltzmann. 3Estas fuerzas se dan lugar en regiones en las que hay una fuerte variación del campo eléctrico local—como ocurre alrededor del punto en el que se concentra un láser—, afectando a cualquier dipolo o partícula dieléctrica. 2 cial prácticamente armónico cerca del punto de focalización. Tanto la posición del mínimo como la constante elástica4de este potencial se pueden controlar con gran precisión [5]. La Termodinámica Estocástica ha permitido extender los conceptos termodinámicas como el calor, el trabajo, o la producción de entropía a la escala mesoscópica, haciendo uso de las herramientas matemáticas de los procesos estocásticos. Al incorporarse el tiempo en las ecuaciones, es posible el control fino de los procesos termodinámicos a través de la variación de ciertos parámetros de control, para sistemas cuya dinámica sea estocástica [1]. Una de las aplicaciones más prometedoras de esta rama de la física es el desarrollo de motores térmicos estocásticos, consistentes en una partícula browniana confinada capaz de realizar un ciclo cerrado, inspirado en alguno de los ciclos clásicos, como el ciclo de Carnot—con dos procesos isotermos y dos adiabáticos.5Estos motores tienen interesantes propiedades en términos energéticos y de eficiencia, que en este caso es una magnitud fluctuante. Además, su construcción podría hacer posible realizar tareas de transporte de partículas o de energía a pequeña escala, algo que supondría notables avances en medicina, e.g. la destrucción selectiva de células cancerígenas [6]. 1.2. Objetivos Este trabajo se ha concebido como una iniciación al campo del control de sistemas estocásticos, que se ha desarrollado en la última década. Nuestro objetivo esencial es controlar la evolución temporal de sistemas con dinámica estocástica. Por simplicidad, consideraremos una partícula browniana confinada por una trampa armónica e inmersa en un fluido en equilibrio. Queremos diseñar protocolos que conecten dos estados de equilibrio de esta partícula, correspondientes a valores distintos de la rigidez de la trampa, en un tiempo finito. Esto es, queremos diseñar lo que se conoce como transformaciones rápidas estado-estado (SST, del inglés "swift state-to-state transfomation"). Estamos interesados en el diseño de estas SST tanto por métodos de ingeniería inversa como por métodos de control óptimo. En particular, dentro del apartado de control óptimo, buscaremos el protocolo SST que minimiza la contribución irreversible al trabajo promedio. 1.3. Movimiento Browniano: Ecuaciones de Langevin y de Fokker-Planck El tratamiento matemático del movimiento browniano se sustenta en la teoría de los procesos estocásticos, en concreto, se modela como un proceso de Markov: el estado del sistema en un ins- 4Depende linealmente de la potencia del láser. 5Adiabáticos: sin intercambios de calor. No en el sentido que se indica más adelante. 3 tante dado sólo depende del estado en un instante inmediatamente anterior, para una cierta escala de tiempo bien elegida. En esta escala de tiempo ocurren muchas colisiones instantáneas e independientes entre sí (y el proceso markoviano es la velocidad), o bien se dan muchas variaciones rápidas de la velocidad en cada instante (y entonces el proceso markoviano sería la posición). Cuando hay una fuerza externa aplicada al sistema que depende de la posición, su evolución está gobernada por la llamada ecuación de Kramers, que describe el proceso markoviano compuesto (posición y velocidad). En el límite sobreamortiguado, el proceso markoviano es la posición de la partícula [2]. Consideremos una partícula «pesada» de masa Minmersa en un fluido en equilibrio a la temperatura T. Por simplicidad, se considera el movimiento unidimensional. La partícula está sometida a un potencial U(x,t)que la confina, que puede depender también del tiempo, como sucede en el caso de control con el que tratamos en este trabajo. La evolución dinámica de la posición es la del movimiento browniano en el régimen sobreamortiguado. Aquí se muestra una forma sencilla de obtener la ecuación que gobierna la evolución de este sistema prototípico fundamental, haciendo una serie de consideraciones plausibles. Con ello se pretende justificar el uso de dicha ecuación, no obstante, no damos la deducción completa pues esto excede los propósitos de este trabajo. Se puede encontrar más información en las referencias [2, 7]. 1.3.1. Ecuación de Langevin En la descripción macroscópica, la ecuación de evolución de la partícula es Mdv dt =−γv−∂U(x,t) ∂x,(1.1) donde ves la velocidad de la partícula y γel coeficiente de fricción del fluido. Cuando el término a la izquierda es despreciable, la ecuación (1.1) se aproxima por6 γ˙x=−∂U(x,t) ∂x,(1.2) que se conoce como el límite sobreamortiguado. En este límite, la velocidad de la partícula sigue instantáneamente a la fuerza impuesta por el potencial U(x,t).La ecuación (1.2) extiende la solución estacionaria, para una fuerza constante, de (1.1).7Por ello, esta solución será cierta siempre que la velocidad inicial—que no tiene por qué verificar (1.2)—ya se haya amortiguado; y la variación relativa de la fuerza −∂xU(x,t),con respecto a la posición o con respecto al tiempo, sea 6Note que se hace un cambio de notación de va ˙x. 7La solución de (1.1) cuando −∂xU(x,t) = F,con Fconstante, es v(t) = v0e−γ Mt+γ−1F1−e−γ Mt, siendo v0=v(t=0),tal que v(t→∞) = γ−1F. 4 insignificante en la escala de tiempo en la cual evoluciona la velocidad.8Esta escala de tiempo es del orden de M/γ. Notemos que en la ecuación (1.2) no aparece ningún término de fluctuaciones. Para incorporarlas, tenemos en cuenta que el fluido está en equilibrio a la temperatura Ty, por tanto, se cumple la relación de fluctuación-disipación. De esta manera, la ecuación con fluctuaciones es la denominada ecuación de Langevin: γ˙x=−∂xU(x,t)+p2γkBTξ(t),(1.3) donde kBes la constante de Boltzmann y ξ(t)es un ruido blanco gaussiano de varianza unidad: ⟨ξ(t)⟩=0,ξ(t)ξ(t′)=δ(t−t′).(1.4) En la imagen de Langevin, ⟨...⟩equivale a promediar sobre una gran multitud de trayectorias. 1.3.2. Ecuación de Fokker-Planck (FP) Una vez obtenida la ecuación de Langevin (1.3), es directo pasar a la descripción de Fokker- Planck mediante el cálculo de los momentos de salto [2], que están definidos como al(x,t) = l´ ım ∆t→0+D(∆x)lEx ∆t,D(∆x)lEx=ˆ+∞ −∞ dx′x′−xlP(x′,t+∆t|x,t),l≥1.(1.5) Aquí, P(x′,t+∆t|x,t)denota la probabilidad de transición entre los estados (x,t)y(x′,t+∆t), en un intervalo de tiempo ∆t.Los momentos de salto son entonces el valor medio de la potencia l-ésima del salto x′−xpor unidad de tiempo. Estos permiten escribir la ecuación de FP ∂tρ(x,t) = ∂x−a1(x,t)ρ(x,t)+ 1 2∂x[a2(x,t)ρ(x,t)],(1.6) donde ρ(x,t)es la función de densidad de probabilidad (FDP). Integrando (1.3), se obtiene que ∆x=−1 γˆt+∆t t dt′∂xU(x,t′)+s2kBT γˆt+∆t t dt′ξ(t′).(1.7) Ahora promediando con x(t)fijo y tomando ∆ttal que ∂xU(x,t)apenas varíe, se deduce que ⟨∆x⟩x(t)=x=−1 γ∂xU(x,t)∆t+s2kBT γˆt+∆t t dt′ξ(t′),(1.8) 8El lado izquierdo de (1.1) debe compensar instantáneamente sobre la escala de tiempo las variaciones de la fuerza. 5 donde por (1.5) y usando (1.4), queda que el primer momento de salto es a1(x,t) = l´ ım ∆t→0+⟨∆x⟩x ∆t=−1 γ∂xU(x,t).(1.9) Utilizando de nuevo (1.7), obtenemos (∆x)2=1 γ2*ˆt+∆t t dt′∂xU(x,t′)ˆt+∆t t dt′′∂xU(x,t′′)+ −2 γs2kBT γ∂xU(x,t)∆tˆt+∆t t dt′ξ(t′) −2 γs2kBT γ∂2 xU(x,t)ˆt+∆t t dt′ˆt+∆t t dt′′x(t′)−xξ(t′′) +2kBT γˆt+∆t t dt′ˆt+∆t t dt′′ξ(t′)ξ(t′′)+o(∆t),(1.10) donde hay términos de orden superior a ∆tcomo se indica. El primer término en (1.10) de (∆x)2 es del orden (∆t)2,por lo que se anula al tomar el límite ∆t→0 para calcular el momento de salto a2.El término cruzado que proviene de (1.7) se desarrolla en serie, esto da lugar a los términos segundo y tercero. Hemos tenido en cuenta que, tanto el término del desarrollo que incluye a∂t∂xU(x,t),como términos superiores del desarrollo, son o(∆t).El segundo término se anula, al igual que en (1.8). Debido a que tanto ξ(t′′)como x(t′)están determinadas para cada realización, sin haber ninguna delta de Dirac en ellas, el tercer término es o(∆t)—es por este argumento que los términos superiores del desarrollo también son o(∆t).Finalmente, sólo queda el cuarto término, acompañado de términos o(∆t).Introduciendo (1.4) en (1.10), el segundo momento de salto—según la definición (1.5)—es a2(x,t) = l´ ım ∆t→0+D(∆x)2Ex ∆t=2kBT γ.(1.11) Siguiendo la misma línea de razonamiento que para obtener el segundo momento de salto, se puede comprobar que D(∆x)lEx=o(∆t)y entonces al(x) = 0,para todo l>2.Esto justifica escribir, usando los momentos de salto (1.9) y (1.11), la ecuación de FP γ∂tρ(x,t) = ∂x{[∂xU(x,t)]ρ(x,t)+kBT∂xρ(x,t)}.(1.12) La ecuación (1.12) se denomina ecuación de Smoluchowski, que es la ecuación de FP en el límite sobreamortiguado. De modo análogo a la ecuación de Langevin (1.3), esta ecuación es vá- 6 de energía libre entre los estados, el tiempo de conexión, etc. [1]. Es especialmente relevante la minimización del trabajo irreversible en un proceso isotermo, un problema que se ha resuelto con bastante generalidad en el caso sobreamortiguado (abordamos formalmente el problema en el Capítulo 4), tanto partiendo del enfoque de Langevin [13] como del enfoque de FP [14]. Aun así, con los métodos actuales, sólo para ciertos potenciales aplicados se pueden obtener protocolos en forma cerrada (técnicas alternativas en forma de protocolos subóptimos, con costes cercanos a los óptimos, también han sido propuestas y comprobadas [15]); en ellos, aparecen discontinuidades en los instantes inicial y final [1]. El caso más simple es el del potencial armónico, para el cual el control óptimo sin ligaduras se obtiene directamente de un principio variacional [16]—esto proviene de que la ecuación de FP, cuando el potencial aplicado es armónico, se puede reducir a una ecuación diferencial ordinaria para la varianza de la FDP (lo demostramos en la sección 2.1). Debido a las restricciones presentes en la aplicación experimental de los protocolos, pueden aparecer ligaduras no honónomas—con forma de desigualdad—en las funciones de control;15 la teoría matemática del control óptimo se encarga de dar respuesta a este problema a través del principio de máximo de Pontryagin, permitiendo encontrar la solución que minimiza el trabajo irreversible, cuando la constante elástica k(t)está restringida entre 0 y un valor km´ax para todo tiempo t. En este caso, aparece un tiempo mínimo de conexión entre estados, que se corresponde con la duración de la evolución para el cambio brusco de k,a uno de los valores extremos. Cuando el protocolo k(t)sin restricciones—que realiza la conexión entre estados—no es posible porque cruza los valores frontera, la solución que proporciona el principio de Pontryagin requiere efectuar un trabajo mayor [11]. Por otro lado, en el problema de minimización del trabajo en el régimen subamortiguado, las principales limitaciones surgen de que el potencial óptimo de conducción depende de la velocidad, cuando lo deseable para su implementación experimental es que se trate de un potencial conservativo de la forma U(x,t)[1]. 15El obstáculo físico más problemático, se relaciona con la incapacidad de crear potenciales que se vuelvan repulsivos en una cierta ventana de tiempo, para procesos rápidos en expansión (aumento del espacio de fase, es decir cuando la constante elástica disminuye). 13 Capítulo 2 Movimiento sobreamortiguado de una partícula browniana en una trampa armónica El sistema prototípico es un coloide inmerso en un fluido en equilibrio a la temperatura T. Recordamos que este coloide se trata de una partícula «pesada» de masa M,sometida a un potencial U(x,t)que la confina adecuadamente. La evolución dinámica de la posición xde la partícula es la del movimiento browniano, en el régimen sobreamortiguado. Así, esta partícula browniana se puede describir o bien por una ecuación de Langevin, o bien por una ecuación de Fokker-Planck. En la introducción explicábamos cómo obtener la ecuación de FP de este sistema para un potencial arbitrario: partíamos primero de la formulación de Langevin y, una vez obtenida la ecuación con fluctuaciones del proceso, pasábamos a la formulación de FP. Esta última ecuación se escribía directamente, haciendo uso de la equivalencia entre ambos enfoques [2]. Aquí mostramos como esta ecuación se simplifica enormemente a una ecuación diferencial ordinaria para la única variable dinámica relevante—la varianza—cuando el potencial aplicado es armónico. En los capítulos siguientes, implementamos en simulaciones numéricas diversos protocolos, deducidos haciendo uso de la técnica de ingeniería inversa, que permiten el control de este sistema. Además del potencial, se podría considerar que la temperatura del baño térmico dependiera del tiempo. Por ejemplo, esto se puede conseguir de modo efectivo mediante un forzado aleatorio del centro de la trampa de confinamiento, lo que permite emular un ruido adicional para la partícula, sin alterar la temperatura del fluido. La fuerza estocástica aplicada debe ser tal que su tiempo de correlación sea mucho menor que el tiempo de relajación intrínseco del sistema, para que pueda ser tratada como ruido blanco. Este procedimiento sólo permite calentar, ya que las varianzas de los dos ruidos se suman. Más sobre el control de la temperatura puede hallarse en las ref. [1, 17]. 14 2.1. Solución de la ecuación de Fokker-Planck para una trampa armónica La ecuación de FP tiene una propiedad muy útil cuando el potencial aplicado es armónico, esto es U(x,t) = 1 2k(t)x2. Consiste en que, cuando la distribución de la posición de la partícula es inicialmente gaussiana, esta distribución se mantiene gaussiana para todo instante. Esto es así independientemente de que el sistema esté «forzado», en el sentido de que la constante elástica k y/o la temperatura Tsean funciones del tiempo. En las siguientes líneas se prueba esta propiedad, que permite derivar una ecuación diferencial ordinaria para la varianza, completamente equivalente a FP para la conexión de estados gaussianos. Partimos de la ecuación (1.12), particularizada para U(x,t) = 1 2k(t)x2: γ∂tρ(x,t) = ∂x{k(t)xρ(x,t)+kBT∂xρ(x,t)}.(2.1) Lo primero es pasar al espacio de momentos mediante la transformación G(q,t) = ˆ+∞ −∞ dxeiqxρ(x,t),(2.2) donde G(q,t)es la función característica. Procediendo de esta manera, (2.1) se convierte en γ∂tG(q,t) = −qk(t)∂qG(q,t)−kBTq2G(q,t),(2.3) de la cual se obtiene γ∂tlnG(q,t) = −qk(t)∂qlnG(q,t)−kBTq2.(2.4) Aquí, lnG(q,t)es la función generatriz de cumulantes, esto es lnG(q,t) = ∞ ∑ n=1 (iq)n n!χn(t),(2.5) siendo χn(t)los cumulantes. Introduciendo este desarrollo en (2.4), e igualando para cada potencia n-ésima de q,se llega a un sistema de infinitas ecuaciones diferenciales ordinarias de primer orden, cuya solución determina la evolución de los cumulantes. Para n≥1,n=2,se tiene una ecuación del tipo γ˙ χn(t) = −nk(t)χn(t),n≥1,n=2,(2.6) mientras que para n=2,se tiene γ˙ χ2(t) = −2k(t)χ2(t)+ 2kBT.(2.7) 15 La solución del subsistema (2.6) es simplemente χn(t) = χn(0)exp−nˆt 0 dt′k(t′) γ,n≥1,n=2.(2.8) Por el momento, todo lo hecho es general para un potencial armónico. Lo siguiente es imponer condiciones iniciales gaussianas, de manera que por (2.2): ρ(x,0) = 1 q2π σ2 0 exp"−(x−µ0)2 2σ2 0#⇐⇒ lnG(q,0) = iqµ0−q2σ2 0 2,(2.9) de donde (al comparar con (2.5)) identificamos χ1(0) = µ0yχ2(0) = σ2 0,siendo nulos el resto de cumulantes. Vemos que bajo esta circunstancia, por (2.8), lnG(q,t)mantiene su forma funcional inicial—χn(t) = 0 para n>2—y, por ende, ρ(x,t)permanece gaussiana ∀t. De hecho, cabe recalcar que cuando se quiere construir un protocolo ESE para una partícula browniana en una trampa armónica, se parte de un estado con µ0=0,pues el estado de equilibrio asociado al potencial U(x) = 1 2kx2cuando hay amortiguamiento tiene valor medio nulo. Esto último se puede ver de la evolución de los cumulantes (2.8) y (2.7), consecuentemente siendo χ2(t)el único que sobrevive para t→∞.Notemos también que en régimen sobreamortiguado, cuanto mayor sea el orden del cumulante, más rápida es la escala de tiempo con la que cae a cero—en particular, el primer cumulante (valor medio) obedece la ecuación sobreamortiguada (1.2). Efectuando un cambio de notación en (2.7), obtenemos la «ecuación de la varianza»,1que contiene toda la información dinámica del problema: d dt x2(t) = −2k(t) γx2(t)+2kBT γ.(2.10) Su solución para un proceso isotermo es x2(t) = "x2(0)+2kBT γˆt 0 exp 2ˆt′ 0 k(t′′) γdt′′!dt′#exp−2ˆt 0 k(t′) γdt′.(2.11) Se sigue de (2.11) que cuando k(t) = kconstante, se tiene x2(t) = kBT k+x2(0)−kBT ke−2k γt,(2.12) así la partícula prácticamente relaja al equilibrio (x2=kBT/k)en un tiempo teq ≃3γ/k. 1La notación presupone que se busca un protocolo SST entre estados de equilibrio, en consecuencia ⟨x⟩(t) = 0. 16 Capítulo 3 Ingeniería inversa para una partícula browniana en una trampa armónica En este capítulo, empleamos la técnica de ingeniería inversa para deducir funciones de control que permiten realizar la conexión, en un tiempo finito arbitrariamente corto, entre dos estados de equilibrio de una partícula browniana en un potencial armónico, en el régimen sobreamortiguado. Primero, se describe el método en esta situación particular, trabajando sobre la ecuación (2.10). En segundo lugar, construimos un protocolo que asegura la conexión, para una evolución polinómica de la varianza. Finalmente, se plantea una forma polinómica más compleja, con el objetivo de encontrar numéricamente un control optimizado—en el sentido de minimización del trabajo— respecto al anterior. Complementamos los desarrollos analíticos con simulaciones numéricas. 3.1. Ingeniería inversa La forma de proceder se basa en imponer un ansatz a la varianza x2(t)y despejar el protocolo de control k(t)de (2.10). De este modo, la conducción viene dada por k(t) = kBT ⟨x2⟩(t)−γ 2 d dt lnx2(t).(3.1) Nuevamente, recalcamos que se hace un control en k(sin embargo, también sería posible en T). Además, el ansatz deberá ser tal que se cumplan las condiciones de contorno: x2(0) = kBT ki ,x2(tf) = kBT kf ,(3.2) d dt x2(t)t=0=0,d dt x2(t)t=tf =0,(3.3) 17 las cuales aseguran que se conectan los estados de equilibrio correspondientes a ki=k(0)y kf=k(tf),1siendo tfla duración del protocolo. Pese a que nos centramos en el movimiento sobreamortiguado en una trampa armónica, por ser el caso más elemental y a su vez el más sencillo de implementar en el desarrollo de aplicaciones, también han sido estudiadas otras situaciones físicas. Por ejemplo, en el régimen subamortiguado el potencial armónico sigue manteniendo el carácter gaussiano de la solución si inicialmente la distribución era gaussiana, salvo que ahora despejar el control k(t)se vuelve una tarea más complicada pues se tiene un sistema de ecuaciones acopladas en el que, aparte de x2, también aparecen v2 y⟨xv⟩como momentos relevantes [1, 2]. Adicionalmente, cabría la posibilidad de conectar estados de equilibrio arbitrarios, no gaussianos.2En el régimen sobreamortiguado, las condiciones que se suelen exigir para obtener solución analítica, hacen que los potenciales de control U(x,t)estén restringidos a aquellos que preservan la forma de la distribución ρ(x,t).3Esto permite extender el método de atajo a potenciales más generales, pero entre estados de la misma familia de potenciales4[18]. No obstante, siempre se puede aplicar un método alternativo que involucra al de ingeniería inversa, en concreto una combinación de protocolos Fast-Forward (FF)5yFast-Backward (FB): tomando como referencia una dinámica conocida,6se puede establecer una «soldadura» de dos protocolos en un estado intermedio ρint(x,t). Esta soldadura conecta efectivamente ambos estados arbitrarios, en un tiempo arbitrariamente corto, bajo un cierto conjunto de potenciales FF [15]. Aun más limitaciones emergen cuando se pretenden conectar estados arbitrarios en régimen subamortiguado, típicamente debido a las suposiciones para hallar soluciones explícitas [1]. 3.2. Conexión ESE polinómica de tercer grado Para conducir al sistema rápidamente entre un estado de equilibrio a temperatura T, con una rigidez kide la trampa, y otro estado de equilibrio, con una rigidez kf, es posible utilizar un polinomio como función de interpolación para la varianza x2(t)y despejar k(t)de (3.1). Con este fin, el polinomio ha de ser tal que cumpla las condiciones de contorno (3.2) y (3.3). Así, como son necesarios cuatro parámetros libres, el menor grado requerido para que el polinomio pueda cum- 1Las condiciones ki=k(0)ykf=k(tf), son equivalentes a (3.3). 2En otras palabras, lnρeq(x)∝−βU(x),con U(x)no armónico. 3Se tiene que ∀t∈(0,tf), lnρ(x,t)∝−βUi(x−µ(t) σ(t)), con µ(t)yσ(t)continuas, tales que lnρ(x,0)∝−βUi(x). 4El mantenimiento de la forma se puede ver como una propiedad similar a la de permanencia gaussiana. 5Consiste en conducir al sistema desde un estado inicial hasta un estado final al que relajaría lentamente, pero sobre una escala de tiempo distorsionada acelerada (no necesariamente en un tiempo finito), a través de la misma sucesión de estados. Cuando lo que se trata es de retornar al sistema al estado de partida, siguiendo la misma sucesión de estados pero invertida, se conoce con el nombre de Fast-Backward. 6Se puede resolver analíticamente FP, dado un potencial U(x,t). 18 plir su propósito es tres, siempre que se busque una función de control k(t)continua.7Podría ser posible una elección de grado inferior si se permitiera una discontinuidad de salto finito para ken los extremos del intervalo [0,tf],esto es interesante y se retoma en el Capítulo 4, en el contexto de protocolos óptimos. La metodología antes descrita es precisamente la seguida en la ref. [12], aunque allí se trabajaba con una variable proporcional a la inversa de la varianza para establecer la conexión. De esta manera, consiguieron ser los primeros en implementar experimentalmente un protocolo ESE, sobre un sistema clásico en contacto con un baño térmico. El sistema se trataba de una esfera micrométrica inmersa en agua, atrapada en un potencial armónico producido por un láser—en el régimen sobreamortiguado—cuya potencia regulaba el valor de la rigidez. Compararon el protocolo ESE con la evolución obtenida al cambiar instantáneamente la rigidez ka su valor final (llamando a este salto protocolo STEP), logrando una relajación por ESE en un tiempo tf100 veces más pequeño que el tiempo de relajación tras el cambio abrupto, i.e. teq ≃3γ/kf.Comprobaron también que la propiedad gaussiana se mantenía cuando el potencial era armónico. Antes de escribir explícitamente el protocolo, conviene introducir la siguiente adimensionalización de las variables: x∗=x p⟨x2⟩i ,x2∗=x2 ⟨x2⟩i ,x2i=kBT ki ;k∗=k ki ;trx =γ kf ,t∗=t trx ;τ=t∗ t∗ f .(3.4) Insertando (3.4) en (3.1), obtenemos k∗(τ) = hx2∗(τ)i−1−k∗ f 2t∗ f d dτlnhx2∗(τ)i.(3.5) La prescripción para la varianza, cumpliendo las condiciones de contorno, es x2∗(τ) = 1+3Λτ2−2Λτ3,(3.6) de modo que x2∗(τ=1) = 1+Λ. Por tanto, el protocolo de control queda k∗(τ) = 1 1+3Λτ2−2Λτ3−3k∗ fΛ(τ−τ2) t∗ f(1+3Λτ2−2Λτ3),(3.7) donde además, se infiere la relación k∗ f=k∗(τ=1) = 1 1+Λ.(3.8) 7Se sigue de (3.1). 19 Previamente a comprobar la fiabilidad de los protocolos en simulaciones numéricas, abordamos los aspectos energéticos de la evolución del sistema y, en particular, el coste que conlleva el batir la escala natural de relajación. Estos aspectos, fundamentados en el marco teórico de las propiedades energéticas estocásticas, han sido presentados brevemente en la sección 1.4 de la introducción. El trabajo medio acumulado ejercido sobre la partícula en la duración del protocolo es ⟨W⟩(t) = ˆt t0 dt′1 2 d dt′k(t′)x2(t′) = 1 2k(t′)x2(t′)t t0−ˆt t0 dt′1 2k(t′)d dt′x2.(3.9) Sabiendo que k(t′)viene dado por (3.1), la ecuación (3.9) se convierte en ⟨W⟩(t) = 1 2k(t′)x2(t′)t t0 +1 2kBTln"x2(t0) ⟨x2⟩(t)#+γ 4ˆt t0 dt′1 ⟨x2⟩d dt′x22 .(3.10) En esta última fórmula se puede identificar ∆F(t)≡⟨∆U⟩(t)−T∆S(t) = 1 2k(t′)x2(t′)t t0 +1 2kBTln x2(t0) ⟨x2⟩(t)!,(3.11) como la variación de energía libre, siendo ⟨∆U⟩(t) = 1 2k(t′)x2(t′)t t0 ,(3.12) ∆S(t) = 1 2kBln"x2(t) ⟨x2⟩(t0)#.(3.13) Además, como el sistema está en equilibrio a la temperatura Ten los instantes inicial y final, entonces k(0)x2(0) = k(tf)x2(tf) = kBT, y el término ⟨∆U⟩(tf) = 1 2k(t)x2(t)tf 0se anula. El trabajo medio acumulado hasta t=tf,para un potencial armónico, queda ⟨W⟩(tf) = 1 2kBTlnkf ki+γ 4ˆtf 0 dt 1 ⟨x2⟩d dt x22 .(3.14) El trabajo medio tiene dos contribuciones: la primera se corresponde con la diferencia de energía libre, y la segunda es estrictamente no negativa y sólo se anula bajo una evolución cuasiestática.8 Por este motivo, se define el trabajo irreversible Wirr ≡⟨W⟩(tf)−∆F(tf) = γ 4ˆtf 0 dt 1 ⟨x2⟩d dt x22 ,(3.15) 8En dicha evolución, tf≫trx ,d dt x2=0. 20 que por la adimensionalización (3.4), se puede escribir como W∗ irr =Wirr kBT=k∗ f t∗ f η,(3.16) con el parámetro ηdeterminado por la fórmula ηhx2∗i=ˆ1 0 dτ1 ⟨x2⟩∗d dτx2∗2 ,(3.17) dependiente del ansatz. El valor del calor medio acumulado se obtiene a partir del primer principio (1.21), por el cual se tiene ⟨Q⟩(t) = ⟨∆U⟩(t)−⟨W⟩(t).(3.18) para un proceso isotermo, con ⟨∆U⟩(t)dada por (3.12). 3.2.1. Simulaciones numéricas En el Apéndice A se encuentra una descripción del método empleado para llevar a cabo las simulaciones. Este se basa en la integración numérica de la ecuación de Langevin (1.3), que para un potencial armónico toma la forma dx dt =−k(t) γx(t)+√2Dξ(t),(3.19) donde D=kBT/γ.La ecuación integrada, en las variables adimensionales introducidas en (3.4), es x∗(τ+∆τ) = x∗(τ)−k∗(τ)x∗(τ)t∗ f k∗ f ∆τ+s2t∗ f k∗ f ∆τN(τ),(3.20) siendo N(τ)una variable aleatoria que sigue una distribución normal, de media nula y varianza unidad,9evaluada en cada instante de tiempo τ. También es posible conocer la FDP de la posición, mediante histogramas construidos en instantes equiespaciados. Una vez se tiene el valor de la FDP en función de la posición para suficientes instantes de tiempo, se puede estimar su valor en tiempos intermedios, por interpolación directa entre las FDP conocidas. Para efectuar las simulaciones, es necesario proporcionar una relación de compresión k∗ f—que da cuenta de la separación entre los estados inicial y final—y un tiempo de duración del protocolo t∗ f. La evolución temporal del parámetro de control k∗(τ)es tal que, partiendo del estado inicial de equilibrio k∗(0) = 1, se alcanza el estado final k∗(τ=1) = k∗ f, en un tiempo t∗=t∗ f. Cuando 9Véase el Apéndice A para entender el origen de esta variable. 21 k∗ fes mayor (menor) que la unidad, estamos comprimiendo (expandiendo) el sistema. A modo de ejemplo ilustrativo, escogemos k∗ f=4 y t∗ f=0,06. Los resultados se muestran en la Figura 3.1. La rigidez del potencial k∗(τ)se modifica conforme a (3.7), hasta un tiempo t∗ f(τ=1)a partir del cual permanece en su valor final k∗ f=4.La simulación prosigue hasta un tiempo t∗=1,5t∗ f para verificar que, efectivamente, la partícula se mantiene en equilibrio para τ≥1. Promediando en cada instante sobre las Ntray =105trayectorias, obtenemos la varianza como función del tiempo. Esta última evoluciona desde un valor unidad hasta el valor σ∗2 f=1/k∗ f=0,25.Hay un acuerdo casi perfecto entre la varianza calculada de este modo y la prescripción teórica (3.6). Para determinar la FDP, se registran valores de x∗en el intervalo x∗∈[−3,3]: es una gaussiana centrada (⟨x⟩=0)en todo instante de tiempo, como corresponde por serlo inicialmente.10 Se observa que la desviación estándar de la FDP, se reduce a la mitad en el equilibrio.Se consigue el equilibrado en un tiempo t∗ f=0,06, que es 50 veces inferior al tiempo de relajación tras el cambio abrupto, de k∗a su valor final k∗ f(STEP).11 Para producir tal reducción de la escala de tiempo, es necesario que la rigidez aumente considerablemente su valor: en nuestro caso, la rigidez máxima durante el control es km´ax =79,5ki,para conectar dos estados tales que kf=4ki.Esto es consistente con los resultados de [12]; allí la rigidez alcanzaba un valor máximo de 37ki,al conectar estados en los que kf=2ki,en un tiempo 100 veces inferior al de STEP. Continuamos con el estudio del balance energético. Vemos que hay un coste, en forma de trabajo irreversible, para la relajación rápida del sistema. Tanto calor como trabajo evolucionan independientemente: el sistema disipa calor durante todo el protocolo, por lo que disminuye monótonamente, mientras el trabajo medio sigue el comportamiento de la rigidez k, aumenta al comprimirse (kaumenta) el sistema, y disminuye al expandirse (kdisminuye) [12]. Al final del protocolo el calor se iguala a menos el trabajo según (1.28). La evolución de ambos puede conocerse evaluando sus expresiones: para el trabajo, se calcula numéricamente (3.10) entre un tiempo t0=0 y un tiempo t≤tf; en cambio, el calor medio se puede calcular a partir del trabajo y la variación media de energía potencial, utilizando la ecuación (3.18). Nosotros nos centramos en el coste, en el trabajo medio requerido, que viene dado por ⟨W⟩∗(t∗ f) = ∆F∗(t∗ f)+W∗ irr (3.21) en variables adimensionales, con ∆F∗(t∗ f) = 1 2ln(k∗ f).(3.22) El trabajo irreversible se obtiene de (3.16), evaluando numéricamente (3.17) para conocer el pará- 10Recuerde (2.8). 11Este tiempo de relajación es teq =3trx (t∗ eq =3). 22 4.1. Caso general: Transporte óptimo de masa En la referencia [13], se demostró que el problema de minimización del calor medio se puede mapear en un problema de transporte óptimo de masa, bajo un campo de velocidad que obedece la ecuación de Burgers. Este resultado se dedujo dentro del formalismo de Langevin. Una vez conocida la solución para el calor, el problema de minimización del trabajo es resuelto inmediatamente, pues dada una conexión isoterma entre dos estados de equilibrio, el trabajo medio se iguala a menos el calor medio por la relación (1.28). Aquí abordamos el problema desde el principio dentro del formalismo de Fokker-Planck. El tratamiento va a ser similar—aunque no idéntico—al que se hace en la referencia [1], donde se vinculan los resultados de la referencia [14] con los obtenidos previamente utilizando la formulación de Langevin [13]. Las expresiones que se utilizan para el trabajo medio realizado sobre el sistema, calor medio y variación de energía potencial media, han sido presentadas en la sección 1.4 de la introducción. No obstante, a continuación deducimos una expresión útil para el trabajo irreversible. El trabajo medio efectuado viene dado, en general, por la expresión ⟨W⟩(t) = ˆt t0 dt′ˆdx∂t′Ux,t′ρx,t′,(4.1) que integrando por partes puede reescribirse como ⟨W⟩(t) = ⟨U⟩(t′) t t0−ˆt t0 dt′ˆdxU x,t′∂t′ρx,t′,(4.2) Ahora, como ρ(x,t′)obedece la ecuación de FP (1.13), escribimos ⟨W⟩(t) = ⟨U⟩(t′) t t0+ˆt t0 dt′ˆdxU x,t′∂xJx,t′,(4.3) donde J(x,t′) = ρ(x,t′)v(x,t′),yv(x,t′)se definió en (1.14). Volviendo a integrar por partes se tiene que ⟨W⟩(t) = ⟨∆U⟩(t)−ˆt t0 dt′ˆdx∂xUx,t′Jx,t′,(4.4) puesto que el término de frontera se anula;1y por (1.14), el trabajo medio es finalmente ⟨W⟩(t) = ⟨∆U⟩(t)+kBTˆt t0 dt′ˆdx J x,t′∂xlnρx,t′+ˆt t0 dt′ˆdx γv2x,t′ρx,t′. (4.5) 1Se tiene que l´ ımx→∞U(x,t′)J(x,t′) = 0. 29 Los dos primeros términos se pueden identificar como la variación de energía libre ∆F(t)≡⟨∆U⟩(t)−T∆S(t),(4.6) siendo ∆S(t)la variación de entropía del colectivo asociado al sistema. El tercer término en (4.5) es la contribución irreversible al trabajo (que ya apareció en (3.15)), que es semidefinida positiva y sólo se anula cuando la evolución es reversible—en dicho caso, v(x,t′) = 0. Operando, se puede identificar con ∆Fde una forma más transparente: integrando por partes el segundo término de la ecuación (4.5), tenemos T kBˆt t0 dt′ˆdx∂t′ρx,t′lnρx,t′=−Tˆt t0 dt′∂t′−kBˆdxρx,t′lnρx,t′≡−T∆S(t), (4.7) en donde ∆S(t)es la dada por la expresión de Gibbs-Shannon. Se ha vuelto a hacer uso de la ecuación de FP.2El trabajo medio acumulado es entonces ⟨W⟩(t) = ∆F(t)+Wirr(t),(4.8) siendo el trabajo irreversible Wirr(t)≡⟨W⟩(t)−∆F(t) = γˆt t0 dt′ˆ+∞ −∞ dxv2x,t′ρx,t′≥0.(4.9) Introduciendo las relaciones (1.28) y (4.6) en la ecuación (4.9) para t=tf,obtenemos3 Wirr(tf) = T∆S−⟨Q⟩ T=T[∆S+∆Sent] = T∆Stot ≥0,(4.10) donde ∆Sent =−⟨Q⟩/Tes la variación de entropía del entorno. Este es, en efecto, el segundo principio dentro del marco teórico de la termodinámica estocástica [9]. Notemos que el trabajo irreversible se invierte en aumentar la entropía total del sistema mas entorno. El problema de minimización del trabajo medio coincide con el de minimización del trabajo irreversible (4.9), para una conexión en el tiempo ∆t=t−t0=tf.En la referencia [14] se hacía uso de cálculo de variaciones. Nosotros, sin embargo, procedemos mediante un método alternativo que involucra el empleo del método de los multiplicadores de Lagrange [1]: sea el campo auxiliar 2Por normalización, ´ρ(x,t)∂tlnρ(x,t)dx =∂t´ρ(x,t)dx =0. 3Recordemos que cuando el proceso es isotermo y conecta estados de equilibrio: ⟨∆U⟩=0. 30 arbitrario ψ(x,t),buscamos el mínimo de Ω[ρ,v,ψ]≡ˆtf 0 dt ˆ+∞ −∞ dxv2ρ+ψ[∂tρ+∂x(ρv)],(4.11) donde se incluye que ρ(x,t)ha de cumplir FP. Aparte de esta, no se considera ninguna restricción adicional al problema. Las ecuaciones de Euler-Lagrange para (4.11) son ∂tψ+1 4(∂xψ)2=0,(4.12) ∂tρ=−∂x(ρv),(4.13) ∂xψ=2v,(4.14) volviéndose a obtener la ecuación de FP en forma de ecuación de continuidad. La ecuación (4.12) es la ecuación de Burgers para el campo ψ(x,t),usando (4.14) podemos escribirla como4 ∂tv+v∂xv=0.(4.15) Esta es la ecuación que aparece en dinámica de fluidos para el campo de velocidad, el cual se encarga del transporte de masa, de acuerdo con la ecuación de continuidad (4.13) [13]. Nosotros, en lugar de trabajar con (4.13), pasamos a la función de distribución acumulada F(x,t) = ˆx −∞ dx′ρx′,t.(4.16) Aplicando esta transformación, (4.13) se convierte en ∂tF+v∂xF=0.(4.17) Para encontrar los campos F(x,t)yv(x,t),es necesario resolver conjuntamente (4.17) y (4.15). Esto lo hacemos recurriendo al método de las características:5una inspección de las ecuaciones permite determinar la pendiente de las curvas características dx dt =v,que es la misma para 4Tomando ∂xen la ecuación para ψ. 5Recordemos que las curvas características de la EDP de primer orden Cx∂xf+Ct∂tf=C,con Cx,CtyCpudiendo depender de la posición y del tiempo, son aquellas tales que la solución de la EDP no puede ser desarrollarada analíticamente. Teniendo en cuenta que d dt f(x,t) = dx dt ∂xf+∂tf,las primeras derivadas parciales del desarrollo de la solución, del problema con condiciones de contorno de Cauchy tales que f(x(t),t) = f0(t),sobre la curva (x(t),t)con f0(t) dada, no se pueden obtener si el determinante del sistema de ecuaciones formado por la EDP y d dt f0(t) = dx dt ∂xf+∂tf, verifica que △= dx dt 1 CxCt =0, De aquí se sigue que las curvas características vienen dadas por la ecuación Ctdx dt −Cx=0. 31 ambas. Además, sobre estas curvas,6tenemos que d dt v=∂tv+v∂xv=0 por (4.15); por lo cual v(x(xi,t),t) = v(xi,0)≡vi(xi),donde xies la coordenada de la curva en el instante t=0.Esto indica que las curvas características son x(xi,t) = vi(xi)t+xi.(4.18) La solución general a (4.15), sabiendo la solución sobre curvas características, es conocida— implícitamente—con el cambio v(vit+xi,t) = vi(xi)⇒ v(x,t) = vi(x−t v(x,t)) ,(4.19) siendo vi(x)una función suave arbitraria. Volviendo a (4.17) y usando (4.18); como (4.17) es una ecuación homogénea, al igual que (4.15), y ambas comparten las mismas curvas características, la solución de (4.17) sobre estas curvas es la misma que la de (4.15), es decir F(x(xi,t),t) = F(xi,0)≡Fi(xi).En t=tftenemos Ff(Θ(xi)) = Fi(xi),(4.20) con Ff(x)≡Fx,tfyΘ(xi)un mapeado tal que, sobre la curva (4.18) dada por cada valor xi, se tiene (xi,0)−→ (Θ(xi),tf).Tanto Ff(x) = ´x −∞dx′ρx′,tfcomo Fi(x) = ´x −∞dx′ρ(x′,0) son conocidas a través de las condiciones de contorno; es más, como ambas son monótonamente crecientes, es posible invertir Ffpara obtener Θ(x) = F−1 f(Fi(x)) .(4.21) Esto permite deducir la función desconocida vi(xi): imponiendo t=tfen (4.18), vi(x)es vi(x) = Θ(x)−x tf .(4.22) Con todo esto, ya está completamente determinada—al menos implícitamente—la solución (4.19) de la ecuación de Burgers. Recordamos que, dado que la estructura de las ecuaciones (4.17) y (4.15) es la misma, la solución de (4.17) tiene el mismo aspecto que (4.19), i.e. F(x,t) = Fi(x−t v(x,t)) ,(4.23) que está determinada por (4.19) una vez se conocen (4.21) y (4.22). En última instancia, por (4.16), 6La notación x(xi,t)se encarga de explicitar cuando nos desplazamos sobre curvas características. 32 derivando (4.23) se tiene ρ(x,t) = ∂xF(x,t),(4.24) lo que completa totalmente el problema de optimización. Si queremos averiguar el control óptimo U(x,t),basta con despejar en la definición del campo de velocidad (1.14),7esto es ∂xU=−γv−kBT∂xlnρ,(4.25) donde tras insertar (4.19) y (4.24) se puede integrar directamente. El mínimo trabajo irreversible Wopt irr se deduce de (4.9). Escribiendo de forma astuta esta última expresión, sobre las curvas características (4.18), Wopt irr se sigue que es Wopt irr =γˆtf 0 dt ˆ+∞ −∞ dxi dx dxi ρ(x(xi,t),t)v2(x(xi,t),t) =γˆtf 0 dt ˆ+∞ −∞ dxi dx dxi ∂xF(x(xi,t),t)Θ(xi)−xi tf2 =γˆtf 0 dt ˆ+∞ −∞ dxi∂xiFi(xi)Θ(xi)−xi tf2 =γˆtf 0 dt ˆ+∞ −∞ dxiρ(xi,0)Θ(xi)−xi tf2 =tfPopt irr ,(4.26) siendo Popt irr =γ t2 fˆ+∞ −∞ dxρ(x,0) [Θ(x)−x]2,(4.27) una potencia constante durante todo el protocolo. Vemos que el trabajo es proporcional a t−1 f, la inversa del tiempo en el que se establece la conexión. La limitación del método se halla en que para calcular el mapeado Θ(x), definido por (4.21), tiene que ser posible invertir la distribución acumulada Fipara componerla con Ff, lo cual sólo se puede hacer en casos sencillos (por ejemplo, en el potencial armónico que consideramos en la siguiente sección). Desde una perspectiva experimental, protocolos sub-óptimos que requieren de unWirr mayor, pero cercanos al umbral inferior y expresables en forma cerrada, pueden ser más que suficiente. Estos protocolos se pueden construir mediante la metodología de «soldadura» discutida en la sección 3.1—su deducción se puede consultar en la referencia [18]. 7Aplicamos ingeniería inversa introduciendo en la ecuación (1.14) la distribución óptima. 33 4.2. Caso armónico Retomando la situación que veníamos tratando en este trabajo, comprobar la fiabilidad del método—de transporte de masa óptimo—para un potencial armónico, es una tarea relativamente simple: los estados por los que pasa el sistema son gaussianos, en concreto, definimos8 ρi(x)≡ρ(x,0) = 1 p2π⟨x2⟩i exp−x2 2⟨x2⟩i,ρf(x)≡ρ(x,tf) = 1 q2π⟨x2⟩f exp"−x2 2⟨x2⟩f#. (4.28) El primer paso es obtener el mapeado Θ(x)utilizando (4.20), i.e. Ff(Θ(x)) ≡ˆΘ(x) −∞ dxf q2π⟨x2⟩f exp"−x2 f 2⟨x2⟩f#=ˆx −∞ dxi p2π⟨x2⟩i exp−x2 i 2⟨x2⟩i≡Fi(x). (4.29) Haciendo un cambio de variable tal que xi=r⟨x2⟩i ⟨x2⟩f xf, ˆΘ(x) −∞ dxf q2π⟨x2⟩f exp"−x2 f 2⟨x2⟩f#=ˆ s⟨x2⟩f ⟨x2⟩i!x −∞ dxf q2π⟨x2⟩f exp"−x2 f 2⟨x2⟩f#,(4.30) como los integrandos son iguales, esto implica Θ(x) =  s⟨x2⟩f ⟨x2⟩i x.(4.31) Note que en este caso, no es posible calcular analíticamente las distribuciones acumuladas y sin embargo, se determina el mapeado.9Insertando (4.31) en (4.22), y haciendo uso de (4.19), el campo de velocidad se extrae de v(x,t) = r⟨x2⟩f ⟨x2⟩i [x−t v(x,t)]−[x−t v(x,t)] tf ,(4.32) 8Recuerde las condiciones de contorno (3.2). 9Esto se puede hacer porque los estados tienen la misma forma. 34 de manera que, despejando, tenemos que el campo de velocidad es v(x,t) = r⟨x2⟩f ⟨x2⟩i−1 tf+ r⟨x2⟩f ⟨x2⟩i−1!t x.(4.33) Conocido v(x,t),la distribución acumulada F(x,t)se obtiene de (4.23). Derivando, llegamos por la expresión (4.24) a ρ(x,t) = ∂xFi(x−t v) = (1−t∂xv)ρi(x−t v),(4.34) y la solución es entonces ρ(x,t) = 1 p2π⟨x2⟩(t)exp−x2 2⟨x2⟩(t),(4.35) siendo x2(t) = x2i 1+ s⟨x2⟩f ⟨x2⟩i−1  t tf  2 (4.36) la evolución óptima para la varianza, que describe completamente la transformación entre estados. A continuación, se presenta una vía alternativa que permite llegar a la misma solución: partiendo de la propiedad de permanencia en estados gaussianos (demostración en la sección 2.1), se escribe la ecuación para la varianza (2.10). Por otro lado, la expresión del trabajo medio (4.1), queda ⟨W⟩(tf) = ˆtf 0 dt 1 2 d dt k(t)ˆ+∞ −∞ dxx2ρ(x,t) = ˆtf 0 dt 1 2 d dt k(t)x2(t).(4.37) Integrando por partes al igual que en, e introduciendo (2.10), el trabajo irreversible es Wirr(tf)≡⟨W⟩(tf)−∆F(tf) = γ 4ˆtf 0 dt 1 ⟨x2⟩d dt x22 (4.38) que ya se obtuvo en la sección 3.2. Tomando variaciones en (4.38), se tiene δWirr =γ 4ˆtf 0 dt ∂L ∂ω −d dt ∂L ∂˙ ωδω ,L(ω,˙ ω) = ˙ ω2 ω,(4.39) con ω≡x2yδω (0) = δω tf=0.Cuando el trabajo irreversible Wirr es mínimo, se cumple que δWirr =0 para cualquier variación arbitraria δω(t)en (0,tf).De modo que, la correspondiente 35 ecuación de Euler-Lagrange de (4.39) es ∂L ∂ω −d dt ∂L ∂˙ ω=0,(4.40) que tiene como integral primera ˙ ω∂L ∂˙ ω−L=c2,(4.41) con cuna constante. Sustituyendo L(ω,˙ ω) = ˙ ω2 ω,se tiene ˙ ω2 ω=c2⇒˙ ω 2√ω=c1,(4.42) con c1=c/2,que se integra en √ω=c1t+c2⇒ω= (c1t+c2)2,(4.43) recordando que ω≡x2,y que las constantes c1yc2se determinan imponiendo las condiciones de contorno (3.2), esto es x2(t) = x2i 1+ s⟨x2⟩f ⟨x2⟩i−1  t tf  2 .(4.44) Así reobtenemos la solución dada por (4.36).10 Previo a encontrar el protocolo adecuado para la rigidez k(t),conviene fijarse en un detalle importante del último desarrollo: no se han impuesto condiciones de contorno a la derivada de x2.De hecho, tras una rápida inspección de (4.44), se comprueba que d dt x2=0 en los tiempos extremos; esto va totalmente en contra de la condición de equilibrio de los estados conectados. La manera de evadir esta circunstancia consiste en permitir discontinuidades finitas en el control, el cual salta desde/hacia el correspondiente valor de equilibrio en el instante inicial/final—de este modo, por la ecuación de la varianza (2.10), se asegura el equilibrio. Estos saltos inicial y final en el control no contribuyen al valor del trabajo irreversible, lo que se puede ver de forma general con la ecuación (4.9)—como ρ(x,t)se mantiene constante en el salto instantáneo, y v(x,t)no diverge, no hay trabajo irreversible suministrado en dichos instantes. La presencia de discontinuidades es una característica robusta [1], no exclusiva del potencial armónico, sino propia de los métodos provistos por la teoría de control óptimo.11 Matemáticamente, esto se entiende si se piensa en el 10Este segundo método sólo es posible porque se puede escribir, a partir de la ecuación FP, una ecuación diferencial ordinaria para la varianza. 11Se exige que la función de control sea continua a trozos. 36 problema variacional: se exigen demasiadas condiciones de contorno para el protocolo que minimiza el trabajo irreversible Wirr ,puesto que la ecuación de Euler-Lagrange (4.40) es de segundo orden. En el Apéndice A de la referencia [11], se demuestra que la solución óptima implica el salto finito, y queda claro que imponer condiciones de contorno adicionales no afecta a la solución en (0,tf)—también se muestra un método de «cirugía» para regularizar estos saltos.12 La potencia irreversible es, particularizando (4.27) para el potencial armónico, e introduciendo el mapeado (4.31): Popt irr =γ t2 fˆ+∞ −∞ dxρ(x,0)x2 s⟨x2⟩f ⟨x2⟩i−1  2 =γ t2 f s⟨x2⟩f ⟨x2⟩i−1  2 x2i,(4.45) y entonces el trabajo irreversible óptimo es Wopt irr =kBTtrx tfkf ki"ski kf−1#2 ,(4.46) donde trx se definió en una de las relaciones (3.4). Haciendo uso de estas relaciones, adimensionalizamos (4.46) como Wopt irr ∗=k∗ f t∗ f"s1 k∗ f−1#2 ,(4.47) donde comparando con (3.16) vemos que el factor ηes ηopt =q1 k∗ f−12 .Esto mismo se puede reobtener introduciendo (4.44) en (4.38), e integrando esta última. 4.2.1. Conexión óptima sin restricciones Bajo la adimensionalización (3.4), y contando con las condiciones de contorno (3.2) se puede reescribir (4.44) como x2∗(τ) =  1+  1 qk∗ f−1 τ  2 ,0≤τ≤1.(4.48) 12Hacer el protocolo continuo en los extremos tiene interés práctico. 37 El protocolo óptimo que hace la conexión se conoce insertando (4.48) en (3.5), este es k∗(τ) = 1−k∗ f t∗ f1+1 √k∗ f−1τ 1 √k∗ f−1 1+1 √k∗ f−1τ2,0<τ<1.(4.49) Como discutimos arriba, en los extremos hay un salto finito tal que k∗(0) = 1 y k∗(τ=1) = k∗ f. 4.2.2. Simulaciones numéricas Recordamos que el procedimiento para simular se detalla en el Apéndice A. La ecuación (3.20) es evaluada en cada instante de tiempo con k∗(τ)dado por (4.49). Los parámetros se corresponden con los de protocolos anteriores: k∗ f=4,t∗ f=0,06.Los resultados se muestran en la Figura 4.1. El parámetro de control obedece (4.49) en el intervalo 0 <τ<1 —en τ=0+,k∗(0+) = 34,3 ; en τ=1−,k∗(1−) = 70,7—con un salto brusco aplicado en las fronteras tal que k∗(0−) = 1 y k∗(τ=1+) = k∗ f.En este caso la rigidez aumenta monótonamente en 0 <τ<1.La simulación continúa hasta el tiempo τ=1,5 para comprobar que el sistema se mantiene en equilibrio cuando τ>1.La evolución de la varianza se determina promediando, en cada instante de tiempo y sobre las Ntray =105trayectorias, las desviaciones de la posición respecto al valor medio. La varianza evoluciona continuamente desde el valor inicial σ∗2 i=1 hasta el valor final de equilibrio σ∗2 f=1/k∗ f=0,25.El acuerdo entre varianza simulada y la curva teórica (4.48) es excelente. La FDP se obtiene registrando valores de x∗en el intervalo [−3,3]: se observa que se mantiene como una gaussiana de valor medio nulo para todo instante de tiempo—se puede identificar (aunque no es fácil de ver), para τ=1 una relajación más brusca que la de los polinomios del Capítulo 3. El equilibrado se alcanza una vez más en un tiempo t∗ f=0,06,que es 50 veces inferior al de relajación por STEP. La diferencia entre (b) de la Figura 4.1 y la correspondiente en la Figura 3.1 o 3.2 no parece grande salvo por la falta de suavidad en los extremos del control. Por otro lado, si nos fijamos en el apartado (a) de las respectivas figuras, claramente los protocolos son muy dispares. El comportamiento del trabajo y del calor a lo largo del protocolo es complejo; como conocer la evolución de estas cantidades no es el objetivo, sino estudiar el coste del proceso, tratamos este comportamiento de forma cualitativa. Según la primera igualdad de la ecuación (3.9), como kes monótonamente creciente en (0,tf),el trabajo es entregado al sistema de forma continua durante el protocolo (compresión). Esto es así salvo en los tiempos extremos, cuando es inyectado instantáneamente en t=0,y cuando es extraído del sistema en t=tf; en ambos tiempos, la energía potencial ⟨∆U⟩(t)dada por (3.12), sufre un cambio inmediato con el respectivo cambio brusco de k,lo que se asocia a un trabajo. Por otro lado, de acuerdo con la expresión (3.18) y usando (4.8), el 38 Apéndice A Simulaciones numéricas con Langevin y determinación de la función de distribución Aquí se explica cómo se han hecho las simulaciones, que consisten en integrar numéricamente la ecuación de Langevin, así como la forma de construir la distribución de la variable posición. Comenzamos partiendo de la ecuación para un potencial armónico dx dt =−k(t) γx(t)+√2Dξ(t),(A.1) donde D=kBT γyξ(t)es un ruido blanco gaussiano, tal que cumple ⟨ξ(t)⟩=0,ξ(t)ξ(t′)=δ(t−t′).(A.2) Para integrar la ecuación, se puede tomar un salto de tiempo ∆tlo suficientemente pequeño, y transformarla en una ecuación en diferencias. ¿Cómo tiene que ser el ∆tpara que la solución así obtenida sea correcta? La respuesta a esto, tras una observación de (A.1), es que debemos mantener la información de la dinámica k(t), por lo que k(t)debe permanecer constante cuando t−→t+∆t. Por ello, elegimos el paso de tiempo de tal manera que k(t+∆t)−k(t) k(t)≲0,01; ∀t.No obstante aparece un segundo problema: el comportamiento discontinuo de ξ(t), lo cual obliga a hacer la integral Z(t;∆t) = ˆt+∆t t dt′ξ(t′).(A.3) Vamos a demostrar que Z(t,∆t)es una variable aleatoria gaussiana, de media nula y varianza ∆t: teniendo en cuenta que el ruido ξ(t)define una variable gaussiana en cada instante de tiempo, Z(t,∆t)es por tanto una variable gaussiana, por ser suma de estas. Adicionalmente, por la primera 45 propiedad de (A.2), su valor medio es nulo: ⟨Z(t;∆t)⟩=ˆt+∆t t dt′ξ(t′)=0.(A.4) Y por la segunda propiedad se tiene que Z2(t;∆t)=ˆt+∆t t dt′ˆt+∆t t dt′′ξ(t′)ξ(t′′)=ˆt+∆t t dt′ˆt+∆t t dt′′δ(t′−t′′) = ∆t,(A.5) por lo que Z(∆t)no depende de ty sigue una distribución normal dada por PZ(z;∆t) = 1 √2π∆texp−z2 2∆t.(A.6) Sabiendo esto, se puede escribir la ecuación integrada de (A.1) como x(t+∆t) = x(t)−k(t) γx(t)∆t+√2D∆tN(t),(A.7) donde N(τ)es una variable aleatoria distribuida normalmente, de media nula y varianza unidad, evaluada en cada instante de tiempo t. Para obtener una trayectoria completa x(t)en el intervalo temporal [0,ttray], se sigue el método de Euler: primero, dividimos el intervalo en un número de pasos Npasos =ttray/∆t; segundo, se impone una condición inicial, que en nuestro caso vendrá dada por la distribución gaussiana PX(x0) = 1 q2πkBT ki exp"−x2 0 2kBT ki#,(A.8) con ki=k(t=0); por último, de forma iterativa, calculamos xl+1=xl−k(tl) γxl∆t+√2D∆tNl,(A.9) donde el subíndice ldenota la magnitud evaluada en el tiempo tl=l∆t. Una vez que se tiene un elevado número de trayectorias Ntray ,se puede determinar la FDP agrupando los datos en histogramas. Para ello, se discretiza el intervalo [xm´ın,xm´ax]en Nbins cajas de anchura ∆xbin , tal que ∆xbin =xm´ax−xm´ın Nbins ; este intervalo se elige de manera que se asegure que más del 99% de los resultados caigan dentro. Cuando un punto cae fuera del intervalo, se contabiliza dentro de las divisiones extremas, esto asegura la normalización. El procedimiento se basa en contabilizar en cada instante de tiempo cuántas partículas, del total Ntray ,quedan dentro de ca- 46 da caja. Por ejemplo, si en el instante tlnos fijamos en la caja xa<x<xb,con ∆xbin =xb−xa suficientemente pequeño, el valor local de la FDP PXxb+xa 2,tlverifica ∆xbin PXxb+xa 2,tl=Nx∈(xa,xb) cuentas (tl) Ntray .(A.10) con Nx∈(xa,xb) cuentas (tl),el número de cuentas en la caja xa<x<xb,en el instante de tiempo tl.Así, se puede conocer aproximadamente la FDP en cada instante de tiempo tl. Esta última metodología, se emplea experimentalmente para conocer la FDP de un colectivo de partículas brownianas, ya sea cuando el sistema consta de una colección de partículas manipuladas en conjunto, o una única partícula que repite muchos ciclos de evolución. Ambos enfoques son válidos para obtener una FDP, y equivalentes cuando no hay interacciones en el caso de un sistema de partículas. 47 Bibliografía [1] D. Guéry-Odelin, C. Jarzynski, C.A. Plata, A. Prados, and E. Trizac, “Driving rapidly while remaining in control: classical shortcuts from Hamiltonian to stochastic dynamics,” Rep. Prog. Phys. 86(3), 035902 (2023). [2] N.G. Van Kampen, Stochastic Processes in Physics and Chemistry, 3rd ed. (Elsevier Science & Technology, 2007). [3] A. Einstein, and R. Brown, Investigations on the Theory of the Brownian Movement, Repr. of Methuen ed. 1926 (Dover Publ, New York, 1967). [4] S. Ciliberto, “Experiments in Stochastic Thermodynamics: Short History and Perspectives,” Phys. Rev. X 7(2), 021051 (2017). [5] K.C. Neuman, and S.M. Block, “Optical trapping,” Review of Scientific Instruments 75(9), 2787–2809 (2004). [6] I.A. Martínez, É. Roldán, L. Dinis, and R.A. Rica, “Colloidal heat engines: a review,” Soft Matter 13(1), 22–36 (2017). [7] J.-L. Barrat, and J.-P. Hansen, Basic Concepts for Simple and Complex Liquids (Cambridge University Press, 2003). [8] U. Seifert, “Stochastic thermodynamics: principles and perspectives,” Eur. Phys. J. B 64(3–4), 423–431 (2008). [9] E. Aurell, K. Gawe¸dzki, C. Mejía-Monasterio, R. Mohayaee, and P. Muratore-Ginanneschi, “Refined Second Law of Thermodynamics for Fast Random Processes,” J Stat Phys 147(3), 487–505 (2012). [10] K. Sekimoto, Stochastic Energetics (Springer Berlin Heidelberg, Berlin, Heidelberg, 2010). [11] C.A. Plata, D. Guéry-Odelin, E. Trizac, and A. Prados, “Optimal work in a harmonic trap with bounded stiffness,” Phys. Rev. E 99(1), 012140 (2019). 48 [12] I.A. Martínez, A. Petrosyan, D. Guéry-Odelin, E. Trizac, and S. Ciliberto, “Engineered swift equilibration of a Brownian particle,” Nature Phys 12(9), 843–846 (2016). [13] E. Aurell, C. Mejía-Monasterio, and P. Muratore-Ginanneschi, “Optimal Protocols and Optimal Transport in Stochastic Thermodynamics,” Phys. Rev. Lett. 106(25), 250601 (2011). [14] Y. Zhang, “Work needed to drive a thermodynamic system between two distributions,” EPL 128(3), 30002 (2020). [15] C.A. Plata, A. Prados, E. Trizac, and D. Guéry-Odelin, “Taming the Time Evolution in Overdamped Systems: Shortcuts Elaborated from Fast-Forward and Time-Reversed Protocols,” Phys. Rev. Lett. 127(19), 190605 (2021). [16] T. Schmiedl, and U. Seifert, “Optimal Finite-Time Processes In Stochastic Thermodynamics,” Phys. Rev. Lett. 98(10), 108301 (2007). [17] A. Patrón, A. Prados, and C.A. Plata, “Thermal brachistochrone for harmonically confined Brownian particles,” Eur. Phys. J. Plus 137(9), 1011 (2022). [18] Material suplementario de [15]. [19] A. Prados, “Optimizing the relaxation route with optimal control,” Phys. Rev. Research 3(2), 023128 (2021). [20] C.A. Plata, D. Guéry-Odelin, E. Trizac, and A. Prados, “Building an irreversible Carnot-like heat engine with an overdamped harmonic oscillator,” J. Stat. Mech. 2020(9), 093207 (2020). [21] C.A. Plata, D. Guéry-Odelin, E. Trizac, and A. Prados, “Finite-time adiabatic processes: Derivation and speed limit,” Phys. Rev. E 101(3), 032129 (2020). 49