Un artículo de métodos numéricos publicado en Physical Review Letters

Los que trabajamos (investigamos) en métodos numéricos y somos físicos de formación tenemos una espinita clavada, publicar en Physical Review Letters. Siempre uno quiere creerse que hace física computacional y que algún día podrá publicar en PRL, pero ojeando los índices de artículos de esta revista cada semana… y sabiendo que publicar en PRL cada día es más difícil… Pero la esperanza es lo último que se pierde y ver un artículo de métodos numéricos en PRL le da a uno una gran alegría. Nunca hubiera creído que el artículo de los físicos Mark K. Transtrum, Benjamin B. Machta y James P. Sethna, de la Universidad de Cornell, Ithaca, New York, titulado “Why are Nonlinear Fits to Data so Challenging?,” que ya ojeé hace meses en ArXiv, y al que estuve a punto de dedicar una entrada en este blog, haya acabado siendo publicado en Physical Review Letters 104: 060201, 12 Feb. 2010. ¡Enhorabuena!

El artículo estudia el ajuste de los parámetros de un modelo a partir de datos experimentales. Este problema se puede resolver fácilmente con Mathematica, Matlab y cualquier otro software científico al uso. Sin embargo, conforme el número de parámetros crece las dificultades se acumulan y los métodos numéricos convencionales que utilizan dichos programas encuentran dificultades de convergencia. Yo recuerdo una vez que estuve ajustando los parámetros de un modelo de la respiración de monos bajo estrés. Matlab (usando la variante del algoritmo de Levenberg-Marquardt desarrollado por Powell) necesitó varios días de trabajo (hace ya unos 10 años) y el resultado al final no me convenció. Ajustar, ajustaba, pero los valores de los parámetros eran poco “razonables” y no acabó de gustarme el resultado. Seguramente el problema tenía múltiples soluciones y había una solución que se me escapaba que ajustaba los datos y parecía “razonable”. La verdad sea dicha, tras cierto tiempo dedicado al asunto, preferí abandonarlo sin publicarlo. Los que sois investigadores sabéis que eso nos pasa muchas veces. Y cuando el tópico es colateral (yo estudio propagación de ondas en medios no lineales dispersivos), uno lo abandona sin remordimiento alguno.

El artículo de Transtrum et al. trata de explicar la razón por la que muchos algoritmos de optimización utilizados para ajustar parámetros de un modelo encuentran soluciones no razonables. Según ellos, es debido a que entran en una región del espacio de parámetros en la que el modelo no responde a cambios en dichos parámetros. La única opción para escapar de dichas regiones “malas” es ajustar a mano los parámetros (es decir, reparametrizar el problema con inteligencia e intuición). Los autores del artículo tratan de explicar este fenómeno aludiendo a la teoría de la interpolación para la aproximación de funciones, que subyace a todos los algoritmos de optimización. Transtrum et al. introducen una noción de curvatura (extrínseca) en el espacio de parámetros del modelo (una idea de geometría diferencial poco explotada en análisis numérico) que permite identificar estas regiones de búsqueda “malas” y la utilizan para mejorar la convergencia del algoritmo de Levenberg-Marquardt, permitiéndole salir de las mismas.

