Interactions between Cosmic Rays and the Atmosphere: Modeling and Practical Applications

INTERNATIONAL DOCTORAL SCHOOL OF THE USC Irma Riádigos Sánchez PhD Thesis INTERACTIONS BETWEEN COSMIC RAYS AND THE ATMOSPHERE: MODELING AND PRACTICAL APPLICATIONS Santiago de Compostela, 2021 Doctoral Programme in Nuclear and Particles Physics TESE DE DOUTORAMENTO INTERACTIONS BETWEEN COSMIC RAYS AND THE ATMOSPHERE: MODELING AND PRACTICAL APPLICATIONS Irma Riádigos Sánchez ESCOLA DE DOUTORAMENTO INTERNACIONAL DA UNIVERSIDADE DE SANTIAGO DE COMPOSTELA PROGRAMA DE DOUTORAMENTO EN FÍSICA NUCLEAR E DE PARTÍCULAS SANTIAGO DE COMPOSTELA AÑO 2021 DECLARACIÓN DO AUTOR/A DA TESE D./Dna. Irma Riádigos Sánchez Título da tese: Interactions between Cosmic Rays and the Atmosphere: Modeling and Practical Applications Presento a miña tese, seguindo o procedemento axeitado ao Regulamento, e declaro que: 1) A tese abarca os resultados da elaboración do meu traballo. 2) De ser o caso, na tese faise referencia ás colaboracións que tivo este traballo. 3) Confirmo que a tese non incorre en ningún tipo de plaxio doutros autores nin de traballos presentados por min para a obtención doutros títulos. 4) A tese é a versión definitiva presentada para a súa defensa e coincide a versión impresa coa presentada en formato electrónico E comprométome a presentar o Compromiso Documental de Supervisión no caso de que o orixinal non estea na Escola. En Santiago de Compostela, 22 de Decembro de 2021. Sinatura electrónica AUTORIZACIÓN DO DIRECTOR / TITOR DA TESE INTERACTIONS BETWEEN COSMIC RAYS AND THE ATMOSPHERE: MODELING AND PRACTICAL APPLICATIONS D. Vicente Pérez Muñuzuri D. Diego González Díaz INFORMA/N: Que a presente tese, correspóndese co traballo realizado por Dna. Irma Riádigos Sánchez, baixo a miña dirección/titorización, e a utorizo a súa presentación , considerando que reúne os r equisitos esixidos no R egulamento de Estudos de Doutoramento da USC, e que como director desta non incorre nas causas de abstención establecidas na Lei 40/2015. De acordo co indicado no Regulamento de Estudos de Doutoramento, declara tamén que a presente tese de doutoramento é idónea para ser defendida en base á modalidade de Monográfica con reprodución de publicacións , nos que a participación da doutoranda foi decisiva para a súa elaboración e as publicacións se axustan ao Plan de Investigación. En Santiago de Compostela, 22 de decembro de 2021 Fdo. Vicente Pérez Muñuzuri Fdo. Diego González Díaz Las cosas podr´ ıan haber sucedido de cualquier otra manera y, sin embargo, sucedieron as´ ı. - Miguel Delibes, en su novela “El Camino” Summary This thesis collects the research work developed in the field of Physics about the different interaction processes between cosmic radiation and the atmosphere. It examines the experimental aspects of the object of study by exploring its practical applications for atmospheric sciences, as well as giving a theoretical viewpoint obtained through numerical modeling as it is necessary to validate results or test new hypotheses. Nowadays, the research of cosmic rays is no longer devoted exclusively to the field of astronomy or astrophysics. During the last few years, a large number of significant results have opened the door to everyday applications of cosmic rays. One of the most remarkable cases was the discovery of a hidden chamber in the Great Pyramid of Giza using muon tomography, a technique similar to radiography that allows discerning the internal structure of dense objects. Other related practical applications also include the use of cosmic rays in volcanology to obtain images of volcanos’ innards or the detection of possible nuclear waste in the transport of cargo containers. All this has been possible thanks to progress in the development of new detectors and measurement techniques in nuclear and particle physics. Cosmic rays, contrary to what their name suggests, are not rays but radiation in the form of high-energy subatomic particles that reach the Earth from outer space in all directions. Their origin is very diverse, the ones with the lowest energies come from the Sun, whereas the most energetic cosmic rays originate in other parts of our Galaxy and even in much more distant places, such as other galaxies. The most energetic cosmic radiation has its origin in the most violent and extreme processes in the Universe, such as supernovae or black holes. When these high-energy particles (also called primary cosmic rays) reach the atmosphere, they immediately interact with air molecules triggering a series of nuclear reactions from which new particles emerge, which in turn repeat the same process giving rise to a cascade of secondary particles that travels through the atmosphere until it reaches the Earth’s surface. Cosmic rays are therefore ubiquitous particles in the atmosphere and, in particular, muons, one of the products created when a primary cosmic ray hits atmospheric nuclei. This kind of particle has a great penetrating power that makes it perfect for the applications mentioned above, especially because muons can pass through objects without damaging them. Moreover, they are of natural origin and do not depend on any artificial source for their generation, so they are available anywhere and any time. In the first part of this dissertation, we consider the use of cosmic rays for atmospheric monitoring, in particular as a tool for measuring the vertical profile of atmospheric temperatures. The idea is to employ cosmic radiation traversing the atmosphere analogously to how weather satellites work. Satellites are able to record the temperature of different layers of the atmosphere by measuring the electromagnetic radiation of different frequencies radiating from it, which depends on the atmospheric state. Thus, we propose to use surface cosmic-ray measurements performed at different observation angles and energies for the same purpose. This idea arises Summary from the discovery of the correlation between variations of the cosmic-ray flux measured both at the Earth’s surface and underground with temperature variations at different heights. The cause lies in the development of cosmic-ray cascades throughout the atmosphere, which depend on the air density along the path, affecting the production and absorption of secondary particles. As a consequence, cosmic-ray rates measured at the surface are not constant and vary over time in correlation with the atmospheric temperature profile and atmospheric pressure. In this work, we analyze the experimental data obtained with a high resolution 2 m2 cosmic-ray detector located at the Faculty of Physics (University of Santiago de Compostela). The objective is to study in detail the variations of cosmic-ray rates measured with the device in order to determine and characterize variations of atmospheric origin in the data. Measurements of correlations between cosmic rays and temperature are usually performed with underground detectors with large areas and volumes or at the surface with detectors of moderate sizes. The former type of device is installed underground to reduce the influence of cosmic-ray radioactivity on the measurements since the ultimate goal of their research is usually related to particle, nuclear, neutrino, or dark matter physics. Regarding ground-based detectors, they have been placed in countless locations around the world. In fact, a Global Muon Detector Network (GMDN) was created to continuously monitor cosmic-ray variations. This global network is used for space weather applications, such as forecasting large geomagnetic storms. Their work also covers the analysis of the atmospheric effect in the observed cosmic-ray data, however, the reason for its characterization is the removal of such effect from the data so that they can observe only the variations associated with space weather. Our challenge is to achieve the opposite, we want to isolate the atmospheric variations present in the data obtained with a small ground-based detector. It should be noted that this is the first time that a multigap timing RPC detector (MtRPC) has been used in the study of such correlations. In addition, we have implemented a program for the simulation of cosmic-ray air showers that allows the introduction of real atmospheric profiles. We simulate the atmospheric effects to corroborate and understand the experimental results. Next, we focus on a practical application by designing a monitoring station aimed to obtain the vertical temperature profile of the atmosphere from cosmic-ray measurements. To date, despite the improved understanding of the influence of the atmosphere on cosmic-ray rates, the technological potential of such a possibility has been barely explored. In this thesis, the limits of the cosmic-ray inversion problem for temperature estimation are examined, presenting a configuration that combines a ground-based and an underground detector placed at an optimal depth. In addition, unlike previous works, we use the angular information. The last part of the thesis is devoted to a slightly different research topic: the influence of cosmic rays on atmospheric processes. In particular, we study how the ionization produced by cosmic rays as they pass through the atmosphere affects cloud formation. At the end of the last century, significant correlations were found between solar activity and global cloud cover. The proposed hypothesis argued that during periods of low solar activity, since the Sun’s magnetic field weakens, cosmic radiation of galactic origin can more easily penetrate the solar system and reach the Earth with a higher flux. Hence, atmospheric ionization increases and through some mechanism that is not yet fully understood is able to favor cloud formation. During peaks of maximum activity, the opposite will occur. The finding of these correlations, if the causality was confirmed, would mean that cosmic rays could play a role in climate variation because any significant alteration in global cloud cover modifies the terrestrial albedo causing changes in global warming or cooling. One of the biggest uncertainties in climate predictions vi IRMA RI´ ADIGOS S´ ANCHEZ is the influence of clouds and how they vary under different conditions. The reason is that contrary to what one might think, the exact mechanisms of cloud formation are hardly known in detail. These uncertainties lead to a large dispersion in the predictions regarding the average temperature increase by the end of the century, which ranges from 1.4◦C to 4.5◦C. In the latter case, we would be talking about an extreme scenario with catastrophic consequences for the climate and human beings. It is therefore vital to forecast temperature change as accurately as possible. Thus, it is essential to go deeper into the study of cloud formation processes. Cloud formation is caused by the presence in the air of small particles (or “cloud seeds”) that act as condensation nuclei for water vapor. These tiny particles are atmospheric aerosols. However, there are many unanswered questions about how aerosols form in the atmosphere and how they affect clouds. Recent studies have found strong correlations between solar activity and aerosol properties, even though other similar analyses were not able to reproduce these results. Thus, this topic is quite controversial and has not yet helped to clarify the link between cosmic rays and clouds. In this thesis, we study the effect of charged aerosols due to cosmic-ray ionization in the condensation and coagulation processes of aerosols. We implement some mathematical models and numerical methods in order to carry out computational simulations of complex atmospheric physics and chemistry. For such purpose, we use the state-of-the-art GEOS-Chem atmospheric simulation program that includes a highly accurate microphysics scheme for the description of aerosol growth processes. Furthermore, the program features the option to vary atmospheric ionization to simulate a period of maximum or minimum solar activity. We launch pairwise simulations to contrast the results of both types of scenarios. The difference between both cases will give us information about the relevance of the change of the cosmic-ray flux due to solar activity in the formation of aerosols. vii Resumen Esta tesis recoge el trabajo de investigaci´ on desarrollado en el campo de la F´ ısica sobre distintos procesos de interacci´ on que se dan entre la radiaci´ on c´ osmica y la atm´ osfera. Para ello, se enfoca el objeto de estudio desde un punto de vista experimental, explorando sus aplicaciones pr´ acticas para las ciencias atmosf´ ericas, y tambi´ en desde un punto de vista m´ as te´ orico a trav´ es de la modelizaci´ on num´ erica, necesario para validar resultados o comprobar nuevas hip´ otesis. En la actualidad, el estudio de los Rayos C´ osmicos ya no se centra exclusivamente en el campo de la astronom´ ıa o astrof´ ısica. Durante los ´ ultimos a˜ nos, se han obtenido una gran cantidad de resultados significativos que han abierto la puerta a la aplicaci´ on de los rayos c´ osmicos en la vida cotidiana. Uno de los casos con m´ as repercusi´ on ha sido el descubrimiento de una c´ amara oculta en la Gran Pir´ amide de Guiza a trav´ es de la tomograf´ ıa de muones, una t´ ecnica similar a la radiograf´ ıa que permite discernir la estructura interna de objetos densos. Entre las aplicaciones pr´ acticas relacionadas tambi´ en cabe mencionar su uso en la vulcanolog´ ıa para obtener im´ agenes del interior de los volcanes, o la detecci´ on de posibles residuos nucleares en el transporte de grandes contenedores. Todo ello ha sido posible gracias al avance en el desarrollo de nuevos detectores y t´ ecnicas de medici´ on en la f´ ısica nuclear y de part´ ıculas. Los rayos c´ osmicos, al contrario de lo que su nombre sugiere, no son rayos, sino que son una radiaci´ on en forma de part´ ıculas subat´ omicas de alta energ´ ıa que llegan a la Tierra procedentes del espacio exterior y en todas direcciones. Su origen es muy diverso, las de menor energ´ ıa proceden del Sol, mientras que las de mayor energ´ ıa se originan en otras partes de nuestra Galaxia e incluso en lugares mucho m´ as distantes, como en otras galaxias. La radiaci´ on c´ osmica m´ as energ´ etica tiene su origen en los procesos m´ as violentos y extremos del Universo, tales como supernovas o agujeros negros. Cuando estas part´ ıculas de gran energ´ ıa (tambi´ en denominados rayos c´ osmicos primarios) entran en contacto con la atm´ osfera, interaccionan inmediatamente con las mol´ eculas del aire desencadenando una serie de reacciones nucleares de las que emergen nuevas part´ ıculas, que a su vez repiten el proceso generando una cascada de part´ ıculas secundarias que viaja por la atm´ osfera hasta llegar a la superficie terrestre. Por lo tanto, los rayos c´ osmicos son part´ ıculas omnipresentes en la atm´ osfera y, en especial, los muones, uno de los productos que se crean cuando un rayo c´ osmico primario impacta con los n´ ucleos atmosf´ ericos. Esta clase de part´ ıculas tiene un gran poder de penetraci´ on que hace que sean perfectos para las aplicaciones mencionadas con anterioridad, sobre todo, porque atraviesan los objetos sin da˜ narlos. Asimismo, son de origen natural y no dependen de ninguna fuente artificial para su generaci´ on, por lo que est´ an disponibles en todo momento y en cualquier parte. En la primera parte de esta tesis, se considera el uso de los muones para la monitorizaci´ on atmosf´ erica, en concreto, como herramienta para medir el perfil de temperatura de la atm´ osfera. La idea es utilizar la radiaci´ on c´ osmica que atraviesa la atm´ osfera de forma an´ aloga al funcionamiento de los sat´ elites de observaci´ on. En su caso, los sat´ elites son capaces de Resumen registrar la temperatura de diferentes capas de la atm´ osfera a trav´ es de la medici´ on de la radiaci´ on electromagn´ etica de diferentes frecuencias que radia de esta y que depende de su estado. De este modo, se propone usar las medidas de rayos c´ osmicos en superficie realizadas a diferentes ´ angulos de observaci´ on y energ´ ıas (equivalente a profundidad bajo tierra) para el mismo prop´ osito. Esta idea surge del descubrimiento de la correlaci´ on de las variaciones de rayos c´ osmicos medidos en superficie y bajo tierra con las variaciones de temperatura a diferentes alturas. Su causa radica en la evoluci´ on de las cascadas de rayos c´ osmicos a lo largo de la atm´ osfera, que dependen de la densidad del aire a su paso, afectando a la producci´ on y absorci´ on de las part´ ıculas secundarias generadas en las diferentes interacciones nucleares. Como consecuencia, las tasas de rayos c´ osmicos medidas en superficie no son constantes y var´ ıan con el tiempo en correlaci´ on con las temperaturas y la presi´ on atmosf´ erica. La clave est´ a en que los muones c´ osmicos de diferentes energ´ ıas se ven afectados de forma diferente por las variaciones de temperatura, siendo los rayos c´ osmicos m´ as energ´ eticos mucho m´ as sensibles a la temperatura de la estratosfera, por ejemplo. Esta peculiaridad se puede aprovechar para construir un modelo que permita resolver el problema inverso y reconstruir la temperatura de la atm´ osfera a partir de la informaci´ on proporcionada por los diferentes “canales” de detecci´ on de rayos c´ osmicos. Los “coeficientes de temperatura” son unos pesos que se calculan de forma te´ orica y que proporcionan informaci´ on sobre c´ omo afecta la variaci´ on de temperatura en cada capa de la atm´ osfera a la tasa medida a un cierto nivel de observaci´ on. Adem´ as, estos coeficientes dependen del ´ angulo de observaci´ on y la energ´ ıa de los muones. A lo largo de la tesis, estos coeficientes juegan un papel fundamental para el desarrollo de las diferentes metodolog´ ıas de an´ alisis. En este trabajo, se analizan, por un lado, los datos experimentales obtenidos con un peque˜ no detector de rayos c´ osmicos de alta resoluci´ ony2m2de superficie situado en la facultad de F´ ısica (Universidad de Santiago de Compostela). El objetivo consiste en analizar en detalle las variaciones de las tasas de rayos c´ osmicos medidas con el dispositivo con el fin de determinar y caracterizar las variaciones de origen atmosf´ erico en los datos. Las medidas de las correlaciones entre rayos c´ osmicos y temperatura se realizan normalmente con detectores bajo tierra que cuentan con grandes ´ areas y vol´ umenes de detecci´ on o en la superficie con detectores de tama˜ nos m´ as reducidos. En relaci´ on al primer tipo de detectores, se instalan bajo tierra para reducir la influencia de la radiactividad de los rayos c´ osmicos en las medidas ya que el objetivo de estas investigaciones est´ a relacionado con la f´ ısica nuclear, de part´ ıculas, de neutrinos o la f´ ısica de la materia oscura. En cuanto a los detectores en superficie, se instalaron una gran cantidad en diversas localizaciones alrededor del mundo. De hecho, la Red Global de Detecci´ on de Muones (en ingl´ es, GMDN) fue creada para monitorizar de forma continuada las variaciones de rayos c´ osmicos. Un detector en superficie es mucho m´ as sensible a las variaciones atmosf´ ericas, adem´ as tambi´ en est´ a afectado por fen´ omenos relacionados con la actividad solar. Esta red global de detectores se usa para aplicaciones de clima espacial, como por ejemplo, para la predicci´ on de grandes tormentas geomagn´ eticas. Tambi´ en analizan el efecto atmosf´ erico presente en los datos medidos de rayos c´ osmicos con el objetivo de eliminarlo y as´ ı poder observar solamente las variaciones asociadas al clima espacial. Por lo tanto, mientras que otros experimentos est´ an interesados en el estudio de las correlaciones atmosf´ ericas para su posterior eliminaci´ on de las medidas, nuestro reto consiste en conseguir lo contrario, aislar las variaciones atmosf´ ericas presentes en las medidas obtenidas con un detector peque˜ no en superficie. Cabe destacar que esta es la primera vez que se consigue reproducir tales correlaciones con un detector basado en c´ amaras de placas resistivas de m´ ultiples bandas x IRMA RI´ ADIGOS S´ ANCHEZ (MtRPC). Esta tecnolog´ ıa consiste en detectores gaseosos de respuesta r´ apida formado por dos placas paralelas cargadas de forma opuesta y hechas con un material altamente resistivo. Entre las dos placas se encuentra el gas, donde interaccionar´ a la part´ ıcula incidente. Esta clase de detectores tiene la caracter´ ıstica de proporcionar una buena resoluci´ on temporal y espacial, muy necesario para un dispositivo de peque˜ no tama˜ no cuyo objetivo es la detecci´ on de rayos c´ osmicos en superficie. Para este trabajo, se mide la correlaci´ on entre las tasas medidas y la presi´ on atmosf´ erica en superficie, que se caracteriza con la intenci´ on de eliminar las variaciones relacionadas con ella de las tasas y poder as´ ı estudiar mejor las variaciones de temperatura, que podr´ ıan quedar apantalladas en caso contrario. El uso de un algoritmo espec´ ıfico para dicha tarea tambi´ en nos permitir´ a eliminar de los datos las variaciones relacionadas con eventos interplanetarios derivados de la actividad solar. Finalmente, se estimar´ an los coeficientes de temperatura espec´ ıficos del detector de forma experimental aplicando una t´ ecnica estad´ ıstica basada en el An´ alisis de Componentes Principales. Esto nos permitir´ a discernir las variaciones de temperatura inherentes en los datos, tanto estacionales como otras variaciones inusuales con periodos m´ as cortos duraci´ on. De forma adicional, se ha desarrollado la implementaci´ on de un programa para la simulaci´ on de cascadas de rayos c´ osmicos en la atm´ osfera que permite introducir datos reales de perfiles atmosf´ ericos. Otro software de simulaci´ on de cascadas disponible (AIRES, CORSIKA,...) tiene la desventaja de usar un modelo de atm´ osfera est´ andar, el cual no permite el estudio de los fen´ onemos que intentamos observar de forma experimental. Es por ello que nuestra intenci´ on consiste en desarrollar un c´ odigo nuevo desde cero para poder validar la parte experimental. Por lo tanto, se han escrito una serie de rutinas simplificadas que permiten simular las cascadas de part´ ıculas en la atm´ osfera. Con esta metodolog´ ıa, intentaremos entender lo que sucede en la atm´ osfera y cu´ ales son los factores relevantes que afectan la evoluci´ on de las cascadas. El prop´ osito final es obtener un programa que nos permita introducir datos reales de la atm´ osfera para simular el flujo de muones que alcanzar´ ıa nuestro detector. En la segunda parte de la tesis, m´ as enfocada en la aplicaci´ on pr´ actica, nos centramos en el estudio del desarrollo de una estaci´ on de monitorizaci´ on para la obtenci´ on de la temperatura de la atm´ osfera a partir de rayos c´ osmicos. Hasta la fecha, a pesar de la mejora en la comprensi´ on de la influencia de la atm´ osfera en las tasas de rayos c´ osmicos, apenas se ha explorado el potencial tecnol´ ogico de dicha aplicaci´ on. En esta tesis, se examinan los l´ ımites del problema de inversi´ on de rayos c´ osmicos para la estimaci´ on de la temperatura, presentando como propuesta una configuraci´ on que combina una estaci´ on de detecci´ on en superficie y otra bajo tierra a una profundidad ´ optima. Adem´ as, a diferencia de anteriores trabajos, usamos la informaci´ on angular de las medidas. La metodolog´ ıa del estudio consiste en simular tasas de rayos c´ osmicos que contengan las variaciones inducidas de temperatura, utilizando como datos de entrada los coeficientes de temperatura te´ oricos junto con perfiles reales de temperatura obtenidos de la base de datos de rean´ alisis del ERA5 (ECMWF). Para obtener una muestra realista de datos, tambi´ en se incluyen las fluctuaciones de origen estad´ ıstico as´ ı como el efecto de la absorci´ on de los muones en la roca. La serie temporal resultante se usa en el problema inverso para obtener el perfil de temperaturas que ser´ a comparado con los datos de temperatura originales. En la resoluci´ on del problema inverso, se analizan diferentes escenarios y los resultados se contrastan con trabajos previos. El objetivo de este estudio es establecer una l´ ınea de trabajo para futuros proyectos experimentales que tengan como meta la construcci´ on de una estaci´ on de sondeo de la atm´ osfera a partir de rayos c´ osmicos. xi Resumen La ´ ultima parte de la tesis est´ a dedicada a una tem´ atica de investigaci´ on ligeramente distinta: la influencia de los rayos c´ osmicos en los procesos atmosf´ ericos. En concreto, estudiamos c´ omo la ionizaci´ on que los rayos c´ osmicos producen a medida que atraviesan la atm´ osfera afecta a la formaci´ on de nubes. A finales del siglo pasado, se observaron correlaciones significativas entre la actividad solar y la cobertura global de nubes. La hip´ otesis planteada razona que en periodos de baja actividad solar, puesto que el campo magn´ etico del Sol se debilita, la radiaci´ on c´ osmica de origen gal´ actico puede penetrar con m´ as facilidad en el Sistem Solar y llegar hasta la Tierra con un flujo mayor. Como consecuencia, la ionizaci´ on atmosf´ erica aumenta y mediante alg´ un mecanismo a´ un desconocido es capaz de favorecer la formaci´ on de nubes. As´ ı, en periodos m´ ınimos de actividad solar, la cobertura de nubes del planeta ser´ a mayor, mientras que en picos de actividad m´ axima ocurrir´ a lo contrario. El hallazgo de estas correlaciones, de confirmarse la causalidad, significar´ ıa que los rayos c´ osmicos estar´ ıan jugando su papel en la variaci´ on del clima. Esto se debe a que cualquier alteraci´ on significativa en la cobertura global de nubes modifica el albedo terrestre provocando cambios en el calentamiento o enfriamiento del planeta. Una de las mayores incertidumbres en las predicciones clim´ aticas es la influencia de las nubes y c´ omo estas var´ ıan bajo diferentes condiciones. El motivo es que al contrario de lo que uno pueda pensar, apenas se conocen con detalle cu´ ales son los mecanismos exactos por los cuales se forman las nubes. Estas incertidumbres provocan una gran dispersi´ on en los resultados de las predicciones del aumento de la temperatura media a finales de siglo, que var´ ıan entre el 1.5◦y 4.5◦C. Cabe destacar que este rango de incertidumbre en las predicciones no se ha conseguido reducir a lo largo de los ´ ultimos a˜ nos. En el segundo caso pronosticado, estar´ ıamos hablando de un escenario extremo con consecuencias catastr´ oficas para el clima y el ser humano. Por lo tanto, resulta vital proyectar con la mayor exactitud posible el cambio en la temperatura. As´ ı que es imprescindible ahondar en el estudio de los procesos de formaci´ on de nubes. La formaci´ on de las nubes se produce gracias a la presencia en el aire de peque˜ nas part´ ıculas (o “semillas de nube”) que act´ uan como n´ ucleos de condensaci´ on para el vapor de agua. Estas min´ usculas part´ ıculas son los aerosoles atmosf´ ericos, que pueden ser de origen natural o antropog´ enico, como por ejemplo, part´ ıculas de sal procedentes del mar, material volc´ anico, polvo del desierto, productos de incendios forestales o la quema de combustibles, etc. Sin ellos, las nubes no existir´ ıan en la Tierra y el clima ser´ ıa radicalmente diferente. Adem´ as, no existir´ ıa la vida. Sin embargo, existen muchas preguntas sin resolver acerca de c´ omo se forman los aerosoles en la atm´ osfera y su efecto en las nubes. En los ´ ultimos 20 a˜ nos, se ha estudiado m´ as a fondo las correlaciones entre la actividad solar y las nubes, incluyendo tambi´ en los aerosoles en el an´ alisis. Mientras que algunos trabajos encontraron correlaciones significativas entre datos de cobertura de nubes y propiedades de los aerosoles, otros an´ alisis similares no fueron capaces de reproducir dichos resultados. As´ ı, este tema resulta bastante controvertido y a´ un no ha ayudado a esclarecer el v´ ınculo entre rayos c´ osmicos y nubes. Sin embargo, un imporante experimento llevado a cabo en el CERN y denominado CLOUD, tiene como objetivo estudiar dentro una gran c´ amara con condiciones atmosf´ ericas controladas la creaci´ on y el crecimiento de los aerosoles. Una de las misiones principales de este proyecto consite en investigar la relaci´ on que existe entre los rayos c´ osmicos y las nubes. Para ello, se utiliza como fuente de radiaci´ on part´ ıculas aceleradas en el sincrotr´ on para emular la ionizaci´ on producida por los rayos c´ osmicos gal´ acticos en la atm´ osfera. Con los primeros resultados publicados en el 2011, este experimento se ha convertido en uno de los xii IRMA RI´ ADIGOS S´ ANCHEZ primeros en demostrar una relaci´ on directa entre radiaci´ on c´ osmica y aerosoles. Aunque los resultados demuestran que los rayos c´ osmicos no juegan un papel fundamental en el cambio clim´ atico, s´ ı son relevantes en una peque˜ na proporci´ on para el proceso de nucleaci´ on los aerosoles. Y se demuestra que bajo ciertas condiciones atmosf´ ericas, pueden favorecer de manera significativa su crecimiento. Por lo tanto, no se debe descartar su efecto. Los modelos emp´ ıricos derivados de este experimento han sido puestos a prueba en modelos inform´ aticos de qu´ ımica atmosf´ erica para simular con precisi´ on los procesos de formaci´ on de los aerosoles y su influencia en la nubosidad. Con ello se ha observado que la ionizaci´ on afecta a los peque˜ nos aerosoles pero estos cambios no son lo suficientemente eficaces como para trasladarse a las part´ ıculas m´ as grandes que forman los n´ ucleos de condensaci´ on. A pesar de todo ello, esto no descarta por completo la influencia de los rayos c´ osmicos en los aerosoles. Otros procesos como la condensaci´ on o la coagulaci´ on de los aerosoles tambi´ en parecen estar afectados por la ionizaci´ on y han sido explorados en menor medida. En la ´ ultima parte de la tesis, vamos a tener en cuenta estos dos ´ ultimos procesos y los implementamos en GEOS-Chem, un programa de simulaci´ on atmosf´ erica que incluye un esquema de microf´ ısica de gran precisi´ on en la descripci´ on de los procesos de crecimiento de los aerosoles. Este modelo es uno de los m´ as completos y accesibles que se pueden encontrar, adem´ as de que est´ a elaborado por cientos de cient´ ıficos alrededor del mundo. Esto hace que sea uno de los m´ as actualizados y complejos de su campo. En concreto, GEOS-Chem es un modelo dise˜ nado para el transporte qu´ ımico tridimensional que permite realizar simulaciones de la composici´ on atmosf´ erica a una escala global o regional. Adem´ as, se puede acoplar con otros modelos clim´ aticos o meteorol´ ogicos, como puede ser el modelo WRF (uno de los m´ as utilizados en el mundo para la predicci´ on regional a corto plazo). Otra de las caracter´ ısticas m´ as destacables es que ya tiene implementado el efecto de los iones generados por la radiaci´ on c´ osmica en el proceso de nucleaci´ on de los aerosoles. La nucleaci´ on es una de las fuentes m´ as importantes de part´ ıculas atmosf´ ericas y consiste en la agregaci´ on de peque˜ nos conglomerados moleculares desde la fase gaseosa y se sabe que este proceso puede ser estimulado por la presencia de iones. En este caso, el modelo GOES-Chem cuenta con una parametrizaci´ on de la nucleaci´ on que depende de varios par´ ametros entre los que se encuentran la tasa de ionizaci´ on atmosf´ erica. La nucleaci´ on es un proceso que involucra a las part´ ıculas m´ as peque˜ nas, por el contrario, la condensaci´ on y coagulaci´ on son los procesos responsables de que los peque˜ nos grupos de part´ ıculas crezcan a tama˜ nos superiores. Si tenemos en cuenta la presencia de iones atmosf´ ericos, estos condensan sobre los aerosoles proporcion´ andoles carga. As´ ı, los aerosoles pueden acumular un variado n´ umero de cargas en su superficie. Como consecuencia, esta distribuci´ on de cargas en los aerosoles va a afectar al proceso de coagulaci´ on, puesto que si dos aerosoles que colisionan tienen cargas de id´ entico signo, aparecer´ a una fuerza de repulsi´ on que inhibir´ a su uni´ on. Por el contrario, si los aerosoles transportan cargas de signo contrario, su coagulaci´ on estar´ a m´ as favorecida. De esta forma, resulta relevante incorporar la distribuci´ on de cargas en el modelo. Sin embargo, introducir de forma expl´ ıcita la distribuci´ on de carga de las part´ ıculas es m´ as complejo y har´ ıa que el c´ alculo fuese extremadamente lento. Esto se debe a que el modelo divide la distribuci´ on de tama˜ no de los aerosoles en cuarenta bins y para cada uno resuelve un par de ecuaciones en cada paso de tiempo, una que calcula el n´ umero de aerosoles en el bin y otra que calcula su masa. As´ ı que el modelo tiene que resolver 80 ecuaciones en cada paso de tiempo y en cada punto del espacio (o malla espacial del modelo). Si suponemos que las part´ ıculas solo pueden transportar una carga negativa o positiva, eso significar´ ıa tener que triplicar los bins para que tambi´ en se tenga en cuenta el n´ umero de part´ ıculas cargadas en xiii Resumo expl´ ıcita de todos os bins de carga. Ademais, tam´ en permite acelerar a velocidade de c´ omputo. O noso obxectivo consite en implementar o c´ alculo das cargas e o seu efecto na coagulaci´ on para ver como a ionizaci´ on afecta ao crecemento dos aerosois. Unha vez implementados os cambios oportunos no c´ odigo, posto que o programa ten como caracter´ ıstica a opci´ on de variar a ionizaci´ on atmosf´ erica inducida polos raios c´ osmicos para simular un per´ ıodo de actividade solar m´ axima ou m´ ınimo, corremos simulaci´ ons a pares para contrastar os resultados das d´ uas situaci´ ons. A diferenza entre ´ ambalas d´ uas daranos informaci´ on acerca da relevancia do cambio do fluxo de raios c´ osmicos debido a actividade solar na formaci´ on dos aerosois. xx Nomenclature and Abbreviations αTEffective Temperature Coefficient αMSS Mass-weighted Temperature Coefficient µ±Muon π±Charged Pion Eth Threshold Energy RcGeomagnetic Cutoff Rigidity Te f f Effective Temperature WTTemperature Coefficient CCN Cloud Condensation Nuclei CME Coronal Mass Ejection CR Cosmic Rays EAS Extensive Air Shower FD Forbush Decrease GCR Galactic Cosmic Rays GeV Gigaelectronvolt GMDN Global Muon Detector Network mwe meter water equivalent PCR Principal Component Regression RMSE Root-mean-square Error RPC Resistive Plate Chamber SSW Sudden Stratospheric Warming UHECR Ultra-high-energy cosmic ray List of Figures 1.1 The cosmic ray energy spectrum . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 AirShowerillustration .............................. 5 1.3 Air Shower simulations from Corsika . . . . . . . . . . . . . . . . . . . . . . 6 1.4 Scheme comparing a multi-directional muon telescope with a muon hodoscope 9 1.5 Scheme of cosmic ray detection methods . . . . . . . . . . . . . . . . . . . . . 11 1.6 Schematic view of a Resistive Plate Chamber . . . . . . . . . . . . . . . . . . 13 1.7 Schematic representation of the atmospheric temperature effect on the muon intensity observed at surface and underground using the products αi·(∆T)j. . 20 1.8 Differential temperature coefficients WTfor vertical direction (θ=0◦) at several thresholdenergies................................. 22 1.9 Illustration of the temperature effect . . . . . . . . . . . . . . . . . . . . . . . 23 1.10 Illustration of the possible link between GCR and cloud-covering . . . . . . . . 26 1.11 Schematic diagram of aerosol-related processes . . . . . . . . . . . . . . . . . 29 2.1 Schematics of the TRAGALDABAS detector . . . . . . . . . . . . . . . . . . 32 2.2 Photo of the TRAGALDABAS detector . . . . . . . . . . . . . . . . . . . . . 34 2.3 Experimental observations of cosmic ray rate for single vertical tracks along with ground-level pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.4 Examples of ERA-INTERIM atmospheric temperature profiles for Santiago de Compostela .................................... 36 2.5 Temperature correlations between atmospheric layers . . . . . . . . . . . . . . 37 2.6 Temperature correlations between atmospheric layers for a region above 40◦. 38 2.7 Example of the barometric coefficient regression fit . . . . . . . . . . . . . . . 39 2.8 Temperature time series of the atmosphere in Santiago de Compostela from October 2015 to January 2017 . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.9 Temperature anomalies for a Sudden Stratospheric Warming event . . . . . . . 41 2.10 (a) Distribution of temperature coefficients obtained with PCR for vertical tracks compared with the theoretical distributions and (b) Slopes obtained via directlinearregression .............................. 43 2.11 (Top) Cosmic ray rate corrected by pressure as observed by the TRAGALDABAS detector and calculated via PCR method compared with the effective temperature and (Bottom) December 2016 SSW event temperature map 46 2.12 Forbush Decrease event observed on 22 June 2015 with the TRAGALDABAS detector ...................................... 48 3.1 Mountain profile along the railroad tunnel under the Mount Tobazo, where the Canfranc Underground Laboratory (LSC) can be found . . . . . . . . . . . . . 50 List of Figures 3.2 Angular distribution of high-energy muons measured at Canfranc Underground Laboratory(LAB2400).............................. 51 3.3 Estimated muon counting statistics as a function of the threshold energy at CanfrancExperiment............................... 52 3.4 Statistic and absolute temperature variations for Ethcosθ=71.13 GeV, an area of 50 m2and a measurement time of 24 h . . . . . . . . . . . . . . . . . . . . . . 53 3.5 Rate variations corresponding to the temperature effect versus data including both temperature and statistical variations. For Ethcosθ=71.13 GeV, A=50 m2 and ∆t=24h.................................... 55 3.6 ∆Te f f compared during a SSW event and a summer period for two different thresholdenergies................................. 55 3.7 Temperature coefficients for vertical direction (θ=0◦) at several threshold energies and including soft muons (E<0.4GeV)................ 56 3.8 Vertical muon intensity and threshold energy as a function of vertical slant depth 58 3.9 Slabschematics.................................. 59 3.10 RMSE of estimated temperatures as a function of the pressure level after cosmic-ray data inversion. Single-channel analysis based on vertical muons using (a) hard and soft component, with the combination of both; (b) using the underground component at different depths . . . . . . . . . . . . . . . . . . . 62 3.11 RMSE of estimated temperatures as a function of the pressure level after cosmic-ray data inversion. Single-channel based on vertical muons combining hard, soft, and underground components. Experimental data from Kohno et al. isincludedforcomparison ............................ 63 3.12 (a) RMSE of estimated temperatures as a function of the pressure level after cosmic-ray data inversion in the multi-channel case, using the hard component. (b) Distribution of temperature coefficients for hard muons at several angles . . 64 3.13 (a) RMSE of estimated temperatures as a function of the pressure level in the multi-channel case, using the soft component. (b) Distribution of temperature coefficients for soft muons, at several angles . . . . . . . . . . . . . . . . . . . 64 3.14 RMSE for the estimated temperatures as a function of the pressure level using underground muons in the multi-channel case. (b) Distribution of temperature coefficients for muons above 10 GeV, at several angles . . . . . . . . . . . . . 65 3.15 (a) RMSE for the estimated temperatures as a function of the pressure level, in the multi-channel analysis. (a) Hard + underground at different depths. (b) Hard + soft + underground at different depths . . . . . . . . . . . . . . . . . . 66 3.16 RMSE between estimated and real temperature for each atmospheric pressure level. Comparison between single-channel and multi-channel analysis along with the results obtained in Miyazaki’s work . . . . . . . . . . . . . . . . . . . 67 3.17 RMSE for the estimated temperatures as a function of atmospheric pressure level, for a three station/multi-channel analysis at an optimal depth. Each line represents a different value of the size of the detector . . . . . . . . . . . . . . 68 3.18 RMSE of the estimation of temperature for each atmospheric pressure level, for a three station/multi-channel analysis at an optimal depth where different levels of Gaussian noise have been aggregated to the rates . . . . . . . . . . . . . . . 69 3.19 Combined temperature coefficients calculated for the atmospheric pressure levels of 50, 150, 200, 500, 850, and 1000 hPa . . . . . . . . . . . . . . . . . . 70 xxiv IRMA RI´ ADIGOS S´ ANCHEZ 3.20 RMSE of the estimation of temperature as a function of atmospheric pressure level for a scenario with uncorrelated atmospheric temperatures . . . . . . . . . 71 3.21 (a) Observed temperature for 2019. (b) Estimated temperature for 2019 using the three station/multi-channel analysis at an optimal depth . . . . . . . . . . . 72 3.22 Difference between observed and estimated temperature for 2019 . . . . . . . . 73 4.1 Schematic flow chart for the air shower simulations . . . . . . . . . . . . . . . 77 4.2 sf ......................................... 81 4.3 Relative variations of muon rates at the surface (Eth =3.2 GeV) as a function of the surface pressure variations . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.4 Distribution of pion-production height as a function of atmospheric height . . . 83 4.5 Simulated maximum height of production as a function of surface pressure variations for pions and muons . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.6 Simulated maximum height of production for muons as a function of height variations in the pressure level of 300 hPa . . . . . . . . . . . . . . . . . . . . 85 4.7 Simulated temperature coefficients as a function of height . . . . . . . . . . . . 85 5.1 CR induced ionization rates in the atmosphere compared between the solar maximumandminimum ............................. 91 5.2 Correction coefficients (Wk,i) calculated with the exact summations together with the optimized computation between particles of size rkand ri....... 96 5.3 Percent change in CN3, CN10, CN40, CN80, and CCN for various atmospheric regions....................................... 97 5.4 (a) Zonal-mean nucleation rates in the X=0.8 simulation for the solar minimum. (b) Percentage change in the nucleation rate between solar maximum andsolar-minimum................................ 98 5.5 (a) Zonal-mean nucleation rates in the standard simulation for the solar minimum. (b) Percentage change in the nucleation rate between solar maximum andsolar-minimum................................ 99 5.6 (a) Zonal-mean nucleation rates in the X=0.7 simulation for the solar minimum. (b) Percentage change in the nucleation rate between solar maximum andsolar-minimum................................ 100 5.7 Percentage change between the solar-minimum and the solar-maximum case of zonal-mean CN3, CN10, CN40, and CN80 concentrations (X=0.8) . . . . . . 100 A.1 October-December 2016 time series of normalized polar geopotential height anomalies from 1000 hPa to 1 hPa . . . . . . . . . . . . . . . . . . . . . . . . 108 B.1 Temperature coefficients obtained from PCA regression applied to simulated dataofsubperiod1 ................................ 110 B.2 Temperature coefficients obtained from PCA regression applied to simulated data of subperiod 2 with different levels of noise . . . . . . . . . . . . . . . . . 111 B.3 Examples of linear regression fits for the predicted temperatures in the multi-channel analysis at an optimal depth . . . . . . . . . . . . . . . . . . . . 113 B.4 Correction coefficients Wk,icalculated for X=0.7 and X=0.9......... 114 B.5 Percentage change between the solar-minimum and the solar-maximum case of zonal-mean CN3, CN10, CN40, and CN80 concentrations (X=0.9) . . . . . . 115 xxv List of Tables 1.1 Principal particle interactions in a cosmic ray air shower . . . . . . . . . . . . 6 1.2 Values of the hadronic interaction and decay lengths for pions as well as decay energies ...................................... 21 2.1 Barometric coefficients for the different sub-periods analyzed. . . . . . . . . . 39 2.2 Values of the temperature coefficients αTand αMSS obtained in this work . . . 45 List of Publications The following list indicates the publications derived from this research that were used to create this thesis. I declare that I am the main author of all these publications, I am properly authorized to use these articles in this context, and they were not used in any other theses. I also declare that there is a non-doctoral co-author (D.G.C.) who participated in the first publication. • I. Ri´ adigos, D. Garc´ ıa-Castro, D. Gonz´ ales-D´ ıaz, and V. P´ erez-Mu˜ nuzuri. Atmospheric temperature effect in secondary cosmic rays observed with a 2 m2ground-based tRPC detector. Earth and Space Science, 7(9), p. e2020EA001131, 2020 Impact Factor (JCR): 2.900 (2020) Quartile (JCR): Q2 (Geosciences, Multidisciplinary) • I. Ri´ adigos, D. Gonz´ alez-D´ ıaz, and V. P´ erez-Mu˜ nuzuri. Revisting the limits of atmospheric temperature retrieval from cosmic-rays measurements. Earth and Space Science, 9, e2021EA001982, 2022. Impact Factor (JCR): 2.900 (2020) Quartile (JCR): Q2 (Geosciences, Multidisciplinary) • I. Ri´ adigos et al. The charge on aerosols from Cosmic Rays and the enhancement of cloud condensation nuclei formation. In preparation, 2022. My contribution to both published publications included the analysis, design, code implementation, development of the methodology, generation of plots and tables, and writing of the original drafts. Co-author D.G.C. was responsible for the correction and preparation of the data as well as the part of the barometric effect analysis reported in the first publication. Chapter 2 and 4 of this thesis are based on the results of the first and second articles, respectively. These articles are open access with permissions to be used in this thesis. The uses-permission for all the results taken, modified, or derived from these publications can be verified under the following link: https://agupubs.onlinelibrary.wiley.com/hub/ journal/23335084/open-access.html. These permissions are included in order to avoid any legal issues related to plagiarism or self-plagiarism. The next list indicates other publications derived from the thesis research, but they were not included as part of the manuscript: Proceedings • J. A. Garz´ on, J. Collazo, J. Cuenca-Garc´ ıa, D. Castro, J. Otero, M. Yermo, J. J. Blanco, T. Kurtukian, A. Morozova, M. A. Pais, A. Blanco, P. Fonte, L. Lopes, G. Kornakov, Chapter 1. INTRODUCTION The broad range of energies stems from the diverse types of sources that originate CR. For this reason, CR are typically categorized into several groups according to their origin. A small part of them (ranging from few tens of keV to several GeV) are generated in the Sun during periods of intense solar flares or caused by interplanetary phenomena (shock waves) associated with coronal mass ejections. The Galactic Cosmic Rays (GCR) represent the bulk of the spectrum. With energies as high as 1016 eV, they are created in supernova explosions, pulsars, double stars, and other objects in the galaxy. However, there are still many open questions about their origins and some aspects remain still unknown [51, 65]. The extragalactic CR are particles with exceptionally high energies (up to 1021 eV) and also the least common and enigmatic, they are referred to as ultrahigh-energy cosmic rays (UHECRs). Very little is known about their origins and what mechanism is accelerating them. The reason is the lack of statistics in the observation of this kind of particles. Above 1016 eV, the flux of CR drops to below one extragalactic particle per square meter per year. Large detectors would be required to collect a sufficient number of them. In any case, some of the candidates that are being considered to explain their origin are pulsars, neutron stars, gamma-ray bursts and jets from active galactic nuclei, to mention a few [99]. We know that some of the extreme phenomena reported above must be responsible for generating the most energetic particles ever measured here on Earth. In fact, the most powerful particle accelerators in the world have been able to accelerate protons to a record energy of 7 TeV, a far cry from the fastest CR [44]. In conclusion, cosmic radiation can provide us with new insights into the nature of the Universe. 1.1.1 Energy Spectrum Cosmic rays reaching the top of the atmosphere are called primary cosmic rays. As we will see later, primary cosmic rays interact with the nuclei of the elements that populate the atmosphere, giving rise to a flux of secondary particles, or secondary cosmic rays. The incoming flux of cosmic rays is not constant over time and depends on several factors. On the one hand, the flux of low-energy particles is modulated by the solar wind, i.e., by the solar activity [165]. Some of these variations are connected to the 11-year solar cycle and show a solid anticorrelation with CR flux. The solar wind is a constant flow of plasma (energetic charged particles) released from the Sun. The heliospheric magnetic field is embedded in it and fills the Solar System like some sort of protecting bubble, preventing lower-energy Galactic Cosmic Rays from entering. Therefore, as the strength of the solar wind weakens during periods of minimum solar activity, the flux of GCR increases. Apart from the modulation effect of the solar activity cycle, there are some sporadic moments when the flux of the low-energy region may also change. Spontaneous powerful solar flares or coronal mass ejections are some of the events that can cause changes in cosmic rays with energies between the range of MeV and GeV. Importantly, the variation of the Earth’s magnetic field plays an important role because it deflects the charged particle fraction whithin cosmic rays along their path towards the Earth’s surface. The ones with the weakest energy will either be reflected into space or will be trapped in the intricate magnetic field lines, prevented from reaching the ground. In spite of this, if the incident particle has enough energy, it will follow a nearly straight path towards the surface. Hence, the cosmic-ray flux has a dependence on latitude, longitude and zenith angle of observation, which is a consequence of the shape of the Earth’s magnetic field. The geomagnetic cutoff rigidity, Rc, is the quantity that defines the rigidity value above which the incoming 2 IRMA RI´ ADIGOS S´ ANCHEZ Figure 1.1: The energy spectrum of the primary cosmic radiation. cosmic-ray particle will have an allowed trajectory. But, it should be noted that this is not a fixed value either, since the magnetosphere also varies with time. Finally, the GCR intensity also suffers variations due to the change in the rate of supernovae explosions in the solar neighborhood. However, these changes occur on long time scales of millennia and it can be assumed to remain fairly constant in this context. All these combined effects that have been mentioned above contribute to the total flux of cosmic rays. The intensity of primary cosmic rays as a function of energy is given in Figure 1.1. This energy spectrum features two prominent transition regions where the slope changes. Its shape is so steep that the flux of particles above 100 GeV is much larger than that above 1011 GeV by sixteen orders of magnitude. One particle per square meter arrive in a year with energies of 107GeV. A small detector flown at the top of the atmosphere with an area of 1 m2 would have to wait one year to measure a particle above that energy. There is an increment of the slope in the interval between 106and 107GeV. This region is generally called the knee. It also exhibits a flattening at higher energies above 109GeV, known as the ankle. The origin of these structures is still unclear, but it is assumed to be related to the different mechanisms of generation of cosmic-ray populations. It is suggested that GCR below the knee are accelerated in the shock waves of supernova remnants (SNR) [86], while particles between the knee and the ankle come from different galactic sources, such as pulsars [28]. The most energetic cosmic rays (UHECRs) with an energy greater than 1 EeV are created outside the Galaxy but, as pointed before, with an uncertain origin [22]. The spectra of primary nucleons can be approximated by an inverse power law in energy for the range between several GeV and tens of EeV: 3 Chapter 1. INTRODUCTION IN∝E−(γ+1)particles m2s sr GeV (1.1) where Eis the energy and γ≈1.7 is the integral spectral index. It is also very common to define α=γ+1 as the differential spectral index. Above the knee the spectrum steepens with an index value of γ∼2. 1.2 Cosmic Rays in the Atmosphere We have previously mentioned that any particle pelting the Earth will encounter the magnetic field, which acts as a natural barrier to the weakest CR. The primary radiation that penetrates this shield is hazardous to organic life forms, and this is one of the reasons why astronauts cannot stay for long periods in outer space. Nevertheless, any particle permeating the magnetic field still has to face the atmosphere, our second natural shield against cosmic radiation. These primary particles will seldom reach the ground, rather they will interact with an atmospheric nuclei, usually in the upper atmosphere [71]. If the primary cosmic ray has enough energy, its collision will produce a large number of new particles and nuclear fragments, triggering a chain of nuclear interactions that are capable of generating up to billions of secondary particles. As implied by the description itself, these are generally called secondary cosmic rays. All these particles compose what is known as an Extensive Air Shower (EAS) and propagate randomly through the atmosphere until some of them eventually reach the ground. In some cases, they can cover a surface area of several hundred square kilometers. These particles can be measured using sophisticated detectors placed on the ground. The secondary particles lose energy due to interactions as they move downwards in the atmosphere. At the first stage of the cascade, the number of particles increases dramatically, reaching a maximum at a height of about ∼20 km (called Pfotzer maximum) [18]. However, the daughter particles will have less energy than their predecessors by conservation of energy. Eventually, they will not be capable of generating new particles and will be absorbed along the way. As a result, only a small fraction of the total amount of particles generated will arrive at the surface. In addition, the number of secondary cosmic rays that are produced in a shower will depend on the available energy of the primary cosmic ray. The more energy it has, the more particles can be created. Considering the energy spectrum of primary CR, those which have the capacity of yielding bigger air showers are the most energetic but also the scarcest ones. Low-energy CR (below a few GeV) will barely have enough energy to initiate a cascade. Regardless of this, it must be kept in mind that the Earth is constantly being hit by lots of CR, and even though only a small part of the secondary particles reach the ground, they are enough to produce a flux of ten particles per second passing through a surface the size of a hand. The secondary component of CR is a part of the natural radioactivity present in the environment in which we live. As can be seen in Figure 1.2, different kinds of interactions occur in a cascade that result in the creation of all types of subatomic particles: muons, pions, kaons, neutrinos, neutrons, electrons, positrons, or gamma rays. After the first collision of the primary CR (usually a proton) with an air nucleus, a lot of mesons are produced. In particle physics, mesons refer to a category of subatomic particles which are distinguished by being composed of a quark and an antiquark. Baryons are the other category that encompasses heavier particles formed by quark 4 IRMA RI´ ADIGOS S´ ANCHEZ Figure 1.2: Illustration of an atmospheric air shower initiated by a primary cosmic-ray particle interacting with an atmospheric nucleus. The cascade of secondary particles is divided into the muonic component (blue), the hadronic component (black), the electromagnetic component (red), and the neutrino component (green), undetectable for practical purposes. triplets, such as the well-known protons and neutrons. Baryons and mesons are both hadrons, the family that includes any particle composed of quarks. Returning to the formation of the air shower, we had mentioned before that in the first interaction a lot of mesons are produced, namely the so-called pions. In addition, kaons (another type of meson, having a “strange” quark inside it) and other baryons can also be created. Pions and kaons are not stable and they are susceptible to decay into other particles rather than interact. Charged pions (π±) decay into elementary particles called muons (µ±) and neutrinos (ν). However, mesons that live longer also have the probability to collide with another atmospheric nucleus before decaying and produce a bunch of new particles. The neutral pions (π0) have a very short lifetime as well and tend to decay rapidly into gamma rays. The latter may in turn create electron-positron pairs through interaction with the field of a nucleus. At the same time, electrons and positrons may produce more gamma rays through bremsstrahlung radiation in the field of a nucleus too. Muons from the hadronic component of the CR shower, which are produced mainly by pions and kaons, are far less interacting and can decay to electrons and positrons [75]. All these phenomena give raise to the so-called electromagnetic component of a CR shower. Table 1.1 summarizes the most relevant particle processes in an air shower. Therefore, it is clear that an air shower is composed of a great variety of particles resulting from the multiple interactions that take place as the cascade develops. For a complete picture, the flux of secondary cosmic rays in the atmosphere is typically divided into three components: • Hadronic component: including protons, neutrons, pions, kaons... • Electromagnetic component: consisting of gamma rays, electrons and positrons. • Muonic component: muons. 5 Chapter 1. INTRODUCTION Interactions Decays p+A→p+n+π0+π±+... π+→µ++νµ π−→µ−+¯ νµ π0→γ+γ π±+A→π0+π±+... µ+→e++νe+¯ νµ µ−→e−+¯ νe+νµ γ→e−+e+ Table 1.1: Most relevant particle processes in a cosmic ray air shower. One of the most relevant components is the one that refers to muons since they represent the largest number of charged particles that reach the Earth’s surface. In fact, the muonic component will cover most of the scope of this dissertation. Figure 1.3 shows some examples of atmospheric air showers obtained using CORSIKA (a Monte Carlo program employed to simulate air showers [57]). The cascades are created by protons of different incident energies: 10 TeV (Fig. 1.3a) and 100 GeV (Fig. 1.3b). It can be seen from the illustrations that the primary cosmic ray with the highest energy is capable of generating much more particles that also cover a wider area when they reach the surface. In addition, it has its first interaction at a height of ∼15 km, whereas the less energetic proton interacts at a much lower altitude (∼8 km). 1.2.1 Muons Most muons are produced in the upper atmosphere at an altitude of about ∼15 km. Compared to their relatives, they are much heavier particles: muons have a mass of 105.7 (a) 10 TeV (b) 100 GeV Figure 1.3: 3D development of air showers triggered in the atmosphere by different incident protons. The XY Z coordinates are given in km. The lines represent the trajectories of muons (blue), hadron particles (black), and gamma rays and electrons/positrons conforming the electromagnetic component (red). The showers were simulated using the Corsika software [57]. 6 IRMA RI´ ADIGOS S´ ANCHEZ MeV/c2, which is approximately 200 times greater than that of the electron (me=0.511 MeV/c2) [75]. Muons are very weakly interacting particles and lose energy by primarily emitting bremsstrahlung radiation as they travel through the atmosphere. This radiation, composed of photons, is produced when a charged particle is deflected by another charged particle, such as an atmospheric nucleus. The amount of energy emitted here is inversely related to the mass of the particle, which accounts for why muons can penetrate far deeper into matter than electrons. By way of illustration, the mean energy of muons created at the site of production is 6 GeV and they lose about 2 GeV before reaching the ground. Muons are relativistic particles moving at nearly the speed of light (∼0.999c) and also feature a very short mean lifetime, τ=2×10−6s. From a classical point of view (distance traveled equals speed multiplied by time), a muon produced at a height of 10 km would only travel 600 m before decaying, implying that muons would never reach the surface. In contrast, it is observed that they do. Then, what is misunderstood? The answer is that relativistc effects have not been taken into consideration. Time dilation, which is involved in those cases, requires the Lorentz factor (Γ≡1/p1−β2, where β=υ/cis the particle velocity relative to the speed of light) to be considered. It indicates how much the temporal characteristics of an object that is moving (in particular its lifetime) change for an independent observer, especially at very high speeds. This factor will have a value close to one for classical speeds but becomes higher than one for relativistic scenarios. For those special cases, which include muon movement through the atmosphere, time dilation makes its mean lifetime to be larger, 1.4×10−4s for the example discussed (muon with an energy of 6 GeV), and therefore the distance traveled would be 42 km instead (relativistic effect in length is l=Γβcτ). This allows them to reach the surface before decaying and even go deep underground. To conclude, muons rain down on every single square centimeter of the Earth’s surface and their average energy at the ground is ∼4 GeV [142]. Their penetrating power makes them a suitable tool for imaging dense and large materials without causing any damage to them. As muons travel through objects, they are absorbed in different amounts depending on the density of the material and the energy of the incident particles. Scientists can compare the flux measured after traversing the obstacle with that expected without it, to reconstruct the inner density (“muon imaging” or muography) [30]. This property has made it possible for muons to be used in different research areas to carry out the most astonishing discoveries. In 2017, a group of archeologists discovered a hidden chamber in Egypt’s Great Pyramid by scanning its interior with the help of cosmic muons passing through it [116]. Furthermore, this technique has also been employed in the area of volcanology to reveal the density profile of a volcano to foresee how an eruption could develop [153]. We will see later how this technique can also be applied to retrieve information from the atmosphere. 1.3 Cosmic Ray Measurements There are different methods for cosmic-ray detection that depend mainly on the component to be studied as well as the part of the spectrum (i.e., range of energies) to be covered. On the one hand, direct detection of primary cosmic rays is possible thanks to particle detectors placed on orbiting spacecrafts or the International Space Station. Along with it, balloon-borne instruments reaching high altitudes are launched for the same purpose. However, this kind of detection only allows measurements of low-energy primary CR. 7 Chapter 1. INTRODUCTION Balloons can also be deployed to register the secondary radiation as a function of altitude as they rise in the atmosphere. Aircraft can be considered for these cases as well. On the other hand, there are several ground-based techniques to measure secondary CR. Firstly, neutron monitors are widely used to monitor variations of neutron rates in the energy range between 500 MeV and 20 GeV [145]. A standard neutron monitor consists of a set of counter tubes made of several layers of materials. The outer layer protects the measurements from external sources of noise and allows neutrons from the cascade to pass through. The incident neutron generates additional neutrons via nuclear interactions in the following layer made of lead. Eventually, when they reach the innermost layer, they can be captured by the nuclei of the gas that fills it, emitting other particles that are easier to detect and transform into electrical pulses. This kind of detector is typically used in astrophysics to monitor the Sun’s activity [115]. In this regard, it is a kind of indirect detection of primary CR. Neutron detectors are also used in the field of hydrology to measure the soil water content over wide areas or reveal mountain snowpacks [45, 95, 90, 105]. Another type of detector most frequently used is that which measures muons. There are plenty of different techniques to detect them and the equipment employed is usually very diverse in design and performance. The most basic detectors may consist of ionization chambers or scintillation counters. In the first case, the particles that enter the detector ionize the gas inside it, creating charges that are collected using an electric field. In the second case, a scintillating material that emits photons in response to ionizing radiation is used. The released photons are then converted into electrons via the photoelectric effect, then they are accelerated to strike a series of dynodes that yield more electrons, and so amplifying the initial signal. The resultant output is a pulse proportional to the energy of the traversing particle [162]. The detectors mentioned above are very versatile but they only measure the integrated flux of CR coming from all directions. Anyhow, better performance can be achieved by installing several of these devices in specific layouts. Muon telescopes assemblies involve several detectors (scintillators for instance) positioned along a straight line from smallest to largest thickness. Identification of the type of particle (i.e., its mass and charge) is possible by measuring the energy deposited in each detector. Such detectors are usually rigid structures that can be rotated in the zenith and azimuth directions. Moreover, depending on the detector layout, they require the particle to have a minimum energy (i.e. threshold energy) to be able to pass through the entire detector assembly and be recorded [52]. Multi-directional telescopes are a more attractive alternative that include angular resolution and allow measurements of CR from different directions (e.g., [14]). Generally, the system configuration consists of two arrays of detectors, one on top of the other. Figure 1.4a shows an example of a multi-directional telescope. This type of layout enables the identification of particle trajectories (tracking) by means of signal coincidences. When a particle hits two detectors in each of the layers within a coincidence window of a few nanoseconds, a signal is assigned to the particle that has crossed them. In CR experiments, the particles involved are relativistic, and hence the required time window to establish a coincidence between a couple of detectors spaced by tens of centimeters will be a few nanoseconds. As the window increases, so does the probability of random coincidences, as well as the possibility of several interesting events occurring inside the same window and being missed. To sum up, when a coincidence is recorded, the trajectory of the particle can be determined. The number of detection directions is constrained by the geometry of the setup and the angular resolution relies both on the size of the individual detectors that compose the arrays and the 8 IRMA RI´ ADIGOS S´ ANCHEZ Figure 1.4: Scheme comparing a multi-directional muon telescope with a muon hodoscope. (a) A simple multi-directional telescope assembling formed by two horizontal layers, each one with 4 individual scintillator detectors. A lead layer is located below the upper layer to absorb low-energy brackground radiation. The estimation of the direction of the incident particle is quite limited by the detector configuration. (b) Hodoscope consisting of 4 layers, each one with an array of 6×6 individual detectors. The particle track can be estimated with high accuracy from the number of individual detectors that are hit when the particle crosses the 4 planes. distance between layers. In this context, detection directions are fixed and the threshold energy of such telescopes is a function of the zenith angle since the particles’ path increases with it. On the other hand, traditional telescopes do not have counting rates sufficiently high and tend to have large sizes in order to retain measurement statistics. Muon hodoscopes took the next step in the development of CR detectors by improving the accuracy and resolution of tracking particles with new approaches. One of the main features of hodoscopes is that they can track charged particles from virtually any direction of the upper hemisphere. This provides a continuous measurement of the angular distribution of the particle flux. The arrangement of a muon hodoscope consists of arrays of many segments (i.e. detectors) located at two or more parallel planes. A track can be inferred from the number of segments that light up a signal (trigger) when they are hit by a particle as it passes through the planes. Obviously, the spatial resolution of these detectors is limited by the segment size. Figure 1.4b illustrates the layout of a muon hodoscope. An example of an hodoscope is URAGAN, which operates at the National Research Nuclear University in Moscow. This detector is the first large area muon hodoscope in the world and is made of four independent modules. Each of them is an assembly of eight planes of small discharge tubes equipped with a two-coordinate system of external readout plates (strips). Every layer contains 320 tubes and the total area covered is 3.5×3.5 m2. The system detection requires the coincidence of signals from at least four of the strips of the detection planes within a time window of 250 ns. In addition, the range of threshold energies goes from 0.2 to 0.6 GeV. With all these features, URAGAN allows high accuracy measurements of the surface muon distribution [16]. In general, the detectors mentioned so far are suitable to study secondary cosmic-ray fluxes at the surface and at low energies. We have previously remarked in Section 1.1.1 that high-energy CR are far rarer and their arrival rate per square kilometer is very low. For this reason, a detector of a gigantic area would be needed to measure at least one of those energetic 9 Chapter 1. INTRODUCTION particles. The assembly of such a device would be nonsensical. Fortunately, there are other ingenious ways to detect these cosmic rays. The study of CR in the high-energy region can only be made in practice by observing the air showers that they produce. As the secondary particles pass through the atmosphere, their interaction with the atmospheric molecules produces different kinds of radiation. On the one hand, they excite the gas molecules, mostly nitrogen, resulting in the emission of visible and ultraviolet radiation. The fluorescence light is produced isotropically and can travel several kilometers through the atmosphere to be detected by an optical telescope (fluorescence detectors) [1]. On the other hand, electrons and positrons travel faster than the speed of light in air and hence emit Cherenkov radiation that is measured by the so-called Cherenkov telescopes. This kind of detector comprises a large segmented mirror that concentrates the Cherenkov radiation it receives towards an array of photomultiplier tubes. Apart from this, electrons and positrons can also emit electromagnetic radiation with frequencies of tens of MHz, one of the reasons being their interaction with the Earth’s magnetic field (synchrotron radiation). These radio signals are pointed sharply downwards relative to the shower development and can be recorded locating antennas at the ground level [89]. Finally, another method to detect the secondary products of highly energetic CR is to deploy big arrays of surface detectors over a wide area (∼100 km2). Such arrays can observe a vast part of the celestial hemisphere. And what’s more, they are often built as hybrid observatories, incorporating other types of complementary detectors (e.g. Cherenkov telescopes and fluorescence detectors) to carry out a more complete study of the cascades. The surface array samples the distribution of charged particles at the ground, whereas the other kind of telescope supplies simultaneous measurements of the longitudinal development and lateral distribution of particles (i.e., particle density over a plane normal to the shower axis). This combination of measurements gives very valuable information about the primary cosmic ray, making it possible to reconstruct its energy, mass composition, and direction of arrival, for example. One of the biggest array experiments for the detection of UHECR is the Pierre Auger Observatory in Argentina. It is located in the vast plain of Pampa Amarilla in Mendoza Province and its array of surface detectors currently consists of more than 1600 water tanks distributed over an area of 3000 km2(about 30 times the size of Paris) accompanied by 27 fluorescence detectors. It was designed for a high statistics study of UHECR and retrieve both their energy and arrival direction. It has been taking data since 2004 [104]. Another interesting example of a big array is the IceCube Neutrino Observatory constructed in Antarctica. The experiment counts with thousands of sensors placed deep in the ice and distributed over a cubic kilometer. The set of detectors can be found at depths between 1450 and 2450 m, and they are based on photomultiplier technology. The observatory also includes the IceTop, a surface array with ∼162 tanks of ice to measure showers of secondary particles. Its construction began in 2005 and, since then, it has been continuously incorporating new improvements [64]. The Antarctic observatory is committed to the search for neutrinos, which are nearly massless particles quite challenging to detect. The most energetic population of these particles originates from the most violent astrophysical sources: gamma-ray bursts, black holes, and neutron stars. Therefore, their study provides information for surveying these fascinating astrophysical phenomena. 10 IRMA RI´ ADIGOS S´ ANCHEZ Figure 1.5: Simplified diagram of cosmic ray showers detection techniques. Muons or neutrons are measured with ground-based particle detectors; the electromagnetic component is measured with another kind of detectors by means of Cherenkov and fluorescence light, or antennas for the radio pulses; deep underground detectors are devoted to high-energy muons measurements; big arrays cover wide areas to study extensive air showers from ultrahigh-energy cosmic rays. The challenge of neutrino detection is that they seldom interact with matter. Nevertheless, when they do interact with the water molecules in the ice, several particles such as muons or electrons are created. These charged particles leave a characteristic signal as they pass through the ice that can be recorded by the sensors. In the same fashion, CR can be studied underground. Hadrons, electrons, and gamma-rays are immediately absorbed by the rock when they reach the ground level given that it is denser than atmospheric air. Contrarily, high-energy muons can penetrate deep underground and be detected at great depths. Considering the energy loss processes of the muon passage through matter (ionization of the medium, bremsstrahlung, etc.), the minimum energy required for a muon at the surface to reach a certain depth Xcan be estimated with the following equation [65]: Eth =εeX/ξ−1(1.2) where εdefines a critical energy that equals 500 GeV for muons in rock and ξ≈2.5·105g/cm2. Therefore, if we are interested in studying muons with energies above a specific threshold, the best way to do it is to place an underground detector at the corresponding depth. An alternative approach would be to incorporate a lead shield of a predetermined thickness in a surface detector that only allows higher energy muons to pass through. Clearly, this technique is only feasible for low threshold energies, since higher energies require a greater thickness of lead. Some of the most relevant underground experiments are located at depths greater than 250 m. MINOS is an experiment at the Soudan Underground Mine State Park in Minnesota 11 Chapter 1. INTRODUCTION underground. Nonetheless, when tested, no significant variation was found. This reinforced the notion that the most important production processes of high-energy muons take place only in the upper parts of the atmosphere. All things considered, we have seen so far that for high-energy muons there is a positive correlation due to the decrease in air density when temperature increases (positive effect), giving as a result that more pions will decay into muons; for low-energy muons on the other hand, a negative correlation is found due to a similar dependency of the muon decay on air density changes (negative effect). Method of weighted temperature As discussed earlier, the real temperature variations are not uniform throughout the atmosphere, and the production of muons or pions cannot be approximated to take place at a single level. As a solution, a weighted or “effective” temperature Te f f can be calculated, which is the equivalent to the temperature that an isothermal atmosphere would have in order to produce the same modulation of the muon intensity as an atmosphere with the actual temperature distribution T(X). An example of an ad hoc expression proposed for the effective temperature was [17] Te f f =T20 +T40 +T80 +T125 +T250 +1 2T500/5.5 (1.6) Here, more weight is being given to upper atmospheric levels (20, 40, 80, 125 and 250 hPa) [17]. The formulation could be improved later by employing sophisticated models for nuclei and meson generation and propagation in the atmosphere. Anyway, this definition of effective temperature is very useful to study the effect of temperature variations in underground detectors. In such case, while the temperature of the troposphere undergoes considerable daily variations, the temperature of the stratosphere remains practically constant (except for occasional abrupt variations). On a seasonal scale, the slow variations of the temperature of the stratosphere and the decrease of the air density will reduce the probability of mesons to interact and a larger fraction of them will decay into muons. As underground rates are largely oblivious to the troposphere conditions, an underground detector will be sensitive to the small seasonal variations in the temperature of the upper atmosphere. Temperature coefficients Over the years, numerous attempts have been made to correlate CR intensity variations with atmospheric variations to obtain the partial temperature coefficients, in lay terms, the distribution of temperature coefficients as a function of height that relates the temperature variation in each layer to the corresponding variation of the total measured rates. However, in practice, it was very difficult to obtain the correct values. In 1986, the temperature effect had already been estimated for different seasons. Admittedly, far from completely solving the problem of the temperature effect, some underground experiments detected several anomalies in their values. Particularly, they had observed a semi-annual modulation in the muonic intensity measured at the Matshushiro station. This semi-annual variation displayed a very striking contrast compared to the annual 18 IRMA RI´ ADIGOS S´ ANCHEZ variations commonly observed at other underground stations closer to the surface (e.g., Misato or Sakashita) [135]. The distribution of temperature coefficients had previously been calculated theoretically by several researchers (e.g., [49]) and satisfactory results had been obtained under certain measurement conditions. Their calculations, however, were based on simplified models of the propagation of nuclei and pions and using a muon production spectrum deduced from surface experiments. These considerations proved to be insufficient. In 1986, Sagisaka recalculated the coefficients both for the barometric and temperature effect, solving the anomalies observed in other stations [135]. In his review, Sagisaka expressed the partial temperature coefficients as the sum of two terms: one representing the negative effect and the other the positive effect. In this way, he attributed the reason for the semi-annual variation found in Matsushiro to the dependence on the muon threshold energy of the atmospheric temperature effect. Furthermore, it was remarked once again that the temperature deviations varied with altitude in a complex way at each season. Proof of this is that the variation of the tropopause’s temperature is nearly opposite in sign to that in the troposphere (antiphase). Considering the partial temperature coefficients introduced by him as αi(h), where i=1,2 corresponds to low and high threshold energies, respectively; and the seasonal variations of the atmospheric temperatures as (∆T(h))j for j=1,2,3,4 being the seasons; the different products αi·(∆T)jcan be estimated. By integrating each of these products from the top of the atmosphere to the surface, the total temperature effect can be obtained for both underground and shallow detectors. Figure 1.7 shows an schematic view of the products as explained in Sagisaka’s work [135]. In the low-energy case (Eth ∼1 GeV), the distribution of coefficients αihas negative values and is nearly independent of the atmospheric height (Fig. 1.7 top). Thus, the contribution of the product α1·(∆T)jin the troposphere to the integral (∆I)jdominates compared to the contribution of the stratosphere (the troposphere has a higher air mass percentage). As a consequence, the temperature effect is opposite to the surface temperature, giving rise to an annual variation with its peak in winter. This can be appreciated in Figure 1.7 where the magnitude of the effect is positive for α1·(∆T)1and α1·(∆T)2, which correspond to the winter and spring months. On the other side, in the case of deep underground detectors, the contribution of the product α2·(∆T)jin both the stratosphere and tropopause sets the seasonal trend. The tropopause temperature peaks towards the end of the spring, much earlier than the troposphere, thus the maximum muon rate occurs during the warmest months. The example shown has been done with mid-latitude temperatures. In the Matsushiro’s location, the temperature variations in the upper atmosphere have a different behaviour than the ones showed in Figure 1.7. In their case, the contribution of the product α2·(∆T)jin both the troposphere and stratosphere to the integral are almost compensated by each other throughout the seasons. The minor differences in the balance between the two contributions yield a semi-annual variation with two maxima, one in summer and the other in winter. This revealed that the aforementioned balance was very sensitive to the height of the tropopause (∼200 hPa), the temperature profile of the stratosphere, and the shape of the distribution of the temperature coefficients as well. Given these points, it should be noted that there exists a competition between the processes of interaction and decay that will determine the evolution of the air showers and thus the muon rates at the ground. This competition mainly concerns mesons because protons, gammas, and electrons do not decay. On the one hand, the decay process depends on the mean lifetime τof the particle. On the other hand, the interaction depends on the amount of traversed matter, that 19 Chapter 1. INTRODUCTION Figure 1.7: Schematic representation of the atmospheric temperature effect on the muon intensity observed at surface and underground as explained in [135]. (Top) αi(i=1,2) are the partial temperature coefficients for two energy thresholds Eth =1 GeV and Eth =100 GeV corresponding to surface and underground locations, respectively. (Left) (∆T)j(j=1,2,3,4) are the seasonal variations of the atmospheric temperature from the yearly average at a mid-latitude location (40◦). The seasons are grouped by: December-January-February (DJF), March-April-May (MAM), June-July-August (JJA), and September-October-November (SON). (Right) The products αi·(∆T)jare plotted for each combination of iand jtogether with the magnitude of the temperature effect, which is shown with a colored circle: the red circles with the plus sign represent the positive and the blue ones with the minus sign the negative. 20 IRMA RI´ ADIGOS S´ ANCHEZ is, it is a function of the density of the medium. By way of illustration, Table 1.2 provides typical values for pion interaction and decay lengths. As can be seen from the table, if the decay length is bigger than the hadronic interaction length, ddec ≫λ, mesons can live long enough to interact and produce other mesons. In contrast, when ddec ≪λ, mesons decay before interacting and generate: gamma-rays in the case of neutral pions; muons and neutrinos in the case of charged pions (contributing to the muonic component of the cascade). Meson Decay Channel Decay Length d=βΓcτ≃Γcτ[cm] Interaction Length at 1012 eV [g/cm2] Edec [eV] π±π±→µ±νµdπ±≃780 ·Γλπ±≃120 7 ·1018 π0π0→γ γ dπ0≃2.5×10−6·Γλπ0≃120 2 ·1010 Table 1.2: Values of the interaction and decay lengths for pions as well as decay energies. Γis the Lorentz factor: Γ=E/(mπc2). The competition is going to depend on the energy of the particle Eand the density ρof the upper atmosphere (above ∼15 km). Hence, the energy Edec at which the decay and interaction compete needs to be calculated. Comparing the interaction length dint with the decay lenght ddec =λ/ρ, we have: ddec =dint ⇒Γcτ=λ ρ⇒Edec mc2cτ=λ ρ(1.7) where the decay energy is given by Edec =λ cτρ mc2(1.8) Now, if E≫Edec, then the particle will live long enough to interact. But if E≪Edec, the particle will decay. Table 1.2 provides the results obtained for a pion with E=1012 eV. In this example, for an air shower initiated by a primary cosmic ray with energy E0, neutral pions will not decay unless E>7·1018 eV, whereas charged pions having E<Edec ∼2·1010 eV will tend to decay. Building on the works of Dorman and Sagisaka, refined calculations have been seen recently, in particular those of Dmitrieva et al. in 2011, leading to a modern formulation of the above problem [47]. According to these studies, these processes can be embedded in a single function WT(Eth,X,h,θ), or DTC (Differential Temperature Coefficient), which provides information on how much an atmospheric layer at an altitude hcontributes to the variations in the flux of muons arriving at a zenith angle θ, with an energy greater than Eth, and at an observation level X. According to Dmitrieva et al., if the atmospheric temperature changes as ∆T(h), the standard muon intensity N0(Eth,X,θ)at a certain observation level Xwill be changed by an amount ∆NT(Eth,X,θ). Therefore, the relative variation of the muon intensity can be written as: ∆NT(Eth,X,θ) N0(Eth,X,θ)=ZX 0 WT(Eth,X,h,θ)∆T(h)dh ≈∑ i WT(Eth,X,hi,θ)∆T(hi)∆hi(1.9) 21 Chapter 1. INTRODUCTION Figure 1.8: Differential temperature coefficients WTfor vertical direction (θ=0◦) at several threshold energies [47]. Figure 1.8 provides some examples of temperature coefficients as a function of atmospheric height for θ=0◦and several values of threshold energies. The shape of the distributions illustrates how muons with high threshold energies are very sensitive to the stratospheric temperatures (h<200 hPa). In simple terms, small changes in the temperature of the stratosphere significantly affect the intensity of high-energy muons reaching the surface. A simplified schematic view of the temperature effect is shown in Figure 1.9 as well. Temperature variations in the atmosphere When estimating the correlation between CR rates and the temperature at a certain atmospheric level (∆N N=αp ∆Tp Tp), the slopes αpof this regression cannot be compared between detectors with similar arrangements but located in different places. The reason is that the atmospheric configuration differs from one place of the planet to another and the variations in ∆Tpcorrelate with the rest of the atmosphere in different ways, impacting the CR variations. For instance, a colder and denser air forms a more compact atmosphere at the poles. Furthermore, these regions are outside of the General Circulation Zone of the atmosphere, so seasonal variations are smoother, with a consistent low-pressure area that weakens towards summer. The polar vortices isolate Antarctica and North Pole from the atmospheric effects of the mid-latitudes (alternation of cyclones-anticylones, etc). Therefore, the atmosphere is preserved in a stable situation which may only be interrupted during some exceptional events called Sudden Stratospheric Warmings (SSWs) (see more information in Appendix A.1) [31]. Also, the tropopause is located at a lower altitude there, around 9 km, while at the equator (with a warmer atmosphere) it is above 17 km. Translated into pressure levels, 225 hPa is the global standard reference level for the tropopause. However, this can be misleading since the tropopause can exist at any position between 100 and 400 hPa over the year [107]. Not only the position of the layers must be taken into account but also the correlations 22 IRMA RI´ ADIGOS S´ ANCHEZ Figure 1.9: Illustration of the temperature effect for an air shower with initial energy E0during winter (left) and (summer). Initial particle and charged mesons are drawn with black lines; red lines are high-energy muons; green lines are low-energy muons; and dashed lines represent electrons and neutrinos. The positive effect is represented by more charged mesons decaying in a less dense atmosphere (right) and the negative effect is seen when more low-energy muons generated lower in the atmosphere can reach the ground level before decaying (left). of temperature variations between different levels. Due to this, comparisons should not be made between the αpof different detectors for the same pressure levels, since they depend on the regional characteristics of the atmosphere. In general, it is convenient to use an effective coefficient αTfor this purpose. The use of such coefficient provides a net value of the total atmospheric effect, thus removing the constraints not only of the differences in pressure levels but also of the diverse correlations between atmospheric layers that occur for different geographical locations. Therefore, the CR variations due to temperature changes can be expressed in terms of the effective temperature coefficient as follows: ∆NT N0 =αT∆Te f f (1.10) and αTand the effective temperature Te f f are defined as αT= n ∑ i=1 WT(Eth,X,hi,θ)∆hi(1.11) Te f f =∑n i=1WT(Eth,X,hi,θ)∆Ti∆hi ∑n i=1WT(Eth,X,hi,θ)∆hi (1.12) This will be addressed in more detail in Chapter 2. 1.5 Influence of Cosmic Rays on the Atmosphere In the previous section, we considered the problem of the influence of the atmosphere on CR rates. However, one can ponder the opposite question of whether or not CR influence the 23 Chapter 1. INTRODUCTION atmosphere. The interaction of CR particles with the atmosphere can lead to a host of interesting effects, several of which can have a big impact on the planetary environment. It is a well-known fact that CR have an influence on atmospheric electric field and, as a consequence, on thunderstorms as well [146]. Moreover, given that CR ionize the atmosphere as they pass through, their connection with the ionosphere is direct. The ionosphere is a layer of the upper atmosphere (from 80 to 1000 km) ionized by high-energy radiation from the Sun and CR. By means of illustration, during the night, without the influence of the Sun, only CR produce the ionization and hence the ionosphere is much less charged at nighttime. Consequently, any variation in the CR flux will undoubtedly affect this part of the atmosphere [169]. CR also influence atmospheric chemistry when they interact with gas particles, triggering a series of physicochemical chain reactions that alter the composition of the atmosphere. For example, the electrons arising from CR ionization can break N2molecules, causing the formation of NOx. This NOxmolecules are potentially dangerous because they can deplete or generate ozone, depending on the conditions. Besides, the change in NOxconcentrations can affect numerous chemical processes by competing with their reactants. Altogether, CR can modify atmospheric chemistry, create an ionosphere, influence atmospheric lightning, produce organic molecules in the atmosphere, destroy stratospheric ozone, and so on. For this reason, when studying the atmospheres of other planets (exoplanets), models that take into account the influence of CR are proposed to investigate their potential habitability [74]. Despite all this, we are more interested in the influence of planetary cloud covering and its possible effect on climate, in the short and long term. This topic will be covered in the last chapter of this thesis. On this occasion, we need to understand the background of this influence and its basic concepts. 1.5.1 The Variability of Solar Activity The main source of energy for the Earth’s surface comes from the Sun. The amount of energy received by a particular location varies over time, especially over the seasons. The reason is the Earth’s motion around the Sun and the tilt of its axis. However, the irradiance, which is the flow of energy radiated by the Sun, could be considered constant. Eventually, it was found that the Sun has cycles of activity and its irradiance can vary too [172]. Solar activity is measured by the number of sunspots on its surface. Namely, sunspots are dark regions that temporarily appear on the Sun and that are colder than their surroundings areas. These regions can be very large, reaching planetary sizes, and are caused by the interaction of the magnetic fields generated by the motion of the Sun’s gases in the outermost layers [144]. These processes generate a lot of activity on its surface, which is known as solar activity. Therefore, sunspots can be considered as an useful indicator of the activity of the star. In the same way as weather on Earth varies from season to season, the Sun’s activity also has its phases, the solar cycles. Granted that solar activity can have an impact here on Earth, scientists need to monitor its status on a regular basis. Sunspots have a lifecycle of the order of weeks and follow the rotational motion of the Sun itself, which lasts 27 days on average. The number of spots in the Sun’s surface has been counted since the beginning of the seventeenth century, by the time the telescope was invented. Thanks to these measurements, it was possible to see that the number of sunspots followed a cycle of about 11 years. The beginning of a 24 IRMA RI´ ADIGOS S´ ANCHEZ cycle is a solar minimum when the Sun presents the least sunspots. Over time, the number of sunspots increases until it reaches a peak, the solar maximum. The cycle ends with the return to the minimum [166]. The irradiation of the Sun changes according to this cycle, although it should be emphasized that this variation is small, in the order of 0.1 %. The effect of this change of the long-term solar irradiance is generally considered to be too small to have any effect on the current climate change [97]. Indeed, they are meant to be negligible only for the last 150 years, approximately, due to the greater magnitude of anthropogenic climate change [186]. On the other hand, large eruptions occur on the Sun, such as solar flares and coronal mass ejections (CMEs), that increase as the Sun approaches the solar maximum. Solar flares happen because the magnetic field lines of the sunspots often twist, cross, and reorganize, causing explosions of energy. These flares release a lot of electromagnetic radiation into space that, if directed towards the Earth, can interfere with radio communications [171]. Solar flares are often followed by a coronal mass ejection. In such a case, a huge amount of plasma (gas of charged particles) from the solar corona with its corresponding strong magnetic field is released. These Sun’s disturbances generate solar storms that travel across the interplanetary space affecting the planets found in their path. Especially, when these charged particles hit the Earth, they are deflected by the Earths’s magnetic field towards the poles, where they interact with the upper layers of the atmosphere producing the auroras. When CMEs are particularly strong, the resulting geomagnetic storms can produce strong induced electric fields which can affect the electrical transmission lines, causing massive power outages. In addition, solar storms can damage satellite electronics and affect terrestrial communications. The whole series of phenomena related to solar activity form what is known as Space Weather. Researchers work hard to improve space weather monitoring and forecasting in order to protect our communications, keep astronauts safe, etc. 1.5.2 Cosmic Rays and the Solar Cycle The intensity and energy spectrum of GCRs is modulated by solar activity. The Sun constantly emits a stream of charged particles called the solar wind. As the solar wind spreads out filling the interplanetary space, it creates an environment of radiation and magnetic fields, forming a giant bubble around the star and planets, known as the heliosphere. This acts as a shield, protecting the planets from galactic cosmic radiation. Additionally, Earth is protected by its own magnetic field, which has several benefits. In Section 1.2, it was mentioned that the magnetosphere is shielding Earth from cosmic radiation. Similarly, it also prevents the atmosphere from being degraded over time by the collision of the solar wind with it. Indeed, planets without magnetic fields, such as Mars and Venus, are exposed to high radiation levels and their atmosphere is gradually vanishing. In the case of Mars, this process is at a very advanced stage, with an atmosphere so rarefied that yields high surface radiation which should be taken into consideration for future space missions [73]. The expansion of solar ejections into interplanetary space produces plasma overdensities at the solar maxima. Therefore, the heliosphere becomes more effective in scattering high-energy GCR entering the solar system. Thus, the flux of GCR arriving at the Earth is reduced, leading to its anticorrelation with the sunspot number [80]. Together with these long-term variations in the flux of GCR, sudden short-term variations can also happen. Forbush Decreases (FDs) are one of the main phenomena responsible for this. 25 Chapter 1. INTRODUCTION A Forbush Decrease is produced by a CME and refers to a sudden decrease in GCR intensity. They are caused by the heliospheric magnetic shock driven ahead of a CME that sweeps away some of the incoming GCR when it reaches the Earth. FDs also vary in amplitude, ranging from 3 to 20 %, and usually last several days. They can be most easily observed at the Earth’s surface with neutron monitors. A FD is characterized by having several phases. The decrease phase starts with the arrival of the shock and may be preceded by a small increase of 2 or 3 %. This preincrease has a very short duration of few hours. The intensity drop is very abrupt and happens in approximately 24 hours. After this stage, a gradual recovery phase occurs in the following days [129]. 1.5.3 Link between cosmic ionization and cloud-covering Figure 1.10: Illustration of the possible link between GCR and cloud-covering: a lower solar activity produces a weaker solar magnetic field that deflects much less GCR and thus more clouds will form due to atmospheric ionization (left). The opposite would yield clearer skies (right). It has been previously pointed out that cosmic ionization of the atmosphere can alter its physical and chemical processes. Because of this, it has been suggested on numerous occasions that solar variability could be correlated with cloud-covering and, therefore, be a contributing factor for climate change by modifying the Earth’s average albedo. However, the effect of the Sun’s activity on the climate is a controversial question, and there is no clear consensus as there are studies that contradict each other. Figure 1.10 shows a sketch of the possible relation between solar activity and cloud cover. At the end of the last century, several papers were published in which they observed cloudiness variations at mid-latitudes associated with Forbush Decreases [127, 128]. In 1997, the first correlation between global cloud cover and CR intensity was reported [151]. The cloud cover seemed to be inversely correlated with the solar activity. Five years of satellite data had been analyzed, in coincidence with a solar minimum, and found variations of 3-4 % in the cloud cover. Apart from this, it was also apparent that the observed variation was larger at higher latitudes than in the Tropics, in agreement with the lower values of the rigidity of the Earth’s magnetic field. Against the significance of this correlation, other scientists argued that it was very difficult to attribute the correlation to cosmic rays because the actual microphysical explanation for such effect was still lacking. Subsequent reassessments of the trends showed divergent results: strong correlations in Forbush Decrease events [157, 78] against weaker or 26 IRMA RI´ ADIGOS S´ ANCHEZ no correlations [147, 100, 158]. Furthermore, some analysis also showed the impact of the methodology used on the differences in results (e.g., [102]). In 1993, Tinsley et al. had already suggested that a possible mechanism to explain the linkage between solar activity and climate could be the atmospheric electricity variations caused by the solar wind and that involve the charging of supercooled water droplets and aerosols (small atmospheric particles) at clouds [156]. Later, Svensmark et al. presented a work where they found correlations not only with cloud cover but also with cloud and aerosol properties during several FDs [148]. These results might be evidence that cloud changes could be driven by changes in aerosols. Yet, the studies did not clarify whether the effect could be the other way around, clouds could be affecting aerosols. Again, a subsequent similar survey found no evidence of significant correlations [32]. Clouds have a strong influence on the Earth’s energy budget, absorbing and reflecting radiation from the Sun and the Earth’s surface. Therefore, small changes can have a big relevance for the climate. However, detailed studies on cloud formation are difficult to carry out because it is not straightforward to replicate cloud formation in laboratories. 1.5.4 The CLOUD Experiment The CLOUD (Cosmics Leaving Outdoors Droplets) project at CERN is a huge experiment that was created to elucidate the possible links between GCR and cloud formation. The experiment focuses on the study of the formation and growth of aerosols that lead to the condensation of cloud droplets. The experiment is performed in a huge chamber (26 m3) equipped with a wide range of sensors to track the evolution of clouds inside. The CERN Proton Synchrotron provides the artificial and adjustable source of “cosmic rays” through a pion beam. The entire setup allows tuning the chamber’s conditions to duplicate those of the real atmosphere, such as temperature, humidity, or levels of ionization. In its more than 10 years of activity, CLOUD has made many striking discoveries. In particular, it was the first experiment to demonstrate a clear relationship between GCR ionization and aerosols formation. In truth, they have found a relatively weak dependence on ion concentrations and that, in the present-day atmosphere, CR intensity cannot meaningfully affect climate via aerosol growth [54]. In spite of this, it has been discovered that ions from GCR can strongly promote the formation rate of biogenic vapors emitted by trees up to a factor of 100. These findings reveal that CR may have played an important role in cloud formation in pre-industrial times [94]. Besides, they demonstrated that nucleation rates could be enhanced thanks to ions under certain conditions. Some of the findings reported are: the enhancement of neutral nucleation caused by cosmic ionization can reach a factor of 15 at the temperatures of the low troposphere; and, in general, ion-induced nucleation is the dominating process in most parts of the troposphere. Nevertheless, it should be taken into account that proving that nucleation of particles is affected by GCR does not imply that this will affect cloud formation. The work of the CLOUD experiment has shed light on our understandings of aerosol nucleation, growth, and their link with clouds and climate. However, they had not yet addressed the connection between GCR and clouds. Alternatively, they have integrated their experimental results in global aerosol models to test the roles of the different processes of atmospheric particle growth and formation [54, 70, 69]. In these investigations, they estimated that ions from GCR are accounting for about half of the nucleation in both the present-day and pre-industrial atmosphere. 27 Chapter 2. Atmospheric Temperature with a High-Resolution Cosmic Ray Detector Figure 2.2: Photo of the TRAGALDABAS detector at the Faculty of Physics of the Univ. of Santiago de Compostela (Spain). 2.3 Input data and processing We report here data from the commissioning phase and early physics run (from October 2015 to January 2017), where only two detector planes, stacked over a height of 120 cm, were used (T2 and T4 in Figure 2.1a). A trigger condition was defined as “at least one fired pad per plane, in time coincidence”. During data analysis, a standard equalization is performed automatically, aimed at the correction of the channel-by-channel variations in the time offsets and signal amplification along the lines of [98]. Provided both charge and time information are stored for each pad, noise signals (displaying zero-charge) can be removed in the next Figure 2.3: Cosmic ray rate for single vertical tracks as observed with the muon telescope (black curve) and ground-level pressure (blue curve). 34 IRMA RI´ ADIGOS S´ ANCHEZ processing step. Finally, “particle tracks” are formed by combinatorially matching the fired pads in both planes with a velocity compatible with the speed of light, within a 3-σtinterval, σt being the time resolution of the detector. This produces the final data sample ready for physics analysis, where any instrumental effects should be greatly minimized. We use in this work a data sub-sample, corresponding to events with a single track (multiplicity M=1), and a zenith angle θlower than 13◦. The former condition means that only cases with one fired pad per plane have been considered. The resulting mean rate is R= 9.05 Hz. The complete data taking period, displayed in Figure 2.3, may be conveniently divided into three phases: • From October 2015 to June 2016. During this first period the detector was run semi-autonomously (several interventions were needed) and problems related to faulty front-end electronics and high voltage instabilities were observed. • From June 2016 to October 2016. Maintenance work was carried out, the faulty electronics modules were replaced, and an online monitor was developed. • From October 2016 to January 2017. The detector run in stable conditions, in a fully autonomous way. Concerning the atmospheric variables, the vertical temperature profiles were retrieved from the European Centre for Medium-Range Weather Forecast (ECMWF) reanalysis, ERA-Interim [42], for the 1979-2017 period, at 37 isobaric levels (1000, 975, 950, 925, 900, 875, 850, 825, 800, 775, 750, 700, 650, 600, 550, 500, 450, 400, 350, 300, 250, 225, 200, 175, 150, 125, 100, 70, 50, 30, 20, 10, 7, 5, 3, 2, 1 hPa), with a horizontal spatial resolution of 0.125◦and a temporal resolution of 6 h. The surface pressure data is provided by a weather station of the Galician Regional MetOffice (MeteoGalicia) located at ∼100 m from the detector. Two exemplary temperature profiles for summer and winter at the Santiago de Compostela location are shown in Figure 2.4. Two important features are revealed from this image. Firstly, we can see how in the wintertime, the atmosphere is more compact and colder, and, as a consequence, the location of the tropopause is lower in the atmosphere (∼300 hPa). Sencondly, the tropopause region (100-300 hPa) is found much cooler in summer than in winter. We will see later how this circumstance becomes relevant. 35 Chapter 2. Atmospheric Temperature with a High-Resolution Cosmic Ray Detector Figure 2.4: Examples of ERA-INTERIM atmospheric temperature profiles for Santiago de Compostela, summer (1 July 2016) and winter (1 January 2017). 2.4 Analysis of Atmospheric Effects As mentioned before, several methods can be used to take into account the temperature effect of secondary cosmic ray particles in the atmosphere [24, 23, 51, 56, 135]. The integral method is one of the most precise [40, 46] but it requires knowing the distribution of the temperature coefficients in the atmosphere, WT(h). These can be theoretically calculated for different threshold energies, zenith and azimuth angles of incidence [47], or extracted from CR data [177]. In this work, we compare both approaches. An experimental determination of WT(h)is not straightforward, and requires special statistical techniques, given the presence of strong correlations between the temperatures of the different atmospheric layers. For illustration, Figure 2.5a shows a scatter plot of the temperatures corresponding to two different layers, from January 2015 to December 2016, and Figure 2.5b shows the pairwise correlation matrix obtained for the temperatures of the different atmospheric layers on top of the detector for the same period. In general, the low stratosphere (∼250-70 hPa) behaves opposite to the troposphere (∼1000-250 hPa) and high stratosphere (∼70-1 hPa). This is because the boundary layer (∼925 hPa) is positively correlated with the rest of the troposphere through convection, while an increase in its temperature will generally result in the low stratosphere cooling down. This is a typical condition observed for latitude regions above 40◦[107]. In the troposphere, heating is associated with convection in the tropics, while the cyclone-anticyclone dynamics in the mid-latitudes is the one forcing the air mixing. Colder regions have a lower tropopause because convection is limited there, for instance, in the polar regions. In general, if the tropopause rises, its temperature decreases. In mid-latitude regions, the cyclonic structures are characterized by a low tropopause, in association with a relatively cold troposphere and a warm lower stratosphere. In contrast, anticyclones uplift the tropopause and relate to a warm troposphere and a cold lower stratosphere. This explains the 36 IRMA RI´ ADIGOS S´ ANCHEZ Figure 2.5: (a) Correlation between temperatures for the atmospheric layers i=125 hPa and j=700 hPa (T0i,jis the mean value of each layer). (b) Pairwise correlations between temperatures of the different pressure levels considered in this analysis for Santiago de Compostela (from January 2015 to December 2016). paradox that tropopause temperatures are lowest where the surface temperatures are highest. To repeat, the extratropical tropopause temperature changes are positively correlated with the lower stratosphere and negatively correlated with those in the troposphere [139]. Figure 2.6 displays the correlations obtained for a high latitude location where the atmospheric conditions are radically different (Novosibirsk, Russia, 55◦N 82◦55’E). On this occasion, the correlations look quite different, especially in the mid-atmosphere. Given these points, any attempt to obtain the temperature coefficients by means of a multivariate regression will result in coefficients whose values do not correspond to the actual values. If explanatory variables of a multiple regression model are strongly correlated, they provide redundant information and violate the condition of non-collinearity required in a least-squares regression. The coefficients will also be highly sensitive to small changes in the model and their sign will be dramatically dependent on the variables considered. In other words, slightly different models might lead to different conclusions. In this way, we would never know the actual effect of each variable. Clearly, any phenomenological model aimed at reliably describing the measured rates needs to start from a sensible set of uncorrelated temperature variables, that need to be obtained beforehand. We adapt for the task the Principal Components Regression (PCR) analysis, which has been successfully used before for this type of studies in [177, 137]. 2.4.1 Barometric Effect Being much subtler, the temperature effect must be analyzed once the pressure effect has been removed. Moreover, in the case of gaseous detectors, the efficiency is a function of the ratio of the applied electric field Eand pressure P(represented by E/Pand dubbed reduced field), so even a high voltage and T-controlled environment is not sufficient to stabilize the detector response completely [108]. The above dependency means that the detector efficiency is anticorrelated with pressure and will add to the barometric effect at ground. Considering the atmospheric effect first, the relative change in the secondary CR rate caused by variations of the 37 Chapter 2. Atmospheric Temperature with a High-Resolution Cosmic Ray Detector Figure 2.6: Pairwise correlations between temperatures of the different pressure levels for a location at a latitude of ∼55◦N. ground-level pressure has an exponential dependence. To first order approximation, it can be expressed through a linear relation: R R0 =eβatm·∆P→∆R R0P≈β·∆P(2.1) where ∆R R0Pis the relative variation of the CR rate due to the pressure effect, R0represents its average value over the period under consideration, ∆P=P−P0is the deviation of the ground-level pressure with respect to its mean value (P0) over the same period, and β= βatm +βdet is the barometric coefficient, with βatm representing the atmospheric effect and βdet the detector contribution. The barometric coefficient was obtained separately for four different sub-periods, that displayed slightly different stability conditions (Table 2.1). An iterative linear fit was performed, with data outside a 2-σinterval removed from the fit (Figure 2.7). Compatible barometric coefficients were obtained, whose mean value was determined to be β=−0.59 ±0.02 %/hPa. This methodology allows us to remove any outliers in the data caused by detector instabilities and occasional space weather effects such as Forbush decreases or interplanetary events. Finally, the barometric effect is removed using: ∆R R0T=∆R R0obs −∆R R0P(2.2) where ∆R R0obs are the experimental CR variations and ∆R R0Tthe remaining variations due to the temperature effect. 2.4.2 Temperature effect Variations of the measured rate of the secondary cosmic component due to the atmospheric temperature effect can be approximated by a linear combination of some temperature 38 IRMA RI´ ADIGOS S´ ANCHEZ Figure 2.7: Example of the linear fit method used to obtain the barometric coefficient for one of the subperiods (3 October to 23 December 2015). The green lines delimit the points left out of a 2-σinterval after an iterative procedure. The blue line is the resulting regression line. Period β[%/hPa] 3 October 2015 - 23 December 2015 -0.607 ±0.003 24 December 2015 - 20 September 2016 -0.602 ±0.002 21 September 2016 - 22 November 2016 -0.583 ±0.003 23 November 2016 - 10 January 2017 -0.568 ±0.002 Table 2.1: Barometric coefficients for the different sub-periods. coefficients and the temperature variations at natmospheric layers [47], as discussed in Section 1.4.2 of the introductory chapter. Again, the corresponding expression is: ∆R R0T= n ∑ i=1 WT(hi)∆Ti∆hi(2.3) where ∆R R0Tare the relative variations due to the temperature effect; WT, given in % K−1atm−1, is the corresponding temperature coefficient for the atmospheric layer iat pressure hi;∆Ti= Ti−T0iare the temperature variations within the same layer with respect to its mean value (T0i), and ∆hi=hi−1−hiis the layer thickness, in atm. 39 Chapter 2. Atmospheric Temperature with a High-Resolution Cosmic Ray Detector Defining kxi=WT(hi)∆hi, equation (2.3) can be rewritten as ∆R R0T= n ∑ i=1 kxi∆Ti(2.4) that we denote formally as y=Xkx(2.5) where yis the vector of the measured relative variations ∆R R0T;Xis the (m×n) data matrix of the temperature variations whose columns are the temperature variations of the ith pressure level and kxrefers to the vector of temperature coefficients, that we want to estimate. As mentioned earlier, the coefficients of this model can not be obtained by an ordinary regression. For our purpose, we decided to use the Principal Component Regression (PCR) technique [91]. This method is applied when a dataset of variables shows multicollinearity, in our case, the temperature variations. The idea is to build new uncorrelated variables (called principal components), maintaining the information conveyed by the original ones, and use them as the new predictors to estimate the unknown regression coefficients of the model. The PCA consists of an orthogonal linear transformation that converts the original variables to a new coordinate system. The principal components (PCs) represent the directions of the data containing the highest variance. So, the first step is standardizing the ∆Timeasurements in X, dividing them by their standard deviations (over the analyzed period). This standardization is needed to prevent the variables with the highest variance from dominating. It causes a change in the notation, too. To keep it simple, we maintain the current notation but taking into account that all the following calculations are based on standardized variables. The principal components are the eigenvectors (directions) obtained from the covariance matrix of Xand sorted by the amount of explained variance. This set of orthogonal vectors forms a new basis in the new coordinate system. The matrix Xcan be transformed using the matrix of eigenvectors, defined as A(n×n), in the following way P=XA (2.6) where Pis now the matrix (m×n) containing the new variables in the new space. We got a set of uncorrelated variables because they were built using orthogonal eigenvectors. As a consequence, a new model can be built using variables P: y=Pkp(2.7) Now, the new set of coefficients kpcan be obtained directly using least-squares regression. Taking into account equation (2.6), we can write y=XAkp(2.8) The regression coefficients kpcan be transformed back into the original space using equation (2.5) and (2.8) kx=Akp(2.9) and multiplying by standard deviations in order to go back to the original scale. The year-to-year variability of the temperature data may affect the determination of the 40 IRMA RI´ ADIGOS S´ ANCHEZ principal components, particularly if exceptional temperature changes took place during the data acquisition period, such as Sudden Stratospheric Warmings (see more information in Appendix A.1). This can be seen in Figures 2.8 and 2.9, where several stratospheric temperature anomalies (warmings and coolings) can be seen during winter periods. Therefore, as a first step, we use a training dataset from a time series of the last 30 years to determine the PCs of the temperature data and avoid the influence of outliers corresponding to exceptional events. Figure 2.8: Temperature time series of the atmosphere in Santiago de Compostela from October 2015 to January 2017. Figure 2.9: Temperature anomaly for March 2016 were an example of a Sudden Stratospheric Warming is observed in the first half of the month. The upper stratosphere warms rapidly in a few days propagating way down into the troposphere in the next weeks. PCR typically uses only a significant subset of all the principal components P′to increase reliability. The components with higher variances are usually selected as the regressor variables for being the most important. No standard method exists for deciding how many components to retain. Anyhow, a good number of components should carry a high percentage of the total 41 Chapter 2. Atmospheric Temperature with a High-Resolution Cosmic Ray Detector variance (>70%). In order to decide the number of PCs to keep, we previously performed an analysis with reconstructed cosmic ray variations using a theoretical distribution of the temperature coefficients as a proxy [47]. These variations represent ideal data (i.e., without noise) only affected by the atmospheric temperature. Then, we apply the PCR method to these data to see how many PCs need to be kept in order to retrieve the original coefficients. To make this study more realistic, we follow the typical procedure of adding extra noise to the original data in three levels (low, medium, and high) to be able to analyze the performance of the technique. It was observed that with two components it is possible to restore the correct values of the coefficients until an acceptable level of noise. Including more components destabilizes the result (further details of the analysis can be found in Appendix B.1.1). Finally, the vector of coefficients k′ pis estimated by regressing the observed vector of cosmic ray data on the selected principal components P′using least-squares regression. So equation (2.7) is reduced to y=P′k′ p(2.10) where P′is now a matrix (m×r) whose columns are the corresponding subset of columns of P (and r<n). Using equation (2.9), k′ pcan be transformed back to the space of the actual temperature variables, providing the regression coefficients kxthat characterize the original model. Also, the relation kxi=WT(hi)∆hiintroduced before is taken into account when converting the estimated regression coefficients to the distribution of temperature coefficients WT, having a dimension %/K·atm. It must be noted that Partial Least Squares (PLS) regression could be an alternative to this technique because it is similar to PCR in that both select components that explain the most variance in the model. The difference is that PLS incorporates the response variable (the CR rate, in this case) into the analysis. One of the main reasons for not using this method is that our set of CR measurements is limited to a period of just two years, which would prevent a robust analysis. 2.5 Results and Discussion Figure 2.10a shows the distribution of temperature coefficients for the secondary cosmic component recorded at sea level for vertical incidence (blue line), where only events with a zenith angle θlower than 13◦were selected. The PCR method was applied to four different sub-periods of the data in order to account for systematic effects, expected to be mostly of instrumental origin at this stage, but also any remaining space weather phenomena having similar timescales to temperature variations. The average value and error bars are obtained from this combined analysis. As mentioned, the detector is placed in the first floor of a two-floor building. Therefore, the composition of the overburden material has been taken into account to estimate the value of the muons threshold energy, Eth ∼0.15 GeV, which is important to compare our coefficients with the theoretical ones (these can be found on the basis of integrations of the distributions describing muon production and propagation in the atmosphere). For illustration, the theoretical distributions for different muon energy thresholds and zenith angle θ=0◦(grey lines in Figure 2.10a) as given in [47] are shown. A good agreement with the ones obtained with PCR is 42 IRMA RI´ ADIGOS S´ ANCHEZ observed for low thresholds, while above 0.75 GeV a systematic deviation appears, specially close to ground level. Figure 2.10: (a) Distribution of temperature coefficients obtained with PCR for vertical tracks and comparison with the theoretical distribution energy thresholds: 0.15, 0.75, and 3.2 GeV. (b) Slopes ˜ WTobtained through a direct linear regression for the same data sample (θ<13◦) Although being compatible, the estimated values in the troposphere (>300 hPa) are systematically above the theoretical ones. This might be a consequence of the method itself but could as well reflect the presence of the soft component in the measured rates, given that it anticipates a positive correlation with the lower layers of the atmosphere [51]. The values also differ at high altitudes, in this case due to the constraint in the selection of the number of principal components in the analysis. The PCR method computes the principal components taking into account the variance in the temperature data. Then, the first components reproduce the general variations in the troposphere and stratosphere (seasonal changes). Variations in the high atmosphere are considerably more complex than in the surface, as illustrated in Figure 2.5. Increasing the selected number of components would help to reduce this effect in an ideal situation. However, the optimum number of selected components in our case is the one that allows to obtain the best results of the coefficients without the solution being destabilized by noise and other instrumental effects. On the other hand, it should be noted that the accuracy of the ECMWF reanalysis is worse at high altitudes (0-200 hPa) due to the lack of data (less satellite/balloon observations, etc), so the temperatures have an inherent source of error that surely increases the difficulties in obtaining the coefficients at those heights. Figure 2.10b shows for illustration the slopes determined before applying the PCR. Each coefficient e WTis obtained by direct regression between the relative variations of the CR intensity and the temperature variations for different layers: ∆R R0T=e WT(hi)∆Ti∆hi(2.11) These coefficients are dominated by the multiple correlations between the atmospheric layers. In particular, the slopes in the troposphere are negative, become positive in the low stratosphere, and return to negative values in the high stratosphere. The distribution of temperature 43 Chapter 3. Revisting the limits of atmospheric temperature retrieval from cosmic-ray measurements An analogy can be drawn between satellite and cosmic-ray measurements when targeting atmospheric temperature retrieval. Weather satellites employ the observations of electromagnetic radiations emitted by the atmosphere, that depend on its state. Temperature-dependent weighting functions, based on competing emission and absorption processes, need to be established beforehand for each atmospheric layer. Besides, they depend on the energy of the measured radiation [60]. The inverse problem of retrieving the atmospheric temperature profile can be solved by combining several energy channels. By the same token, CR measurements at different angles and energies (ground/underground) may be used for the same purpose. The following section will introduce a real case of data inversion, which will serve as an illustration to understand the constraints that must be taken into account when solving the inverse problem. Subsequently, we will proceed to characterize and evaluate in detail the methodology for retrieving the atmospheric temperatures. 3.2 A Case Study of Inverse Problem: Canfranc (LSC) Flux and angular distribution of high-energy muons have been measured underground at Canfranc Underground Laboratory (LSC), located under Mount Tobazo (1980 m) in the Aragonese Pyrenees [161]. Figure 3.1 shows the characteristic mountain profile along the railway tunnel under Tobazo, where the laboratory is located. To measure the CR flux, a muon monitor based on an array of scintillators arranged in 3 planes has been installed [101]. This station has the capability to track the muon trajectories and thus extract the angular distributions. It is worth mentioning the relatively small size of the detector, with an active area of 0.95 m2. Figure 3.2 shows the angular distribution of the muon flux measured in the LAB2400 as a function of the zenith and azimuth angle. The experiment used a sample recorded between October 2015 and March 2018. The asymmetry seen in the figure corresponds to the profile of the mountain above the laboratory. Muons are absorbed as they pass through the rock, so the more slant depth traversed, the less intensity detected in that particular direction. As a matter of fact, the maximum intensity observed corresponds to the direction of the Rioseta valley. These intensity values correspond to the average obtained during the period of Figure 3.1: Mountain cross section along the railroad tunnel that joins Spain and France under Mount Tobazo (1980 m). The LSC can be found at ∼800 m under Tobazo. 50 IRMA RI´ ADIGOS S´ ANCHEZ Figure 3.2: Muon intensity measured at Cranfanc Underground Laboratory (LAB2400) as a function of the zenith and azimuth angle (θ,φ). The maxima observed around θ=40◦and φ=150◦corresponds to the direction of the Rioseta valley [161]. measurement. Now, we would like to estimate the corresponding temperature effect in this data. To do this, we first need to know the threshold energies of the observations. There is an empirical relation that correlates the muon intensity from a given direction with the slant depth [106, 9]: I(X)≈IuX0 Xη e−X X0(3.1) where Iu=2.15 ±0.08 ×10−6cm−2s−1sr−1,η=1.93+0.20 −0.12 and X0=1155+60 −30 mwe (meter water equivalent) [7]. At the same time, we have seen in the Introduction (Section 1.3) the equation that gives the minimum energy, namely threshold energy, required for a muon at the surface to reach a depth X(Eq. 1.2). Therefore, we can estimate the depth for the different directions of the data displayed in Fig. 3.2 using equation 3.1 and then calculate the threshold energy by means of equation 1.2. Once we obtain the threshold energies for each angular bin in Fig. 3.2, we can plot a histogram grouping values of similar threshold energies to evaluate the statistics of the measurements. For the purpose of the analysis, we transform intensity rates to number of counts, N: N=I·A·∆t·∆Ω (≡R0∆t)(3.2) where Iis the intensity in units of cm−2s−1sr−1,Ais the area of the detector, ∆tis the measurement time, and ∆Ω corresponds to the solid angle of the angular bin. The statistical fluctuations associated to Nare then given by √N. Assuming a measurement time of 24 h and including the detector size, we obtain the number of muons recorded for each angular bin, that we group with Eth ·cosθ, as shown in Figure 3.3. The bars represent the errors associated with 51 Chapter 3. Revisting the limits of atmospheric temperature retrieval from cosmic-ray measurements Figure 3.3: Estimated number of muons for each energy threshold Eth ·cosθfor a time step of 24 h at the LSC site. The error bars represent the standard error for Poisson counting. counting. Figure 3.3 gives us an idea of the daily counting values of the detector as well as the features of the measurements. First of all, it is apparent from this plot that Canfranc has the potential to measure a wide range of threshold energies, from ∼70 GeV to ∼1000 GeV. It should be recalled that this is a consequence of the unique profile of the mountain. Secondly, the counting values are about 20 on average but there is a significant difference between the maximum and the minimum counting. Furthermore, these estimated values are relatively small and thus have statistical errors of ∼20 %. Now, we are interested in determining the temperature-induced variations at Canfranc location. Using the values of Eth ·cosθ, we can compute the corresponding effective temperature using the analytic expression given for high energies [6]: Te f f ≃∑N n=0∆XnT(Xn)(Wπ n+WK n) ∑N n=0∆Xn(Wπ n+WK n)(3.3) where Wπ,K(X)are weights, i.e. temperature coefficients, depending on pressure level Xnand Ethcosθ: Wπ,K≃(1−X/Λ′ π,K)2e−X/Λπ,KA1 π,K γ+(γ+1)B1 π,KK(X)(⟨Ethcosθ⟩/επ,K)2(3.4) The parameters A1 π,Kand B1 π,Krefer to meson production and attenuation in the atmosphere, 1/Λ′ π,K≡1/ΛN−1/Λπ,Krelates the attenuation lengths of the primary cosmic rays, pions and kaons, which are ΛN,Λπand ΛK, respectively. The temperatures values at different pressure levels included in the expression 3.3 are obtained from ERA5 reanalysis database for the Canfranc location [81]. The point is that we have to take into account the different variations that affect the real CR measurements if we want to use the data to access the temperature profile through the inverse problem. As we have seen in the previous chapter (Chapter 2), muon intensities are subject to variations of different origins. Here, we only assume those due to the temperature effect 52 IRMA RI´ ADIGOS S´ ANCHEZ Figure 3.4: The black line represents the statistic variations over time corresponding to the muon intensity of the threshold energy 71.13 GeV, for an area of 50 m2and a measurement time of 24 h; while the red line corresponds to the absolute temperature variations related to the same threshold energy. Part two and one of eq. 3.5, respectively. (systematic) and those referring to the statistical fluctuations of the measurement (error): ∆R R0obs =∆R R0T +∆R R0err =αT ∆Te f f Te f f0 +∆R R0err (3.5) In this equation, it is mandatory to have low statistic variations in order to measure with good accuracy the temperature variations. If statistic fluctuations are higher than temperature variations, we won’t be able to obtain any temperature information. From equation 3.2, two ways of dismissing statistic variations, i.e. increase counting, can be deduced: increase the measurement time or the size of the detector. Obviously, ∆Ω can be increased as well, but in the case of Canfranc, the relation between ∆Ω and the observation threshold is fixed by the mountain topology. Figure 3.4 presents the variations corresponding to Eth =71.13 GeV for time steps of 24 hours and assuming that the LSC detector has an area of 50 m2. In this scenario, the mean counting is 502.7, which has an error of 4.46 %. The black line represents the statistic variation as a function of time. It can be appreciated that whereas temperature variations (red line) are below 2 %, statistic fluctuations are so much higher (>4 %). In the original case (A=0.95 m2), the mean counting for the same threshold energy is 9.55 with an error of 30 % (see Fig. 3.3), making it even harder to retrieve any information of the temperature. It would be necessary to reach particle countings higher than ∼40000 (with an associated error of 0.5 %, similar to MINOS experiment [6]) to be able to obtain any valuable information. The only remaining option would be to increase the measurement time to more than two months, which is clearly not usable for temperature monitoring. As a consequence, we have shown that it would not be possible to observe the temperature variations in the real data with the current configuration of the LSC. The above can be written more formally as follows. Assuming a certain rate Rmeasured in a time ∆t, then the counting is N=R·∆tand the standard deviation ∆N=√N=√R∆t. Therefore, 53 Chapter 3. Revisting the limits of atmospheric temperature retrieval from cosmic-ray measurements ∆N N0 =1 √R∆t(3.6) If we want to estimate the temperature variations associated with the rates, we use the well-known relation ∆R/R=αT∆Te f f /Te f f0. So, it must be fulfilled that: ∆N N0 ≲∆R R0 =αT ∆Te f f Te f f0 (3.7) that can be rewritten using equation 3.6: 1 √R∆t≲αT ∆Te f f Te f f0 (3.8) and we get ∆Te f f Te f f02 ∆t≥1 Rα2 T (3.9) Considering that temperature variations are of the order of ∆Te f f Te f f0∼10−2−10−3, then ∆Te f f Te f f02 ∼10−4−10−6(3.10) From equation 3.9 and assuming that ∆t∼hours and αT≲1, we can obtain the value of the rates that we should have: R∼102muons/s ∼105−106muons/h Therefore, this is the limit to being able to retrieve temperature variations. Figure 3.5 shows how the data would look with real variations (gray line) for the case of 50 m2, where it can already be seen that the noise prevents discerning the temperature variations that are shown, by comparison, with the red line. To give a more quantitative idea of the dissimilarity, the correlation coefficient between both time series represented in the figure is R=0.11. Last but not least, it should also be noted that despite the potential of the LSC location, which allows accessing a wide range of threshold energies, and apparently to different effective temperatures, the reality is rather different. Figure 3.6 compares the effective temperatures corresponding to the maximum and minimum threshold energies of Canfranc data for two different periods, a “standard” atmospheric period and a Sudden Stratospheric Warming event. Both variations are practically identical and reflect the variations of the stratospheric temperature because the distributions of temperature coefficients associated with both threshold energies peak in the stratosphere (see the high threshold energies in Figure 3.7). As a result, both effective temperatures are supplying information about the same area of the atmosphere. And this is not very useful since our objective is to obtain the atmospheric temperatures at different levels. So far, we have seen that the temperature coefficients give us information about how the observed rates vary for each degree of temperature change at a given atmospheric level h. Moreover, we have seen that the coefficients are the sum of two terms: the positive effect 54 IRMA RI´ ADIGOS S´ ANCHEZ Figure 3.5: Rate variations corresponding to the temperature effect (red line) compared with data that includes both temperature and statistical variations (grey line). For Ethcosθ=71.13 GeV, A=50 m2and ∆t=24 h. Figure 3.6: (Left) ∆Te f f for two different threshold energies (71.13 GeV and 1066.91 GeV) during a SSW event. (Right) Same but for a “standard” atmospheric situation, such as summer. The black dashed line represents the variation of temperature at 50 hPa. 55 Chapter 3. Revisting the limits of atmospheric temperature retrieval from cosmic-ray measurements Figure 3.7: Temperature coefficients for vertical direction (θ=0◦) at several threshold energies corresponding to different underground depths [135], shown in black on the right axis. Coefficients for vertical muons observed at ground, tagged respectively by passage or absorption in 10 cm-lead and labelled as “hard” (Eth =0.4 GeV) and “soft” (E<0.4 GeV) are shown by continuous lines (black and red, respectively). The latter have been obtained from [51], and have axis on the left with the same units as the one on the right. related to mesons and the negative effect associated with muons. The sign of the coefficient gives information about the net value of the total effect produced by each layer, which depends on Eth. In the case of low threshold energies, the negative effect dominates, which means that for every temperature increase in that layer, the measured change in rates will have the opposite sign. However, as the threshold energy increases, the positive effect prevails. The latter effect becomes prominent around 100 hPa (∼15 km) where the peak of meson production takes place. This can be appreciated in Figure 3.7. As seen also in this Figure, selection by threshold energy (equivalently, depth) is the most natural way to separate the positive and negative effects in the weights, allowing a priori a higher sensitivity to the behaviour of the different atmospheric layers. One of the first attempts to estimate the atmospheric temperature profile using CR measurements with multiple detectors was carried out by [96]. By combining measurements at ground for soft muons (stopped in 10 cm of plastic scintillator), hard muons (passing through 10 cm of lead) and underground muons (80 mwe) performed with 1m2-area detectors they claimed a daily accuracy in the range 2-2.5 K for the atmospheric regions corresponding to 100, 500, and 900 hPa, over a period of about half a year. The detector suite was accompanied by measurements from a neutron monitor, in order to independently identify solar or interplanetary events that could bias the temperature estimate. Despite the sophistication of the approach, neither the depth of the detector was optimized nor angular information was used. In the 56 IRMA RI´ ADIGOS S´ ANCHEZ light of those very promising results, it is surprising for us that such natural extensions were not pursued. In fact, more recent studies have endeavored to obtain the temperature for more atmospheric layers by using a single detector and measurements at different angles (for example, [176, 178, 175]). Clearly one early limitation was detector complexity, as measurements were done over large areas by resorting to large scintillator tiles. However, multi-directional muon detectors (sometimes called hodoscopes) are nowadays common-place and available as part of the Global Muon Detector Network (GMDN) [133], for instance. Furthermore, the revival of the fields of muon tomography [126] and muography [116] has led to the adaptation of new technologies from particle physics, including for instance the development of extruded plastic scintillator [124], micropattern gaseous detectors [67, 109], as well as classic [15] and timing [173] resistive plate chambers (RPCs), just to name a few. They all offer affordable ways to cover large areas at high angular resolution. As an example, we have already seen in Chapter 2 that the effective atmospheric temperature has been measured at ground with a 2 m2timing RPC station, the first time that this technology, capable of time resolutions down to 50-60 ps and precise angular reconstruction on areas of several m2[170, 27], has been used for the task. In view of these powerful technological assets, the existence of new detailed calculations of the atmospheric coefficients, as well as the latest generation of accurate temperature data from the European Centre for Medium-Range Weather Forecast (ECMWF), reassessing the technological potential of cosmic rays for atmospheric temperature forecast seems very timely if not imperative. 3.3 Methods 3.3.1 Temperature effect Even for a perfect detector, cosmic-ray rates are subject to variations of diverse origins: those due to changes in the solar activity and in the atmosphere thermodynamic state are the most important. The former act as a potential systematic bias to the atmospheric temperature estimate, and in the remainder of this work we will assume implicitly that they can be isolated and eliminated. Although a natural way to perform this task is through the complementary use of neutron detectors (highly insensitive to atmospheric temperature variations), underground muon detectors (as the one proposed in text) may be sufficient, as solar and interplanetary events have less influence at underground depths. This occurs because interplanetary phenomena affect low-energy primary cosmic rays (in the range of a few MeV and GeV), which in turn are responsible for originating the low-energy secondary muons at sea level [65, 51]. Along these lines, changes in the measured CR rates induced by temperature variations can be approximated by the expression already presented in previous chapters (see equation 1.9). 3.3.2 Cosmic-ray intensity simulation For a mid-latitude location around 40◦, the intensity of vertical muons at ground is typically Ig≈70 m−2s−1sr−1[72, 76], the value used hereafter. In this case, the angular dependence has been parameterized as dN/dcos(θ)∝(cosθ)2for either soft (<0.4 GeV) or hard (>0.4 GeV) muons (e.g., [131]). Given the very high statistics for any angular bin, the particular choice of this distribution does not influence the temperature retrieval for muons reconstructed at ground level. 57 Chapter 3. Revisting the limits of atmospheric temperature retrieval from cosmic-ray measurements Figure 3.8: Parameterizations used in this work: vertical muon intensity (left axis) and threshold energy (right axis) as a function of vertical slant depth. A typical soil density of ρ=2.7 g/cm3has been assumed. In the case of underground measurements, expression 3.1 has shown a simple relation between muon intensity and slant depth, X. Since this parameterization is not accurate for shallow depths (<20 m) it has been complemented here by the one presented in [29], which is obtained from an approximation of the surface muon spectrum together with muon range tables. Figure 3.8 shows the muon intensity (in units of m−2s−1sr−1) as a function of slant depth using the aforementioned parameterizations for an average soil density of 2.7 g/cm3. For the simulation of underground muons we have assumed in the following an isotropic distribution impinging on a homogeneous soil slab, with rates and threshold energies for each angle obtained from eqs. 3.1 and 1.2. The “slab” denomination is because for each zenith angle θk, the corresponding depth ρkis calculated with ρk=ρ0/cosθ(see Figure 3.9) For a given detector with a fixed detection area, the number of counts measured over a period of time and solid angle is given by equation 3.2. For simplicity, the detection efficiency and acceptance have been assumed to be angle-independent and close to 1. The statistical fluctuations associated to Nare then given by √N. On account of that, the variations of CR rates can be expressed as previously indicated by equation 3.5, where (∆R/R0)obs are the experimental CR rate variations, (∆R/R0)Tthe changes due to the temperature effect, and (∆R/R0)err the associated statistical fluctuations. Specially underground, the detector area and measuring time become critical variables. Based on a preliminary analysis, and practical considerations, we set for a detector area of 4 m2, and a time interval of 6 h, although the impact of these choices in our analysis is evaluated at the end of the work. In sum, the procedure followed for the simulation of CR variations over a specific period of time can be sketched as: • For underground detectors, CR intensity and threshold energy are estimated using Eq. 3.1 and 1.2, respectively, for a certain slab thickness over the detector (depth). The approximation of [29] is used to calculate intensities for depths shallower than 20 m. Isotropic emission at ground level is assumed, for muons reaching underground. • For surface detectors, a vertical CR intensity of Ig∼70 m−2s−1sr−1is assumed, following 58 IRMA RI´ ADIGOS S´ ANCHEZ Figure 3.9: Relations of depths for a certain zenith angle θkin the slab assumption. ρ0is the vertical slant depth. an angular distribution like dN/dcos(θ)∝cos(θ)2. For the soft component, the value of the intensity is assumed to be about a fraction 0.4 of the hard component [51]. • Counting rates, N, are calculated for a fixed detector size, time interval and solid angle as indicated in Eq. 3.2. • Variations in CR rates due to the temperature effect are calculated with Eq. 1.9 using the temperature time series from ERA5 and the temperature coefficients for the corresponding energy, Eth, and angle, θ(linearly interpolated when needed). A more detailed discussion of the estimation of these coefficients is presented below in Section 3.3.4. • Poisson noise with a mean value of Nwas added. 3.3.3 Temperature data Vertical profiles of atmospheric temperatures were retrieved from ECMWF reanalysis, using the ERA5 dataset which offers 37 isobaric levels (1000, 975, 950, 925, 900, 875, 850, 825, 800, 775, 750, 700, 650, 600, 550, 500, 450, 400, 350, 300, 250, 225, 200, 175, 150, 125, 100, 70, 50, 30, 20, 10, 7, 5, 3, 2, and 1 hPa), with a horizontal spatial resolution of 0.25◦ and a temporal resolution of 6 h [81]. A mid-latitude location at 40◦was chosen, in this case corresponding to Santiago de Compostela (Spain). 3.3.4 Temperature coefficients The distributions of temperature coefficients (WT) have been calculated before by several authors. Dorman supplied the most extensive calculations for a wide variety of threshold energies and zenith angles [50]. These estimates were later re-evaluated by Sagisaka and Dmitrieva [135, 47]. The former provided coefficients for various combinations of threshold energy and zenith angle whereas the latter introduced up-to-date parameters in the calculations to give a vast database of coefficients, focused on threshold energies for surface hodoscopes. 59 Chapter 3. Revisting the limits of atmospheric temperature retrieval from cosmic-ray measurements Figure 3.15: RMSE between estimated and real temperature for each atmospheric pressure level, in the multi-channel analysis. (a) Hard + underground muons at different depths. (b) Hard + soft + underground muons at different depths. Depths smaller than 18 m and larger than 21 m are excluded because they do not exhibit any discernible difference compared to those. selected a depth of 55 m.w.e, which is equivalent to a threshold energy of ∼11 GeV or, in other words, to a soil-thickness of 20 m. However, the early temperature coefficients there assumed are much more peaked than the ones used here, and would correspond to a threshold of 50 GeV (80 m-depth) if resorting to more modern estimates as those shown in Figure 3.7. The performance of the two methods is very different too, as can be appreciated in Figure 3.16 where the RMSE from the three-station/single-channel analysis as proposed in Miyazaki and Wada (black dashed line) is compared with the present one (blue line). Despite the assumed depth is the same in both cases, the difference in the assumed weights makes all the difference, allowing us to establish the relevance of the 20 m-depth as an actual optimum for atmospheric studies. The improvement is even more apparent when considering a multi-channel analysis (orange line): up to a factor of two or more can be gained in critical atmospheric regions like the tropopause and stratosphere compared to earlier simulation work. The reconstruction reaches a best value of 0.8 K at 850 hPa and a worse one around 2.2 K for the tropopause region and up to 50 hPa, becoming the temperature intrinsically inaccessible above the 10 hPa layer. Moreover, Miyazaki and Wada estimated RMSE-values between 1 and 3 K when disregarding statistical noise from counting, for seven pressure levels between 1000 and 100 hPa. Interpreting that as the intrinsic limit to the inversion problem and comparing to the analogous result in our analysis (red-dashed line in Figure 3.16), it can be concluded that a three-station/multi-channel analysis with optimized depth provides an overall improvement of around a factor 2 also in that situation. 66 IRMA RI´ ADIGOS S´ ANCHEZ Figure 3.16: RMSE between estimated and real temperature for each atmospheric pressure level. Observations performed in 6 h time intervals and for 2 m ×2 m detection area. Black dashed: three-station/single-channel analysis using weights from [114]. Blue: three-station/single-channel analysis at an optimized underground depth. Orange: three-station/multi-channel analysis at optimized underground depth. Dashed red: three-station/multi-channel analysis at optimized underground depth, neglecting statistical fluctuations in particle counting. The grey dashed line represents the intrinsic spread of the temperatures for each atmospheric layer. 3.5 Discussion Within the relatively simple inversion algorithm proposed in this work, the width of the optimal-depth plateau in Figure 3.15 exhibits a non-trivial dependence with the chosen angular binning. For depths around the optimal one, this is exacerbated since the overall correlation between rates and temperature (eq. 2.12) changes sign (hence, it vanishes) as a function of the angular channel. It does so in a way that is both abrupt and critically dependent on the precise shapes of the weights. While this fact suggests that a finer resolution is desirable (1 degree is technically possible without great effort), we opted to leave such a study outside this work. The reason is two-fold: i) the binning becomes too thin compared to the four angular bins available for our underground coefficients in [135], and so, in the absence of new calculations of those, the present simulation work would depend largely on the interpolation method and ii) the ×10 increase in the number of fitting parameters would require of a more dedicated optimization study than intended here. In our case we have relied on the python package Statsmodels, and the Ordinary Least Squares method included in it, without constraints in the fitting parameters (some examples of regression plots can be found in Figure B.3 of Appendix B.2.1). The results were little sensitive to the method chosen or how the regression was conditioned (initial values, parameter range, function tolerance and linear constraints between variables). Studies performed with a mildly increased binning (5 degrees) show indeed that the results of the regression within the 19-20 m plateau become much more stable. Once the existence of an optimal depth-plateau has been established, it is important to understand how the performance of the inversion algorithm depends on the size of the detection area and the presence of systematic errors in particle counting (that we have simulated getting random samples from a Gaussian distribution with the width being a certain percentage of the average rate and subsequently adding them to the data series). To make the latter more 67 Chapter 3. Revisting the limits of atmospheric temperature retrieval from cosmic-ray measurements realistic, we follow the typical procedure of adding extra noise to the data in different levels (low, medium, and high) to analyze the performance of the technique. Figure 3.17 seems to show a preference towards areas at the scale of few m2, as there are just marginal gains compared to the case of infinite area (a condition well fulfilled for practical purposes above 300 m2), whereas performance deteriorates very perceptibly below 1 m2. A larger detector will be more resilient against systematic variations in counting too, with a 4 m2 detector able to tolerate up to ∼0.3% additional fluctuations in counting, assuming they are uncorrelated at every time step (Figure 3.18). This poses a very stringent requirement for the detection system, whose overall efficiency should be kept stable within these values. From this point of view, plastic detectors coupled to photon sensors represent a most natural choice, although a gaseous detector with redundant layers could become more affordable/practical at the expense of a larger design complexity. Results presented here are difficult to interpret from an atmospheric physics perspective, so in order to get a better grasp of how the minimization process works, a close examination of what we define here as “combined temperature coefficients” will show to be useful. For that we rewrite equation 3.11 taking into account expression 1.9: ∆ˆ Ti= nst ∑ k=1 nch ∑ j=1 cjki nl ∑ p=1 WT(Eth,jk,θj,hp)∆Tp∆hp!k (3.17) Rearranging the terms for the same pressure level pwe obtain: ∆ˆ Ti= nl ∑ p=1 nst ∑ k=1 nch ∑ j=1 cjkiWT(Eth,jk,θj,hp)!∆Tp∆hp(3.18) from where we define the combined temperature coefficient as: Wpi ≡ nst ∑ k=1 nch ∑ j=1 cjkiWT(Eth,jk,θj,hp)(3.19) Figure 3.17: RMSE between estimated and real temperature for each atmospheric pressure level, for a three station/multi-channel analysis at an optimal depth. Each line represents a different value of the size of the detectors. The grey dashed line represents the standard deviation of the temperatures for each atmospheric layer. 68 IRMA RI´ ADIGOS S´ ANCHEZ Figure 3.18: RMSE of the estimation of temperature for each atmospheric pressure level, for a three station/multi-channel analysis at an optimal depth. Different levels of Gaussian noise have been aggregated to the rates: 0.1, 0.5, 1, 4, and 10%. The grey dashed line represents the standard deviation of the temperatures for each atmospheric layer. When solving the inverse problem for the temperature of a certain pressure layer, one may expect that the coefficients cjki of the regression should in principle have values such that in eq. 3.18 the combined coefficients are able to enhance the p-term of the same temperature layer (i), minimizing the contribution from the rest. As the shapes of the temperature coefficients are not flexible enough to accommodate this condition for any arbitrary layer, the correlation between atmospheric layers becomes an essential ingredient. Figure 3.19 shows some examples of the combined coefficients obtained for the retrieval of the temperature at 50, 150, 200, 500, 850, and 1000 hPa for different depths of interest: 15, 19, 20, 21, and 100 m. For the highest and lowest atmospheric layers (1000 and 50 hPa, represented at bottom right and top left in Fig. 3.19 botton right, respectively), the combined coefficients have indeed the highest value on that layer. The presence of a significant contribution from the other atmopsheric layers limits temperature reconstruction here. Indeed, the effect of the width of the combined coefficients can be seen in the reconstruction of the 850 hPa layer (chosen since it corresponds to the layer with the lowest RMSE in the analysis) that is very similar to the 1000 hPa one. On the other hand, the combined coefficients for the case of 200 and 500 hPa (middle subplots in Fig. 3.19), especially at optimal depths, rely strongly on the correlation between layers. In the case of 200 hPa for 15, 21, and 100 m detector depths, the regression has chosen to give positive weight to the stratospheric levels and negative one to the tropospheric ones. As these temperatures are correlated, this approach minimizes both contributions enhancing the ones for the intermediate layers. The combined coefficients for 19 and 20 m use the atmosphere information differently, which is a consequence of the better use of the underground component included in the regression. That is, for 20 m (Eth ∼10 GeV), the troposphere exhibits positive combined coefficients and even softly peaks in the corresponding 200 hPa layer. As the temperatures at 200 hPa are still correlated with the tropospheric ones a better sensitivity ensues, inherited from the weights’ shapes at underground depths (Figure 3.14b). Another interesting example is the situation for the 500 hPa level. The depths of 15, 21, and 100 m prefer to rely on coefficients from layers above 200 hPa to predict the temperature. 69 Chapter 3. Revisting the limits of atmospheric temperature retrieval from cosmic-ray measurements Figure 3.19: Combined temperature coefficients calculated for the atmospheric pressure levels of 50, 150, 200, 500, 850, and 1000 hPa in the three-station/multi-channel configuration. The curves represent different depths of the underground station: 15, 19, 20, 21, and 100 m. However, a better accuracy can be achieved if relying on coefficients from layers below 800 hPa (depths 19-20 m), as this makes use of the positive correlation between troposphere and stratosphere and thus it balances both contributions to the temperature of the 500 hPa layers. Even more extreme is the case of the 150 hPa level (top right in Fig. 3.18), which is strongly anticorrelated with the high stratosphere and troposphere. The regression is succesfully able to use this information for the depths of 10 and 20 m, even allowing the coefficients to peak in the nearby region around 100 hPa (19 m depth). By comparison, the regression performed for other depths attempts desperately to give weight to all the stratospheric levels down to 200 hPa, in order to achieve some sensitivity to the 150 hPa one. As the high stratosphere is anticorrelated with the 150 hPa layer, the strategy is bound to fail, producing the large peak-structure on the RMSE observed in this region for non-optimized depths, and that has been presented abundantly throughout the text. As a consequence of the above observations, it becomes clear that the cosmic-ray method for atmospheric temperature retrieval involves a suitable weighting of the information from the entire atmosphere, for each layer whose temperature is being resolved. In particular, the fact that the variations in the lowest part of the stratosphere and around the tropopause (∼100-200 hPa) are partly anticorrelated with the variations in the troposphere and upper stratosphere 70 IRMA RI´ ADIGOS S´ ANCHEZ Figure 3.20: RMSE of the estimation of temperature for each atmospheric pressure level, for a three-station/multi-channel analysis. The simulated temperature variations have been taken from data but introduced in uncorrelated fashion, for each 6 h interval. σ(T) is the standard deviation of temperatures. (∼70-1 hPa) determine the shapes of the combined coefficients for intermediate atmospheric layers. Simulating a custom atmosphere with uncorrelated temperature variations allows to better illustrate the situation. In this case the total spread of each layer is kept, and assigned randomly at each time interval. With 37 layers, and in light of the broad Wp’s obtained above, this should represent a very harsh situation for the linear regression, as shown in fact in Figure 3.20. It is relevant to see that the weights in the region 19-20 m still offer the highest sensitivity, confirming that a higher variability of shapes of the temperature coefficients distribution in that region plays a part in the overall performance of the regression. The dominant contribution, however, must come from the long-distance correlations between different layers, and the ability of the weights to accommodate them, as indicated by Figure 3.19. Overall, our study shows that the introduction of different angles into the analysis of the inverse problem helps improving the estimates of the atmospheric temperature profile dramatically. A temperature predictability (RMSE) ranging from 0.8-1K in the low troposphere to 2.2 K in the 50 hPa level was obtained, deteriorating for higher stratospheric layers. For a temperature monitoring emplacement dedicated to improve atmospheric forecasts, this range of heights is more than enough since the most interesting atmospheric phenomena occur in the tropospheric layer. However, there is a unique phenomenon that takes place in the stratosphere and is attracting a lot of atmospheric scientists due to its capacity to modify the weather at the surface. We are referring to Sudden Stratospheric Warmings (SSW). A SSW is an event that occurs in polar vortices when the stratospheric temperature suffers an abrupt increase in a short period of time. In such events, the vortex may collapse, releasing cold air towards lower latitudes that could impact the surface weather. This situation is more likely to happen in the northern hemisphere (see more information in Appendix A.1). In many cases, the monitoring of the temperature at 10 hPa (for which a modest 5 K RMSE was obtained in this analysis) is useful when a major event occurs. Still, the observation of any other lower stratospheric level provides very valuable information about these events. In such a case, the setup proposed in this work would be at least complementary to this kind of research. For illustration, a comparison between the input temperatures and the ones estimated with the inversion method proposed in 71 Chapter 3. Revisting the limits of atmospheric temperature retrieval from cosmic-ray measurements text is shown in Figure 3.21. The difference between the observed and estimated temperatures is also included in Figure 3.22. The bias calculated for these predictions is zero for the pressure levels presented. However, it can be seen that the upper atmosphere displays the most significant errors, which correspond to the occurrence of the aforementioned SSW events. The model loses detail when capturing these events, as we have mentioned above, yet it is able to reproduce its influence on lower levels (between 20 and 200 hPa). Figure 3.21: (a) Observed temperature for 2019. (b) Estimated temperature for 2019 using the detector configuration and data analysis discussed in text. 72 IRMA RI´ ADIGOS S´ ANCHEZ Figure 3.22: Difference between observed and estimated temperature for 2019 using the detector configuration and data analysis discussed in the text. 3.6 Conclusions The purpose of the present study was to determine the best configuration of a monitoring station integrated by cosmic-ray telescopes. One of the most significant findings to emerge from this study is the possibility to retrieve the temperature of the atmosphere from the surface with good accuracy up to a considerable height (∼20 km). An implication of this is that atmospheric muon detectors can be used in scientific research beyond the field of astrophysics. Returning to the question posed at the beginning of this chapter, it is now possible to state that multi-directional telescopes can enhance the estimate of atmospheric temperatures. Our results have shown that one can achieve a high degree of accuracy with error margins between 0.8 and 2.2 K that could be improved in practice applying more advanced statistical techniques such as those employed in satellite observations [59, 134]. For 4 m2-scale detectors, this performance requires of an outstanding detector stability (below 0.3% on its counting efficiency). While detector inefficiencies to minimum ionizing particles down to 0.1% are not alien to particle physics instrumentation, the requirement will pose significant constraints on the chosen technology and detector design. On the other hand, the current work was limited by the use of simulated CR data and further work needs to be done in the area of experimentation to evaluate the actual limits of the estimates. In spite of this, our findings establish several courses of action for future research. A good line for future work would be to contrast the retrieved temperatures from real cosmic-ray data against temperature data from balloon measurements. Additionally, we have found evidence of an optimal depth to place one of the detectors. Depths of 19-20 m seem feasible and affordable, with no need to go deep into a mountain as is done for other CR research. Such an underground station could be easily located in dams, tunnels or subway stations for instance. Moreover, continuous measurements of vertical temperature provided in this way would no doubt complement satellite measurements, as this technology would be much more affordable and easy to assemble and maintain. 73 Chapter 4 Simulating the atmospheric effects with a dummy model Abstract: In this chapter we present the atmospheric effects analyzed by means of a cosmic-ray air shower simulation code. This tool serves as an additional testing to corroborate and better understand the correlations between cosmic-ray variations at the surface and atmospheric properties. 4.1 Introduction Chapter 2 has paved the way for introducing the influence of the atmosphere on the observations of CR at the surface. We have seen how tricky it can be to identify the individual effects in the total rates, especially when the atmosphere exhibits a complex behavior. Over the past decades, a number of works have studied in detail the correlations between CR observations and atmospheric variations [e.g., 120, 23, 184, 118]. However, it can be difficult to make a correct interpretation of these relationships when there are so many factors at play. Numerical simulations of EAS can provide a deeper insight into the CR variations, especially because the influencing factors can be controlled. For example, the primary CR spectrum can be fixed, which makes it possible to exclude variations caused by space weather. 4.2 Simulation Overview Numerous programs are available for the realistic simulation of EAS [65, 82]. Some of the most widely used are CORSIKA and AIRES [57, 138]. These programs consist of a series of routines and subroutines that track the paths in the atmosphere followed by the particles created in the cascade produced by a high-energy primary cosmic ray. Figure 1.3 in the introduction section shows a pair of examples of 3D air showers simulated with CORSIKA. As we have seen so far, the Earth’s atmosphere is the medium in which air showers propagate and, therefore, their evolution strongly depends on its state (density, temperature, etc). Thus, the program must incorporate a realistic model of the atmosphere so that the simulations are as accurate as possible. All programs that simulate EAS use what is known as US Standard Atmosphere as a model [13]. It is based on experimental data and constitutes Chapter 4. Simulating the atmospheric effects with a dummy model the temperatures calculated using equations 4.14 and 4.15: p=p0e−gMh RT (4.16) where p0is the surface pressure. Finally, the atmospheric density and depth can be calculated using equations 4.10 and 4.13. 4.3 Simulation Setup We focus the simulation study on uderstanding the observed atmospheric effects by means of the muon/meson balance in the atmosphere. The pressure and temperature effects are analyzed separately. In each simulation run (i.e., each atmospheric profile), 104primary protons are launched down to the surface. This number has been carefully selected so that the statistical fluctuations of the output data are low enough to allow observation of the atmospheric effects. Then, the energy of the primary particle is selected from the probability distribution 4.1 whithin the interval ranging from 100 to 104GeV. The recorded outputs are the energy of muons reaching the surface and the production height of pions and muons. For the surface muon energies, we only save the ones having an energy above ∼4 GeV. 4.4 Results 4.4.1 Altitude Effects To study the pressure effect, we have modeled the seasonal variations of the atmosphere during 2015 and analyzed their effects on the muon flux variations at the surface. We have averaged the real data to obtain the atmospheric profiles for the different seasons. For the purpose of this exercise, it is not necessary to analyze the whole time series, day by day. On the one hand, it would be computationally very expensive and, on the other hand, considering only seasonal variations makes it easier to study the barometric effect. Figure 4.3 shows the simulated results for the barometric effect. In agreement with what we have obtained experimentally, the variations of muon flux are anticorrelated with the surface pressure variations. We have mentioned several times throughout this dissertation that mesons produced in the cascades can either interact or decay into muons. Therefore, we analyze the production heights of pions and muons to understand the reason of the observed seasonal behaviour. Figure 4.4 displays one example of the distribution of the pion production height. As can be appreciated, our simulation predicts a maximum shower production level at an altitude of 7 km. Although it should be noted that in our case it appears at larger depths than the real one (about 12-15 km). However, we must remember that we are employing a simplified 1D-model. The particles of the cascade are moving vertically downwards in the atmosphere, traversing less amount of atmosphere than if they were traveling with a certain angle of inclination, as in real life, and therefore in our case they travel down deeper in the atmosphere. The exact value of the peak is obtained by fitting the data with a Gaussian function in the atmospheric region around the maximum (blue line in Fig. 4.4). In the exampled shown, a value of 7.273 km has been obtained, which corresponds to altitudes in the 300 hPa region. 82 IRMA RI´ ADIGOS S´ ANCHEZ Figure 4.3: Relative variation of muon rates at the surface (Eth =3.2 GeV) as a function of the surface pressure variations. Figure 4.4: Distribution of pion-production height as a function of the atmospheric height. Height zero represents the surface level. The blue line represents the gaussian fit to the data in order to obtain the maximum height of production. 83 Chapter 4. Simulating the atmospheric effects with a dummy model Figure 4.5: Simulated maximum height of production as a function of surface pressure variations for pions (a) and muons (b). The maximum of production should shift when the surface pressure varies. To analyze its displacement, we fit all the simulated distributions to a Gaussian function in the same way as indicated above. The results are shown in Figure 4.5a as a function of the variations in surface pressure. Note that maximum height of production appears higher for lower surface pressure values. Figure 4.5b looks at the muon production maximum heights. First of all, the values of the peaks are similar to those of the pions, although slightly lower (∼200 m). This proves that pions are decaying almost immediately after creation, just as expected. As a consequence, the muon production peak shifts in the same way with pressure variations. Further inspection of the results shows that the peak variation is only a few hundred meters between different seasons (raging from 7200 to 7700 m in the case of pions). This is due to the fact that the surface pressure variations between seasons for the year we are analyzing are very small. The height of the maximum depends on the altitude of interaction of the primary particles. The interaction probability of a proton depends on its energy and the amount of air mass traversed, i.e., atmospheric pressure. If its fate is to interact after traversing 300 hPa of mass, the only thing that changes is the altitude where that specific pressure level is located. As a result, the production height of pions and muons changes according to this. To test this hypothesis, we plot the maximum peak of the muon production versus height variations for the 300 hPa pressure level. The results can be seen in Figure 4.6. The lowest altitudes of the 300 hPa pressure level coincide with the lowest values of the production peak. 4.4.2 Temperature correlations We simulate a period of one year with steps of 25 days to obtain a small sample of cosmic-rays flux. A full year with 6-hour steps would take a long time, but a shorter period is enough to corroborate the experimental results. The simulated rates are correlated with the temperature variations at different heights to obtain the regression slopes e WT, similar to those obtained in Chapter 2 in Figure 2.10b. Figure 4.7 shows the results of the regression of the simulated data, indicating a good agreement with the expected values based on the theoretical weights WTby Dmitrieva et al. (assuming a threshold energy of 3.2 GeV). The fact that we are 84 IRMA RI´ ADIGOS S´ ANCHEZ Figure 4.6: Simulated maximum height of production for muons as a function of height variations in the pressure level of 300 hPa. Figure 4.7: Slopes e WTobtained through a direct linear regression with simulated cosmic ray data (yellow line) compared to the expected values (blue line). able to reproduce the theoretical results confirms that the atmospheric effects seen in the real data are indeed dominated by the absorption/production processes of the secondary particles in the atmosphere, mainly muons and pions. As already indicated in the discussion of Figure 2.10, when doing a direct regression the correlations between atmospheric layers also play a role, and indeed the strong variations in the region 50-250 hPa are a footprint of the tropopause dynamics. 85 Chapter 5 Modeling the Influence of Cosmic Rays on the Atmosphere Abstract: The flux of cosmic rays in the atmosphere has been reported to correlate with cloud and aerosol properties. Several mechanisms have been proposed and tested to explain this effect, leading to the conclusion that the induced effects were minor. However, these studies did not disprove the link between cosmic rays and clouds (i.e., climate). Since then, some different mechanisms that could be relevant to aerosol growth have been postulated. In this chapter, we use a global chemistry transport model to include the effects of charging on the microphysical development of aerosols. We will compare the variations of cloud condensation nuclei (CCN) concentrations between the solar maximum and solar minimum. This study aims to discover the complex relationship between GCR and aerosols. *This chapter includes content from the following article: I. Ri´ adigos et al.. The Charge of aerosols from Cosmic Rays and the enhancement of cloud condensation nuclei formation. In preparation, 2022. 5.1 Introduction In the previous chapters, we have covered the subject of the atmospheric effect on the CR flux measured near the surface. However, as was mentioned in Chapter 1, the charged population of CR can affect those atmospheric processes where the ionization, electric field, or particle charges play an important role. The charging of the atmospheric aerosols and the subsequent creation of CCN, as pointed out in Chapter 1, are among those processes. Thus, in this last chapter, we address the opposite question to the one that has been covered so far in the previous chapters, i.e. how CR can affect atmospheric conditions, in particular those related to cloud formation. It should be also mentioned that while in Chapters 2-4 we were focused mainly on the muon component of secondary CR, in this chapter the entire CR flux is considered to calculate atmospheric ionization rates. Besides, the barometric/temperature effect previously analyzed is direct, however, the CR-cloud effect that we are going to study is considered to be an indirect effect, because there is no linear proportionality between changes in CR flux and the variations of CCN. One of the most debated aspects of aerosols is the role that Galactic Cosmic Rays play in their growth. We have seen in the introductory part of this thesis (Section 1.5.3) Chapter 5. Modeling the Influence of Cosmic Rays on the Atmosphere that numerous studies have reported strong correlations between the cosmic-ray flux and aerosol-cloud properties (e.g., [148]). To date, literature has emerged that offers contradictory findings about this issue [102]. In addition to this, it has not been easy to find the corresponding process to explain such connection. Numerous mechanisms have been proposed but have failed when testing their relevance to CCN formation. One of the most promising hypotheses has been the so-called “aerosol clear-sky mechanism”, which involves the nucleation process. Nucleation is a process affecting the finest atmospheric particles by which they aggregate to form small clusters, giving birth to the smallest aerosols. It is one of the most important processes in the creation of aerosols in the atmosphere. From here, condensation and coagulation are responsible for these small clusters to grow to CCN sizes (≳100 nm). The former causes the growth of aerosols through condensation of vapors (generally sulfuric acid and low-volatility organics), and the latter refers to the attachment of two colliding aerosols to form a larger one. It is well known that the presence of small ions created from the ionization of atmospheric particles by GCR can enhance the nucleation rates. If the newly created particles do not get lost along the way, they can grow to CCN sizes. One of the possible ways of being lost is through coagulation with existing CCN particles. The balance between coagulation losses and growth to CCN will determine how much GCR induced-ionization can influence the number of CCN and eventually cloud cover. These are the fundamentals of the clear-sky mechanism. During the past decade, Jeff Pierce has presented numerous reports and collaborated on several studies to shed light on the relevance of this mechanism in the atmosphere [122, 143, 121]. The state-of-the-art CLOUD experiment at CERN has made significant contributions to our understanding of the relation between aerosol nucleation and GCR. The CLOUD collaboration was the first to experimentally demonstrate the impact of GCR on nucleation rates and has provided parameterizations for this relationship that can be easily incorporated into atmospheric models for further analysis [54]. However, they have determined that the impact of GCR on nucleation (i.e., the clear-sky mechanism) is not sufficient to explain the correlations found between cosmic-ray flux and clouds. These results have opened the door for proposing new mechanisms that can explain this elusive link. Recently, another group of researchers has explored theoretically and experimentally the possibility that cosmic rays may enhance the condensation rates of aerosols [149]. The proposed hypothesis states that an increase in ionization results in faster aerosol growth by condensation, which prevents them from being lost by coagulating with existing particles. The point here is that the charge that aerosols acquire when ions condense on them has traditionally been neglected. This has been assumed because the flux from neutral molecules (such as sulfuric acid) to aerosols by condensation is much higher when compared to the mass flux from the ions. By illustration, the typical ratio between them is 10−3. Svensmark et al. argue that this small ion flux should not be underestimated. So far, only the condensation of neutral molecules in aerosol growth has been taken into account. Svensmark et al. also consider that aerosols can be charged by the condensation of ions and at the same time neutral gas can condensate onto charged aerosols. They develop a model where all these interactions are taken into account. This model considers that aerosols can be positively, negatively, or neutrally charged and that the condensible gas contains positive and negative ions from cosmic rays as well. They define βas the interaction or attachment coefficient (in units of m3/s) between gas molecules and aerosols. The key is that βhas different values depending on whether the particles in play are charged or uncharged. If the electrostatic interactions between charged particles are taken into account when computing this coefficient, it is observed that for small particles its value 88 IRMA RI´ ADIGOS S´ ANCHEZ is greater than the original coefficient of neutral particles. As a consequence, small charged particles would condense faster than neutral ones. In their work, Svensmark et al. present experimental results showing that the presence of ions seems to support their hypothesis under some atmospheric conditions. However, when estimating the interaction coefficients they make too many simplifications, such as considering a constant temperature or setting the mass of the neutral gas to a value of 100 AMU. This is only a small representation of the atmospheric conditions. To test the real impact of this mechanism, the proposed scheme should be incorporated into an atmospheric model. In 2020, Svensmark et al. presented a numerical approach for this mechanism that could be implemented in atmospheric models [152]. However, numerous challenges arise when it comes to implementing this model. On the one hand, the calculation of the interaction coefficients for charged particles requires a lot of computational resources. This implies that tables should be previously created with the values of the coefficients as a function of different variables, which could then be used to interpolate the corresponding values in the model. As we will see later in this chapter, there is a much faster and accurate approach to do this calculation. On the other hand, it has been shown that in the real atmosphere particles are capable of acquiring multiple elementary charges [93, 168]. So this must be taken into account and turns out to be one of the objectives of this work. To take into account that aerosols can accumulate a large number of charges on their surface, their charge distribution must be estimated. The charge distribution affects the coagulation process since if two colliding aerosols carry charges of identical signs, a repulsive force will appear and inhibit their union. In contrast, if the aerosols carry charges of opposite signs, their coagulation will be enhanced. Thus, it becomes relevant to incorporate the charge distribution in the model but, as we will see, doing it explicitly is almost impossible and computationally very expensive. However, there is an approach that can be adopted for this purpose and will be explored in the following. Improving our knowledge of aerosol growth is important to better understand changes in clouds. Any change in the global cloud cover modifies the terrestrial albedo by increasing or decreasing the warming effect on climate. One of the biggest unknown factors in climate prediction is how clouds vary under different conditions. As we mentioned, one of the reasons is that the exact mechanisms of CCN formation are hardly known in detail. These uncertainties lead to a large dispersion in the climate predictions regarding the average temperature increase for the following decades. Therefore, it is crucial to go deeper into the study of cloud formation processes as climate change is one of the greatest concerns of our generation. In short, numerical models that describe in detail all aerosol microphysical processes are very demanding from a computational point of view. Expanding upon previous works, we want to include the effects of atmospheric charging from CR in aerosol growth to test and understand the relationship between cosmic radiation and aerosols. The most important processes for aerosol growth are nucleation, condensation, and coagulation. So far, only the former has been parameterized and included in global aerosol models taking into account the effect of ions in the process (ion-induced nucleation) [181, 54, 69]. However, multiple investigations have proven that this mechanism alone is not strong enough to produce substantial changes in the final CCN concentrations (e.g., [143]). Therefore, our work is devoted to introduce the charging effect of ions into the condensation and coagulation processes using a global 3-D atmospheric chemistry model called GEOS-Chem. This model has been designed to simulate atmospheric composition on a global and regional scale. Besides, it can be coupled with other 89 Chapter 5. Modeling the Influence of Cosmic Rays on the Atmosphere climatic or meteorological models, such as WRF (one of the most widely used in the world for short-term regional forecasting). GEOS-Chem is one of the most complete and accessible models that can be found, in addition to the fact that it is developed by hundreds of scientists around the world. Therefore, it is one of the most updated and complex models in its field. Another outstanding feature is that it has already implemented the effect of the ions generated by GCR in the nucleation process. Hence, it is one of the most suitable models to carry out our study. In the following sections, we present the methodology used to achieve our objectives and we will show that the effect of charging on the aerosol processes can be quite relevant and further research should be carried out in this direction. 5.2 Atmospheric Simulations with GEOS-Chem We use the global chemical transport model GEOS-Chem v12.1.0 (https://zenodo. org/record/1553349) with a horizontal resolution of 4◦×5◦and 47 vertical layers (from the surface up to 0.01 hPa). The model is driven by assimilated meteorological data from MERRA-2 reanalysis (https://gmao.gsfc.nasa.gov/reanalysis/MERRA-2/). This 3D model includes two aerosol microphysics schemes: TOMAS and APM. In this work, we use the TOMAS (TwO-Moment Aerosol Sectional) package [4, 160]. An advantage of using TOMAS is that it provides a higher resolution for all chemical species, especially for small sizes, which is very relevant to our study. This microphysics model simulates two independent moments (number and mass) of the aerosol size distribution for a number of discrete size bins: Nk=Zxk+1 xk nk(x)dx (5.1) Mk=Zxk+1 xk xnk(x)dx (5.2) where Nkand Mkare the total number and mass of aerosol in the kbin, nk(x)is the number of particles with masses between x+dx, and xkis the lowest limit of the kbin. In addition, the package incorporates modules for computing nucleation, condensation, and coagulation. We use the version of TOMAS40 that includes 30 bins logarithmically spaced ranging from 10 nm to 10 µm to represent the aerosol diameters, plus ten additional sub-10nm bins with a lower limit of 1 nm. The latter provides a high resolution for small particles, which is necessary for the simulations we want to carry out. Particularly, GEOS-Chem includes the ion-mediated nucleation (IMN) mechanism [180, 181] to calculate the nucleation rates taking into consideration the influence of atmospheric ions. The IMN depends on five key parameters: sulfuric acid concentration, temperature, relative humidity, ionization rate, and surface area of preexisting particles. The nucleation rates as a function of these parameters are extracted from a look-up table that covers a wide range of atmospheric conditions in order to be used as an input for GEOS-Chem. The global atmospheric ion rates due to CR are calculated following the model given by Usoskin and Kovaltsov [167]. The contribution of radioactive materials from soil to ionization rates is also included. The Usoskin and Kovaltsov approach consists of a very simple numerical model which computes the cosmic induced ionization in the atmosphere from the surface up to the stratosphere, all over the world. The model is parameterized by the modulation potential φ 90 IRMA RI´ ADIGOS S´ ANCHEZ (given in units of GV), which is used to easily calculate the variations in the induced ionization caused by the Sun’s activity. This parameter is utilized to determine the energy spectrum of GCR at the Earth’s orbit, which is modulated by the solar variations. It takes a typical value of 1 GV at solar maximum and 0.4 GV for the solar minimum. The smaller the value, the more CR enter the atmosphere. Figure 5.1: (a) Zonal-mean CR induced ionization (in units of ion-pairs cm−3s−1) in the atmosphere in the solar maximum (φ=1 GV). (b) Zonal-mean percent change in the induced ionization between the solar minimum (φ=0.4 GV) and solar maximum. [167] Figure 5.1a shows the atmospheric ionization rates for the solar-maximum case calculated using the method described in [167]. The rates are averaged over longitude in order to examine the regional differences. As can be seen, the ionization rates in the solar maximum are generally higher in the upper troposphere than near the surface. Furthermore, they are also higher towards the poles because magnetic rigidity is smaller, i.e., the Earth’s magnetic field shields much less. Figure 5.1b compares the ion-pair formation rates from CR between the solar minimum and the solar maximum. The changes between both situations reflect how the polar regions are most susceptible to CR changes. Apart from this, it should be remarked that the strongest FDs can cause changes in the atmospheric ionization rates similar to the changes between a solar maximum and solar minimum. Thus, a comparative study between solar peaks can be used to estimate how the changes would look like in a FD event [143]. 5.3 Approach to simulate charge distributions 5.3.1 Condensation Ions produced in the atmosphere by CR can charge aerosols through diffusion charging, which refers to the attachment of ions from the environmental background to the surface of particles [85, 174, 130]. We formulate the interactions governing the temporal dependences of ions and aerosols. The terms for the temporal changes in ion concentrations nqare: 91 Chapter 5. Modeling the Influence of Cosmic Rays on the Atmosphere Figure 5.4: (a) Zonal-mean nucleation rates in the X=0.8 simulation for the solar minimum. (b) Percentage change in nucleation rate between the solar-minimum and the solar-maximum simulations for X=0.8. Redish values refer to faster nucleation during the solar minimum case. size increases. The reason was already explained in previous works (see [143]). The Standard case is characterized by including only the effect of cosmic-rays ionization on the nucleation process, previously referred to as the IMN mechanism (see Section 5.2) and already included in GEOS-Chem. An increase in ionization causes more particles to have the capability to grow to larger sizes. The new small particles created compete for condensable material and start to grow more slowly, taking longer to reach 40 nm and 80 nm. Slower growth rates lead to an increase in the coagulation sink. As a consequence, the increase produced in the concentrations of smaller particles delays in reaching larger particles sizes. That is why CN80 and CCN show the least significant improvements. Surprisingly, this is not the case for X=0.8. Although there is a general decrease in the enhancement from CN3 to CCN, the magnitude of the change is maintained from 40 nm to CCN sizes. Furthermore, CCN concentrations display a slight increase when compared to CN80. For X=0.8, it is also noteworthy that whereas the change is twice as large as the Standard case for CN3, CN10, and CN40 concentrations, the change is more than three times greater for CN80 and CCN. One of the possible explanations is that diffusion charging is having a significant impact on small aerosols. For such a case, the enhancement or inhibition of small paticles with other size ranges is very dependent on the mobility value (as can be seen comparing the enhancement factors from Figure 5.2 and the ones included in Appendix B.4). As a result, it seems that fewer particles are lost by the coagulation sink and more can survive to CCN sizes. On the contrary, when the mobility ratio is set to 0.7 and 0.9 this effect disappears. In these cases, it could be that the coagulation is so inhibited that aerosols cannot grow by this pathway and compete for the condensable gases, slowing down their growth. In the Standard case, there is no major difference between the free troposphere and low troposphere changes. However, when including the charging effect in the simulations, it seems like the free troposphere region presents higher increases for smaller particles and reverses the tendency in CN80 and CCN concentrations. This difference is more pronounced in the case of X=0.8. The free troposphere has lower concentrations of CN40 and CN80 than the boundary layer and hence it might be more sensitive to variations in CN40 and CN80 concentrations. We can take a look at the nucleation rates to have a glimpse of what is happening. Figure 5.4a shows the zonal-mean nucleation rates for the solar minimum case when X= 98 IRMA RI´ ADIGOS S´ ANCHEZ Figure 5.5: (a) Zonal-mean nucleation rates in the standard case for the solar minimum. (b) Percentage change in nucleation rate between the solar-minimum and the solar-maximum simulations for the standard simulations. 0.8. As it can be seen, the highest nucleation values are found above 600 hPa in the mid-high latitude regions. At the same time, Figure 5.4b displays the percentage change in the nucleation rates between the solar minimum and the solar maximum. This percent change in the mid-high latitude regions is 5-15 %, with some peaks with even higher values. Furthermore, the overall percent change shows positive values in practically all regions. The contrast is striking when compared to the standard situation (Figure 5.5b). Here, the spatial distribution of the percent variations is quite different. On the one hand, almost all zonal locations show an increase of 1-6 % with respect to the solar maximum (i.e., higher cosmic-ray intensity). It should be noted that these values are in agreement with those reported in previous studies [181, 143]. On the other hand, the highest difference is seen in the mid-latitudes of the northern hemisphere where the nucleation rates have also the largest values (see Fig. 5.5a). The latter makes sense since the period of the simulation coincides with wintertime in that hemisphere. It is possible that the location of the increase regions might change with other periods. However, in the X=0.8 case, the nucleation rates show the same spatial distribution (Fig. 5.4a) but this does not result in an increase in the northern hemisphere only, rather it can be seen in both (Fig. 5.4b). This is evidence that the charged coagulation is taking effect. The smallest particles are less negatively charged than larger ones. Therefore, the coagulation between the smallest positive particles and the negatively charged larger bins is enhanced, leading to the reduction of the concentrations of the large particles. This could be the reason why in Figure 5.3 the CN80 and CCN concentrations show fewer percent changes in the free troposphere case (X=0.8) than in the lower troposphere. For illustration, Figure 5.6 shows the results when the mobility ratio is X=0.7. In such case, the percent differences display the same features as in the standard simulations. The X=0.9 case is not included because it displays similar nucleation values. Figure 5.7 shows the percent change in the zonal-mean CN values between the solar-minimum and the solar-maximum simulations when the mobility ratio is X=0.8 ( changes for X=0.9 can be found in Appendix B.3.2). In general, the changes in the concentrations (CN3, CN10, CN80, and CCN) increase during the solar cycle at nearly all zonal regions. Moreover, the biggest changes take place in the same mid-high latitude regions as in Figure 5.4b. Apart from this, changes in CN3 and CN10 are more pronounced in the 400-600 hPa region but the tropical upper troposphere also shows an increase, around 5 %. 99 Chapter 5. Modeling the Influence of Cosmic Rays on the Atmosphere Figure 5.6: (a) Zonal-mean nucleation rates in the X=0.7 case for the solar minimum. (b) Percentage change in nucleation rate between the solar-minimum and the solar-maximum simulations for X=0.7. Figure 5.7: Percentage change between the solar-minimum and the solar-maximum case of zonal-mean CN3, CN10, CN80, and CCN concentrations (X=0.8). 100 IRMA RI´ ADIGOS S´ ANCHEZ However, CN80 and CCN show a zonal sensitivity quite different because the largest variations are shifted to higher latitudes and appear in the lower troposphere. We observe that changes in the concentrations are a bit less than those found in the nucleation rates. In the mid-high latitudes, Figure 5.4b depicts changes in the nucleations between 10-20 %, whereas it slightly drops to 10-15 % in CN3. The situation for CN80 and CCN changes radically, but it should be noted that the average variations in CCN are 2 % with local changes as high as 10 %. Performing the analysis by regions, the CCN concentrations show the following increases: 4 % for the polar regions, 2.7 % for mid-latitude regions, and 0.82 % for the tropical areas. This would corroborate the theory that the highest changes in ionization rates at polar areas, as seen in Figure 5.1, would result in larger aerosol changes in those regions. So far, we have determined the response of CCN to changes in cosmic rays. We have just seen above that the changes in ion formation rates do not lead to similar changes in nucleation rates, and in turn changes in nucleation do not cause the same changes in CCN. The reason is that the processes involving the growth of aerosols are complex and compete with each other, giving rise to feedbacks that can enhance or dampen cosmic ray changes. For example, small particles can take two paths: coagulate with larger particles to reach CCN sizes or grow through condensation to form a new CCN. Previously, it has been proved that when only the effect of cosmic rays on nucleation was taken into account, the processes competed in such a way that the enhancement was depressed under certain conditions [143]. In this work, we are bringing more variables into play that make things more complicated. So, it is difficult to have a clear picture of what is happening but we have been able to draw a general description of the results that give us a lot of information. To illustrate this challenge, one can take a look at the correction factors such as those shown before in Figure 5.2, which are supposed to give information on what is happening with the coagulation of particles of different sizes. However, when compared to corrections factors for the other mobility cases (included in Fig. B.4 of the Appendix), one can appreciate that the differences are in the smallest particles, whose values change drastically from Wk,i<1 to Wk,i>1 depending on X. One possible way to understand which particle sizes are determining the fate of the aerosols would be to consider simulations with only the charge distribution on small particles or the opposite situation where only large particles carry charges. This remains pending as future work. Now, the question is whether changes in CCN that we have been reported in this work can lead to similar changes in cloud albedo or cloud cover. Some other dampened mechanisms may exist that cause CCN to be lost before cloud cover changes. However, it has been stated that in order to produce changes in cloud cover over the solar cycle, one needs changes higher than ∼1 % in CCN. In our case, these requirements are met when X=0.8. This opens up the possibility that the processes we have implemented may be the missing link between cosmic rays and clouds. However, we should be cautious with this statement, as there is still a lot of work to be done. For instance, future work will focus on increasing the period of the simulations to ensure that the observed changes hold over time and at different epochs of the year. Furthermore, we can observe that the value of the ion mobility ratio has a dramatic impact on the results. Here, we have assumed that it has a fixed value throughout the atmosphere, but the truth is that in real life this does not have to be the case. We know that Xhas values less than 1 in the atmosphere because negative ions have higher mobility than positive ions [84]. However, the sensitivity of the results to this ratio may be the reason behind the discrepancies 101 Chapter 5. Modeling the Influence of Cosmic Rays on the Atmosphere in the correlations between solar activity and clouds. 5.6 Conclusions This part of the thesis has addressed the problem of the influence of cosmic rays on the growth of atmospheric aerosols. Through the methodology previously proposed in other works [93, 168], we have been able to incorporate for the first time the effect of the diffusion charging in the microphysical evolution of atmospheric particles of a state-of-the-art atmospheric model, GEOS-Chem. Thus, the simulations performed with this model have provided a more realistic and detailed view of the indirect effect of cosmic rays on the final concentrations of CCN. The evidence from this study suggests that the ionization induced by cosmic radiation in the atmosphere may favor the growth of small particles to CCN sizes under certain conditions (X=0.8). We observed that changes in CCN concentrations between the solar maximum and solar minimum (2-10 %) may become significantly relevant for cloud formation. This work will serve as a base for future studies and further research is needed to estimate this effect more accurately. 102 GENERAL CONCLUSIONS The impetus for the work discussed in this thesis was to explore the implications of cosmic rays for atmospheric physics. In particular, we have used cosmic rays as a tool to analyze the properties of the atmosphere and showed how the appropriate technology has the potential to deliver the atmospheric profile with good accuracy. The second major finding to emerge from this dissertation is that cosmic rays may be more relevant than previously thought for aerosol growth and cloud formation. This observation was supported by the simulation results and agrees with the hypothesis posed at the beginning of the study. In general, the following conclusions can be drawn from the present work: • In the analysis of correlations between atmospheric variables and cosmic ray measurements, we have commissioned and calibrated a small-size 2 m2multigap timing RPC detector devoted to the detailed study of cosmic rays at ground level and performed the first analysis of the atmospheric temperature effect with this kind of technology. By studying a data sample of about one year, it has been possible to estimate the distribution of temperature coefficients (WT(h)), showing that the contribution of the hard component is dominant and in good agreement with theoretical expectations. We have seen how the presence of strong correlations among the different atmospheric layers precludes the use of conventional regression methods. A Principal Component Regression (PCR), considering the first two components, is sufficient to capture at least 77% of the temperature variability, giving a good description of the WT(h)and the global slope parameter αTexp =−0.279 ±0.051 %/K (compared to a theoretical value of αTtheor =−0.319 %/K). This results in an anticorrelation with the effective atmospheric temperature, which allows to clearly identify its seasonal cycles as well as short-term exceptional events (such as the tropospheric consequences of a Sudden Stratospheric Warming) through measurements performed at ground level. • We have developed a 1D-model Monte Carlo tool to simulate cosmic-ray air showers in real atmospheric scenarios. With this tool, we have been able to estimate the atmospheric effects, and probed them to be in qualitative agreement with a theoretical treatment based on weight coefficients WT(h)as well as data. This study has gone some way towards enhancing our understanding of the measured atmospheric effects by further inspecting the influence of the atmospheric attributes that would be impossible to accomplish with observational data. • Another study was undertaken to evaluate and establish the limits for obtaining the atmospheric temperature through cosmic-ray data. In that way, we have found the best configuration for a monitoring station integrated by cosmic-ray telescopes. Besides, the General Conclusions results have shown that it could be possible to retrieve the temperature of the atmosphere from the surface with good accuracy and at several atmospheric layers. An implication of this is the possibility that atmospheric muon detectors can be employed in scientific research beyond the field of astrophysics. Therefore, we have fulfilled one of the main objectives of this thesis which was to exploit the potential of cosmic rays for the development of practical applications for everyday life. Returning to the question posed at the beginning of the work, it is now possible to state that multi-directional telescopes can improve the estimate of atmospheric temperatures. Multiple analyses have revealed that one can achieve a high degree of accuracy with error margins between 0.8 and 2.2 K up to ∼20 km. It was also shown that it is possible to follow strong temperature variations that happen in the low stratosphere, like those taking place during Sudden Stratospheric Warmings. • With the two previous results, we pave the way for continuous monitoring of cosmic-ray variations using stations equipped with muon telescopes, that would allow real-time atmospheric temperatures to be retrieved for their use in meteorological monitoring or climate quality data records. This kind of station could be part of the global observing system along with the existing weather stations, satellites, and balloon measurements. In addition, it is important to highlight that this could be achieved with more affordable and accessible technology, and much easier to maintain than other alternatives. • In the study of the influence of cosmic rays on the atmosphere, we have found a possible link between cosmic rays and clouds. The enhancement that we have identified in the CCN concentrations assists in our understanding of the role of charged particles in aerosol growth. Previous works have focused on the study of processes affecting small aerosols, such as nucleation [181, 54, 149], but the results presented in this thesis seem to indicate that other processes such as charging coagulation are of great relevance. Furthermore, the results are in agreement with the theory and consistent with previous works, which increases their robustness. However, we should be aware that there is still room for improvement. It is still too early to claim that we have found the key to the missing link between cosmic rays and clouds. But it is certainly a step forward in understanding the relationship. Furthermore, whether or not this proves to be true in the future, it is undoubtedly a process that appears to have some degree of relevance to aerosol growth and should be accounted for in atmospheric models just as nucleation was once incorporated. Applicability and Future Perspectives The discoveries presented in this thesis have many important implications for future practice. On the one hand, the analysis of the temperature retrieval suggests that several courses of action can be pursued and there is still much room for improvement. For instance, the accuracy in temperature retrieval can be improved in practice by employing more advanced statistical techniques such as those employed in satellite remote sensing [59, 134]. Additionally, for 4 m2-scale detectors, the performance achieved in the analysis requires outstanding detector stability (below 0.3% on its counting efficiency). While detector inefficiencies to minimum 104 IRMA RI´ ADIGOS S´ ANCHEZ down to 0.1% are not alien to particle physics instrumentation, the requirement will pose significant constraints on the chosen technology and detector design. Regarding the part of the correction of the variations that are not of atmospheric origin, some experiments have stressed the relevance of complementary neutron detectors to remove those effects related to the primary cosmic ray fluctuations. However, we believe that this would not be necessary as the underground detector could be playing a similar role. Besides, the results of Chapter 2 show that it would not be very complicated to avoid such effects with some specific statistical methods. Anyway, we observe that there is a definite need for real measurements performed with a realistic setup. The installation of a functional station integrated by a couple of detectors, one ground-based and another underground, would be the ultimate test to see how realistic our proposal is. Furthermore, the temperature estimates obtained could be directly contrasted with balloon-sounding measurements, which would allow determining how accurate these estimates would be. In the case of the air shower simulations, a straightforward step would be to implement the changes in a complete simulation program, such as AIRES, which already has some built-in functionality to change the atmospheric model. In this way, more realistic simulations can be achieved, and the temperature effect on the soft component could be examined as well, for instance. Moreover, apart from estimating the distribution of the temperature coefficients experimentally or numerically as up to know, one could obtain them through Monte Carlo simulations, without requiring of any simplifying assumption. Finally, regarding the part of the GEOS-Chem simulations, a new world of possibilities opens up to investigate the consequences of the implementations performed. For example, the CLOUD experiment has shown that cosmic rays could be relevant in pre-industrial times. In particular, they have found that cosmic rays strongly enhance the production of pure biogenic particles by a factor of 10-100 compared to particles without the influence of ions, suggesting that cosmic radiation could have been more relevant for cloud formation in pre-industrial times than in today’s polluted atmosphere [94]. Other studies also argue that 1/3 of the warming in the last century was induced by changes in cosmic rays [35]. This could be related to the biogenic particle nucleation since the amount of atmospheric pollution was less than in the present atmosphere. Future work could be to run simulations for past periods with the anthropogenic emissions off and the charged coagulation in order to study the effect on CCN. Another alternative to test the effects of the charging on coagulation would be to simulate volcanic periods, where coarse aerosols are ejected into the atmosphere. This kind of simulation would serve as a sensitivity test for the implemented approach since it is assumed that for larger particles and high concentrations, the impact of charging can be significant. 105 Appendix A Atmospheric Dynamics Abstract: We include here some basic concepts about atmospheric dynamics. In particular, the Sudden Stratospheric Warming events are described more in detail in order to have a better understanding of what they are. A.1 Sudden Stratospheric Warming The stratosphere is the layer located immediately above the troposphere. The top of the stratosphere occurs around 50 km. As its name suggests, it is an atmospheric layer stratified into other layers, with the cooler ones located lower in the stratosphere. Indeed, the temperature profile of the stratosphere is characterized by the fact that the temperature increases with height, in contrast with the temperature of the troposphere, which decreases with altitude. The increase in temperature stems from the presence of ozone which absorbs ultraviolet radiation from the Sun. As a result of the temperature stratification of the whole layer, vertical mixing and convection are much rarer than horizontal mixing. Therefore, the layers of air are quite stable [140]. The troposphere is the layer of the atmosphere to which meteorologists pay most attention because it is where most weather phenomena take place. However, interactions between the troposphere and the stratosphere are also closely monitored, as they can deeply impact the weather down at the surface. One of these phenomena is the Sudden Stratospheric Warming (SSW), which has the ability to alter atmospheric patterns. The stratospheric polar vortex is a large and persistent low-pressure region located in the Poles. It weakens towards summer and gains intensity in winter. SSW events are very common and occur when the polar vortex starts to weaken in late winter. A SSW refers to a rapid and large warming in the stratosphere over a short period of time, usually a couple of days [31]. A normal stratosphere exhibits a polar vortex rotating counterclockwise (similar to cyclones) with very low temperatures. When a SSW takes place, the vortex gets weaker and can split in two or rearrange out of its usual position over the Pole. In the most extreme events, polar winds may even reverse and begin to rotate clockwise. A SSW is the result of a Rossby wave, also known as a planetary wave, lifting from the troposphere to the stratosphere, carrying a warmer air mass into the upper atmosphere. Atmospheric Rossby waves form primarily as a result of the land’s orography. For example, the great mountain systems of the Northern Hemisphere, such as the Himalayas or the Alps, can cause the dominant westerly winds in the mid-latitudes to ripple when they encounter these