Los lectores de este blog que no hayan oido hablar del algoritmo de Levenberg-Marquardt y que nunca hayan optimizado los parámetros de un modelo a partir de datos experimentales seguramente no sabrán lo importante que es esta tarea numérica en la práctica diaria de los físicos experimentales, biólogos, economistas, etc. El mejor ajuste suele ser de tipo mínimocuadrático, se minimiza cierta función coste definida positiva, normalmente la suma de los cuadrados de las desviaciones entre el modelo y los datos experimentales. Muchos modelos tienen parámetros que dependen de forma no lineal, con lo que este proceso de minimización requiere métodos numéricos para la resolución de sistemas no lineales de ecuaciones. Encontrar la solución óptima es difícil porque el espacio de parámetros (posibles modelos) está repleto de mínimos locales que no corresponden al mejor ajuste posible (el mínimo global). Más aún, la función coste puede ser muy sensible a ciertos valores de los parámetros y muy poco sensible a otros. La idea de Transtrum et al. es que en el espacio de los parámetros hay una hipersuperficie (variedad diferencial) en la que se encuentra la solución óptima que puede pasar desapercibida al algoritmo de optimización porque tiene un grosor muy delgado para que la pueda resolver correctamente, le llaman “hipercinta” (hiper-ribbon). Para poder detectar la presencia de este tipo de “hipercintas” proponen el uso de una noción de curvatura basada en asumir que en cada paso el método numérico recorre una geodésica en el espacio de parámetros a cierta “velocidad” y “aceleración” geodésicas. La curvatura se calcula gracias a estas segundas derivadas y permite acelerar la convergencia de un método numérico (lo aplican sólo al algoritmo de Levenberg-Marquardt, LM).

¿Qué speedup obtienen con el nuevo algoritmo? Para el ejemplo que incluyen en el artículo (un problema de biología de sistemas basado en una red metabólica modelada con leyes de Michaelis-Menten) el algoritmo LM acelerado obtiene la respuesta óptima el 65% de las veces, cuando el algoritmo LM sin acelerar la obtiene sólo el 33%. ¡Un gran logro! Bueno, no tanto, calcular la solución en 2 horas en lugar de en 4 horas no parece un gran avance científico. Pero … me corroe la envidia.

¿Enviaré mi próximo artículo numérico a PRL a ver si cuela? No sé, no sé, me lo pensaré… pero tengo que venderlo bien ya que la respuesta que espero es que Physical Review Letters publica física no métodos numéricos. Bueno… casi siempre.

Publicado en PRL: Nuevo límite de exclusión combinado (CDF+DZero) para la masa del bosón de Higgs en el Tevatrón del Fermilab

En este blog ya nos hicimos eco de esta noticia en Nuevo límite de exclusión combinado (CDF+DZero) para la masa del bosón de Higgs en el Tevatrón del Fermilab, 26 Enero 2010. Ahora la noticia es que los tres artículos ya citados se han publicado en Physical Review Letters y son de acceso gratuito: [1] T. Aaltonen et al. (CDF Collaboration) “Inclusive Search for Standard Model Higgs Boson Production in the WW Decay Channel Using the CDF II Detector,” Phys. Rev. Lett. 104: 061803, February 12, 2010; [2] V. M. Abazov et al. (DØ Collaboration) “Search for Higgs Boson Production in Dilepton and Missing Energy Final States with 5.4  fb-1 of pp̅ Collisions at sqrt[s] =1.96  TeV,” Phys. Rev. Lett. 104: 061804, February 12, 2010; y [3] T. Aaltonen et al. (CDF and DØ Collaborations), “Combination of Tevatron Searches for the Standard Model Higgs Boson in the W+W Decay Mode,” Phys. Rev. Lett. 104: 061802, February 12, 2010.

Además, se ha publicado un comentario en Physics al respecto, faltaría más. En concreto, Klaus Mönig (DESY), “First bounds on the Higgs boson from hadron colliders,” Physics 3: 14, Feb. 12, 2010 [DOI, aún no funciona]. Los dos experimentos en activo del Tevatrón en el Fermilab, Chicago, han publicado los últimos resultados de su búsqueda directa del bosón de Higgs del Modelo Estándar, tanto ambas colaboraciones independientemente, como ambas uniendo sus datos. Permitidme de nuevo que os hable del bosón de Higgs usando como guía el artículo de Mönig.

En la física de partículas elementales, las interacciones entre las partículas son descritas por simetrías tipo gauge que determinan como interactúan las partículas entre sí gracias al intercambio de partículas mediadoras (portadoras) de la interacción (fuerza). Por ejemplo, la fuerza electromagnética entre partículas cargadas eléctricamente está mediada por los fotones (las partículas de luz). Las teorías gauge del Modelo Estándar tienen un gran problema: requieren que las partículas mediadoras de fuerza sean de masa nula (como el fotón). Sin embargo, sabemos que no es así ya que las partículas mediadoras de la fuerza débil (los bosones vectoriales W y Z) tienen masa en reposo no nula (por eso la fuerza débil tiene un alcance muy corto). Para resolver este problema se propuso a mediados de los 1960 el mecanismo de Higgs (recientemente se han propuesto otras alternativas). Según este mecanismo, la simetría gauge está rota, no se cumple, a baja energía. Se produce una transición de fase parecida a la congelación del agua líquida en hielo. El agua tiene simetría rotacional completa, alrededor de un punto es igual en todas direcciones, sin embargo el hielo tiene simetría cristalina hexagonal. La isotropía del agua se rompe en el hielo, que presenta un grupo de simetría más reducido. Por debajo de cierta temperatura (cero grados a presión atmosférica) el agua consume energía y se transforma en hielo. Lo mismo ocurre con la teoría electrodébil.

En física de partículas elementales la ruptura de la simetría electrodébil se produce gracias a un campo de energía (llamado campo de Higgs) que se acopla a los bosones vectoriales electrodébiles, dotando de masa a las bosones W y Z, pero no al fotón. En este proceso el campo de energía se consume, pero no completamente. Parte de este campo de energía no es consumido en el fotón (por eso permanece a baja energía sin masa). La teoría predice que dicha energía permanece a baja energía en forma de una partícula aún no descubierta denominada bosón (escalar) de Higgs. Todavía no se ha obsevado experimentalmente ninguna partícula de tipo bosón escalar. Son las partículas cuya descripción matemática es más sencilla, pero el modelo estándar nos da mucha libertad sobre sus propiedades. Podría haber un sólo bosón de Higgs neutro, que podría ser o no ser igual a su antipartícula, o también podría haber dos bosones de Higgs cargados, uno la antipartícula del otro, o incluso  muchos más. Sólo experimentalmente se podrán conocer los detalles del mecanismo de Higgs (si es que es correcto). Entre las pocas cosas que se saben hay que destacar que para que la existencia del bosón de Higgs sea compatible con el comportamiento de la fuerza débil mediada por los bosones W y Z medido en test de alta precisión del Modelo Estándar es necesario que su masa sea menor que 1 TeV/c². Poco más sabemos sobre esta partícula aún no descubierta.

Para entender el mecanismo exacto de la generación de masa en el Modelo Estándar es necesario encontrar el bosón de Higgs o bien excluir su existencia y recurrir a un mecanismo alternativo al mecanismo de Higgs. La búsqueda del bosón de Higgs se inició en los 1990 en el LEP, un colisionador de electrones y positrones en el CERN cerca de Ginebra, Suiza. Se pensaba que el bosón de Higgs tendría una masa inferior a 100 GeV/c² (similar o un poco por encima de la masa de los bosones vectoriales W y Z). En el LEP el mecanismo de producción más probable de un bosón de Higgs se basaba en la producción conjunta con un bosón Z, cuya masa es de unos 90 GeV/c². El bosón de Higgs no fue descubierto en el LEP2, que lograba colisiones con una energía en el centro de masas de hasta 205 GeV. Lo único que logró fue excluir una masa para el bosón de Higgs por debajo de 114,5 GeV/c². Los datos del LEP2 y los datos recabados en el Tevatrón del Fermilab, Chicago, parecen indicar que el Modelo Estándar le da preferencia a un bosón de Higgs de baja masa, menor de 200 GeV/c².

En los colisionadores de hadrones, como el Tevatrón (donde colisionan protones contra antiprotones) o el LHC del CERN (colisiones entre protones), también se pueden producir bosones de Higgs. Según el Modelo Estándar una pequeña fracción se produce en la misma forma que en el LEP, por la aniquilación de un par quark-antiquark produciendo un Higgs junto a un bosón W o Z. Sin embargo, la mayor cantidad de bosones de Higgs se produce por la aniquilación (colisión) de dos gluones, las partículas portadoras de la fuerza fuerte. Esto puede parecer extraño, ya que los gluones, como el fotón, no tienen masa (en reposo) y por tanto no se acoplan con el Higgs (que dota de masa en reposo a las partículas). Sin embargo, el principio de incertidumbre permite la producción de partículas intermedias con masa, que sí se acoplan al Higgs y permiten su producción en una colisión entre gluones. En el Modelo Estándar estas partículas son generalmente los quarks arriba (up), sin embargo, también podrían contribuir nuevas partículas aún desconocidas.

La tasa de producción del bosón de Higgs en el Tevatrón es mayor que en LEP ya que la energía que desarrolla en las colisiones protón-antiprotón es mucho mayor, algo menor de 2 TeV. Sin embargo, debido a la enorme complejidad de las colisiones entre hadrones (que son partículas compuestas de quarks y gluones, los llamados partones) el ruido de fondo (los resultados de colisiones que no nos interesan) es extremadamente alto y muchas de las colisiones que producen Higgs no son detectadas en los detectores de los experimentos del Tevatrón (DZero y CDF). Si el bosón de Higgs es más ligero que unos 130 GeV/c²,  se descompone principalmente en un par de quarks bottom (fondo) abajo (down) lo que se confunde fácilmente con el ruido (aparecen millones de quarks bottom abajo en las colisiones protón-antiprotón debido a la fuerza fuerte). A mayores masas, la desintegración más probable de un Higgs es en dos bosones W (mucho más fácil de discriminar en el Tevatrón). Los dos experimentos del Tevatrón, DZero CDF están estudiando ambos tipos de desintegración aunque sus resultados están sesgados hacia el segundo caso (donde su poder de discriminación es mucho mayor). Por ello, los nuevos resultados publicados en Physical Review Letters por ambos detectores se centran fundamentalmente en la desintegración de un Higgs en dos bosones W.

CDF y DZero han estudiado las desintegraciones de un Higgs cuyo estado final contiene dos bosones W, dos bosones W y un bosón Z, o tres bosones W. Para poder discriminar estas desintegraciones han considerado sólo los casos en los que dos de los bosones W se desintegran en un par de leptones, un electrón (muón o tauón) y un neutrino (que se observa como una pérdida de energía en el producto final), o en el caso de que aparezca un bosón Z, su desintegración en dos leptones. La complejidad de las colisiones protón-antiprotón obliga a utilizar sofisticadas técnicas de análisis para poder discriminar este tipo de señal dentro del ruido de fondo. Tanto CDF como DZero utilizan técnicas basadas en redes de neuronas artificiales que han sido entrenadas gracias a simulaciones de Montecarlo. Aún así, el tratamiento estadístico de los resultados obtenidos de las colisiones seleccionadas por las redes de neuronas artificiales requiere un análisis de máxima verosimili tud que tenga en cuenta todas las incertidumbres asociadas a los errores sistemáticos esperados según los resultados de las simulaciones de Montecarlo de las colisiones. El análisis es tan complejo que cada año que pasa las técnicas de análisis mejoran significativamente (sin que hayan cambiado los detectores experimentales).

Para todas las masas posibles del Higgs consideradas en el análisis experimental, la sección transversal (la tasa de producción de Higgs) es compatible con cero, es decir, no hay evidencia de que se hayan producido bosones de Higgs. Los datos permiten excluir cierta región de masas para el bosón de Higgs mediante un proceso que compara la sección transversal esperada según el Modelo Estándar y el límite experimental de dicha sección transversal observado en los experimentos. Este cociente (llamado R) permite excluir un intervalo de masas para el bosón de Higgs si su valor es menor que la unidad. La combinación de los resultados experimentales de CDF y DZero mediante complejas técnicas de análisis estadístico sólo ha permitido excluir el bosón de Higgs en un pequeño rango de masas alrededor de 160 GeV/c². Este límite de exclusión es válido sólo si se asume que el Modelo Estándar es correcto. Sin embargo, la mayoría de las extensiones del Modelo Estándar predicen una sección transversal para la producción del Higgs mayor, por lo que el límite de exclusión observado se cree que es pesimista y el intervalo calculado con otras extensiones del Modelo Estándar sería mayor.

El ajuste estadístico desarrollado por las colaboraciones CDF y DZero ha tenido en cuenta todo los datos conocidos sobre la masa del bosón de Higgs, incluyendo los datos de precisión de la teoría electrodébil, así como los límites de exclusión publicados con anterioridad. Como muestra la figura que abre esta entrada, los nuevos datos de exclusión indican que en el Modelo Estándar el rango de masas más probable para el bosón de Higgs está entre 115 y 150 GeV/c², con preferencia a masas bajas, cercanas a 115 GeV/c². Esta es la región en la que el Tevatrón (y el LHC) tiene más dificultades a la hora de encontrar un bosón de Higgs, caso de que exista en dicha región de masas. El Tevatrón funcionará sólo durante dos años (todavía no se ha asignado financiación para que funcione a partir de enero de 2012 y muchos creen que dicha financiación no será aprobada). ¿Podrá encontrar el Tevatrón al bosón de Higgs en los dos años que le restan de vida? Nadie lo sabe, lo que sí está claro es que el rango de masas de exclusión del Higgs se incrementará y quizás llegue a alcanzar la región de 120 a 160 GeV/c². Si el Tevatrón no encuentra el Higgs, la búsqueda a partir de 2012 estará en manos del LHC del CERN, que a máxima energía (14 TeV), habrá que esperar como mínimo a 2013 para ello, será capaz de descartar un bosón de Higgs con una masa menor de 1 TeV/c². Peter Higgs debe estar ansioso tras 40 años de espera. Lo mejor para él es que con toda seguridad la búsqueda terminará en la próxima década. Con 70 años (cumplidos en mayo de 2009) será testigo de excepción de su descubrimiento o de su refutación.

Publicado en Science: El nivel del mar en Mallorca hace 81000 años era 1 metro más alto del actual

Los espeleotemas (estalactitas, estalagmitas y otras formaciones minerales en cuevas) son testigos envidiables del cambio climático y en las cuevas marinas del nivel del mar. Un grupo de investigadores mallorquines, europeos y norteamericanos han estudiado el nivel del mar en el oeste del Mediterráneo gracias a los espeleotemas en cuevas costeras de la isla española de Mallorca. El estudio, publicado en Science, demuestra que el nivel de mar era un 1 metro más alto hace 81 000 años, lo que significa que la cobertura de hielo en dicha época era similar a la actual. Esto contradice los estudios previos basados en los ciclos de glaciaciones, ya que la glaciación más reciente, la wisconsiense o glaciación de Würm, tuvo que tener un final mucho más brusco y rápido de lo que se pensaba, hace unos 100 000 años. El nuevo estudio ha encontrado variaciones del nivel del mar de hasta 2 metros por siglo en los últimos 100 000 años, mucho más de lo que se pensaba. Además, el nuevo estudio parece apoyar la versión más simple de la teoría de Milankovitch para explicar el comportamiento cíclico de las glaciaciones. El artículo técnico es Jeffrey A. Dorale, Bogdan P. Onac, Joan J. Fornós, Joaquin Ginés, Angel Ginés, Paola Tuccimei, David W. Peate, “Sea-Level Highstand 81,000 Years Ago in Mallorca,” Science 327: 860-863, 12 February 2010. Nos comentan dicho artículo R. L. Edwards, “Ice Age Rhythms,” Science 327: 790-791, 12 February 2010, y Phil Berardelli, “Can Sea Level Rise and Fall With Lightning Speed?,” ScienceNOW, 11 February 2010.