UNIVERSIDAD NACIONAL AUTÓNOMA DE MÉXICO
POSGRADO EN ASTROFÍSICA
ASTROFÍSICA OBSERVACIONAL
STUDY OF THE STELLAR PHYSICS THROUGH INFRARED
INTERFEROMETRIC IMAGING
TESIS
QUE PARA OPTAR POR EL GRADO DE:
DOCTOR EN CIENCIAS (ASTROFÍSICA)
PRESENTA:
JESÚS ABEL ROSALES GUZMÁN
TUTOR PRINCIPAL
DR. JOEL SÁNCHEZ BERMÚDEZ
(INSTITUTO DE ASTRONOMÍA, UNAM - CU)
MIEMBROS DEL COMITÉ TUTOR
DRA. IRENE ANTONIA CRUZ GONZÁLEZ ESPINOSA
(INSTITUTO DE ASTRONOMÍA, UNAM - CU)
DR. JULIO CÉSAR RAMÍREZ VELEZ
(INSTITUTO DE ASTRONOMÍA, UNAM - ENSENADA)
Ciudad Universitaria, Ciudad de México, Abril de 2024
 
UNAM – Dirección General de Bibliotecas 
Tesis Digitales 
Restricciones de uso 
  
DERECHOS RESERVADOS © 
PROHIBIDA SU REPRODUCCIÓN TOTAL O PARCIAL 
  
Todo el material contenido en esta tesis esta protegido por la Ley Federal 
del Derecho de Autor (LFDA) de los Estados Unidos Mexicanos (México). 
El uso de imágenes, fragmentos de videos, y demás material que sea 
objeto de protección de los derechos de autor, será exclusivamente para 
fines educativos e informativos y deberá citar la fuente donde la obtuvo 
mencionando el autor o autores. Cualquier uso distinto como el lucro, 
reproducción, edición o modificación, será perseguido y sancionado por el 
respectivo titular de los Derechos de Autor. 
 
  
 
Agradecimientos
***
Sin duda alguna, el doctorado ha sido uno de los viajes más largos, dif́ıciles e, irónicamente,
gratificantes en mi vida, tanto personal como académica. Por ello, deseo expresar mi
profundo agradecimiento a todas aquellas personas que compartieron sus conocimientos
y experiencia conmigo en este camino.
En primer lugar quiero agradecer a mi asesor Joel Sánchez por haberme guiado a lo
largo de este proceso, por esas discusiones que me llevaron a entender las cosas que se me
complicaban y por alentarme siempre a desarrollar nuevas ideas para poner en práctica.
Agradezco también a los miembros de mi comité tutor, Irene Cruz y Julio Ramirez por
haberme mostrado siempre su apoyo en los momentos dif́ıciles y guiarme siempre hacia
adelante en cada reunión.
Agradezco a mis sinodales Julio Ramirez, Gisela Ort́ız, Jesús Toalá, Laurence Sabin
y Aina Palau por sus valiosos comentarios y correcciones que ayudaron a mejorar en
gran manera el contenido que se encuentra a lo largo de éstas páginas.
Finally, I would like to express my heartfelt gratitude to Claudia Paladini, whose
invaluable guidance has been instrumental throughout my doctoral journey. Not only
has she shared her wealth of knowledge with me, but she has also extended her generosity
by inviting me to collaborate with her in Chile. Thank you very much Claudia for all
your kindness.
Yo Jesús Abel Rosales Guzmán agradezco al proyecto CONAHCyT “Ciencia de Fron-
tera” CF-2019/263975 y al proyecto UNAM PAPIIT IA 105023 por el apoyo brindado
para la realización de esta tesis. Agradezco también a CONAHCYT por el apoyo fi-
nanciero otorgado a través de la beca doctoral con número 760678.
I
II
Resumen
***
¿Qué mecanismos dan lugar a la pérdida de masa y a la formación de estructuras
asimétricas de gas y polvo en estrellas AGB (Asymptotic Giant Branch, por sus siglas
en inglés)? ¿Qué tan importante es la composición qúımica del gas y polvo circundante
a las estrellas para la formación de vientos? ¿Cómo funciona la convección en estrellas
diferentes al Sol? Al d́ıa de hoy, éstas son preguntas que no tienen una respuesta conc-
reta. Aunque las estrellas AGB son de las principales contribuyentes al enriquecimiento
del medio interestelar, la manera en la que el material procesado durante la evolución de
la estrella es devuelto al medio interestelar es un fenómeno que no se encuentra caracter-
izado completamente. En particular las estrellas AGB de tipo M (o ricas en O, Caṕıtulo
1) son de las menos entendidas debido a la enorme complejidad qúımica existente en las
partes más cercanas a la estrella, donde ocurren las primeras interacciones que dan lugar
a la subsecuente pérdida de masa.
Una de las dificultades encontradas al tratar de observar las regiones más internas
de estrellas AGB es la alta resolución angular necesaria para su estudio. Por ejemplo,
la envolvente circunestelar de polvo de una estrella AGB, observada a 11 µm, tiene un
tamaño angular aproximado de 50 milisegundos de arco (mas) equivalentes a 50 UA a una
distancia de 1 kilopárcsec (kpc). Para poder resolver un objeto con este tamaño angular
necesitaŕıamos un telescopio de al menos 45 m de diámetro. Si nos vamos a regiones
todav́ıa más internas con tamaños de unos pocos mas, necesitaŕıamos telescopios cada
vez más grandes, con un costo que se dispara de manera exponencial. Es entonces que
el uso de interferómetros (Caṕıtulo 2) adquiere una enorme importancia para estudiar
estructuras con tamaños angulares tan pequeños como algunos mas con el fin de brindar
pistas que nos ayuden a entender que ocurre en el intrincado y complejo entorno de las
estrellas AGB.
III
IV
Esta tesis se centra en el estudio de estrellas AGB con interferometŕıa infrarroja. De
manera más concreta, presento un análisis de las partes más internas (1–3 R∗) de la
estrella Mira R Car utilizando datos interferométricos adquiridos con los instrumentos
PIONIER y GRAVITY del VLTI (Very Large Telescope Interferometer, por sus siglas
en inglés) en las bandas H (1.6 - 1.99 µm) y K (2.0 - 2.4 µm) en el infrarrojo cercano.
El trabajo de tesis está dividido en 5 caṕıtulos que abarcan tres partes principales.
Los primeros dos caṕıtulos ofrecen una introducción a las estrellas AGB y a la interfer-
ometŕıa infrarroja, respectivamente. En el Caṕıtulo 3 titulado Interferometric analy-
sis of the M-type Mira star R Car with PIONIER-VLTI exploramos la superficie
estelar de R Car en la banda H en dos épocas diferentes separadas por alrededor de 7
años. La resolución angular alcanzada nos permitió detectar una serie de estructuras
brillantes en la superficie de la estrella que cambian de posición y tamaño en el tiempo.
Fuimos capaces de determinar el tamaño de estas estructuras y realizamos compara-
ciones de tamaño con extrapolaciones de diferentes modelos de convección que fueron
desarrollados para estrellas gigantes, enanas y enanas blancas. Aunque a la fecha no ex-
isten modelos espećıficos de convección para cada AGB, no es del todo incorrecto pensar
que los procesos que gobiernan la convección deben ser similares en todas las estrellas
a lo largo del diagrama Hertzsprung–Russell (HR). Es por ello que pudimos comparar
nuestros resultados con extrapolaciones de modelos diseñados para otro tipo de estrel-
las. Derivado de estas comparaciones encontramos que las estructuras observadas en
la superficie de R Car son consistentes en tamaño con las extrapolaciones usadas, lo
cual sugiere que tales estructuras brillantes pueden estar asociadas a células convecti-
vas. Poder observar y caracterizar tales estructuras estelares y sus cambios temporales
nos abre la puerta para brindar nueva y mejor evidencia observacional que nos per-
mita entender mejor los modelos actuales de convección, los cuales, a pesar de ser muy
abundantes, carecen de un soporte observacional en el cual basarse.
El mecanismo más aceptado para la pérdida de masa en la fase de AGB es un
proceso de 2 pasos, (i) las pulsaciones de la estrella elevan el gas a una distancia a
partir de la cual el polvo puede comenzar a condensarse y (ii) la radiación de la estrella
acelera al gas y polvo hasta escapar del potencial gravitacional. Aunque este esquema
funciona, algo que no se entiende del todo bien es ¿Qué tipos de granos de polvo son
eficientes para transportar dicha radiación, ya sea mediante absorción o dispersión? Para
arrojar luz sobre este problema es indispensable poder observar las capas de gas más
V
cercanas a la superficie de la estrella, por lo que, además de una buena resolución angular,
necesitamos una ventana donde el polvo sea lo más transparente posible. En las estrellas
AGB de tipo M se cumple el segundo requisito observando en la banda K, razón por la
cual en el Caṕıtulo 4 titulado Interferometric Analysis of the Mira star R Car
with GRAVITY-VLTI llevamos a cabo la caracterización de las capas gaseosas más
próximas a la superficie estelar en esta banda en dos épocas diferentes separadas por
alrededor de 1 mes.
La importancia de este estudio radica en que las condiciones f́ısicas de temperatura
y densidad en las que se encuentra el gas, principalmente el CO, son compatibles con las
regiones donde se origina el polvo que da lugar a los vientos y a la pérdida de masa en
este tipo de estrellas. En el análisis de estos datos encontramos que la distribución de las
capas de gas que rodean a la estrella es asimétrica y cambia de forma en periodos cortos
de tiempo. Al ser capaces de estimar la temperatura, tamaño y densidad de las capas
de gas más internas, pudimos determinar de manera observacional si los tipos de granos
de polvo propuestos como ”wind-drivers” pueden o no ser eficientes para transportar la
radiación de la estrella.
Una manera de comprobar si las afirmaciones hechas en el Caṕıtulo 4 sobre los
posibles ”wind-drivers” que se encuentran distribuidos en diferentes capas alrededor de
la estrella es mediante la observación y caracterización de sus transiciones moleculares.
Dado que éstas no pueden ser detectadas en las bandas H y K, para comprobar su
existencia en las cercańıas de la estrella seŕıa necesario observar en las bandas L (3.2
- 3.9 µm), M (4.5 - 5.0 µm) y N (8 - 13 µm) principalmente, ya que es ah́ı donde se
espera que ocurran tales transiciones. En el Caṕıtulo 5 se plantea como trabajo a futuro
la caracterización de las capas más externas de gas y polvo, observadas en las bandas
mencionadas previamente.
Esta tesis tiene como objetivo cient́ıfico inmediato brindar, por primera vez, un es-
tudio detallado desde adentro hacia afuera de las diferentes regiones de una AGB rica
en O para brindar evidencias observacionales y mejorar los modelos actuales de con-
vección y de pérdida de masa, además de motivar los estudios observacionales periodicos
y detallados con técnicas de alta resolución angular.
VI
Producción Cient́ıfica
***
Art́ıculos de primer autor
Mi trabajo doctoral tuvo como resultado la escritura de dos art́ıculos originales, los
cuales listo a continuación. El primero de ellos ya se encuentra publicado y el segundo
está en proceso de revisión por pares. Ambos fueron enviados a la revista Astronomy &
Astrophysics.
1. Imaging the innermost gaseous layers of the M-type Mira star R Car
with GRAVITY-VLTI
A. Rosales-Guzmán, J. Sanchez-Bermudez, C. Paladini, A. Alberdi, W. Brandner, E.
Cannon, G. González-Torà, X. Haubois, Th. Henning, P. Kervella, M. Montarges, G.
Perrin, R. Schödel, M. Wittkowski
Publicado en 2023 (A&A 674, A62)
2. A new dimension in the variability of AGB stars: convection patterns
size changes with pulsation
A. Rosales-Guzmán, J. Sanchez-Bermudez, C. Paladini, B. Freytag, M. Wittkowski,
A. Alberdi, F. Baron, J.-P. Berger, A. Chiavassa, S. Höfner, A. Jorissen, P. Kervella,
J.-B. Le Bouquin, P. Marigo, M. Montargès, M. Trabucchi, S. Tsvetkova, R. Schödel, S.
Van Eck
Enviado a A&A y en proceso de revisión por pares.
VII
VIII
Art́ıculos en colaboración
Además de los art́ıculos previamente mencionados, durante mi doctorado tuve la opor-
tunidad de colaborar en dos art́ıculos y un proceeding. Estas contribuciones las listo a
continuación junto con una descripción de mi participación en los proyectos.
1. The dusty circumstellar environment of Betelgeuse during the Great
Dimming as seen by VLTI/MATISSE.
E. Cannon, M. Montargès, A. de Koter, A. Matter, J. Sanchez-Bermudez, R. Norris, C.
Paladini, L. Decin, H. Sana, J.O. Sundqvist, E. Lagadec, P. Kervella, A. Chiavassa, A.
K. Dupree, G. Perrin, P. Scicluna, P. Stee, S. Kraus, W. Danchi, B. Lopez, F. Millour,
J. Drevon, P. Cruzalèbes, P. Berio, S. Robbe-Dubois, and A. Rosales-Guzmán
Publicado en 2023 (A&A, 675, A46)
En este art́ıculo, liderado por Emily Cannon, se da una explicación a la disminución
del brillo de la estrella gigante roja Betelgeuse, ocurrida a finales de 2019 e inicios
de 2020, utilizando datos interferométricos en la banda N del infrarrojo medio. Mi
participación en este proyecto consistió en revisar si era posible reconstruir una imágen
de la estrella a partir de los datos interferométricos. Dada la calidad de estos datos no
fue posible reconstruir la imagen pero se dieron condiciones para lograr esto con futuras
observaciones.
2. An impressionist view of V Hydrae.
L. Planquart, C. Paladini, A. Jorissen, A. Escorza, E. Pantin, J. Drevon, B. Aringer,
F. Baron, A. Chiavassa, P. Cruzalèbes, W. Danchi, E. De Beck, M.A.T. Groenewegen,
S. , Höfner , J. Hron, T. Khouri, B. Lopez, F. Lykou, M. Montarges, N. Nardetto, K.
Ohnaka, H. Olofsson, G. Rau, A. Rosales-Guzmán, J. Sanchez-Bermudez, P. Scicluna,
L. Siess, F. Th’evenin, S. Van Eck, W.H.T. Vlemmings, G. Weigelt, M. Wittkowski
En segunda revisión por pares.
En este art́ıculo, liderado por Lèa Planquart, se explora la posibilidad de que la
estrella AGB V Hydrae tenga una compañera binaria en su región de formación de
polvo, a menos de 10 R∗. Se exploran además los posibles efectos de esta binariedad en
IX
la pérdida de masa y en la forma de la envolvente circunestelar de V Hydrae. Para llevar
acabo los objetivos mencionados, se analizaron datos interferométricos de V Hydrae
en las bandas L, M y N en el infrarrojo medio. Estos datos fueron tomados con el
instrumento MATISSE del VLTI (Ver Cap. 2 para una explicación del instrumento).
Como coautor de este art́ıculo, he tenido la oprtunidad de revisar el draft y las siguientes
revisiones del referee y darle a Léa sugerencias para mejorar la discusión.
Proceedings
1. Compressed sensing for infrared interferometric imaging
Joel Sanchez-Bermudez, A. Rosales-Guzmán, Héctor Morales, Antxon Alberdi, Anand
Sivaramakrishnan, Rainer Schödel, and Jörg-Uwe Pott
Publicado en 2020 (Proceedings Volume 11446, Optical and Infrared In-
terferometry and Imaging VII; 114461O)
En este proceeding, liderado por Joel Sánchez, se explora la técnica de compressed
sensing para la reconstrucción de imágenes utilizando datos interferométricos tomados
con el telescopio James Webb. Como coautor de este trabajo, obtuve las observables
interferométricas y realicé las gráficas que se presentan a lo largo del escrito.
X
Contents
***
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XVII
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XIX
List of Acronyms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XXIII
Motivation and outlook of this thesis 1
1 AGB stars 5
1.1 Evolution from the MS to the early-AGB phase at M = 1 M⊙ . . . . . . . 6
1.2 Evolution in the AGB phase . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.2.1 Stellar Variability . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.3 Mass loss during the AGB phase . . . . . . . . . . . . . . . . . . . . . . . 15
1.3.1 Early indications of mass loss . . . . . . . . . . . . . . . . . . . . . 15
1.3.2 Mass loss from a theoretical point of view . . . . . . . . . . . . . . 16
1.3.3 The role of convection in the mass-loss processes . . . . . . . . . . 21
2 Interferometry 31
2.1 Principles of interferometry . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.1.1 Visibility amplitude . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.1.2 Visibility phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.1.3 The effect of the atmosphere on the interferometric observations . 35
2.2 Interferometric facilities around the World . . . . . . . . . . . . . . . . . . 39
2.2.1 CHARA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
2.2.2 VLTI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
2.3 Image reconstruction in optical interferometry . . . . . . . . . . . . . . . . 45
2.3.1 BSMEM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
XI
XII CONTENTS
2.3.2 SQUEEZE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
2.4 Studies of AGB stars through IR interferometric observations . . . . . . . 49
2.4.1 Model fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
2.4.2 Image reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . 50
2.4.3 Model Atmospheres . . . . . . . . . . . . . . . . . . . . . . . . . . 50
2.4.4 Radiative transfer modeling . . . . . . . . . . . . . . . . . . . . . . 51
2.4.5 The importance of observational evidence . . . . . . . . . . . . . . 51
3 Interferometric analysis of the M-type Mira star R Car with PIONIER-
VLTI 57
3.1 Our source: R Car . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.2 Introduction to the problem of convection in AGB stars . . . . . . . . . . 59
3.3 Observations and data reduction . . . . . . . . . . . . . . . . . . . . . . . 60
3.4 Analysis and results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
3.4.1 Diameter estimation: parametric uniform-disk model . . . . . . . 61
3.4.2 Asymmetries in the stellar disk: regularized image reconstruction . 64
3.4.3 Power spectrum analysis . . . . . . . . . . . . . . . . . . . . . . . . 65
3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
3.5.1 Sizes of the stellar disk and of the surface structures . . . . . . . . 68
3.5.2 Convective pattern size(s) and their correlation with the stellar
parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
3.6 Appendix: Additional figures . . . . . . . . . . . . . . . . . . . . . . . . . 72
4 Interferometric analysis of the M-type Mira star R Car with GRAVITY-
VLTI 75
4.1 The problem of dust and wind formation in AGB stars . . . . . . . . . . . 75
4.2 Observations and data reduction . . . . . . . . . . . . . . . . . . . . . . . 78
4.2.1 GRAVITY observations and interferometric data reduction . . . . 78
4.2.2 GRAVITY spectrum calibration . . . . . . . . . . . . . . . . . . . 78
4.3 Analysis and results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.3.1 Geometrical model . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.3.2 Spherically thin layer: MOLsphere model . . . . . . . . . . . . . . 85
4.3.3 Image reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . 88
4.3.4 Pure-line CO images . . . . . . . . . . . . . . . . . . . . . . . . . . 96
CONTENTS XIII
4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
4.4.1 CO column densities from the single-layer model . . . . . . . . . . 98
4.4.2 Observational restrictions to the wind driving candidates for M-
type Mira stars . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
4.4.3 The role of CO asymmetries on mass-loss estimates . . . . . . . . . 102
4.5 Appendix: Additional figures . . . . . . . . . . . . . . . . . . . . . . . . . 103
5 Conclusions and future work 111
5.1 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
5.2 Prospects for the following decade: the ngVLA . . . . . . . . . . . . . . . 119
5.3 Closing remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
A Python codes 123
XIV CONTENTS
List of Figures
***
1.1 Overview of an AGB star . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2 Evolution of a 1 M⊙ star . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3 Examples of AGB stars with their CSE . . . . . . . . . . . . . . . . . . . 10
1.4 The light curve of o Ceti . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.5 Curves of constant condensation distance as a function of the power law
coefficient and the value of the condensation temperature . . . . . . . . . 20
1.6 Best-reconstructed images of the stellar surface of the S-type AGB star
π1 Gru . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
1.7 Characteristic granule size of the convection cells in π1 Gru . . . . . . . . 24
1.8 An illustration of the movement of a convective blob across the medium . 27
2.1 Fringe pattern on a detector’s screen . . . . . . . . . . . . . . . . . . . . . 33
2.2 Example of a u−v plane achievable with GRAVITY using three Telescope
arrays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.3 Comparison between single-dish and interferometric observables for two
objects with different sizes. . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.4 Simulated images of three models with their corresponding visibility curves 36
2.5 Visibility curves of two binary systems . . . . . . . . . . . . . . . . . . . . 37
2.6 An illustration to define the closure phases . . . . . . . . . . . . . . . . . 38
2.7 Footprints of modern O/IR interferometers . . . . . . . . . . . . . . . . . 39
2.8 The CHARA array and its instruments . . . . . . . . . . . . . . . . . . . 41
2.9 The Paranal Observatory . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.10 The three second-generation instruments at the VLTI . . . . . . . . . . . 46
XV
XVI LIST OF FIGURES
2.11 Reconstruction of the interferometric imaging contest data set of 2004. . . 50
2.12 Overview of the number of AGB stars studied with interferometric high-
resolution observations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.1 Visual light curves of R Car for our 2014 and 2020 data . . . . . . . . . . 64
3.2 PIONIER u-v coverages for our 2014 and 2020 data . . . . . . . . . . . . 65
3.3 Best reconstructed images for our 2014 and 2020 datasets . . . . . . . . . 66
3.4 Radial averaged power spectra of the reconstructed images . . . . . . . . 67
3.5 Logarithmic scale of the convective patterns (xg) versus effective temper-
ature (Teff) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
3.6 Comparison between the observables extracted from the 2014 mean re-
constructed image and the data . . . . . . . . . . . . . . . . . . . . . . . . 73
3.7 Comparison between the observables extracted from the 2020 mean re-
constructed image and the data . . . . . . . . . . . . . . . . . . . . . . . . 73
4.1 Visual light curve of R Car for our 2018 data . . . . . . . . . . . . . . . . 80
4.2 R Car calibrated spectra and V2 . . . . . . . . . . . . . . . . . . . . . . . 81
4.3 GRAVITY calibrated interferometric observables . . . . . . . . . . . . . . 83
4.4 Best-fit UD models for the V2 in the pseudo-continuum versus spatial
frequencies for both epochs. . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.5 Best-fit ΘUD as a function of wavelength for the two analyzed epochs . . . 85
4.6 Illustration of single-layer model . . . . . . . . . . . . . . . . . . . . . . . 86
4.7 Best fit to spectra with the single-layer model for the two CO band heads
and the two epochs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.8 Best fit to V2 from the best single-layer models across the first CO band
head for our two epochs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
4.9 Best fit to V2 from the best single-layer models across the second CO
band head for our two epochs. . . . . . . . . . . . . . . . . . . . . . . . . . 89
4.10 Best pseudo-continuum reconstructed images for the January and Febru-
ary epochs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
4.11 Best reconstructed images per spectral channel for the January and Febru-
ary epochs across the first CO band head . . . . . . . . . . . . . . . . . . 93
4.12 Best reconstructed images per spectral channel for the January and Febru-
ary epochs across the second CO band head. . . . . . . . . . . . . . . . . 94
LIST OF FIGURES XVII
4.17 Pure-line CO images of the first CO band head for the January and Febru-
ary epochs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
4.13 Principal components at the first CO band head for the January and
February epochs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
4.14 Principal components for the January and February epochs at the second
CO band head . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
4.15 Representative example of CPs versus spatial frequencies of the data set
that corresponds to the first CO band head . . . . . . . . . . . . . . . . . 106
4.16 Two-dimensional photocenter offsets in mas . . . . . . . . . . . . . . . . . 106
4.18 Condensation temperature as a function of the dust absorption coefficient
for a constant levitation radius Rl = 2.43 R∗ . . . . . . . . . . . . . . . . 107
4.19 Condensation temperature as a function of the dust absorption coefficient
for a constant levitation radius obtained directly from the pure-line CO
images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
4.20 Fit to observed V2 and CPs corresponding to the best reconstructed im-
ages across the pseudo-continuum for the January and February epochs . 108
4.21 Fit to observed V2 and CPs from the best reconstructed images across
the first CO band head for the January and February epochs . . . . . . . 109
4.22 Fit to observed V2 and CPs from the best reconstructed images across
the second CO band head for the January and February epochs. . . . . . 110
5.1 ISO-SWS spectra of three AGB stars . . . . . . . . . . . . . . . . . . . . . 115
5.2 The typical spectrum of an AGB star in the infrared from 0.5 to 25 µm . 116
5.3 Calibrated V2 and CPs of R Car in the L−band (3.0 - 4.0 µm) as a
function of the spatial frequency . . . . . . . . . . . . . . . . . . . . . . . 118
5.4 Calibrated V2 and CPs of R Car in the M band (4.1 - 4.7 µm) as a
function of the spatial frequency . . . . . . . . . . . . . . . . . . . . . . . 118
5.5 Reconstruction tests performed on simulated data to test the capabilities
of the ngVLA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
5.6 uv-coverage of the simulated ngVLA observations . . . . . . . . . . . . . . 121
XVIII LIST OF FIGURES
List of Tables
***
1.1 Early Discoveries of Intrinsic Stellar Variables . . . . . . . . . . . . . . . . 13
2.1 Angular resolution of telescopes with different diameters for three different
band-passes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.2 High-resolution observations of AGB stars . . . . . . . . . . . . . . . . . . 54
2.3 High-resolution observations of AGB stars (Continued) . . . . . . . . . . . 55
3.1 Fundamental astrophyscal parameters of R Car . . . . . . . . . . . . . . . 59
3.2 Log of the observations for the 2014 interferometric data . . . . . . . . . . 62
3.3 Log of the observations for the 2020 interferometric data . . . . . . . . . . 63
3.4 Logarithm of the characteristic size of the patterns observed on the surface
R Car. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.1 Log of the observations for our GRAVITY R Car data . . . . . . . . . . . 79
4.2 Best-fit parameters of the uniform disk geometrical model. . . . . . . . . . 82
4.3 Parameters of the single-layer model for the two CO band heads and our
two epochs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
4.4 Parameters of the single-layer model for other AGBs. . . . . . . . . . . . . 100
XIX
XX LIST OF TABLES
List of Acronyms
***
AAVSO American Association of Variable Star Observers
AFGL Air Force Geophysics Laboratory
AGB asymptotic giant branch
AGN active galactic nuclei
ALMA Atacama Large Millimeter/submillimeter Array
ATs Auxiliary Telescopes
AU Astronomical Unit
BCD Beam Commuting Device
BSMEM BiSpectrum Maximum Entropy Method
CEV Cumulative Explained Variance
CHARA Center for High Angular Resolution Astronomy
CP Closure Phase
CSE circumstellar envelope
DIT Detector Integration Time
E-ELT European Extremely Large Telescope
XXI
XXII List of Acronyms
ESO European Southern Observatory
FITS Flexible Image Transport System
FWHM Full Width Half Maximum
HBB hot-bottom-burning
HR Hertzsprung–Russell
IR infrared
IRAM Institut de Radioastronomie Millimétrique
ISM Interstellar medium
JMMC Jean-Marie Mariotti Center
LPV Long Period Variables
mas milliarcseconds
MATISSE Multi AperTure mid-Infrared SpectroScopic Experiment
MCMC Markov Chain Monte Carlo
MIR mid infrared
MLT Mixing length theory
MS Main Sequence
NDIT Number of Detector InTegrations
ngVLA next-generation Very Large Array
NIR near infrared
OPD optical path difference
PIONIER Precision Integrated Optics Near Infrared ExpeRiment
List of Acronyms XXIII
PNe Planetary nebulae
RGB red giant branch
RSG red super giant
SED Spectral Energy Distribution
TMSS Two-Micron Sky Survey
UTs Unit Telescopes
VLA Karl G. Jansky Very Large Array
VLTI Very Large Telescope Interferometer
WDs white dwarfs
Motivation and outlook of this
thesis
***
Although asymptotic giant branch stars are among the main contributors to the enrich-
ment of the interstellar medium, our understanding of the processes that lead to the
formation of winds and the subsequent mass loss is still limited. While current theoreti-
cal models have provided us with many insights of processes such as convection or wind
formation in these stars, observational evidence to test these models is still scarce.
One of the main reasons for this is the angular resolution required to characterize
phenomena such as those mentioned above, which occur in the regions closest to the
stellar surface. This is where interferometry becomes a crucial tool to provide the obser-
vational evidence needed to improve our understanding of these astrophysical sources.
We will see in the following chapters that interferometric studies of asymptotic giant
branch stars are not very abundant, but each study offers different results that allow us
to increasingly understand the processes taking place in these objects with better detail.
Motivated by the lack of observational evidence to confront the abundant mass-loss
and convection models, this thesis focuses on the study of asymptotic giant branch stars
using infrared interferometry. More specifically, I present a detailed analysis of the
innermost regions (1–3 R∗) of the M-type Mira star R Car using interferometric data
obtained with the Precision Integrated Optics Near Infrared ExpeRiment and GRAVITY
instruments of the Very Large Telescope Interferometer in the H (1.6 - 1.99 µm) and
K (2.0 - 2.4 µm) bands in the near infrared. The analysis of these data cover different
wavelengths (multi-band analysis) and also different phases of the pulsation cycle of the
star (multi-epoch analysis).
1
2 List of Acronyms
This thesis is structured as follows: In Chapter 1, I present an introduction about
asymptotic giant branch stars in general and then I move to more specific topics: con-
vection and the first interactions that lead to the subsequent mass loss in these types of
stars. In Chapter 2 I introduce the basic concepts of infrared interferometry; the most
common techniques used to interpret interferometric data; and the facilities around the
world that can provide these kind of observations. Both chapters are aimed also to em-
phasize the need for high-angular resolution observations to shed light into the complex
phenomena occurring in the innermost regions of asymptotic giant branch stars.
In Chapter 3, I introduce the M-type Mira star R Car, which is the subject of study
in my thesis, then, I present an analysis of its surface using interferometric data in
the H band, obtained with the instrument Precision Integrated Optics Near Infrared
ExpeRiment of the Very Large Telescope Interferometer. The quality of these data and
the angular resolution achieved, allowed us to detect bright structures on the surface
of the star, whose size and position change over time. We compared these results with
extrapolations from various convection models developed for different types of stars.
Although specific convection models for asymptotic giant branch stars are absent, it is
not entirely incorrect to think that the processes governing convection should be similar
across all stars along the Hertzsprung–Russell diagram. Our results suggest that the
observed bright structures on the surface of R Car are consistent with extrapolated sizes,
potentially associated with convective cells. This observational evidence enhances our
understanding of convection models, which, again, currently lack observational support
despite their abundance.
The most widely accepted mechanism for mass loss in the asymptotic giant branch
phase involves a two-step process: (i) stellar pulsations elevate gas to a distance where
dust can begin to condense, and (ii) stellar radiation accelerates the gas and dust, allow-
ing them to escape the gravitational potential. While this scheme works, an aspect not
fully understood is the efficiency of different dust grains in transporting radiation, either
through absorption or scattering. Shedding light on this issue requires the observation
of gas layers closest to the stellar surface. In addition to high angular resolution, a
spectral window with dust transparency is essential. This requirement is met in M-type
asymptotic giant branch stars by observing in the K band. Hence, in Chapter 4 we
characterize the gas layers closest to the stellar surface of R Car in this band during two
separate epochs spaced approximately 1 month apart.
List of Acronyms 3
The importance of this study lies in the fact that the physical conditions of tempera-
ture and density in which the gas, primarily CO, is found, are consistent with the regions
where the dust originates, giving rise to winds and mass loss in such stars. With the
analysis of these data, we found that the distribution of the gaseous layers surrounding
the star is asymmetric and undergoes shape variations over short periods. By success-
fully estimating the temperature, size, and density of the innermost gas layers, we were
able to determine some types of dust grains as wind-drivers for M-type asymptotic giant
branch stars.
Finally, in Chapter 5, I present the conclusions and future work, which includes the
ongoing analysis of interferometric data for R Car in the L, M , and N bands in the
mid infrared. This study aims to detect and characterize the innermost dusty layers
around the star. I close this chapter presenting the future of interferometry for the
next decade with the construction of the next-generation Very Large Array and the
subsequent enhancements to the Very Large Telescope Interferometer.
4 List of Acronyms
1C
h
ap
te
r
AGB stars
***
This chapter gives an introduction to the asymptotic giant branch (AGB) phase, the
different types of AGB stars, and their different mass-loss mechanisms. The contents
in this chapter are based principally on the book by Habing & Olofsson (2013), on the
review by Höfner & Olofsson (2018), and the article by Bladh & Höfner (2012). Other
references will be cited throughout the chapter.
The AGB phase constitutes one of the the final stages in the life cycle of stars
characterized by low-to-intermediate masses within the range of 0.5 - 8 M⊙. Following
their departure from the Main Sequence (MS), these stars undergo an initial red giant
branch (RGB) phase, transitioning from the subsequent central He -burning stage to
harbor an inert C-O core where eventually the electron-degeneracy pressure dominates.
At this stage, nuclear burning within the stellar core ceases to occur. Subsequently,
the stars move through a second RSG phase, denoted as the AGB phase, sustained by
the alternating burning of H and He in thin shells around the inert C-O core. This
marks the concluding phase in their evolutionary trajectory before they transform into
Planetary nebulae (PNe) and then white dwarfs (WDs).
An AGB star comprises distinct components: (i) a compact, very hot, and dense
core, (ii) a large, hot, and comparatively less dense stellar envelope, (iii) a rarified, warm
atmosphere, and (iv) a very large, extremely diluted, and cool circumstellar envelope
(CSE) (see Fig. 1.1). So, the size scales span over more than 10 orders of magnitude,
the density encompasses 30 orders of magnitude, and the temperature spans 7 orders of
magnitude. To fully grasp the AGB star’s evolution, it is crucial to explore the intricate
dynamics of the diverse physical and chemical processes unfolding in its various segments.
The following sections offer an exploration of the processes taking place as a 1 M⊙ star
5
6 CHAPTER 1. AGB STARS
embarks from its MS to the AGB phase. This discussion aims to elucidate fundamental
attributes of such stars, shedding light on their preceding evolutionary stages.
Figure 1.1: Overview of an AGB star (not drawn to scale). Some important physical/-
chemical processes are indicated. Image taken from Höfner & Olofsson (2018).
1.1 Evolution from the MS to the early-AGB phase at M
= 1 M⊙
The important features and schematic evolution of a 1 M⊙ star are presented in Fig.
1.2. For an easy identification in the figure, those features are indicated with numbers
from 1 to 15. I will refer to these points along the section.
The initial stage in the life of a newly formed 1 M⊙ star is characterized by the onset
of nuclear reactions in the core, marking the main sequence. During this phase, the only
source of stellar energy is the burning of H in the core. The star’s structure undergoes
modifications primarily due to alterations in its chemical composition resulting from
nuclear reactions. This phase has a duration of 109 years, constituting the longest part
of the star’s life and represents a stable period where the star maintains a balance
between the gravitational force pulling it inward and the outward pressure generated by
H fusion in the core.
As the amount of H in the core decreases, the star will slowly move upwards in
the HR diagram, almost along the main sequence (points 1-4). The corresponding H
abundance profiles are shown in inset (a) of Fig. 1.2. The star will experiment slight
increases in brightness and temperature, with a relatively stable radius. Then, it will
1.1. EVOLUTION FROM THE MS TO THE EARLY-AGB PHASE AT M = 1 M⊙ 7
move to the right in the HR diagram, as H in the core nears its end. Eventually, the
core is almost pure degenerated He . H will continue to burn in a thick shell around the
core. Inset (b) shows the advance of the H-burning shell during this evolution.
After the H exhaustion at the core, the star transitions from the main sequence to
a state in which H starts burning in a shell surrounding the He-core, causing the core’s
mass to increase. As a consequence, the envelope of the star expands, moving it to the
right in the HR diagram (points 5-7). Upon reaching the Hayashi limit (about point
7), convection extends inward (in mass) from the surface, and the convective envelope
penetrates into the region where partial H burning has occurred in earlier evolution (see
inset (c)). This phenomenon is known as first dredge-up and it enriches the surface with
4He and products of CN (primarily 14N and 13C), which are mixed to the surface (point
8). The growing core eventually becomes degenerate, leading to a rise in temperature
as the star ascends the giant branch until it reaches 100 million K, which is enough for
He to burn into C in the triple-alpha process. Inset (d) shows how neutrino losses from
the center cause the temperature maximum to move outward.
The He burning begins across the central region, causing a sudden temperature
increase. Since the degenerate core cannot expand as a response to the temperature
increase, the rate of nuclear reactions accelerates instead. At some point, the increasing
temperature will remove the degeneracy of the gas in the core, and a violent expansion
will occur. Shortly after the He ignition, an explosion occurs: the He-core flash or
He-flash (point 9).
Post-He flash, the star moves to the horizontal branch, burning 4He in a convective
core and H in a shell (10-13). The position on the horizontal branch depends on the
envelope mass, which is influenced by mass loss during the He-flash. While the luminosity
does not vary much along the horizontal branch, the effective temperatures are higher
for stars with less mass in the envelope.
After the He exhaustion (point 14) the star will ascend the giant branch for the second
time. In the degenerated core there are two energy sources: (i) The He-burning shell
above the C-O core and (ii) the H-burning shell above the He-shell. A deep convective
shell is located above both the H and He shells.
8 CHAPTER 1. AGB STARS
Figure 1.2: Schematic evolution of a 1 M⊙ star. Image taken from Habing & Olofsson
(2013).
1.2. EVOLUTION IN THE AGB PHASE 9
1.2 Evolution in the AGB phase
The core and stellar envelope: Following the conversion of He in the core into C
and O , the He-burning zone shifts outward, causing the C-O core to contract to a size of
approximately 103 km and a temperature of 108 K, akin to the characteristics of a white
dwarf. This core contraction, coupled with the expansion of the envelope, results in a
rapid increase in luminosity, marking the onset of the E-AGB (early asymptotic giant
branch) phase. During this phase, He-burning in a shell becomes the primary source of
energy. The stellar envelope expands to about 108 km, and the surface temperature is
around 3000 K. At a certain point in the early evolution, the envelope becomes pulsa-
tionally unstable, exhibiting a characteristic time scale on the order of a hundred days
(see Sec. 1.2.1).
Double-shell burning and thermal pulses: When a star reaches the luminosity
of approximately the tip of the RGB, denoted byMbol = −4 mag or 3000 L⊙, a new phase
starts. The star becomes capable of burning both He and H in shells. Periodically, the
thin He-layer surrounding the core undergoes rapid burning into C , which is primarily
added to the core’s mass. This brief episode, known as a thermal pulse or He-shell
flash (point 15), induces a moderate luminosity modulation in the star. In the intervals
between thermal pulses, the star resumes burning H . This intermittent phase of He and
H burning is called the TP-AGB phase.
The dredge-up: During a thermal pulse, convection reaches layers undergoing
nuclear processing, transporting enriched material to the stellar surface. This process
results in the accumulation of products from nuclear burning, particularly ”new” C . In
some cases, this enrichment persists, causing the C abundance to surpass that of O .
Consequently, the atmosphere undergoes a transition from O-rich (with C/O ∼ 0.5, also
known as M-type AGB stars) to C-rich conditions (C/O > 1, also known as C-type AGB
stars)1. This represents the second significant dredge-up for a low-mass star and the third
for an intermediate-mass star. However, in certain cases, this transition is prevented by
the hot-bottom-burning (HBB), particularly for AGB stars with initial masses above
a specific limit (around ∼4 M⊙ for solar composition, less at lower metallicity). The
HBB inhibits the formation of C stars by causing CNO-burning to destroy C produced
1This transition may involve the formation of S-stars, covering spectral types MS, S, SC, with 0.5 ≤
C/O ≤ 1 (Höfner & Olofsson 2018).
10 CHAPTER 1. AGB STARS
during He -burning, transforming it into N . Consequently, the most massive AGB stars
are expected to retain their M-type classification throughout their evolution.
Figure 1.3: Examples of the CSE for different AGB stars. Top-left panel: Image of
R Scl (Maercker et al. 2012), top-right panel: Image of U Ant (Kerschbaum et al.
2017). Both images were obtained with the Atacama Large Millimeter/submillimeter
Array (ALMA). Bottom panel: Image of IRC+10216, obtained with the Institut de
Radioastronomie Millimétrique (IRAM) 30 m telescope (Cernicharo et al. 2015). The
different fields of view of these images range from 100 to 200 arcsec.
The atmosphere: The outer regions of the stellar envelope reach temperatures
low enough to facilitate the formation of molecules, and their presence significantly
influences the spectral characteristics of stars. This influence can vary considerably,
particularly between C-type and M-type stars, at least within the near infrared (NIR)
range. Additionally, the structure of the stellar atmosphere is affected to some extent.
Pulsations of the star deposit mechanical energy into the weakly bound outer regions,
causing the atmosphere to extend significantly beyond what a hydrostatic model would
1.2. EVOLUTION IN THE AGB PHASE 11
predict. This extension leads to the development of shocks, and at sufficiently high
altitudes, grain condensation occurs in the post-shock gas.
The CSE: An AGB star will eventually lose mass at a considerable rate in the form
of a slow (∼ 15 km s−1) wind, and this forms a CSE of escaping gas and particles.
The gas is to a large extent molecular with H2 dominating. The molecules are located
in envelopes of different sizes (being CO and H2 the largest ones by far), or shells
of different radii and widths, depending on their resistance to UV photons and their
formation routes (see Fig. 1.3 for examples of AGB stars with CSE). The dust grains
survive much further from the star than do the molecules. The molecular setup, as well
as the grain properties, are a strong function of the C/O ratio. At some point, the CSE
merges into the surrounding Interstellar medium (ISM). In this region the temperature is
below 10 K, and the particle density is less than 10 cm−3. This may occur at a distance
of 1013 km or more from the star, and it defines the outer limit of the AGB star.
Termination of the AGB: Ultimately, the rate of ejection of the matter is higher
than the growth rate of the core, and from this point on, the mass-loss process (see Sec.
1.3) determines the evolution. Eventually, all the matter around the core is dispersed
into space, and this marks the endpoint of AGB evolution, and the beginning of the
post-AGB evolution. In practice, this starts when the stellar envelope mass becomes
less than about 0.01 M⊙. For a very short time (on the order of a few thousand years)
the star becomes a post-AGB star, and then, for a time of a few tens of thousands of
years, it becomes a Planetary Nebula. The nebula is expanding rapidly; it fades, and
only the stellar core remains as a WD. Figure 1.1 shows a schematic view of the different
parts of an AGB with the information of relevant processes occurring.
1.2.1 Stellar Variability
In the previous Section, I mentioned that at some point during the early evolution of an
AGB star, its envelope becomes pulsationally unstable with a characteristic time scale
on the order of hundred days. Here, I will describe how this variability was first observed
and how this discovery influenced the posterior studies of stellar variability. Then, since
the object of the study of this thesis is a Mira star, I will provide some characteristic
properties of this kind of stars.
Like most astronomical discoveries, the first observations of stellar variability were
12 CHAPTER 1. AGB STARS
accidental. Fabricius2 observed, in 1596, that the star o Ceti exhibited a range of
variability from being invisible to reaching second magnitude (see Fig.1.4). As a tribute,
he named the star Mira, meaning the wonderful. This discovery provided further evidence
supporting the notion, advocated strongly by Tycho Brahe after the appearance of the
new star in Cassiopeia in 1572, that the heavens are not static.
Figure 1.4: The visual light curve of o Ceti extracted from the American Association of
Variable Star Observers (AAVSO) database.
A significant milestone occurred in the 18th century when Goodricke made a ground-
breaking discovery. In 1782, he announced the variations of β Persei, also known as
Algol. These variations were found to be periodic and stable. To explain these observa-
tions, Goodricke proposed a model of an eclipsing binary star system, where the stars
are closely situated and unresolved. This was the first explanation for stellar variability
and one that served as a stimulus to search for other binary star systems. The stim-
ulus was such that, in 1784, Goodricke discovered the variability of δ Cephei (δ Cep),
and in the same year, Piggot discovered the variability of η Aquilae (η Aql). Since the
periods of these stars were found to be similar to that of Algol, a natural choice was to
apply the eclipse model to explain their variability. Nowadays, both δ Cep and η Aql are
2David Fabricius (9 March 1564 – 7 May 1617) Born in Esens, Frisia (now Germany), he was a
Lutheran pastor serving small communities in East Frisia, Germany. In addition to his pastoral duties
he studied Astronomy. The study of sunspots, as well as the first variable star recorded by the modern
Western world, are his two main notable contributions to the field of Astronomy
1.2. EVOLUTION IN THE AGB PHASE 13
known to be pulsating variables unrelated to the Algol class of binaries, but the general
explanation of the variations had to wait for more than 150 years after this discovery.
To summarize, Table 1.1 shows a list of the earliest discoveries of what are now known
to be intrinsically variable stars.
Star name Discoverer Year
o Ceti Fabricius 1596
ψ Leonis Montanari 1667
κ Sagitarii Halley 1676
χ Cygni Kirch 1687
δ Cephei Goodricke 1784
η Aquilae Piggot 1784
R CrB Piggot 1795
α Herculis Wm. Herschel 1796
Table 1.1: Early Discoveries of Intrinsic Stellar Variables. Table adapted from (Shore
2003).
During the next century, the search for stellar variability grew in importance for
Astronomy. The work of Argelander at Bamberg, the establishment of the Bonner
Durchmusterung star atlas and the international efforts to construct the Carte du Ciel,
a photographic all-sky atlas, motivated the discovery and cataloging of variable stars.
The most important discovery made during the 19th Century was connected with
observations of variable stars in the Large and Small Magellanic Clouds, with most of
the work being performed at Harvard by H. Levitt. Around 1905, she realized that these
stars were generally of some type as the variable δ Cephei. The extraordinary discovery
was the relationship between the period and the apparent brightness of the stars, the
so-called period-luminosity relation.
Together with the observational discoveries, theoretical work on the cause of vari-
ability continued. Eddington suggested the pulsation model for stellar variability and
derived the linearized, adiabatic pulsation equation (see Eddington 1918, 1919). Subse-
quent works by Cowling (1934) and Rosseland (1949) further clarified the instabilities.
This epoch culminated with the encyclopedic review article by Ledoux & Walraven
14 CHAPTER 1. AGB STARS
(1958), which still serves as an important reference for all discussions of stellar variabil-
ity.
Cox (1963) highlighted the importance of stellar convection zones in driving pulsa-
tion and producing the radial velocity variations in Cepheids. The first machine simula-
tions of nonlinear, nonadiabatic pulsation models were accomplished in the mid-1960s.
Around the same time, the first large-scale grids of stellar interior models were produced.
Schwarzschild & Härm (1965) discovered the thermal instability of double shell sources.
Mira stars
As one may suspect at this time, intrinsic stellar variability can be found in all types
of stars, from the MS to the RGB phase; even WDs have variability. Depending on
the nature of the pulsations, the stars can be classified between radial and non-radial
pulsators. Each category has sub-categories defined according to pulsation periods,
luminosities, temperatures, etc. In the category of radial pulsators exist the Long Period
Variables (LPV) with pulsation periods between 100 - 1000 days. Mira variables fall in
this sub-category. Named after the discovery of o Ceti, Miras are emission-line cool
giants generically arising from a population of low to intermediate mass (1 - 2 M⊙)
stars. Miras typically have luminosities around 103 L⊙ and temperatures less than 4000
K. These stars have some of the largest visual amplitudes observed among pulsating
variables, ranging from 2.5 up to 11 mag, but generally considerably smaller in the
infrared (at a wavelength of few microns, the amplitude is less than one magnitude).
The periods range between 80 and 1000 days, although many of the longest period
systems are quite irregular in their light curve behavior.
Miras are perhaps the best-studied examples of non-linear pulsators, having ampli-
tudes large enough to produce the ejection of matter and non-adiabatic behavior of the
envelopes. One interesting aspect of these stars is that their pulsation periods are suf-
ficiently long, and their amplitudes large enough, that time-dependent convection may
be important in modeling their internal structure. It is worth noting that considerable
theoretical attention is currently being paid to these stars (see e.g, Chiavassa et al. 2024;
Rees et al. 2024).
1.3. MASS LOSS DURING THE AGB PHASE 15
1.3 Mass loss during the AGB phase
I mentioned in previous sections that during the AGB phase the stars will lose most of
their mass in the form of a slow and dense wind. As a consequence, all the processed
material will be returned to the ISM, contributing to its enrichment. Understanding the
wind formation and mass loss mechanisms in AGB stars is then of prime importance to
study the cycle of matter in the Universe. In this section, I will describe our current
knowledge of the mass-loss processes that take place during the AGB phase.
1.3.1 Early indications of mass loss
The first evidence of (substantial) post-main-sequence mass loss came through obser-
vations of the binary system α Her by Deutsch (1956). Circumstellar absorption lines
were seen in the spectrum of the companion star, suggesting circumstellar gas as far
away as 105 R⊙ from de M giant. With an estimated outflow velocity of 10 km s−1,
the gas velocity is well above the escape velocity at this point. The estimated mass-loss
rate was 3 × 10−8 M⊙ yr−1. In a study of similar systems, Reimers (1975) derived an
empirical relation between mass-loss rate and stellar parameters (luminosity, mass, and
radius; Eq. 1.1 ) that bears his name:
dM
dt
= 4 × 10−13 (L/L⊙)(R/R⊙)
(M/M⊙)
M⊙yr−1, (1.1)
where L is the luminosity, M the mass, and R the stellar radius.
There were also indirect indications. Auer & Woolf (1965) found that the Hyades
cluster contains about a dozen WDs; each one must have a mass below the Chan-
drasekhar limit of 1.4 M⊙. But, the Hyades cluster has stars with a mass of 2 M⊙ that
are still on the MS. Consequently, the WDs progenitors must have lost at least 0.6 M⊙
during their post-MS evolution. The advent of infrared observations led to considerable
progress, e.g., the Two-Micron Sky Survey (TMSS) by Neugebauer & Leighton (1969)
and the Air Force Geophysics Laboratory (AFGL) catalog of Price & Walker (1976).
It was shown that late M-giants often have an excess of radiation at infrared wave-
lengths, an effect of dust in a CSE. In some objects, the dust optical depth is so high
that most of the stellar energy is emitted at wavelengths longer than 2 µm. Gillett
et al. (1968) and Woolf & Ney (1969) attributed an emission band around 10 µm to
small silicate particles and Gehrz & Woolf (1971) estimated circumstellar dust masses
16 CHAPTER 1. AGB STARS
and mass-loss rates of M giants. Similar results were found for Carbon stars (Woolf &
Ney 1969). The appearance of different dust compositions around M- and C-type giants
was explained by Gilman (1969) to be a consequence of the high binding energy of the
CO molecule, which prevents the less abundant of the two elements from forming other
molecules and solids.
Wickramasinghe et al. (1966) showed that light pressure on graphite grains can push
both the grains and the surrounding gas out of the stellar gravitational field due to
momentum transfer between dust and gas. Circumstellar OH (in maser emission) was
first detected in 1968 (Wilson & Barrett 1968), and circumstellar CO rotational lines
were observed a few years later towards the carbon star IRC+10216 (CW Leo) (Solomon
et al. 1971). The first detailed physical/chemical model of a CSE was presented by
Goldreich & Scoville (1976). The OH masers were explained as due to infrared pumping
of OH molecules in a shell where they formed through the dissociation of H2O by the
interstellar radiation field.
In the decades following these early findings, significant progress has been made
thanks to both space- and ground-based instrumentation (see e.g., Höfner & Olofsson
2018). These developments have enabled the study of spectral energy distributions
with a wide wavelength coverage, high spectral resolution data on circumstellar atom-
ic/molecular lines, and structural information on CSEs via interferometry of nearby
objects, complemented by large-scale surveys and photometric monitoring of large sam-
ples of stars. In recent years, high-angular resolution techniques, covering wavelengths
from the visual to the radio regime, have started to give direct insights into the dy-
namical atmospheres and wind acceleration regions of AGB stars which is essential for
understanding the mass-loss mechanisms.
1.3.2 Mass loss from a theoretical point of view
Given the high luminosities of AGB stars, it is intuitive to consider radiation pressure
as a driving force for their winds. To initiate a wind, the outward force exerted by
radiation must surpass the inward force from gravitational potential. Radiative acceler-
ation, involving the transfer of momentum from photons to matter through absorption
and scattering, depends critically on the total radiative cross section per mass of stellar
material. Estimating the required opacity to initiate a stellar wind via radiative pressure
involves comparing the strength of radiative acceleration relative to the inward gravi-
1.3. MASS LOSS DURING THE AGB PHASE 17
tational acceleration acting on a material layer. This comparison can be expressed as
follows:
arad
agrav
=
〈κ〉L∗
4πcGM∗
, (1.2)
where 〈κ〉 is the flux mean opacity (total cross-section per unit mass), G is the gravi-
tational constant, M∗ and L∗ are the mass and luminosity of the star in consideration,
respectively. The previous equation defines three regimes; (i) where radiative acceler-
ation dominates over gravitational attraction (arad/agrav > 1), (ii) where gravitational
attraction dominates over radiative acceleration (arad/agrav < 1), and (iii) where both
forces are equal in magnitude. This situation corresponds to a critical value of the
flux-mean opacity, given by:
〈κ〉crit =
4πcGM∗
L∗
≈ 2.6
(M∗
M⊙
)( L∗
5000L⊙
)−1
[cm2 g−1]. (1.3)
It is worth noting that wind driving might be a cumulative process, which means that
radiative pressure is acting on material that is temporarily moving upwards on a ballistic
trajectory due to the effect of shock waves. In that case, a value of arad/agrav that is above
one may be sufficient to lift material out of the potential well of the star. Nevertheless,
〈κ〉crit from Eq. 1.3 provides us with a good estimate for the flux-mean opacity required
to drive and outflow.
The previous discussion raises one important question: What are the sources of
opacity around AGB stars? In the atmospheres of cool giant stars, most of the gas is
in molecular form, and certain species may condense into solid particles in the upper
atmospheric layers. In principle, both gas and dust opacities contribute to the radiative
pressure. In practice, however, the flux-mean opacity will be dominated by species
that are more efficient at absorbing or scattering photons in the NIR wavelength region
around the maximum stellar flux.
Opacity of molecules
The rich spectrum of vibro-rotational band heads3 from abundant molecules coincides
with the stellar flux maximum in the NIR wavelength region, and their influence on ra-
3In the context of molecular spectroscopy, a vibro-rotational spectrum refers to the spectrum of a
molecule that arises from its vibrational and rotational motion. In such spectra, a ”band head” refers
to the most intense or prominent peak within a spectral band.
18 CHAPTER 1. AGB STARS
diative pressure has been explored in the literature (see, e.g., Jorgensen & Johnson 1992;
Elitzur et al. 1989). Despite molecules playing a significant role in shaping the thermal
structures of AGB stellar atmospheres, their collective flux-mean opacity remains insuf-
ficient to impact atmospheric and wind dynamics significantly (see, e.g., Liljegren et al.
2016). Hence, the observed mass loss rates in outflows cannot be solely attributed to
molecules.
Opacity of dust grains
In contrast to molecules, dust grains are good candidates for driving winds of AGB stars,
as they can be very efficient at absorbing and/or scattering the stellar photons, depend-
ing on their chemical composition and size. Certain materials, which are ingredients
of circumstellar dust around AGB stars (e.g., amorphous carbon, or magnesium-iron
silicates) are also efficient absorbers in the NIR range around the stellar flux maximum.
Moreover, if the abundance of dust particles with the right properties is large enough in
the outer atmospheric layers, their collective opacity can easily exceed 〈κ〉crit and trigger
an outflow.
Dust species that are commonly associated with wind driving (amorphous carbon or
silicates) will typically condense at temperatures between 1000 - 1500 K, or below. We
can estimate the temperature, Td(r), of such dust species as a function of the distance
from the center of the star if we assume that (i) the condition of radiative equilibrium sets
the grain temperature, (ii) the radiation field follows a Planck function corresponding
to an effective temperature T∗, geometrically diluted with distance, and (iii) the grain
absorption coefficient can be approximated with a power law of the form κabs αλ
−p in the
relevant wavelength region around the stellar flux maximum, where λ is the wavelength
and p is a material-dependent constant (see, e.g., Bladh & Höfner 2012, for an example
of such fits). Td(r) is then given by:
Td(r) ≈ T∗
( 2r
R∗
)−
2
4 + p , (1.4)
where R∗ is the stellar radius. On the other hand, we can define Rc as the minimum
distance at which dust grains of a given material can form (without evaporating due to
radiative heating), i.e., where the grain temperature is equal to the dust condensation
temperature of the material. According to Eq. 1.4, when setting Td(Rc) = Tc we obtain
1.3. MASS LOSS DURING THE AGB PHASE 19
for Rc:
Rc
R∗
≈ 0.5
(Tc
T∗
)−
4 + p
2 . (1.5)
The absorption coefficient’s wavelength dependence, denoted by the value of p, results
in distinct reactions of various condensates to the proximity of a star. A grain material
with a positive p tends to heat up when exposed to the stellar radiation field. It excels
in absorbing radiation, being more effective in absorption than emission, consequently
pushing the condensation distance farther away from the star. Conversely, a grain ma-
terial with a negative p exhibits the opposite behavior. Therefore, the survivability of a
particular grain type close to a star is dictated by the combination of its condensation
temperature and the slope of the absorption coefficient, p. An example of this trend is
shown by Bladh & Höfner (2012) (Fig. 1.5), where the authors plotted curves of constant
condensation distances as a function of p and Tc using Eq. 1.5.
To gain more feeling about condensation distances, let’s take two examples. For
Fe-bearing silicates with a condensation temperature of 1100 K and p ≈ 2.3, Rc ≈ 10
R∗ for a star with a temperature of 2800 K. The corresponding value for Fe-free silicates
(Tc = 1100 K, p ≈ -0.9), present in the atmospheres of O-rich AGB stars, is Rc ≈ 2
R∗ for a star with an effective temperature of 2800 K. As we see, for Fe-free silicates,
the gas condenses into dust at ∼ 2 R∗, but if the gas is initially located on the surface
of the star, how does it reach these distances? What is the mechanism responsible for
lifting the gas from the surface of the star beyond the condensation radius? The answer
is: pulsation.
The most accepted theory establishes that when the stars pulsate, the induced shock
waves levitate the gas in the atmosphere to a distance known as levitation distance Rl,
which is defined as follows:
Rl
R∗
=
R0
R∗
[
1 − R0
R∗
(
u0
uesc
)2
]−1
, (1.6)
where u0 is the initial velocity of the gas at a distance, R0, and uesc = (2GM∗/R∗)1/2
is the escape velocity. The combination of pulsation-induced shock waves and radiation
pressure on dust grains is the most promising mechanism to explain atmospheric levita-
tion and wind driving. In this scenario, the mechanisms for driving dust differ depending
on the dust properties. I will discuss these differences in the following section.
20 CHAPTER 1. AGB STARS
Figure 1.5: Curves of constant condensation distance (Rc/R∗ = 2, 4, and 10) as a
function of the power law coefficient p and the value of the condensation temperature
Tc, using Eq. 1.5 for two stars, one with T∗ = 2800 K (solid lines) and the other with
T∗ = 2500 K (dashed lines). The filled red circles mark where the selected set of dust
species are situated in the p/Tc-plane. The dust species in parenthesis have uncertain
or interpolated optical data in the NIR. The blue dashed-dotted line marks a limit that
indicates whether certain dust particles are more efficient at absorbing or emitting the
radiation from the star. The image and its description were taken from (Bladh & Höfner
2012).
The mass loss scenario for C- and M-type AGB stars
Theoretical and observational studies support that C-type AGB stars form amorphous
C grains with large opacities that can be accelerated by radiation pressure (see e.g.,
Gautschy-Loidl et al. 2004; Sacuto et al. 2011; Rau et al. 2015, 2017). The mass-loss
rates, wind velocities, spectral energy distributions, and photometric variations resulting
1.3. MASS LOSS DURING THE AGB PHASE 21
from the time-dependent atmosphere and wind models are in good agreement with obser-
vations (see, e.g., Nowotny et al. 2011, 2013; Eriksson et al. 2014). However, due to the
strong absorption by C -dust, the stellar photospheres and their inner atmospheres may
be severely obscured, making detailed comparisons with spatially resolved observations
difficult (Sacuto et al. 2011; Wittkowski et al. 2017).
For M-type AGB stars, the high chemical complexity around them makes identifying
the dust grains driving the dusty winds more complex (see Fig. 1.1). Höfner (2008)
suggested a scenario where the radiative pressure that triggers the outflows is caused
by photon scattering on Fe-free silicate grains. These silicates are highly transparent at
visual and NIR wavelengths, resulting in significantly less radiative heating and smaller
condensation distances than for Fe-bearing silicates. In order to create sufficient radia-
tive pressure by scattering, the grains need to be typically of a size comparable to the
wavelengths where the stellar flux reaches its maximum, which means that the grain size
has to be in the range 0.1 - 1 µm. Recent observational studies report the presence of
dust grains with radii of about 0.1 - 0.5 µm in AGB stars at distances below 2 - 3 stellar
radii (see, e.g., Norris et al. 2012; Ohnaka et al. 2016; Ohnaka et al. 2017). Despite these
results, the evidence to support or discard these models is still scarce, and even when
the theory agrees with these observations, more evidence is needed in order to better
constrain the mass loss models of O-rich AGB stars in different environments.
1.3.3 The role of convection in the mass-loss processes
Convection is a complex process in stellar astrophysics with important physical impli-
cations, such as heat transfer, mixing, and generation of magnetic fields, but also in
pulsation and winds on evolved stars (see, e.g., Montargès et al. 2021; Ohnaka et al.
2016; Ohnaka et al. 2017). These effects, in turn, significantly impact the structure
and evolution of stars, particularly during the AGB phase (Kupka 2004). In fact, in the
mass-loss scenario described in the previous section, non-stationary convective flows also
contribute to the generation of the pulsation-induced shock waves responsible for lifting
the gas beyond the condensation distance (Eq. 1.5). Furthermore, asymmetries observed
in AGB stars and other evolved stars could be caused by the underlying irregular global
convective flows even if they are not visible in the stellar surface (Freytag et al. 2017).
Given its importance, there is a necessity of increasing our current knowledge of
stellar convection both from the theoretical and observational sides. The first source
22 CHAPTER 1. AGB STARS
in which the convection models are calibrated is our Sun. According to Schwarzschild
(1975), the Sun has around 2 million convective cells on its surface with typical sizes
around 1000 - 2000 km across. Extrapolating this result would imply that on the surface
of giant and supergiant stars, there should be only a few large (several thousand larger
than those of the Sun) convective cells.
Studying and characterizing stellar convection has its difficulties, both from theo-
retical and observational points of view. From the observational point of view, one of
the main challenges of studying convection in evolved stars is that their atmosphere is
composed of gas and dust, which obscures the surface and makes it challenging to ob-
serve the convective patterns (Wittkowski et al. 2017). Furthermore, the resolution of
current telescopes is limited, which restricts our ability to observe small-scale features
on the surface of stars. Nevertheless, with increasing angular resolution, interferometric
measurements at IR wavelengths offer a significant advantage for observing the surface
of evolved stars. In fact, Paladini et al. (2018) used the Precision Integrated Optics Near
Infrared ExpeRiment (PIONIER) at the Very Large Telescope Interferometer (VLTI) to
observe the surface of the AGB star π1 Gru, and presented reconstructed images that
reveal the presence of complex patterns that were associated with giant convection cells
(see Fig. 1.6). Those authors also characterized these structures by obtaining their
contrast and size and compared these quantities with existing theoretical models that
relate the convective structure with the stellar parameters (see Fig.1.7). Their results
provided an observational test that served as motivation for studying the surfaces of
other evolved stars with high-angular resolution techniques.
From the theoretical point of view, realistic modeling of the underlying physical pro-
cesses involved in convection is a notoriously tricky problem, and theoretical studies
found in the literature are usually restricted to 1D simulations (see Sec. 1.3.3). Still,
these dynamical 1D atmospheres and wind models have been used successfully to ex-
plore the basic dust-driven mass-loss processes, relying on a parameterized description
of sub-photospheric velocities due to radial pulsations (the so-called piston models; see,
e.g., Höfner 2008; Jeong et al. 2003). The combination of star-in-a-box models such as
CO5BOLD (Freytag et al. 2003) with wind models such as Darwin (Höfner et al. 2016)
can bring new descriptions of the asymmetric distributions of dust and gas observed in
evolved stars. Those models can be confronted with high-angular resolution observations
and analyses of the surface and the innermost CSE’s of AGB stars.
1.3. MASS LOSS DURING THE AGB PHASE 23
Figure 1.6: Best-reconstructed images of the stellar surface of the S-type AGB star π1
Gru using interferometric data from the instrument PIONIER/VLTI in the H−band.
These images were obtained using two different image reconstruction algorithms:
SQUEEZE (upper panels) and MiRa (lower panels). The images reveal the presence
of complex structures associated with giant convection cells. The white circle at the
bottom right of the top panel in a, corresponds to the angular resolution of the obser-
vations (λ/2B ∼ 2 mas). This image was taken from (Paladini et al. 2018).
1D convection and the Mixing Length Theory (MLT)
The following subsections are based on A Review of the Mixing Length Theory of Con-
vection in 1D Stellar Modeling by Joyce & Tayar (2023) and are intended to provide
a general description on how is the convection described in a 1D model, its limitations
and the different approaches that can be employed to avoid those limitations.
Understanding and simulating energy transport in stars poses a complex challenge,
especially due to the intricate nature of convection within stellar interiors. Direct obser-
vations struggle to unveil the details of convection, introducing significant uncertainties
into stellar models. A reliable stellar evolution code must not only capture the di-
verse behaviors of stars over enormous ranges of temperature, density, and pressure but
also span evolutionary timescales, from fractions of milliseconds for nuclear processes to
billions of years. Considering spatial scales ranging from nuclear cross-sections to thou-
24 CHAPTER 1. AGB STARS
Figure 1.7: The characteristic granule size xg obtained from the PIONIER images of π1
Gru (triangle; log(g) = 0.44) is quite different from the solar value (circle; log(g) = 4.4).
As it can be seen, the extrapolations of theoretical predictions (dotted line, dashed line,
and solid line) to lower effective temperatures (Teff = 3200 K, log(g) = 0.44) are in
good agreement with the observations. The image and its description were taken from
Paladini et al. (2018).
sands of solar radii, simultaneous inclusion of all scales in evolutionary codes becomes
computationally unfeasible. To address this, especially in the context of convection —a
3D, turbulent, non-linear, and time-dependent process— parametrizations are essential.
The Mixing length theory (MLT) serves as a practical solution by offering a static 1D
approach for the convenience of stellar modeling in such situations.
The MLT provides a simplified description of the bulk movement of fluids in analogy
to molecular heat transfer. It considers the packets within a convection zone as discrete
”parcels” of fluid, each with local uniform physical characteristics. These parcels are
in pressure equilibrium with their surroundings but not in thermal equilibrium. Con-
sequently, hotter parcels will move towards cooler regions, and cooler parcels will move
towards hotter regions. The buoyancy force, driven by the under-density of hot parcels,
1.3. MASS LOSS DURING THE AGB PHASE 25
causes them to rise and expand, while cooler ones sink and compress. The characteristic
distance a parcel can travel before losing its locally homogeneous characteristics is associ-
ated to its mean-free path, measured in terms of the pressure scale height d ln(P )/d ln(T ),
of the stratified fluid.
In this section, I will formulate the mixing length theory by reproducing the standard
derivation of the amount of flux carried by convection. This derivation is based on a
framework presented in a review of the Mixing Length Theory of Convection in 1D Stellar
Modeling by Joyce & Tayar (2023) which, at the same time, is based on derivations by
Cox & Giuli’s principles of stellar structure (Cox 1968), Kippenhan and Wieger’s Stellar
Structure and Evolution (Kippenhahn et al. 1990), and Cassisi & Salari’s Stars and
Stellar Populations (Salaris & Cassisi 2005).
We first introduce the concept of pressure scale height Hp as a measure of the distance
over which the total pressure, P = Pgas + Prad, changes by a factor 1/e. The following
definition is valid under the assumption of hydrostatic equilibrium:
−d lnP
dr
≡ ρg
P
=
1
Hp
. (1.7)
In the MLT the ”parcel of fluid” is considered to be in pressure, but not in thermal
equilibrium with its surroundings. We next define, at at given radius r, (1) the average
(ambient) temperature gradient of the fluid, with respect to the pressure of all matter,
and (2) a temperature gradient of the fluid parcel itself, also taken with respect to the
total pressure. We define (1) by:
∇ ≡ d lnT
d lnP
, (1.8)
and (2) by
∇parcel ≡ d lnTparcel
d lnP
, (1.9)
where T and Tparcel are the average (ambient) and the parcel temperature, respectively.
Fig 1.8 shows a schematic representation of the previous scenario. The difference between
the parcel’s temperature and the ambient temperature at some radial shift ∆r, in terms
of a first-order Taylor expansion, can be expressed as
∆T (∆r) = Tparcel(r + ∆r) − T (r + ∆r) ≃ ∆r
[
dTparcel
dr
− dT
dr
]
, (1.10)
Assuming that the temperature change over ∆r is small, T ≃ Tparcel, and so
∆r
[
dTparcel
dr
− dT
dr
]
→ ∆rT
[
−d lnT
dr
−
(
−d lnTparcel
dr
)]
. (1.11)
26 CHAPTER 1. AGB STARS
Using the chain rule and the definition from Eq. 1.7, we can rewrite
d lnT
d lnP
=
d lnT
dr
dr
d lnP
=
d lnT
dr
(−Hp). (1.12)
By substituting Eqs. 1.8 and 1.9 into Eq. 1.11 we obtain
∆T (∆r) = ∆r
T
Hp
(∇ − ∇parcel). (1.13)
We note here that the moving parcel can exchange heat with its environment, and
so ∇parcel is a function of the rate at which this exchange takes place. However, the
assumption that the parcel does not change heat with its surroundings, i.e., that it is
moving adiabatically is a reasonable and simplifying assumption for stellar conditions.
In this case,
∇parcel → ∇ad =
d lnT
d lnP
∣
∣
∣
ad
. (1.14)
In the derivation of the MLT, it is important to note that the previous equation (Eq.
1.14) is not interchangeable with Eq. 1.8. The adiabatic temperature gradient, ∇ad,
and the average or ambient temperature gradient, ∇, are not, in general, the same. We
also do not know the definition of ∇ in a convective region. Substituting ∇parcel for ∇ad
and combining with Eq. 1.13 yields
∆T (∆r) = ∆r
T
Hp
(∇ − ∇ad). (1.15)
Now, the convective flux transported by a parcel of fluid moving upwards with velocity
v over some distance λ is given by
Fconv =
1
2
ρvcP
[(
dT
dr
)
ad
−
(
dT
dr
)]
λ, (1.16)
where ρ is density, cP is specific heat, and the factor of 1/2 emerges from the assumption
that half of the material in a given layer is rising while half is falling. Given the following
relation
dT
dr
= − T
HP
d lnT
d lnP
= − T
HP
∇ (1.17)
and Eq. 1.13, we have
Fconv =
1
2
ρvcP
[(
− T
HP
∇ad
)
−
(
− T
HP
∇
)]
λ. (1.18)
Finally, by rearranging the previous equation, we have
Fconv =
1
2
ρvcPT
λ
HP
(∇T − ∇ad) , (1.19)
1.3. MASS LOSS DURING THE AGB PHASE 27
with
αMLT ≡ λ
HP
. (1.20)
The previous definition is that of the mixing length parameter. This is a dimensionless
parameter, characterizing the ”distance”, measured in terms of the pressure scale height
(Eq. 1.7), that a parcel of convective material can travel.
This αMLT can be thought of in many ways, including as the convective mean-free
path, as a measure of the convective efficiency, or as a pseudo-physical quantity that
captures the change in entropy from the base to the top of the convection zone. As
can be seen, Eq. 1.19 is determined according to two things: (i) the difference in the
temperature gradients and (ii) the value of αMLT. In other words, αMLT modulates the
suppression or enhancement of the surface convective flux so that larger values mean
that more flux is carried by convection.
Figure 1.8: An illustration of the movement of a convective blob across the medium.
The distance a convective parcel can travel is measured in multiples of the pressure
scale height, Hp. A larger mixing length implies that the parcel travels over a larger
pressure differential before mixing with its surroundings, which corresponds to more
efficient transport of the flux by convection. Image taken from Joyce & Tayar (2023).
28 CHAPTER 1. AGB STARS
Limitations of the MLT
Although MLT has proven to be useful in describing convection, sometimes it requires
to make naive and outright incorrect physical assumptions. In this section, I will briefly
describe some of these assumptions and contrast them with the real physical phenomena
taking place inside the stars.
Advective processes involve the transportation of material or energy through the
bulk motion of fluids. Models employing a diffusive approximation represent particle-like
fluid parcels diffusing through the region to redistribute heat. However, in reality, the
notion of a homogeneous convective material unit sustained over a considerable distance
is invalid. Convection operates through a continuum of constantly-changing up-flow and
down-flow channels.
The MLT supposes rigid convective boundaries. However, in the real life, con-
vective boundaries are permeable, flexible and are subject by the inertia of convective
motions carrying convective plumes across the convective-radiative boundary. In these
regions, the local mixing occurs and MLT cannot account for this.
The MLT introduces the simplification of strictly vertical paths due to its 1D for-
mulation. This limitation arises because physical convection involves shearing, fragmen-
tation, reorientation, and deletion of flow channels, features that a formulation relying
solely on radial displacements cannot capture.
Time dependency in the MLT is not taken into account. This means that the
formulation cannot capture any physics happening faster than the convective turnover
time. In standard treatments, convective regions are assumed to be instantaneously
mixed over one evolutionary time step.
Flow channels. In physical convection, there is an asymmetry between upflows
(hot material rising) and downflows (cool material falling). The MLT falls to account
for this effect.
The mixing length parameter, αMLT, in many MLT formulations is often consid-
ered a free, numerical parameter with loose connections to any physical property of the
star. However, it is linked to the entropy jump between zones of efficient and less efficient
convection, specifically, the asymptotically adiabatic portion of the convection zone and
the top of the convective envelope. Although characteristics like entropy jump, density
gradient, and convective envelope properties are not directly observable, asteroseismol-
ogy provides an indirect method for investigating these attributes. Direct calibration
1.3. MASS LOSS DURING THE AGB PHASE 29
of αMLT is crucial, and observations of our nearest star offer the most straightforward
means for this calibration.
Extensions to 3D
Scientific efforts in utilizing 3D simulations of convection are on the rise, aiming to re-
fine the representation of mixing length in a manner that aligns more realistically with
convection physics. Recent advancements in 3D simulations enable consideration of not
only local and static mixing length but also global derivatives, asymmetries between up-
flows and downflows, larger-scale star properties, and transverse versus radial differences.
These simulations generally indicate a slight variation in the mixing length concerning
the star’s composition, luminosity, and surface gravity, likely attributed to the physics
of turbulent dissipation (Trampedach et al. 2014; Magic et al. 2015; Sonoi et al. 2018).
In terms of trends observed in 3D simulations between αMLT and global properties
such as metallicity, they often align directionally with trends implied by 1D stellar
model calibrations to match observational data (Bonaca et al. 2012; Creevey et al. 2015;
Tayar et al. 2017; Joyce & Chaboyer 2018a,b; Viani et al. 2018). However, quantitative
agreement is rare. In some cases, 1D-to-observational calibrations suggest a need for a
αMLT variation about ten times larger than indicated by 3D simulations (see e.g., Joyce
& Chaboyer 2018a; Trampedach et al. 2014). This points to a potential need for further
refinement in three-dimensional simulations or additional corrections to the physics of
1D stellar evolution calculations impacting the inferred mixing length.
30 CHAPTER 1. AGB STARS
2C
h
ap
te
r
Interferometry
***
This chapter gives an introduction to the terminology and the basic concepts of opti-
cal/infrared long baseline interferometry. It is based mainly on the books by Lawson
(2000), Monnier (2003), and Buscher (2015).
2.1 Principles of interferometry
A telescope’s angular resolution (∆θ) is its capacity to discern between two separate
close objects in the sky. It depends on the wavelength of interest, λ, and on the size of
the telescope, usually quantified by the diameter D, of its primary mirror. The Rayleigh
criterion establishes the diffraction-limited resolution of a circular telescope as:
∆θtel = 1.22
λ
D
rad. (2.1)
According to this relation, a telescope cannot distinguish objects separated by an angular
distance smaller than ∆θtel. Table 2.1 shows typical angular resolutions in milliarcsec-
onds for three wavelength regimes, considering Eq. 2.1.
Visible NIR MIR
0.4 - 1 µm 1 - 5 µm 10 - 20 µm
D = 8 m 12.5 - 31.5 mas 31.5 - 157 mas 315 - 630 mas
D = 100 m 1 - 2.5 mas 2.5 - 12.5 mas 25 - 50 mas
Table 2.1: Angular resolution of telescopes with different diameters (8 and 100m) for
three different band-passes.
31
32 CHAPTER 2. INTERFEROMETRY
Let’s take an example to illustrate the angular resolution. AGB stars are very bright
in the infrared (see Chapter 1); therefore, NIR or mid infrared (MIR) wavelength ob-
servations are well suited to study these kinds of objects. If the typical size of the dust
shell at 11 µm is ∼ 50 mas (50 AU at adistance of 1 kpc), we would need a telescope of
approximately 45 m to resolve it. However, this is more than the size of the European
Extremely Large Telescope (E-ELT; 38 m) that the European Southern Observatory
(ESO) is currently building.
Another example includes active galactic nuclei (AGN), whose nucleus is about one
mas in the MIR (5 µm). To properly resolve this part of the AGN, a 100 m telescope
would be needed. However, such kinds of instruments are too expensive, so astronomers
have to adopt a different technique to resolve astronomical objects: interferometry.
In this technique, the light of a source is collected through two (or more) spatially
separated apertures and then is combined coherently. Instead of observing the brightness
distribution of the object directly, we observe an interference pattern (see Fig. 2.1). The
interferometric observable derived from these fringes is the so-called complex visibility, V .
Using this technique, we gain angular resolution at the cost of having a lower sensitivity
and no direct imaging.
The angular resolution achieved by a two-station interferometer, considering the full
constructive interference, is:
∆θint =
λ
B
rad, (2.2)
where B is the projected separation in the plane of the sky between the two telescopes,
called the baseline length. This criterion can be relaxed when considering a partial-
constructive interference up to a limit of ∆θint >
λ
2B
, where the interference is completely
destructive (Monnier 2003).
The link between the intensity distribution of a source projected in the plane of the
sky and the complex visibility is given by the Van-Cittert-Zernicke theorem:
The complex visibility is the Fourier transform of the source intensity dis-
tribution on the sky at the spatial frequency corresponding to the projected
baseline on the sky per observing wavelength
This theorem can be expressed by the following equation:
I(α, β) =
∫ +∞
−∞
∫ +∞
−∞
V (u, v)e2πi(αu+βv)dudv, (2.3)
2.1. PRINCIPLES OF INTERFEROMETRY 33
Figure 2.1: (a) Beams of light from a source arriving at angles of ±β on a detector.
(b) The fringe pattern on the detector. (c) A one-dimensional cut through the fringe
pattern. Image taken from Buscher (2015).
where I is the intensity distribution of the observed source, (α, β) are the angular coor-
dinates of the source in the sky (in radians), V is the so-called complex visibility, and
u, v are the spatial frequencies that describe the brightness distribution. The spatial
frequencies are related to the baseline, B, through the following expressions:
u =
Bα
λ
; v =
Bβ
λ
. (2.4)
In interferometry, the term u − v plane is used frequently. This is the complex plane
sampled with the interferometric observations. For instance, a two-aperture telescope
configuration will give one point in the u − v plane, meaning one visibility. Ideally, in
order to reconstruct an image of a given source, one should fill the u − v plane with
several interferometric observations (see, e.g., Fig 2.2).
2.1.1 Visibility amplitude
The visibility is a measurement of the contrast between the fringes in the interference
pattern (Fig. 2.1 and Eq. 2.5). Usually, this quantity is normalized, meaning that
it takes values between 0 and 1. For instance, an unresolved point source will have
34 CHAPTER 2. INTERFEROMETRY
Figure 2.2: Example of the u-v plane achievable with GRAVITY using three Telescope
arrays: short (Bmax ∼ 30 m), medium (Bmax ∼ 100 m) and large (Bmax ∼ 140 m), each
color corresponds to a different wavelength between 2 and 2.4 µm (K band).
high contrast fringes (see the upper row in Fig. 2.3) and a normalized visibility of 1.
For a spatially resolved star, light from across the stellar surface combines incoherently,
causing a decrease in the fringe contrast and, therefore, a drop in the visibility amplitude.
Thus, the normalized visibility amplitude for a resolved source will be less than 1. The
bigger the star, the lower the fringe amplitude (see Fig. 2.3 for an illustration of this
effect). By measuring the visibility, one can determine a star’s size, shape, and surface
features. This can be done, for example, through modeling the brightness distribution
of the source (see Fig. 2.4).
V =
Imax − Imin
Imax + Imin
. (2.5)
2.1.2 Visibility phase
The phase is related to the fringe location with respect to a reference (called phase
center or phase reference). This quantity carries information about the geometry, i.e.,
symmetry of the object. For instance, an object is called asymmetric when the shape
2.1. PRINCIPLES OF INTERFEROMETRY 35
Figure 2.3: Comparison between single-dish and interferometric observables for two
objects with different sizes. Left: two simulated stars of different sizes. Center: the same
objects as seen by a single-dish telescope. Right: the interference patterns produced by
the two objects. It is worth noting that the contrast between fringes is smaller for more
extended objects. This image was taken from the ESO web page.
is not centrally-symmetric; an ellipse is symmetric, but an ellipse with a spot inside is
not 1. Talking about astrophysical objects, a binary whose components have the same
brightness is symmetric, but if the brightness is different, then the binary is asymmetric.
These asymmetries in the binary can also be observed in the visibility curves (see Fig
2.5).
2.1.3 The effect of the atmosphere on the interferometric observations
When performing optical interferometry from the ground, one must remember that quan-
tities like the phase of the complex visibility are very sensitive to atmospheric turbulence.
In fact, many of the critical design parameters of an interferometer, such as site selection,
operating wavelength, aperture size, and inclusion of adaptive optics, are driven by the
boundary conditions set by the atmosphere. The sensitivity of interferometers, and the
1An ellipse is centrally symmetric because it has an axis of symmetry passing through its center. If
we add a spot away from the center of the ellipse, the symmetry breaks because there is no longer a
corresponding point on the opposite side of the ellipse that mirror’s the spot position.
36 CHAPTER 2. INTERFEROMETRY
Figure 2.4: Top panel: Simulated images of three models (From left to right: Uniform,
Gaussian, and Limb-darkened disk). All the models have a diameter (or Full Width
Half Maximum (FWHM) in the case of the Gaussian model) of 15 mas. Bottom panel:
Visibility curves for the three models versus baseline length at 2 µm (see the labels on
the plot). The ”first lobe” of the visibility curve (before the visibility amplitude reaches
the first zero) contains information about the size and shape of the star. The second
and higher-order lobes of the visibility curve also encode information about the surface
features of the star, including limb-darkening, starspots, and convection cells.
precision of astrometric measurements, depend strongly on the seeing (Lawson 2000).
When taking interferometric observations, the atmospheric path length above each
telescope generates phase shifts, η, in the recorded fringes. These phase shifts vary in
proportion to the phase shift differences between the two telescopes forming a baseline.
Hence, the phase of the measured visibilities, arg[V (u, v)], is affected as follows:
arg[V (u, v)] = Φreal + (ηj − ηi). (2.6)
The existence of phase delays, η, modifies the visibility phase so that it can no longer
be associated with the Fourier phase of the source’s brightness distribution. A new
2.1. PRINCIPLES OF INTERFEROMETRY 37
Figure 2.5: Visibility curves of two binary systems. The separation between the peaks in
the visibility curve provides a measurement of the binary separation, while the minimum
visibility reflects the flux ratio (F2/F1) between the components.
interferometric observable is then introduced to correct for this effect: the Closure Phase
(CP) (Baldwin et al. 1986; Jennison 1958). This quantity is obtained by multiplying the
visibility functions of three different baselines that form a closed triangle (see Fig. 2.6).
By doing this, the telescope-dependent errors are canceled:
Φb12,meas = Φb12,real + (η1 − η2)
Φb23,meas = Φb23,real + (η2 − η3)
Φb31,meas = Φb31,real + (η3 − η1)
Φb12,meas+Φb23,meas + Φb31,meas = Φb12,real + Φb23,real + Φb31,real. (2.7)
For instance, a point-symmetric object will have CPs equal to zero, while variations
in the closure phases between −180 and 180 degrees will be observed for asymmetric
objects. Given an interferometric array of N elements (telescopes), the total number of
independent visibilities is equal to N(N − 1)/2, and the number of independent closure
phases is equal to (N − 1)(N − 2)/2. Since we need three phases to construct a CP, the
number of CPs is smaller than the number of interferometric phases.
Apart from the atmospheric fluctuations, the detector noise can also affect the sen-
38 CHAPTER 2. INTERFEROMETRY
Figure 2.6: Phase shifts in the recorded fringes are introduced by the atmospheric tur-
bulence. These phase shifts vary in proportion to the phase shift differences between the
two telescopes that form a baseline and cancel out when we calculate the closure phase.
This image was taken from the University of Arizona.
sitivity of an interferometer. In the case of optical interferometry, this noise is produced
by random arrivals of photons so, it follows a Poisson statistics. Variations of the av-
erage number of photons Np arriving to the detector at a given time are given by:
∆Np =
√
Np. Thus, the total number of photons recorded in a given interval of time is
expressed as:
Np = η
Pν∆ν
hν
τ, (2.8)
where η is the total efficiency of the optical transmission system, Pν is the power spectral
density, τ is the time interval, ν is the frequency of the observation, h is the Planck
constant and ∆ν the bandwidth used. Since the coherence time is very small at IR
wavelengths compared to radio (e.g., at 2.2 µm: τc = 25 ms and at 1 mm: τc = 14 s),
the number of photons recorded is considerably lower in a single exposure. To mitigate
this problem, a series of short exposures are taken and the final detected signal, over
which the fringe contrast (Eq. 2.5) is measured, corresponds to the average of the
total power in the recorded exposures. Since the signal-to-noise-ratio (SNR) for optical
interferometers scales as NV 2, the total power of the visibility V 2 is a better observable
than the amplitude V .
2.2. INTERFEROMETRIC FACILITIES AROUND THE WORLD 39
2.2 Interferometric facilities around the World
This section is intended to give an overview of the facilities around the world that can
provide IR interferometric observations. Figure 2.7 shows the two big interferometers
and others that are currently offered to the international community. In particular, I
will only describe the CHARA and VLTI (Fig. 2.9) and give more emphasis to the VLTI
instruments GRAVITY, PIONIER, and MATISSE (Fig. 2.10). I will briefly describe
the reduction and calibration processes of these instruments, as the analyses reported in
this thesis are based on observations taken with these three instruments. Finally, most
of the information has been adopted from the websites of the instruments as well as the
instrument’s user manuals.
Figure 2.7: This sketch shows the footprints of modern O/IR interferometers, with
telescopes and baselines to scale. The telescopes for LBTI (purple), CHARA (blue), and
VLTI-UTs (red large circles) are fixed, while the NPOI (green) and VLTI-ATs (red small
circles) telescopes are mobile and can be re-positioned. Image taken from Eisenhauer
et al. (2023).
40 CHAPTER 2. INTERFEROMETRY
2.2.1 CHARA
The Center for High Angular Resolution Astronomy (CHARA)2 is an optical interfero-
metric array of six telescopes located in Mount Wilson, California. The array consists
of six 1 m telescopes that are arranged in a Y-shape configuration providing 15 base-
lines ranging from 34 to 331 m (see Fig. 2.8). With the maximum baseline (331 m),
CHARA can resolve objects as small as 0.2 mas at visible wavelengths and 0.5 mas in
the near-infrared.
This instrument is particularly suited for stellar astrophysics, where it is used to
measure the diameters and temperatures of stars, image features such as spots and flares
on their surfaces (e.g., Norris et al. 2021), and map the orbits of close binary companions
(e.g., Anugu et al. 2023). Other projects range from detecting other planetary systems
(e.g., Martin et al. 2023), imaging stars in the process of formation (e.g., Labdon et al.
2023), and studies of bright transient phenomena like novae (e.g., Schaefer et al. 2014).
2.2.2 VLTI
The VLTI 3 (Lena 1979; Beckers et al. 1990; Schöller 2007) is an interferometric array
consisting of the four VLT 8 m Unit Telescopes (UTs) or the four movable 1.8 m Auxiliary
Telescopes (ATs). It provides milli-arcsec angular resolution in low and intermediate
(R=4000) spectral resolution at NIR and MIR wavelengths.
The four 8.2-m UTs and the four 1.8-m ATs are the light-collecting elements of
the VLTI. The UTs are set on fixed locations while the ATs can be relocated on more
than 10 different stations covering a baseline range between 30 and 130 m. All of the
VLTI instruments recombine the light from four telescopes simultaneously. After the
light beams have passed through a complex system of mirrors and the light paths have
been equalized by the delay line system, the light re-combination is performed by the
PIONIER and GRAVITY instruments in the near-infrared and by MATISSE in the
mid-infrared part of the spectrum.
Due to its unique characteristics, the VLTI has become a very attractive means
for scientific research on various objects like young pre-main sequence stars and their
protoplanetary disks (see e.g., Ibrahim et al. 2023), post-main sequence mass-losing stars
(Rosales-Guzmán et al. 2023), binary objects and their orbits (see e.g., Chauvin et al.
2https://www.chara.gsu.edu/public/tour-overview
3https://www.eso.org/sci/facilities/paranal/telescopes/vlti.html
2.2. INTERFEROMETRIC FACILITIES AROUND THE WORLD 41
Figure 2.8: The CHARA array telescopes and its instruments (see the labels in the
image). The image was taken from the CHARA web page.
2023), and extragalactic objects such as active galactic nuclei (see e.g., Isbell et al. 2023).
PIONIER
PIONIER (Le Bouquin et al. 2011) is a four-telescope interferometric combiner at the
VLTI operating in the H band (1.6 - 1.99 µm). It has a limited spectral capability of one
spectral channel (bandwidth ∼ 0.3 µm) or six spectral channels (spectral resolution R ∼
40). It measures simultaneously six squared visibilities and four closure phase values in
order to characterize the spatial structure of a source at an angular scale as small as one
mas at the longest baselines available. Its main characteristics are high efficiency making
it suitable for interferometric imaging, and high accuracy of the measured quantities
making it suitable for high contrast observations.
A typical PIONIER observation is composed of:
• A telescope preset
42 CHAPTER 2. INTERFEROMETRY
Figure 2.9: The Paranal Observatory telescopes. The instruments listed in italics are
not yet in operation. The year in brackets indicates the start of operations. Image taken
from the ESO webpage.
• An instrumental setup
• An instrumental acquisition which requires
– nexp exposures of nscan interferograms (typical 5 exp or 100 scans)
– One dark exposure
– A kappa matrix (measuring each beam individually = 4 exposures)
The observation of calibration stars is interleaved with the observations of the science
star, all with the same instrumental setup. For instance, a typical sequence would be
CAL1-SCI-CAL2-SCI-CAL1.
GRAVITY
GRAVITY (Abuter et al. 2017) is a four-telescope beam combiner operating in K band
(1.9 - 2.4 µm). The instrument provides a spatial resolution λ/2B = 1.58 mas for a
2.2. INTERFEROMETRIC FACILITIES AROUND THE WORLD 43
baseline of 130 m at a wavelength λ = 2µm. It comes with multiple subsystems designed
to accurately adjust and manage the incoming wavefront from the astronomical source
as it travels from the telescopes to the beam combiner. These sub-systems include a
fringe tracker, an IR wavefront sensor, a polarization control system, a pupil-guide, and
a field-guide system.
GRAVITY forms spectrally dispersed interference fringes of either a single astro-
nomical source (single-field) or of two sources simultaneously (dual-field and dual-field
wide). The interferometric quantities provided by this instrument are the absolute and
differential visibility, spectral differential phase, and closure phase. These quantities
allow reconstructing images at a spatial resolution of two mas (when chosing the maxi-
mum baselines), as well as time-resolved differential astrometry per spectral bin at the
accuracy of few tens of µas.
A full observation with GRAVITY is composed of different steps. The first one is the
instrument calibration, necessary to know the characteristics of the detector, the beam
combiner’s transfer matrix (pixel to visibility matrix), and some characteristics of the
optics (spectrograph calibration, fiber dispersion, polarization). Once the instrument is
calibrated, the observations of a science target can be performed either in single-field or
dual-field mode.
Closest to the science observation, a calibrator star in single field mode, or a binary
calibrator star in dual field mode has to be observed. Knowing the diameter of the
calibrator is necessary to calibrate the two transfer functions of the instrument: the one
of the fringe tracker channel and the one of the science spectrometer channel. In the
dual field mode, two sets of observations are acquired by swapping the two targets in
order to calibrate the zero point of the metrology. The pipeline reduces the data coming
from the calibrator observation templates and generates visibilities and phases.
The calibration process can be started once the visibilities of the calibrator and
science targets are recorded. During this process, the transfer function is computed
from Eqs. 2.9 and 2.10, which is used to calibrate the science visibilities according to
Eqs. 2.11 and 2.12.
V 2
inst = (V 2
measured/V
2
expected)CAL, (2.9)
CPinst = (CPmeasured − CPexpected)CAL, (2.10)
V 2
calib = (V 2
measured)SCI/V
2
inst, (2.11)
CPcalib = (CPmeasured)SCI − CPinst, (2.12)
44 CHAPTER 2. INTERFEROMETRY
where CAL stands for the calibrator star and SCI for the science target.
MATISSE
The Multi AperTure mid-Infrared SpectroScopic Experiment (MATISSE) (Lopez et al.
2022) is a four-telescope beam combiner operating in the L (3 - 4 µm), M (4.5 - 5 µm),
and N (8 - 13 µm) bands. The instrument provides spatial resolution λ/2B = 2.25
mas, and λ/2B = 7.25 mas for a baseline of 140 m and for wavelengths λ = 3µm and
λ = 10µm respectively. It can provide data separated in many wavelengths at different
spectral resolutions, low, medium, high, and very high (R ∼ 30, 500, 1000, and 3500) in
the L−M bands, and low and high (R ∼30, 220) in the N band.
In each of the L, M , and N bands, MATISSE uses a multiaxial all-in-one beam
combination, that is the four separated telescope beams are focused onto the detector
by a common lens and produce six dispersed fringe patterns that are mixed into one single
focal spot (an airy disc), representing the so-called interferometric channel. In the L and
M bands, the total incoming flux is split between the interferometric channel and the four
photometric channels (SiPhot mode). This mode allows us to simultaneously perform
the source photometry and thus to properly calibrate the visibilities. In the N band,
the source photometry is obtained after the fringes recording by closing alternatively the
MATISSE shutters to feed the four photometric channels. This allows to improve the
sensitivity by sending 100% of the total incoming flux in the different channels (HiSens
mode).
To account for the thermal background emission, which is significant in the mid-
infrared, two additional modulation steps have to be added to calibrate the data: the
first one is chopping modulation which consists of observing the target and an empty
area of the sky alternatively. The second modulation is an optical path difference (OPD)
modulation, aimed at removing the low-frequency peak contamination, containing the
contribution from the thermal background, to the fringe peaks (in the Fourier transform
of the interferograms). The instrument’s different modes (telescope chopping, detector
integrations, OPD modulation) are synchronized all together to ensure minimum over-
heads when observing. Finally, a beam modulation (or beam commutation) is applied
to the observations by swapping two-by-two beams in a four-step sequence. This fea-
ture is aimed at removing detection artifacts in the phases measured by the instrument
(differential phases and closure phases). These artifacts are instrument phase residuals
2.3. IMAGE RECONSTRUCTION IN OPTICAL INTERFEROMETRY 45
located in the optical train between the Beam Commuting Device (BCD) and the detec-
tor (Lopez et al. 2022). The BCD is a module used to calibrate closure and differential
phases.
A MATISSE observational sequence is composed of twelve one-minute-long exposures
with simultaneous interferometric and photometric data in the L and M bands and four
interferometric exposures followed by eight photometric exposures in the N band. The
first four exposures are taken without chopping while the following eight exposures are
taken with chopping. Each exposure is taken in one of the four configurations of the
BCD. Overall this leads to four reduced OIFITS24 (Duvert et al. 2017) files per pointed
target (one for each BCD configuration) containing, per exposure, six dispersed squared
visibilities, differential visibilities phases, and three independent dispersed closure phases
(out of four in total). Then, the interferometric calibration is applied to the data: the
targets used as calibrators (CAL) are corrected from their expected observables given the
predicted value of their angular diameters, and each science target (SCI) is corrected
from the transfer function estimated. The calibrated OIFITS2 files of the four BCD
configurations are merged to produce a single calibrated OIFITS2 file for each CAL-SCI
pair.
2.3 Image reconstruction in optical interferometry
As mentioned in previous sections, one of the caveats of performing interferometric
observations is that we are not able to obtain a direct image of an observed object, but
instead, we have its visibilities and phases (closure phases). Hence, interferometric data
are challenging to interpret, and we must rely on different techniques to understand
our observations. Perhaps, the most illustrative method of interpreting interferometric
information is to perform image reconstruction. However, owing to the sparse u − v
coverage, the image reconstruction process is ill-posed as there are more unknowns (pixels
in the image) than measurements. To sort out these effects, the problem of image
reconstruction is approached using a regularized minimization of the form:
x+ = argminxf(x), (2.13)
4The Optical Interferometry exchange format (OIFITS) is a standard for exchanging calibrated data
from optical (visible/infrared) interferometers. The standard is based on the Flexible Image Transport
System (FITS).
46 CHAPTER 2. INTERFEROMETRY
(a) (b)
(c)
Figure 2.10: The three second-generation instruments at the VLTI: (a) PIONIER, (b)
GRAVITY, and (c) MATISSE. The images were taken from the respective ESO instru-
ment web pages.
were x+ represents the solution we are looking for and
f(x) = fdata(x) + µfprior, (2.14)
where the likelihood term fdata(x) measures the discrepancy between the model and the
available data, while the regularization term fprior(x) measures the discrepancy with the
prior information. The hyperparameter µ is used to adjust the relative weight of the
constraints set by the measurements and the ones set by the priors (Renard et al. 2011;
Sanchez-Bermudez et al. 2018).
There are many image reconstruction packages. The most commonly used are:
BSMEM (Baron & Young 2008), SPARCO (Kluska et al. 2014), MiRA (Thiébaut 2009)
and SQUEEZE (Baron et al. 2010). They mainly differ on how they minimize the
2.3. IMAGE RECONSTRUCTION IN OPTICAL INTERFEROMETRY 47
function f (Eq. 2.14). The reconstructed images reported in this thesis were obtained
using BSMEM (which uses a gradient descent algorithm) and SQUEEZE (which uses a
Markov Chain Montecarlo; MCMC). In the following sections, I will describe the capa-
bilities of each one of these two algorithms. These sections are based on the article by
Sanchez-Bermudez et al. (2016b), available at the Jean-Marie Mariotti Center (JMMC)
page.
2.3.1 BSMEM
The BiSpectrum Maximum Entropy Method (BSMEM) was introduced in 1992 to
demonstrate image reconstruction from aperture synthesis data, and it has been ex-
tensively enhanced since then (Baron & Young 2008). The algorithm employs a fully
Bayesian approach to the inverse problem of finding the most probable image given the
evidence, using the Maximum Entropy approach to maximize the posterior probability
of an image. BSMEM uses a trust region method with non-linear conjugate gradient
steps to minimize the sum of the log[likelihood] (χ2) of the data given the image and a
regularization term expressed as the Guy-Skilling entropy (Gull & Skilling 1999):
fprior(x) = Σn[xnlog(xn/x̄) − xn + x̄n], (2.15)
with x being the default image, which is the one that would be recovered in the absence
of any data. The likelihood term used by BSMEM assumes independent Gaussian noise
statistics for the amplitude and phase of the measured data. The optimization engine
is MemSys created by Maximum Data Consultants Ltd., which implements the strategy
proposed by Skilling & Bryan (1984) and automatically finds the most likely value of
the hyperparameter µ that weights fprior with respect to fdata.
The user can choose the default image (or model image) x between a Gaussian,
a Lorentzian, a uniform disk, or a delta function centered in the field of view, which
conveniently fixes the location of the reconstructed object. It is worth noting that the
bispectra and power spectra are invariant to translations. We will come back to this
point in Chapter 4.
2.3.2 SQUEEZE
The primary engine of SQUEEZE is a MCMC algorithm. In interferometry, MCMC was
first successfully implemented in image reconstruction in the radio by Sutton & Wandelt
48 CHAPTER 2. INTERFEROMETRY
(2006) then in the optical/IR by the MACIM software (Ireland et al. 2006). Like its
predecessors, SQUEEZE uses the simulated annealing technique to converge to a global
minimum of the posterior probability. An N×N pixel image is described by a collection
of M elements, which may be pixel fluxes, wavelet coefficients, or model parameters.
Defined steps, such as changing the flux of elements or moving them, are potentially
taken on these elements in order to modify the image.
Regularization functions in SQUEEZE
Compared to BSMEM, SQUEEZE allows different regularization functions, either us-
ing them individually or combined. Some of the regularization functions allowed by
SQUEEZE are the following (Sanchez-Bermudez et al. 2016b):
• Entropy: With this regularizer, one ensures that the image that best fits the data
is the one with the minimum information. It is implemented in SQUEEZE as the
sum of the natural logarithm of the Γ function for each one of the non-zero pixels
in the image:
R(x) = Σ log |Γ(x)|, ∀x 6= 0 (2.16)
• Dark Energy: This counts all the zero elements in the image. The function
implemented in SQUEEZE also adds values to the edges of the image. The equation
of this regularizer is the following one:
R(x) = Σx, ∀x = 0 (2.17)
• Uniform disk: corresponds to the squared sum of the Lp-norm, with p=0.5, of
the local gradient of the pixels (xl) in the image. The SQUEEZE implementation
ignores the image borders. The general formula of this regularizer is:
R(x) = (Σ
√
||∇x||l,s)2, (2.18)
with ||∇x||l,s =
√
(xl+1,s − xl,s)2 + (xl,s+1 − xl,s)2 (2.19)
• Total variation: Instead of quantifying the intensity distribution in the image,
this regularizer quantifies its spatial gradient distribution, helping us to describe
uniform areas in the sought image with steep but localized changes (Strong & Chan
2003).
2.4. STUDIES OF AGB STARS THROUGH IR INTERFEROMETRIC OBSERVATIONS49
R(x) = (Σ||∇x||l,s)2, (2.20)
with ||∇x||l,s as defined by Eq. 2.19.
• L0 sparsity norm: If an image is reconstructed using this regularizer, it will be
forced to be the sparsest solution. Unlike other algorithms, SQUEEZE has the
advantage of being able to optimize non-convex functions like the L0 norm. The
mathematical formula of this regularizer is:
R(x) = ΣL0 (2.21)
with L0 = 1, ∀x > 1
• L2-norm: This regularization function is part of the ℓ-p norm functions (with p
= 1.5, 2, and 3), which tend to produce smooth images as they reduce the variance
in the pixels. The expression of the L2 norm is the following:
R(x) = Σx2 (2.22)
• Transpectral L2 regularization: This regularization function is used for poly-
chromatic image reconstructions. The function computes the L2-norm of the pixels
across the different spectral channels in the data sets. The expression of this reg-
ularizer is:
R(x) = Σn
√
ΣlΣs(xλ0
l,s)2 + (xλ1
l,s)2 + ...+ (xλn
l,s )2 (2.23)
Figure 2.11 shows a reconstructed image using BSMEM and SQUEEZE. The (sim-
ulated) image corresponds to the ”2004 interferometric imaging contest” and is always
used for testing the algorithms once installed.
2.4 Studies of AGB stars through IR interferometric ob-
servations
In this section, I will give a brief review of some of the more common analyses that
are performed on interferometric data. These analyses have provided important clues
for our understanding of the physics around AGB stars and include: Model fitting,
image reconstruction, model atmospheres and radiative transfer modeling. The works
referenced in the text can be found also in Tables 2.2 and 2.3.
50 CHAPTER 2. INTERFEROMETRY
Figure 2.11: Reconstruction of the interferometric imaging contest data set of 2004. The
”true” image is shown in the left panel, while the SQUEEZE and BSMEM reconstructed
images are shown in the center and right panels, respectively. As can be seen, the
BSMEM image is smoother due to the maximum entropy regularizer employed. The
image was taken from Baron et al. (2010).
2.4.1 Model fitting
Early in this chapter, I mentioned that the estimation of the angular size of observed
objects relies on geometric modeling of the V2. This method is widely utilized by various
authors for estimating the angular size of their sources, as evidenced in works such as
Ireland et al. (2004b), Lacour et al. (2009), Rau et al. (2015), Ohnaka et al. (2016),
Wittkowski et al. (2016), Ohnaka et al. (2019), Chiavassa et al. (2020), Perrin et al.
(2020), and Rosales-Guzmán et al. (2023).
2.4.2 Image reconstruction
When the interferometric observables are good enough, it is possible to perform an image
reconstruction. Various algorithms can be applied for this task, and the quality of the
reconstructed images is contingent upon factors such as data quality, the selection of
regularizers, among others.
2.4.3 Model Atmospheres
The atmosphere of an AGB star is intricate, revealing evidence of diverse atoms and
molecules through observed spectral transitions. Estimating the physical conditions re-
sponsible for these transitions involves atmospheric modeling. Models like PHOENIX
2.4. STUDIES OF AGB STARS THROUGH IR INTERFEROMETRIC OBSERVATIONS51
(Hauschildt & Baron 1999) can generate synthetic spectra and reproduce several spec-
tral transitions at different wavelength ranges for comparison with observed spectra
and interferometric observables (V2), allowing the determination of fundamental stel-
lar parameters such as temperatures, gravities, size in both continuum and molecular
transitions, optical thickness, and more. Atmospheric modeling has been widely used to
determine properties of AGBs (see e.g. Wittkowski et al. 2008, 2011, 2016; Rau et al.
2017; Wittkowski et al. 2018). However, it should be noted that these models are spheri-
cally symmetric, so they can only be used with V2 data, meaning that information about
asymmetries contained in the CPs is not taken into account.
2.4.4 Radiative transfer modeling
In radiative transfer modeling, a star with a temperature (Teff) and a radius (R∗) emit-
ting photons in a specified distribution (usually a blackbody distribution) is simulated.
The photons emitted at the star’s surface are free to travel through the CSE, and this
interaction results in an intensity Iν that can be calculated and then compared with
observational data. In this type of model, the measured intensity will depend on the
composition and distribution of the dust, making it possible to determine dust proper-
ties such as temperature, distribution, and optical thickness as a function of wavelength,
among other properties. These models provide the user with a Spectral Energy Distribu-
tion (SED), but it is also possible to obtain synthetic images, from which interferometric
observables can be extracted and compared with observational data (see e.g., Karovicova
et al. 2011, 2013; Rosales-Guzmán et al. 2023).
2.4.5 The importance of observational evidence
Recent high-angular resolution interferometric observations of nearby AGB stars, con-
ducted in both visual and infrared wavelengths, have unveiled intricate, non-spherical
distributions of gas and dust in the close CSE (see, e.g., Wittkowski et al. 2017; Paladini
et al. 2018). Temporal monitoring further indicates changes in atmospheric morphology
and grain sizes over weeks and months, as demonstrated by studies like Ohnaka et al.
(2017) and Norris et al. (2012). Spherically symmetric atmosphere and wind models
prove to be inadequate for investigating such dynamic phenomena.
The development of a predictive mass loss theory requires realistic models encom-
passing stellar convection and pulsation, both of which significantly impact the stellar
52 CHAPTER 2. INTERFEROMETRY
atmosphere and wind. However, a well-elaborated mass loss theory remains ineffective
without observations for comparison or to establish boundary conditions for these mod-
els. Hence, there is a growing interest in providing analyses of the innermost parts of
AGB stars through high-angular resolution observations. These observations often in-
volve multi-epoch and multi-band analyses of various sources, enabling the monitoring
of gas and dust distributions over time and the tracking of chromatic changes across
different wavelength bands. In addition to traditional analyses of high-resolution in-
terferometric data, chromatic image reconstruction has emerged as a crucial tool for
comprehending the morphology of AGB stars at various scales—from the surface to the
CSE (see, e.g., Wittkowski et al. 2017; Paladini et al. 2018; Rosales-Guzmán et al. 2023).
In Sec. 2.3, I highlighted the numerous challenges associated with reconstructing
images from interferometric data. When these challenges are mixed with the intrinsic
variable nature of AGB stars, conducting multi-epoch and multi-band analyses through
image reconstruction becomes a titanic task. To illustrate this fact, I performed a search
into the ESO Telescope Bibliography webpage5 with the objective of finding the bulk of
interferometric studies performed on AGB stars since 2004. Tables 2.2 and 2.3, and Figs.
2.12(a) to 2.12(c) show those findings. In all the plots, I present the effective temperature
(Teff) of the stars versus their log luminosity in solar luminosities log(L/L⊙).
An overview of the total number of AGB stars (including M-, S- and C-types) with
interferometric analyses is presented in Fig. 2.12(a). Figure 2.12(b) contains a sub-
sample including only M-type AGB stars. I selected the stars considering whether they
were analyzed through image reconstruction (either in combination with model fitting
or other analysis techniques) or analyzed through model fitting only. Finally, in Fig.
2.12(c), I show all the M-type AGB stars that include a multi-epoch and multi-band
image reconstruction.
The sample size significantly decreases as more stringent requirements are imposed
on the analyses conducted on the complete sample of stars. Ultimately, only a handful
of sources undergo multi-epoch and multi-band interferometric image reconstructions.
This highlights the need to expand the number of interferometric observations of nearby
AGB stars. Such efforts are indispensable for acquiring additional information that will
enhance our understanding of these intricate phenomena.
5https://telbib.eso.org
2.4. STUDIES OF AGB STARS THROUGH IR INTERFEROMETRIC OBSERVATIONS53
(a) (b)
(c)
Figure 2.12: Overview of the total number of AGB stars studied with interferometric
high-resolution observations. The panel a) shows both M- and C-type AGB stars (see
the labels in the plot). The panel b) shows a sub-sample containing only M-type AGBs.
The colors indicate whether image reconstruction was performed or not. Finally, the
panel c) shows only M-type AGB stars where the main analysis is conducted through
image reconstruction. The colors indicate whether a multi-epoch and/or multi-band
analysis was performed. Note here that the source analysed in this thesis is labeled in
the plot. In all the plots we present the effective temperature, Teff versus log(L/L⊙) for
the different sources.
54 CHAPTER 2. INTERFEROMETRY
N
a
m
e
S
p
.
T
y
p
e
T
eff
[K
]
lo
g
(L
/
L
⊙
)
B
an
d
s
A
n
al
y
se
s
E
p
o
ch
s
In
st
ru
m
en
t
R
ef
er
en
ce
A
Q
S
gr
C
-t
y
p
e
26
0
0
3.
7
N
M
o
d
A
tm
1
M
ID
I
R
au
et
al
.
(2
01
7)
C
h
i
C
y
g
S
-t
y
p
e
26
0
0
3.
5
H
M
o
d
F
it
,
P
I
4
IO
T
A
L
ac
ou
r
et
al
.
(2
00
9)
C
L
L
ac
M
-t
y
p
e
32
93
3.
6
9
H
Im
R
ec
,
M
o
d
F
it
1
C
H
A
R
A
C
h
ia
va
ss
a
et
al
.
(2
02
0
)
IK
T
au
ri
M
-t
y
p
e
22
9
4
3
.9
2
K
P
ol
ar
Im
,
M
o
d
F
it
1
A
M
B
E
R
O
h
n
ak
a
et
al
.
(2
01
6)
O
C
et
M
-t
y
p
e
24
5
0
3
.5
7
K
1D
,
3D
M
o
d
A
tm
2
A
M
B
E
R
W
it
tk
ow
sk
i
et
al
.
(2
01
6
)
O
C
et
M
-t
y
p
e
2
45
0
3.
5
7
H
Im
R
ec
,
M
o
d
F
it
1
IO
T
A
P
er
ri
n
et
al
.
(2
02
0
)
P
i
G
ru
S
-t
y
p
e
32
00
3
.9
H
Im
R
ec
1
P
IO
N
IE
R
P
al
ad
in
i
et
al
.
(2
01
8)
R
A
q
r
M
-t
y
p
e
22
50
3.
63
K
1D
,
3D
M
o
d
A
tm
1
A
M
B
E
R
W
it
tk
ow
sk
i
et
al
.
(2
01
6
)
R
C
a
r
M
-t
y
p
e
2
80
0
3
.6
2
H
Im
R
ec
1
P
IO
N
IE
R
M
on
n
ie
r
(2
00
3)
R
C
ar
M
-t
y
p
e
28
0
0
3.
62
H
Im
R
ec
2
P
IO
N
IE
R
R
os
al
es
et
al
.,
(i
n
p
re
p
)
R
C
ar
M
-t
y
p
e
28
0
0
3.
62
K
Im
R
ec
,
M
o
d
F
it
1
G
ra
v
it
y
R
os
al
es
-G
u
zm
án
et
al
.
(2
02
3)
R
C
a
r
M
-t
y
p
e
2
80
0
3
.6
2
V
is
M
o
d
F
it
1
S
U
S
I
Ir
el
an
d
et
al
.
(2
00
4b
)
R
C
n
c
M
-t
y
p
e
23
90
3.
6
6
K
1D
,
3D
M
o
d
A
tm
3
A
M
B
E
R
W
it
tk
ow
sk
i
et
al
.
(2
01
6)
R
C
n
c
M
-t
y
p
e
23
9
0
3
.6
6
N
R
ad
T
ra
n
sf
M
o
d
2
M
ID
I
K
ar
ov
ic
ov
a
et
al
.
(2
01
3
)
R
D
or
M
-t
y
p
e
27
10
3.
6
4
K
Im
R
ec
,
M
o
d
F
it
1
A
M
B
E
R
O
h
n
ak
a
et
al
.
(2
01
9)
R
L
eo
M
-t
y
p
e
25
7
0
3.
6
K
1D
,
3D
M
o
d
A
tm
1
A
M
B
E
R
W
it
tk
ow
sk
i
et
al
.
(2
01
6)
R
L
ep
C
-t
y
p
e
2
80
0
3
.8
5
N
M
o
d
A
tm
1
M
ID
I
R
au
et
al
.
(2
01
7
)
R
P
eg
M
-t
y
p
e
24
8
0
3
.6
2
K
M
o
d
A
tm
4
G
ra
v
it
y
W
it
tk
ow
sk
i
et
al
.
(2
01
8)
Table 2.2: High-resolution observations of AGB stars. All the references were retrieved
from a search conducted on the Telbib web page and include studies from 2004 to 2023.
The letters stand for ModAtm = Model atmospheres, ImRec = Image reconstruction,
ModFit = Model fitting, PI = Parametric imaging, RadTransMod = Radiative transfer
modeling, and PolarIm = Polarimetric imaging.
2.4. STUDIES OF AGB STARS THROUGH IR INTERFEROMETRIC OBSERVATIONS55
N
am
e
S
p
.
T
y
p
e
T
ef
f
[K
]
lo
g
(L
/
L
su
n
)
B
a
n
d
s
A
n
al
y
se
s
E
p
o
ch
s
In
st
ru
m
en
t
R
ef
er
en
ce
R
R
A
q
l
M
-t
y
p
e
2
12
7
3.
5
3
N
R
ad
T
ra
n
sf
M
o
d
13
M
ID
I
K
ar
ov
ic
ov
a
et
al
.
(2
01
1)
R
S
cl
C
-t
y
p
e
2
64
0
3.
74
H
,
K
Im
R
ec
,
M
o
d
F
it
3
A
M
B
E
R
,
P
IO
N
IE
R
W
it
tk
ow
sk
i
et
al
.
(2
01
6)
R
S
cl
C
-t
y
p
e
27
00
3
.8
5
K
,
N
S
ed
M
o
d
el
in
g
4
M
ID
I,
V
IN
C
I
S
ac
u
to
et
al
.
(2
01
1
)
R
V
ol
C
-t
y
p
e
28
00
3.
85
N
M
o
d
A
tm
1
M
ID
I
R
au
et
al
.
(2
01
7)
R
U
V
ir
C
-t
y
p
e
30
00
3.
7
8
N
M
o
d
F
it
1
M
ID
I
R
au
et
al
.
(2
01
5
)
S
O
ri
M
-t
y
p
e
2
62
7
3.
79
J
,
H
,
K
M
o
d
A
tm
,
M
o
d
F
it
1
A
M
B
E
R
W
it
tk
ow
sk
i
et
al
.
(2
00
8)
S
O
ri
M
-t
y
p
e
26
27
3
.7
9
N
R
ad
T
ra
n
sf
M
o
d
13
M
ID
I
K
ar
ov
ic
ov
a
et
al
.
(2
01
3
)
T
L
ep
M
-t
y
p
e
28
00
N
A
J
,
H
Im
R
ec
,
M
o
d
F
it
1
A
M
B
E
R
L
e
B
ou
q
u
in
et
al
.
(2
00
9)
U
H
y
a
C
-t
y
p
e
32
0
0
3.
55
N
M
o
d
A
tm
1
M
ID
I
R
au
et
al
.
(2
01
7)
U
O
ri
M
-t
y
p
e
26
41
3.
85
H
Im
R
ec
1
IO
T
A
P
lu
zh
n
ik
et
al
.
(2
00
9)
W
H
y
a
M
-t
y
p
e
25
00
3.
73
K
P
ol
ar
Im
,M
o
d
F
it
1
A
M
B
E
R
O
h
n
ak
a
et
al
.
(2
01
6)
W
H
y
a
M
-t
y
p
e
25
00
3.
73
K
M
o
d
F
it
2
A
M
B
E
R
H
ad
ja
ra
et
al
.
(2
01
9
)
W
v
el
M
-t
y
p
e
27
0
0
3.
8
8
K
1D
,
3D
M
o
d
A
tm
2
A
M
B
E
R
W
it
tk
ow
sk
i
et
al
.
(2
01
6)
X
H
y
a
M
-t
y
p
e
27
00
3.
58
K
1D
,
3D
M
o
d
A
tm
2
A
M
B
E
R
W
it
tk
ow
sk
i
et
al
.
(2
01
6)
X
H
y
a
M
-t
y
p
e
27
00
3.
58
J
,
H
,
K
Im
ag
e,
m
o
d
el
in
g
of
V
2
1
A
M
B
E
R
H
au
b
oi
s
et
al
.
(2
01
5)
X
T
rA
C
-t
y
p
e
26
00
4.
0
0
N
M
o
d
A
tm
1
M
ID
I
R
au
et
al
.
(2
01
7
)
Y
P
av
C
-t
y
p
e
32
00
3.
55
N
M
o
d
A
tm
1
M
ID
I
R
au
et
al
.
(2
01
7)
Table 2.3: High-resolution observations of AGB stars (Continued).
56 CHAPTER 2. INTERFEROMETRY
3C
h
ap
te
r
Interferometric analysis of the
M-type Mira star R Car with
PIONIER-VLTI
***
This chapter is based on the paper:
A new dimension in the variability of AGB stars: convection patterns size
changes with pulsation recently sent to the A&A journal and which is in the review
process.
A. Rosales-Guzmán, J. Sanchez-Bermudez, C. Paladini, B. Freytag, M. Wittkowski,
A. Alberdi, F. Baron, J.-P. Berger, A. Chiavassa, S. Höfner, A. Jorissen, P. Kervella,
J.-B. Le Bouquin, P. Marigo, M. Montargès, M. Trabucchi, S. Tsvetkova, R. Schödel, S.
Van Eck
3.1 Our source: R Car
In this thesis, we analyze interferometric data of a prototypical M-type Mira star: R
Car. This source exhibits variations in the V band from +3.9 to +10.5 mag. The
proximity and brightness in the NIR and MIR of this source make it an ideal candidate
for observation with infrared interferometry. In the following section I will provide
relevant fundamental parameters of this star and a brief description on how they were
obtained. All the parameters are gathered in Table 3.1.
57
58 CHAPTER 3. H BAND INTERFEROMETRIC ANALYSIS OF R CAR
One of the first estimates of the angular size of R Car was carried out by Ireland
et al. (2004b). These authors obtained interferometric observations of the star in the V
band (between 700-900 nm), using the Masked Aperture-Plane Interference Telescope
(MAPPIT) system at the Anglo-Austrian Telescope (AAT). By fitting the V2 using a
Gaussian function, the authors derived an average FWHM of 20 mas in this band, with
a size increase to 40 mas around 712 nm, corresponding to the location of a TiO band.
McDonald et al. (2012) derived an effective temperature, Teff = 2800 K, luminosity,
L = 4164 L⊙, and period P = 314 days for R Car by comparing BT-Settl model atmo-
spheres (Allard et al. 2003) to SEDs created from different surveys such as Hipparcos
Tycho, Sloan Digital Sky Survey (SDSS), among others.
The radius of a pulsating star can be used to estimate the stellar mass when the
characteristic period, the pulsation constant, Q, is obtained from theoretical considera-
tions (Kamezaki et al. 2012). For instance, Xiong & Deng (2007) performed a detailed
theoretical study on the radial pulsation of red giants and provided the period ratio and
the characteristic period, Q for various models of giants. According to the authors, the
relation between the period ratio and the pulsation constant has the form:
QF = −0.178(Q1O/QF) + 0.171, (3.1)
where QF and Q1O are the pulsation constants for the fundamental and the 1st overtone
modes, respectively. The relation between the period and and the mean density of a
pulsating star can be written in the following form (e.g., Cox 2017):
P (ρ/ρ⊙)1/2 = Q, (3.2)
where P is the period of the star, ρ its mean density, ρ⊙ = 1.41g cm−3 is the mean
density of the sun, and Q is the pulsation constant 1. Using Eqs. 3.1 and 3.2, Takeuti
et al. (2013), estimated a mass M = 0.87 M⊙ for R Car.
Groenewegen et al. (1999) found a mass loss rate < 1.6 × 10−9 for R Car trough
model fitting of the SED in the NIR and from the observed CO (3-2) line intensities.
The observations that lead to these estimations were carried out with the IRAM 30 m
telescope in September 1996.
1Strictly speaking, Q is not constant. However, compared with other quantities that appear in
pulsation theory, it remains in a relatively narrow range, making it useful for deriving the mean density
of stars based on the pulsation period (Takeuti et al. 2013; Cox 2017).
3.2. INTRODUCTION TO THE PROBLEM OF CONVECTION IN AGB STARS 59
Parameter Value Ref.
Distance [Pc] 182 ± 16 Brown et al. (2021)
Period [Days] 314 ± 1 days Lebzelter et al. (2005); McDonald et al. (2012)
M [M⊙] 0.87 Takeuti et al. (2013)
Teff [K] 2800 McDonald et al. (2012)
Luminosity [L⊙] 4164 McDonald et al. (2012)
Ṁ [M⊙ yr−1] < 1.6 × 10−9 Groenewegen et al. (1999)
V band FWHM 15-20 mas Ireland et al. (2004b)
H, K band magnitudes -0.77, -1.23 Ducati (2002)
Table 3.1: Fundamental astrophysical parameters of R Car.
The photosphere of R Car was imaged with PIONIER, as part of the ”Interferometric
Imaging Contest” (Monnier et al. 2014). The data were obtained in January 2014 at
the pulsation phase2 φ = 0.78 ± 0.02. The reconstructed image shows a disk with an
angular diameter of about 10 mas (∼1.9 AU), and a couple of bright spots interpreted
as convective cells. No further analysis has been carried out on those images until now.
3.2 Introduction to the problem of convection in AGB
stars
The late phases of stellar evolution are largely affected by mass loss due to convection
and pulsation processes, which in turn, significantly affect the structure and evolution
of stars (Kupka 2004). In the case of convection, despite its importance, our current
knowledge is still limited. Observationally, we mainly rely on studying the convective
cells on the surface of the Sun. Schwarzschild (1975) suggested that the convection
length scale is proportional to the pressure scale height of a stellar atmosphere. Hence,
while the Sun has a couple of million cells on its surface, AGB stars will have (up to) a
few tens of them. However, testing this hypothesis directly on resolved stellar surfaces
has not been possible until very recently.
2Throughout the text we will refer to the pulsation phases of R Car obtained from the Visual light-
curve of the object (see Fig.3.1), and they are defined as φ = (t − t0)/T where t is the observing epoch,
t0 is a reference epoch, and T the period of the star.
60 CHAPTER 3. H BAND INTERFEROMETRIC ANALYSIS OF R CAR
High-angular resolution techniques are necessary to resolve those structures and IR
interferometers played an important role in the last decade measuring the diameters of
stellar disks and later on resolving smaller scale structures, which were related to con-
vection (see e.g., Le Bouquin et al. 2011; Wittkowski et al. 2017; Perrin et al. 2020). One
of the main challenges of studying convection in AGB stars is that their atmospheres are
composed of molecular gases (e.g., H2O, CO, HCN, etc.) and dust, obscuring the pho-
tosphere and making it difficult to observe the convective patterns. This is particularly
important for carbon-rich AGB stars as shown in the H band images by Wittkowski
et al. (2017). Oxygen-rich AGB stars are therefore better-suited targets to study con-
vection because the dust is mostly transparent at the H (λ0 = 1.76 µm) and K (λ0 =
2.2 µm) bands, thus, allowing us to probe the photosphere observationally.
Paladini et al. (2018) used the PIONIER (Precision Integrated-Optics Near-infrared
Imaging ExpeRiment; Le Bouquin et al. 2011) instrument to image in the H band the
stellar disk of the AGB star π1 Gru. The interferometric images revealed the presence
of bright spots that were associated with convective patterns. Those authors charac-
terized their size to compare it, for the first time, with predictions from mixing-length
theoretical models (e. g., Freytag et al. 1997; Trampedach et al. 2013; Tremblay et al.
2013). This reveals that the observational size of the convective structures agrees with
a scale proportional to ∼10 times the pressure scale height of the atmosphere. Climent
et al. (2020) and Norris et al. (2021) conducted similar comparisons for a couple of red
supergiant (RSG) stars. However, the estimated size of the convective patterns in RSG
stars did not follow the theoretical predictions extrapolated from the Sun’s convection.
In this work, we analyze the 2014 PIONIER R Car data together with more recent
data obtained in 2020. Our goal is to further characterize the structures on R Car’s
surface to connect them with the pulsation of the star. Sect. 3.3 describes our observa-
tions and data reduction process. Sect. 3.4 presents our image analyses and results. In
Sect. 3.5, we discuss our results in terms of the physical parameters of the star and the
angular scales obtained for the convective elements. Finally, we present our conclusions
in Chapter 5.
3.3 Observations and data reduction
Two datasets of observations were considered for our analysis. Both of them were taken
using the 4-telescope beam combiner PIONIER. The first one consists of archival data
3.4. ANALYSIS AND RESULTS 61
from February of 2014; the second one corresponds to observations taken between Jan-
uary and March of 20203 (see Tables 3.2 and 3.3). The observations were taken at
pulsation phases φ = 0.78 ± 0.02 and φ = 1.01 ± 0.10 for 2014 and 2020, respectively.
All data were taken with the 4 Auxiliary Telescopes (ATs). The baselines used
reached a minimum and maximum resolution of 23.5 and 1.2 mas (at λ0=1.67 µm and
in the super-resolution regime of λ/2Bmax; being B the baseline length), respectively.
The observations provided 6 spectral channels across the H band between 1.516 and
1.760 µm. The scientific data of R Car were interleaved with interferometric calibrators
(i.e., point-like sources with similar brightness as the science target and within a few
degrees of distance) to estimate the atmospheric and instrumental response function (the
calibrators are listed in Tables 3.2 and 3.3).
The calibrated interferometric observables, squared visibilities (V2), and closure
phases (CPs), were obtained using the pndrs package (Le Bouquin et al. 2011). The
final calibrated observables are the average of five consecutive five-minute object obser-
vations. Considering the calibration error in the observables, the V2 and CPs reached
a precision of σV2 ∼ 0.0025 and σCP ∼ 1 deg, respectively. The 2014 calibrated data
were obtained from the archive ”Optical Interferometry Data Base” (OiDB) supported
by the JMMC. The interferometric observables and the uv-coverage(s) of both epochs
are included in Fig 3.2.
3.4 Analysis and results
3.4.1 Diameter estimation: parametric uniform-disk model
To obtain the H band angular size of R Car, we applied the geometrical model of a
uniform disk (UD) to the V2 data. The visibility function of this model is defined by
the following equation (Berger & Segransan 2007):
VUD(u, v) = 2F r
J1(πρΘUD)
πρΘUD
, (3.3)
where J1 is the first-order Bessel function; ρ =
√
u2 + v2 where u and v are the spatial
frequencies sampled by the interferometric observations; ΘUD is the angular diameter of
3More data were obtained for the star in 2019, however such data cover only one array configuration
and are too far from the 2020 data in terms of pulsation phase to be combined with.
62 CHAPTER 3. H BAND INTERFEROMETRIC ANALYSIS OF R CAR
Table 3.2: Log of the observations for the 2014 data. HD 80603 (Diameter = 0.94 mas)
and HD 81502 (Diameter = 1.23 mas) were selected as calibrator stars.
No. Date Target Sta. Conf. τmin
0 - τmax
0 [ms] N. Obs.
HD 80603
1 2014-01-22 R Car D0-A1-C1-B2 2.18 - 8.8 98
HD 81502
HD 80603
2 2014-01-24 R Car D0-H0-G1 3.77 - 5.03 12
HD 81502
HD 80603
3 2014-01-27 R Car D0-G1-H0-I1 0.78 - 4.31 80
HD 81502
HD 80603
4 2014-01-28 R Car D0-G1-H0-I1 1.83 - 14.92 82
HD 81502
HD 80603
5 2014-01-29 R Car A1-G1-I1 2.58 - 10.33 52
HD 81502
HD 80603
6 2014-01-30 R Car A1-G1-J3 1.97 - 11.88 56
HD 81502
HD 80603
7 2014-02-02 R Car K0-A1-G1-J3 3.85 - 8.47 18
HD 81502
3.4. ANALYSIS AND RESULTS 63
Table 3.3: Log of the observations for the 2020 data. HD 80404 (Diameter = 1.866 mas)
and HD 71878 (Diameter = 2.949 mas) were selected as calibrator stars.
No. Date Target Sta. Conf. τmin
0 - τmax
0 [ms] N. Obs.
HD71878
1 2020-01-20 R Car A0-G1-J2-K0 17.0 - 21.1 10
HD 80404
HD71878
2 2020-01-21 R Car A0-G1-J2-K0 9.3 - 15.2 10
HD 80404
HD71878
3 2020-02-05 R Car A0-B2-D0-C1 3.2 - 3.8 20
HD 80404
HD71878
4 2020-02-14 R Car A0-B2-D0-C1 5.3 - 7.5 10
HD 80404
HD71878
5 2020-02-20 R Car K0-G2-D0-J3 9.0 - 14.5 10
HD 80404
HD71878
6 2020-03-20 R Car A0-B2-D0-C1 5.4 - 5.8 10
HD 80404
the uniform disk profile; and F r is the scaling factor that accounts for the over-resolved
flux in the observations. To fit the data, we used a MCMC algorithm based on the
Python package emcee (Foreman-Mackey et al. 2013). We let 250 independent chains
evolve for 1000 steps using the data of each one of the spectral bins independently.
Finally, we averaged the best-fit UD diameters. This gave us an estimate of Θ2014
UD =
10.23 ± 0.05 mas (1.86 ± 0.01 AU) for the 2014 epoch and Θ2020
UD = 13.77 ± 0.14 mas
(2.51 ± 0.23 AU) for the 2020 epoch. For instance, the diameter of the Sun is 9.35 ×
10−3 AU so, the derived sizes are approximately 270 times the diameter of the Sun!
64 CHAPTER 3. H BAND INTERFEROMETRIC ANALYSIS OF R CAR
Figure 3.1: Visual light curves of R Car obtained from the AAVSO database as func-
tion of the pulsation phase around the dates of our interferometric observations. The
blue stars and the black circles correspond to the data of the 2014 and 2020 epochs, re-
spectively. The vertical-red shaded area indicates the phase of the 2014 interferometric
observations, while the green line and the green, orange, and purple areas correspond to
the 2020 ones at different station configurations: short, long, and astrometric, respec-
tively (see Tables 3.2 and 3.3).
3.4.2 Asymmetries in the stellar disk: regularized image reconstruc-
tion
To characterize the asymmetric structure of R Car observed in the large CP variations,
we reconstructed aperture-synthesis images. For the 2014 image, we used the same setup
used by J. Sanchez-Bermudez and described in Monnier et al. (2014). Since this image
was already compared with other methods and reconstructions, confirming the validity
of the structures observed, we did not made any further changes to the imaging setup.
For the 2020 data, we recovered images using a combination of SQUEEZE (Baron
et al. 2010) and BSMEM (Baron & Young 2008). We recovered the images using a
pixel sampling of 0.3 mas/pixel, with a pixel grid of 256 × 256 pixels (which enclo-
sures a field-of-view of ∼ 77 mas). First, we ran SQUEEZE using 17000 iterations
over 1000 independent chains. The images were recovered independently at each of the
six wavelength channels. Two regularization functions were used: Total Variation, and
Compactness with hyperparameter values of one and twenty, respectively (see Sanchez-
3.4. ANALYSIS AND RESULTS 65
(a) (b)
Figure 3.2: PIONIER u-v coverages of the data sets used in this study. The left panel
corresponds to the 2014 data set, and the right one corresponds to the 2020 epoch. The
colors in the panels show the different wavelengths covered with the band-pass of the
instrument (see the labels on the plot).
Bermudez et al. 2018, for a more complete description of each regularizer). With this
setup, we ensure the convergence of most of the chains with reduced χ2 < 4. We aver-
aged the chains that converged into a single image. Finally, that mean image was used
as starting point for BSMEM (which uses entropy as regularizer), to smooth it without
increasing the residuals in the fitting of the observables. BSMEM converged in less than
100 iterations. Figure 3.3 shows our mean best-fit images from the 2014 and 2020 data
and Figs. 3.6 and 3.7 show the comparison between the observables extracted from the
reconstructed images and the data.
3.4.3 Power spectrum analysis
The H band images of R Car’s surface shown in Fig. 3.3 reveal the presence of several
large bright spots. The comparison between the two images shows that those features
change with time. We performed the analysis of the images’ power spectral density (PSD)
to estimate the characteristic size of the structures (Paladini et al. 2018). The first step is
to calculate the squared modulus of the Fourier Transform of the reconstructed images.
The power spectrum of an image shows the flux distribution per spatial scale or spatial
frequency, the smallest frequency being the one that encloses all the normalized intensity
66 CHAPTER 3. H BAND INTERFEROMETRIC ANALYSIS OF R CAR
(a) (b)
Figure 3.3: Best reconstructed images for the 2014 (left) and 2020 (right) datasets.
The white circle located in the right corner of the images corresponds to the maximum
resolution of the interferometer.
in the image; and the highest frequency is the one that maps the normalized intensity
of the most compact features. To remove the stellar disk contribution from the PSD,
we set the images’ pixel values below 70% and 90% of the intensity’s peak (for the 2014
and 2020 epochs, respectively) to the median flux of the pixels above the mentioned
threshold values. To characterize the weight of the power, P(k), at a given frequency, k,
we computed the momentum of the power-spectrum PSD(k) = k × P(k)/Ptot.
In order to estimate the loci of the spatial frequencies that dominate the PSDs, we
averaged them radially. Fig. 3.4 shows the resulting radial averaged PSDs for both
epochs. We fitted a Gaussian to the PSDs’ peaks to determine the dominant spatial
frequencies that trace the characteristic sizes of the structures in the surface of our object.
For the Gaussian fitting, we used the Python tool lmfit (Newville et al. 2016). Since
we only have a single image from the 2014 epoch, the reported uncertainties correspond
to the standard deviation of the Gaussian fit. For the 2020 data, the model fitting was
performed individually per spectral channel and, the reported errors correspond to the
mean of the standard deviation of the best-fit peak positions across the wavelengths.
The Gaussian center was limited within the range of -0.6 < log(k) < -0.2 for the 2014
epoch and -0.7 < log(k) < -0.3 for the 2020 epoch, respectively. The measured maxima
in the PSDs trace sizes of 1.82 ± 0.16 mas (4.94 ×1010 ± 4.89 ×109 m) and 2.51 ± 0.32
3.4. ANALYSIS AND RESULTS 67
Figure 3.4: Radial averaged power spectra of the reconstructed images. The red and
blue lines show the logarithm of the radial averaged power spectrum versus the logarithm
of the spatial frequencies (in cycles/mas) of the 2014 and 2020 data, respectively. We
applied a vertical shift of 0.5 to the 2020 powerspectra for better visualization. The gray
shaded area corresponds to spatial frequencies higher than the best resolution element
of the interferometer.
mas (6.82×1010 ± 9.59 ×109m) for the 2014 and 2020 epochs, respectively.
The derived sizes are bigger than the effective resolution elements of both epochs.
However, since the PSD of our images can be affected by the interferometer resolution
and the image reconstruction process, we quantify their impact on the determination
of the measured sizes of the structures. For that purpose, we simulated two models
with the following setups: (i) in the first one, we created two Uniform-Disk spots with a
diameter equal to the maximum angular resolution of the interferometer (1.2 mas); (ii)
in the second one, we generated two Uniform-Disk spots with a diameter equal to the
size of the surface structures derived from the 2020 image (1.8 mas).
We extracted the interferometric observables from these two models using the inter-
ferometric u-v coverage from the 2020 epoch. We used this epoch because those data
are sparser than the 2014 one, hence, we expect bigger u-v plane systematics on the
reconstructed images. We performed the reconstructions following the same process as
the original data. After convolving our images with the resolution of the interferometer,
68 CHAPTER 3. H BAND INTERFEROMETRIC ANALYSIS OF R CAR
and extracting the characteristic (dominant) size of the structures in the images using
the PSD analysis, we found that the sizes of the reconstructed spots were quite similar
to the original diameters used in the models: 1.30 ± 0.01 mas and 1.81 ± 0.03 mas for
the first and second simulated reconstructed images, respectively. The biggest effect is
observed in the reconstructed image of the first model, with a difference ∆ ∼ 0.1 mas,
when the simulated spot diameter is at the maximum resolution of the interferometer.
However, when the simulated spot diameter equals the minimum sizes of the structures
measured, the size is properly recovered with the PSD analysis. Hence, with those tests,
we confirmed that the effect of the interferometer resolution and u-v coverage do not
appear to have a significant effect on the measured sizes of the bright spots in the surface
of R Car with the used analysis.
3.5 Discussion
3.5.1 Sizes of the stellar disk and of the surface structures
The epochs analyzed were obtained on a time baseline larger than the pulsation period
of the star, hence, are part of different pulsation cycles. However, because of the stable
nature of the light-curve it is reasonable to assume that the star’s properties have not
changed fundamentally between the two observing epochs. In fact, the photometry
of the star is remarkably similar (see Figure. 3.1). This assumption will underlie our
interpretation and for the rest of the text we will speak in terms of pulsation phase
rather than years of observation.
The PSD analysis shows that the diameter of the star is smaller in phase 0.78,
and larger after the photometric maximum at phase 1.01. Ireland et al. (2005) made
predictions of diameter variations with pulsation for the broad H band. The model
predicts a shift between the size variation and the light-curve. The phase shift depends
on the band of observation. In H band, the object is expected to be smaller around phase
0.8, and larger before the minimum of the pulsation phase. The computed diameters
agree with this trend.
It is even more remarkable that the surface patterns are changing in size between
the two pulsation phases. On average, they are smaller when the star is smaller, and
larger when the star is larger. While this effect is expected from basic physics, to our
knowledge, this is the first time that it is shown for the same star using images and with
3.5. DISCUSSION 69
a quantitative analysis such as the PSD. These recovered PIONIER-VLTI images are,
thus, opening a new dimension in the study of the pulsation of stars. After several years
when optical interferometry measured diameter changes through the pulsation of the
AGB stars, the time and technology seems mature to measure the change in diameter
of surface structures throughout the pulsation cycle.
3.5.2 Convective pattern size(s) and their correlation with the stellar
parameters
Paladini et al. (2018) used interferometric observations to compare the structures on
the surface of the S-type AGB star π1 Gru with theoretical predictions derived from the
mixing-length theory of convection (Ulrich 1970). Those authors used empirical relations
of convective cell sizes obtained from 2D and 3D convective models calibrated with Solar
cell sizes and applied to different stellar types (Freytag et al. 1997; Trampedach et al.
2013; Tremblay et al. 2013). These authors found that the convective structures on the
surface of π1 Gru were consistent with the extrapolations of those convection models,
suggesting that the convective process could be similar for different stellar types across
the HR diagram.
To uncover the nature of the observed bright structures on the surface of R Car, we
conducted the same experiment and compared our measurements described in Sect. 3.4.3
with the extrapolations of the models used by Paladini et al. (2018). First, according to
Freytag et al. (1997) (FR97 from now on), the horizontal granulation size, xg, is related
to the characteristic pressure scale height, Hp, by the following expression:
xg = 10Hp , (3.4)
where Hp = RgasTeff/µ g, being Teff the effective temperature of the star; µ the mean-
molecular weight; g is the surface gravity of the star, and Rgas is the ideal gas constant.
Eq. (3.4) could be rewritten using the following logarithmic form:
log(xg,Freytag) = log(Teff) − log(µ) − log(g) + 0.92 . (3.5)
We calculated the predicted convective cell size of R Car using Teff = 2800 K (Mc-
Donald et al. 2012), Rgas = 8.314 × 107 erg K−1 mol−1 and µ = 1.3 g mol−1 for a
non-ionized solar mixture with 70% Hydrogen, 28% Helium and 2% heavier elements
70 CHAPTER 3. H BAND INTERFEROMETRIC ANALYSIS OF R CAR
(Carroll & Ostlie 2017). To estimate log(g), we used the mass reported in Takeuti et al.
(2013) of 0.87 M⊙. Considering the radii derived in Sect. 3.4.1, we estimated log(g)
values of -0.241 and -0.498 for the 2014 and 2020 epochs, respectively4.
Second, Trampedach et al. (2013) (TRA13 from now on) define the following expres-
sion to calculate the cell size:
log(xg,Trampedach) = (1.321 ± 0.004) log(Teff)
− (1.0970 ± 0.0003) log(g)
+ (0.031 ± 0.036).
(3.6)
The same Teff and log(g) values used in Eq. 3.5 were employed for this prescription.
Third, using the CIFIST grid, Ludwig et al. (2009) and Tremblay et al. (2013) (TRE13
from now on) derived the following equation for xg:
log(xg,Tremblay) = 1.75 log[Teff − 300 log(g)]
− log(g) + 0.05[Fe/H] − 1.87 .
(3.7)
In the latter formula we assumed for R Car solar metallicity [Fe/H]=0. The log-
arithm of the previously derived values of xg is reported in Table 3.4. Figure 3.5 also
displays the logarithm of the derived xg for the two reconstructed images (i.e. at dif-
ferent pulsation phases) as well as the three prescriptions. For comparison, we plot the
value of xg reported by Paladini et al. (2018) for π1Gru. In angular scales, the xg from
the models vary from 1 to 5 mas, being the lower and upper boundaries set by the
TRE13 and TRA13 prescriptions, respectively. The xg derived from our reconstructed
images have sizes of 1.82 ± 0.16 mas and 2.51 ± 0.32 mas for the 2014 and 2020 data,
respectively. This poses them in between the boundaries predicted by the models. Our
highly precise measurements are only possible due to our interferometric data. It is im-
portant to highlight that the derived sizes are larger than the minimum resolution of our
interferometer (1.3 mas), including the effects introduced by the image reconstruction
process.
4Notice that we also expect a change in temperature as a function of the radius for the two phases.
This effect should not be large given that the maximum variability in the H band is around 0.5 mag
(Ireland et al. 2004a). However, with the actual data, we could not quantify this change, so we used the
same effective temperature for the two data sets.
3.5. DISCUSSION 71
Table 3.4: Logarithm of the characteristic size of the patterns observed on the surface
R Car.
Model
2014
φ ∼ 0.78
2020
φ ∼ 1.01
log(xg,Freytag) 10.49 10.75
log(xg,Trampedach) 10.85 11.13
log(xg,Tremblay) 10.42 10.70
log(xg,PSD) 10.69 ± 0.04 10.83 ± 0.06
Our estimations of the characteristic size of the bright structures on the surface of
R Car also follow the predicted pattern from different pulsation phases, being smaller
when the star is contracted and larger when the star expands. Despite this agreement
with theoretical predictions, recent 3D convection models for AGB stars suggest that
convection acts in a more complex way than the FR97, TRA13, and TRE13 scale rela-
tionships suggest (see e.g., Colom i Bernadich 2020). For example, lower surface gravity
implies low densities both in the atmosphere and in the outer layers of the convection
zone. This effect leads to large convective velocities, which imply much more violent
convection when compared with, for example, the Sun. Then, large convective velocities
cause large amounts of material to overshoot above the top of the convective-unstable
layers. This causes the material to cool down and become optically thick due to molec-
ular opacities. In this case, what could be observed as bright structures on the surface
of the star are not the characteristic cell sizes but the material through the gaps created
between the elevated optically thick material.
This effect appears to be more important in more massive objects. For example,
Chiavassa et al. (2011) already suggested that the radiative energy transport in RSGs
is more efficient than in the Sun. RSGs are prone to higher pressure fluctuations than
AGBs. As a consequence, two kinds of surface structures might appear in those kinds
of stars: (i) small-scale, short-living convective features whose sizes are consistent with
the models used in Paladini et al. (2018) and in this work, but that the interferometers
can not resolve due to the (on average) more distant location of RSGs (as highlighted by
Climent et al. 2020) and; (ii) large-scale, long-living features whose sizes are consistent
with modified scale relationships predictions for which turbulent motions are included
(Norris et al. 2021; Chiavassa et al. 2011).
72 CHAPTER 3. H BAND INTERFEROMETRIC ANALYSIS OF R CAR
Still, the fact that for R Car and π1 Gru the measured bright spot sizes are consistent
with the extrapolations from the FR97, TRA13, and TRA13 scale relationships suggests
that we are observing the typical convective cell sizes, probably, because both AGB stars
have similar surface gravities and progenitor masses close to solar. Thus, their convective
processes are more similar to the ones of the Sun, which appears not to be the case for
more massive objects.
Figure 3.5: Logarithmic scale of the convective patterns (xg) versus effective temperature
(Teff). The different lines correspond to the predictions of convective cell sizes from
Freytag et al. (1997), Trampedach et al. (2013) and Tremblay et al. (2013) (see the
labels on the plot). The blue and red circles correspond to the characteristic scale of
the surface structures obtained directly from the reconstructed images of R Car, and the
green triangle shows the logarithmic value of xg reported by Paladini et al. (2018) for
π1Gru.
3.6 Appendix: Additional figures
In this section, we show the comparison between the synthetic and observed V2 and CPs
for our 2014 and 2020 epochs (see Sect. 3.4.2 ). In all figures, the red dots correspond
to the synthetic V2 and CPs extracted from the reconstructed images, while the data
are shown with grey dots.
3.6. APPENDIX: ADDITIONAL FIGURES 73
Figure 3.6: Comparison between the observables extracted from the 2014 mean re-
constructed image and the V2 (left panel) and CPs (right panel) data versus spatial
frequencies. The red dots correspond to the synthetic V2 and CPs extracted from the
reconstructed images, while the data are shown with gray dots.
Figure 3.7: Comparison between the observables from the 2020 reconstructed image to
the V2 (left panel) and CPs (right panel) data versus spatial frequencies. The red dots
correspond to the synthetic V2 and CPs extracted from the reconstructed images, while
the data are shown with gray dots.
74 CHAPTER 3. H BAND INTERFEROMETRIC ANALYSIS OF R CAR
4C
h
ap
te
r
Interferometric analysis of the
M-type Mira star R Car with
GRAVITY-VLTI
***
This chapter is based on the paper:
Imaging the innermost gaseous layers of the M-type Mira star R Car with
GRAVITY-VLTI
A. Rosales-Guzmán, J. Sanchez-Bermudez, C. Paladini, A. Alberdi, W. Brandner,
E. Cannon, G. González-Torà, X. Haubois, Th. Henning, P. Kervella, M. Montarges, G.
Perrin, R. Schödel, M. Wittkowski, 2023, A&A 674, A62
4.1 The problem of dust and wind formation in AGB stars
When1 low-to-intermediate-mass stars reach the AGB phase they are characterized by an
inert C and O nucleus surrounded by He and H burning shells (Höfner 2008). The early-
AGB phase begins with He-shell burning and subsequently evolves to a thermal pulse
1The calibrated GRAVITY data used in this work are available at the Optical Interferometric
Database of the Jean-Marie Mariotti Center (JMMC). Examples of the Python code used for optimizing
the uniform disk model described in Sec. 4.3.1, the MOLsphere model in Sect. 4.3.2, the bootstrapping
method described in Sect. 4.3.3, and the PCA analysis described in Sect. 4.3.3 are available to the
reader in this Github repository and in the appendix of this thesis.
75
76 CHAPTER 4. K BAND INTERFEROMETRIC ANALYSIS OF R CAR
AGB phase where the H-shell burns. At the beginning of this phase, the atmosphere
of the star is dominated by O-rich molecules (such as CO, TiO, SiO, H2O, and VO;
Wittkowski & Paladini 2014), but eventually (and depending on the mass of the star)
the third dredge-up takes place, transporting the C formed in the He burning shell to the
atmosphere. Depending on the efficiency of this dredge-up, the star will be of M-type
(C/O ∼ 0.5), C-type (C/O > 1), or an intermediate star (0.5 < C/O < 1) (Höfner &
Olofsson 2018).
Models of dust formation in AGBs suggest that pulsations are not sufficient to form
the observed outflows and mass-loss rates. However, it is believed that there is a momen-
tum transfer due to propagating shock-waves in the atmosphere that lift gas creating
density-enhanced layers around the star where the nucleation of dust can occur (Jeong
et al. 2003).
The models suggest that the mechanisms for driving dust differ depending on the
dust composition. C-type AGBs form amorphous carbon grains with large opacities that
can be accelerated by radiation pressure (Gautschy-Loidl et al. 2004; Sacuto et al. 2011;
Rau et al. 2015, 2017). For M-type AGB stars, the high chemical complexity around
them makes identifying the dust grains driving the dusty winds more complex. Theo-
retical models of the formation of dust-driven winds in M-type AGBs suggest that the
distance(s) from the photosphere to the innermost molecular layers (located at a few
stellar radii) determine how far the shock waves can levitate the gas. This is strongly
coupled with the condensation region that the dust of a given chemistry must be formed
to produce the observed winds (Nowotny et al. 2005a; Bladh & Höfner 2012; Liljegren
et al. 2016). Therefore, high-angular-resolution observations of the innermost molec-
ular layers are particularly important to constrain the gas-dust coupling and the dust
chemistry (see, e.g., Norris et al. 2012).
During the last decade, due to the enhanced spectral and spatial resolution made
available, near-infrared interferometry has become an important technique to map the
innermost structure of AGBs’ wind in order to link the observed morphology with the
mass-loss models. Most of the analyses based on interferometric data consist of model
fitting the observables directly (squared visibilities -V2- and closure phases -CPs-) with
parametric models of the star’s atmosphere and/or more sophisticated radiative transfer
and hydrodynamic simulations with codes like CODEX (Ireland et al. 2008, 2011) or
CO5BOLD (Freytag et al. 2002, 2003). Some examples of these studies include Ireland
4.1. THE PROBLEM OF DUST AND WIND FORMATION IN AGB STARS 77
et al. (2004b), Wittkowski et al. (2008, 2011, 2016), and Haubois et al. (2015). However,
due to the complexity of the environment, it is still difficult to fully reproduce the
asymmetries traced by the CPs of the interferometric data with those models. Therefore,
complementary techniques, such as aperture-synthesis image reconstruction can be used
to better depict the inner circumstellar layers of AGB stars. Examples of interferometric
imaging applied to AGB stars include the analysis presented by Wittkowski et al. (2017)
on the C-rich AGB star R Sculptoris using the H band (λ0 ∼ 1.76 µm) beam combiner
PIONIER (Le Bouquin et al. 2011) at the Very Large Telescope Interferometer (VLTI).
In that study, the reconstructed images are consistent with a dynamic atmosphere and
a stellar disk dominated by giant convection cells, which has a strong impact on the
distribution of the molecular layers and the formation of dust at a few stellar radii.
A further example includes interferometric imaging applied to the long period vari-
able (LPV) Mira using the Interferometric Optical Telescope Array (IOTA) in the H
band, as reported by Perrin et al. (2020). In that study, a large localized absorbing
patch is observed in the images. This patch is interpreted as a local dust formation
zone, and it is consistent with the combined effect of pulsation and convection in the
mass-loss process of AGBs (see also Paladini et al. 2018). These examples highlight the
importance of interferometric imaging to confront the observed morphologies with the-
oretical models of dust-driven winds and mass-loss processes in AGBs. However, there
is still limited observational information, in particular for the case of M-type AGBs.
In this work, we aim to characterize the innermost environment of R Car using in-
frared interferometry in order to link its morphology with models of the formation of
dust-driven winds. For that purpose, we present new K band (λ0 = 2.2µm) interfer-
ometric data obtained in 2018 with the instrument GRAVITY (Abuter et al. 2017) at
the VLTI. These data were analyzed to depict the innermost structure of R Car at two
phases within one pulsation cycle. Our analyses include: (i) parametric model fitting
of the squared visibilities (V2); (ii) image reconstruction of both the pseudo-continuum
and across the first (λ0 ∼ 2.29 µm) and second ( λ0 ∼ 2.32 µm) CO band heads; and
(iii) an estimation of the CO layer(s) physical parameters such as temperature, size, and
optical thickness by using a single-layer model. This chapter is structured as follows:
observations and data reduction are presented in Sec. 4.2; in Section 4.3, we describe
the results of our geometrical model fitting, our image reconstruction procedure, and
our single-layer model. Our results are discussed in Section 4.4. Finally, we present our
78 CHAPTER 4. K BAND INTERFEROMETRIC ANALYSIS OF R CAR
conclusions in Chapter 5.
4.2 Observations and data reduction
4.2.1 GRAVITY observations and interferometric data reduction
R Car was observed in January and February2, 2018 for a total of eight hours, with
the beam-combiner GRAVITY. The data were collected using six baselines with the
Auxiliary Telescopes (ATs) configuration A0-B2-C1-D0, which provides minimum and
maximum baselines of 7.67 m (equivalent to an angular resolution of ∆θ ∼ 59 mas in this
wavelength range) and 32 m (∆θ ∼ 15mas), respectively. The data cover a wavelength
range between 2.0 µm and 2.4 µm with a spectral resolution of R∼ 4000. The complete
log of observations is reported in Table 4.1.
We reduced and calibrated our data using the latest release of the GRAVITY pipeline
(version: 1.5.0). Detailed descriptions of the recipes are available in the ESO pipeline
manual and in (Lapeyrere et al. 2014). The same standard data reduction procedures
were applied to the scientific target and the calibrator star (HD 80404; which has an
angular size of 1.8 mas) to correct for the DARK, FLAT, and BAD calibrations, as well as
to extract the interferometric observables. The calibrated V2 of the target were obtained
following the standard procedure described in the ESO Gravity pipeline manual, which
consists of dividing the raw observables by the ones of the calibrator star. Calibrated
CPs were obtained by subtracting the raw CPs of the calibrator from the ones of the
source. Due to the observed changes in the observables at similar baselines, the data
were split into two epochs (January and February, with equivalent pulsation phases
of 0.56 and 0.66, respectively, see Fig. 4.1). The calibrated R Car data files can be
downloaded directly from the Optical Interferometric Data Base (OiDB) archive hosted
at the Jean-Marie Mariotti Center3.
4.2.2 GRAVITY spectrum calibration
GRAVITY data provide us with K band spectra (R∼4000) of the target and calibra-
tor. For the posterior analysis of the CO band heads observed in the R Car spectrum,
we calibrated and normalized it to remove the telluric contribution from the Earth’s
2We also have one observation night in March but we combined it with our February data.
3https://oidb.jmmc.fr/index.html
4.2. OBSERVATIONS AND DATA REDUCTION 79
Table 4.1: Log of the observations for our GRAVITY R Car data. HD80404 was selected
as a calibrator star. DIT and NDIT are the Detector Integration Time in seconds and
Number of Detector InTegrations, respectively. The coherence time (τ0) was taken from
the Paranal Astronomical Site Monitoring (ASM).
No. Date
Time
(UT)
Target DIT NDIT Seeing (”) τ0 (ms)
1 2018-01-27 05:30:21 R Car 3 100 0.76 5.4
2 05:42:30 R Car 3 100 0.63 5.5
3 06:25:34 R Car 3 100 0.81 5.2
4 06:37:45 R Car 3 100 0.89 6.0
5 06:51:46 HD80404 30 10 0.76 6.1
6 07:03:09 HD80404 30 10 0.99 5.8
7 07:29:58 R Car 3 100 0.73 6.3
8 07:42:04 R Car 3 100 0.67 5.0
9 07:56:18 HD80404 30 10 0.85 5.0
10 08:07:34 HD80404 30 10 0.90 4.3
11 08:22:08 R Car 3 100 0.50 4.5
12 08:34:15 R Car 3 100 0.58 4.4
13 2018-02-22 05:48:54 R Car 3 100 0.47 12.3
14 06:01:00 R Car 3 100 0.63 9.8
15 06:17:41 HD80404 30 10 0.52 8.6
16 06:29:02 HD80404 30 10 0.65 11.8
17 2018-02-25 05:35:44 R Car 3 100 0.33 11.6
18 05:47:51 R Car 3 100 0.49 10.8
19 06:02:42 HD80404 30 10 0.45 7.7
20 06:14:04 HD80404 30 10 0.55 6.2
21 2018-03-12 00:51:56 R Car 3 100 0.79 4.3
22 01:04:06 R Car 3 100 0.74 5.4
23 01:10:39 R Car 3 100 0.76 5.4
24 01:25:50 HD80404 30 10 0.74 6.0
25 01:37:13 HD80404 30 10 0.78 6.4
80 CHAPTER 4. K BAND INTERFEROMETRIC ANALYSIS OF R CAR
Figure 4.1: Visual light curve of R Car based on AAVSO database (blue points). The
red line represents the best fit of a sine function to the curve which results in a period of
P = 314 ± 1 days. The dashed vertical lines and the colored diamonds show the epochs
of our 2018 GRAVITY observations (see the labels on the plot).
atmosphere using the following procedure. First, we averaged the R Car and calibra-
tor spectra to obtain the high signal-to-noise spectra FRCar
avg and FHD 80404
avg , respectively.
To remove the telluric lines from the averaged R Car spectrum, we used the following
expression:
FRCar
λ =
FRCar
avg
FHD 80404
avg
× Fmodels
avg , (4.1)
where Fmodels
avg was calculated by averaging 15 theoretical spectra from the theoretical
spectra web server 4(Allard et al. 2010) of our calibrator spectral type A7 Ib. Second,
we flattened and normalized FRCar
λ (in the range of 2.21 - 2.34 µm) using a third-order
polynomial function. Finally, to ensure that the CO vibro-rotational transitions are
centered at their expected wavelengths, we applied a wavelength correction (equivalent
to a Doppler velocity of 17 kms−1) for the local standard of rest (LSR), taking into
account the proper motion and the distance of the source. The normalized R Car
spectrum of each epoch is presented in Fig. 4.2.
4http://svo2.cab.inta-csic.es/theory/newov2/
4.3. ANALYSIS AND RESULTS 81
Figure 4.2: R Car calibrated spectra and V2. Top panel shows the normalized calibrated
spectra of the GRAVITY data. The blue spectrum corresponds to the January epoch,
while the red one corresponds to the February one (see the labels on the plot). The
middle and bottom panels show plots of the V2 versus wavelength for the two epochs
observed. In these panels, the colors correspond to different baseline configurations.
4.3 Analysis and results
4.3.1 Geometrical model
To obtain the angular size of R Car across the K band, we applied a geometrical model
of a uniform disk (UD) to the V2 data. The visibility function of this model is given by
the following equation (Berger & Segransan 2007):
VUD(u, v) = 2(Fr)
J1(πρΘUD)
πρΘUD
, (4.2)
where ρ =
√
u2 + v2, u and v are the spatial frequencies sampled by the interferometric
observations, J1 is the first-order Bessel function, and ΘUD and Fr are the angular diam-
eter of the uniform disk profile and a scaling factor that accounts for the over-resolved
82 CHAPTER 4. K BAND INTERFEROMETRIC ANALYSIS OF R CAR
Table 4.2: Best-fit parameters of the uniform disk geometrical model.
λ[µm] ΘUD [mas] Fr
January February January February
2.016 20.44+0.17
−0.18 19.00+0.13
−0.13 0.808+0.010
−0.009 0.842+0.009
−0.008
2.069 19.71+0.12
−0.13 18.14+0.09
−0.10 0.910+0.008
−0.008 0.908+0.008
−0.007
2.122 18.37+0.08
−0.09 16.59+0.07
−0.07 0.940+0.006
−0.006 0.932+0.006
−0.006
2.175 17.36+0.07
−0.06 15.39+0.06
−0.06 0.954+0.006
−0.006 0.942+0.006
−0.006
2.228 16.67+0.06
−0.06 14.84+0.05
−0.05 0.964+0.005
−0.005 0.948+0.005
−0.006
2.281 16.97+0.06
−0.06 15.41+0.06
−0.05 0.945+0.005
−0.005 0.931+0.006
−0.005
2.333 19.12+0.08
−0.08 18.18+0.07
−0.08 0.932+0.006
−0.006 0.923+0.007
−0.007
2.386 20.83+0.10
−0.10 20.18+0.10
−0.10 0.928+0.006
−0.007 0.915+0.007
−0.007
flux in the observations, respectively. Prior to modeling the pseudo-continuum data, we
averaged 200 spectral bins into one single measurement. This allowed us to obtain eight
spectral channels across the 2.0 µm - 2.4 µm range. Considering the calibration error in
the averaged quantities, the final precision of the V2 is σV2 ∼ 0.01 and the precision of
the CPs σCP is ∼ 1 deg. Figure 4.3 shows the calibrated V2 and CPs as well as the u-v
coverage of the two epochs of observation for the eight spectral channels probed.
To fit the data, we used a Markov chain Monte Carlo (MCMC) algorithm based on
the Python package emcee (Foreman-Mackey et al. 2013). We let 250 walkers evolve for
150 steps using the data of each spectral bin independently. Table 4.2 shows the best fit
of ΘUD and Fr for each one of the observed epochs. Figure 4.4 shows the V2 data and
their corresponding best-fit models. Figure 4.5 shows the evolution of the best-fit ΘUD
versus wavelength of the two epochs of observation. From this figure, we observe two
significant changes in the pseudo-continuum structure of R Car.
First, the size of the target changes with wavelength, decreasing as we move toward
the center of the bandpass (λ ∼ 2.23 µm) and then increasing again. This trend is
expected due to the presence of molecules such as H2O at ∼ 2 µm and CO at ∼ 2.29-
2.33 µm lying above the stellar surface (see, e.g., Wittkowski et al. 2008). The smallest
ΘUD value (at ∼ 2.23 µm) corresponds to the pseudo-continuum-emitting region due to
4.3. ANALYSIS AND RESULTS 83
(a)
(b)
Figure 4.3: GRAVITY calibrated interferometric observables. The left plot shows the
V2 and CPs plotted versus spatial frequencies for the data sets collected and the right
one shows the corresponding u-v coverage. Panel (a) corresponds to the January epoch,
while panel (b) corresponds to the February one. The colors in the plots show the
different wavelengths (see the labels on the plot).
the photosphere. This size is the best estimate of the stellar disk diameter in the K
band. Ireland et al. (2004b) derived a wavelength-dependent Gaussian full width at half
maximum (FWHM) between 15 and 20 mas at 700-920 nm. Our best-fit UD diameter
at ∼ 2.23 µm indicates that the star appears to be larger near the V band than in the K
84 CHAPTER 4. K BAND INTERFEROMETRIC ANALYSIS OF R CAR
band. This size difference is due to the presence of dust layers obscuring the stellar-disk
in the V band, but that are transparent in the K band, allowing us to have a closer look
at the continuum-emitting region.
Figure 4.4: V2 of the pseudo-continuum versus spatial frequencies for both epochs. The
lines indicate the best-fit UD models. Each color corresponds to a different wavelength
(see the label on the plot). We highlight that the spatial frequencies increase in opposite
directions depending on the epoch. This allowed us to have a better visual comparison
of the trends present in the V2 data between the two epochs.
Second, the disk diameter is larger during the January epoch by up to 12% compared
with the February epoch. From the light curve reported in Fig. 4.1, it can be observed
that the V magnitude of the star increased in the February epoch, when the diameter of
the star is smaller. This behavior comes with an expected increment of the temperature
during February, which makes the star look brighter at the cost of it having a smaller
radius (see, e.g., Wittkowski et al. 2012, 2018; Haubois et al. 2015). A rough estimate
of the temperature difference (TJ/TF) between the observed epochs could be obtained
using the Steffan-Boltzmann law and the magnitude difference (|∆Vmag| = 1.8) from
Fig.4.1, resulting in TJanuary/TFebruary ∼ 0.62.
4.3. ANALYSIS AND RESULTS 85
Figure 4.5: Best-fit ΘUD as a function of wavelength for the two analyzed epochs. Each
color corresponds to a different epoch (see the labels on the plot).
4.3.2 Spherically thin layer: MOLsphere model
The GRAVITY spectra show prominent CO band heads (see Fig. 4.2). The data show
drops in the V2 of the CO 2-0 (λ0 ∼2.2946 µm) and CO 3-1 (λ0 ∼2.3240 µm) vibro-
rotational transitions compared with the pseudo-continuum, which indicates that the CO
line-emitting region is more extended (see, e.g., Wittkowski et al. 2012). The innermost
CO layers are important to trace the dust nucleation region since most of the carbon or
oxygen present in the atmosphere is concentrated in the stable bound of this molecule
(Höfner & Olofsson 2018).
To characterize the regions where these transitions originate, we employed a single-
layer model (called MOLsphere) that has proven to be successful in explaining the ab-
sorption of the CO molecule that appears in the visibilities of evolved stars (Perrin et al.
2004, 2005; Montargès et al. 2014; Rodŕıguez-Coira et al. 2021). This model consists of
a stellar disk with a compact layer around it. The star is modeled by a stellar surface of
radius R∗ that emits as a black body at a temperature of T∗. It is surrounded by a com-
pact spherical layer with a radius of RL that absorbs the radiation emitted by the star
and re-emits it like a black body. The MOLsphere is characterized by its temperature,
86 CHAPTER 4. K BAND INTERFEROMETRIC ANALYSIS OF R CAR
TL, radius, RL, and its optical depth, τλ. The region between the stellar photosphere
and the layer is assumed to be empty (see Fig. 4.6 for a sketch of the model). The
analytical expression of the model is given by
Ir
λ =

























Bλ(T∗)e(−τλ/cosθ)
+Bλ(TL)[1 − e(−τλ/cosθ)] if r ≤ R∗
Bλ(TL)[1 − e(−2τλ/cosθ)] if R∗ < r ≤ RL
0 otherwise,
, (4.3)
where Ir
λ = Ir
λ(T∗,TL,R∗,RL, τλ); T∗ and TL are the temperatures of the photosphere
and of the CO layer, respectively; R∗ and RL are the angular radius of the star and the
layer, respectively; τλ is the optical depth of the molecular layer at wavelength λ; Bλ(T)
is the Planck function (at wavelength λ and temperature T); and β is the angle between
the radius vector and the line-of-sight so that cosβ =
√
1 − (r/RL)2.
Figure 4.6: Illustration of single-layer model. The yellow disk represents the star, and
the orange ring represents the layer. The parameters in the image are as described in
Sect. 4.3.2. This image was adapted from Rodŕıguez-Coira et al. (2021).
For the fitting, we adopted a T∗=2800 K and R∗=5 mas (Monnier et al. 2014;
McDonald et al. 2012). We used the disk estimation reported in the H band because
it is considerably less affected by molecular contribution, in contrast with the K band
where molecules such as CO, H2O, and OH, among others (see, e.g., Wittkowski et al.
4.3. ANALYSIS AND RESULTS 87
2018; Paladini 2011), are present over the entire band. The unknown parameters are
then TL, RL, and τL. We performed a two-step procedure to estimate these parameters.
As a first step, for each pair of TL and RL, we performed a least-squares minimization
to find the best-fit value of τL that reproduces the spectrum FRCar
λ . The next step
consisted of using an MCMC method based on the Python library emcee to estimate the
best combinations of TL, RL, and τL that reproduce the V2.
Figure 4.7: Best fit to spectra with the single-layer model for the two CO band heads
and epochs. The gray dots correspond to the data points and the solid-lines correspond
to the best-fit model (see the labels on the plot).
To account for the over-resolved flux when modeling the V2 data, we added a scaling
factor of Fr. This value was estimated simultaneously with the other parameters in
the model. Table 4.3 shows the corresponding results for the good quality spectral
channels across the two CO band heads. Figures 4.7 to 4.8 4.9 show the best fit to R
Car’s spectra and V2 obtained with the single-layer model. From the best-fit models,
we derive a temperature range between 1300 and 1500 K, which is consistent with the
expected temperature values to produce the CO vibro-rotational transitions observed
in the K band (Geballe et al. 2007). The values of the optical depth, τL, support that
the layer is optically thick across all the spectral channels, with higher τL toward the
centers of the band heads. In Section 4.4, we discuss the implications of these results in
the context of the wind-driving dust candidates that can coexist with the CO.
88 CHAPTER 4. K BAND INTERFEROMETRIC ANALYSIS OF R CAR
Figure 4.8: Best fit to V2 from the best single-layer models across the first CO band
head for our two epochs. The colored dots and the shaded regions correspond to the
best-fit models and their corresponding error-bars. Each color correspond to a different
wavelength (see labels on the plot). The data are shown with gray dots. We highlight
that the spatial frequencies increase in opposite directions depending on the epoch.
This allowed us to have a better visual comparison of the trends present in the V2 data
between the two epochs. Also, for better visualization, we displaced the V2 vertically.
4.3.3 Image reconstruction
Regularized reconstructed images
The geometrical model described in Sect. 4.3.1 and the MOLsphere model in Sect.
4.3.2 show that the source is symmetric. However, the CPs observed in our data differ
considerably from zero or 180 degrees, which indicates that the source is asymmetric (see
Fig. 4.3). So, in order to depict the complex nature of the source across the K band and
at the CO band heads, we reconstructed aperture-synthesis images using the software
BSMEM (Baron & Young 2008). This code uses an iterative regularized minimization
algorithm of the following form:
f(x) = fdata(x) + µfprior(x) , (4.4)
4.3. ANALYSIS AND RESULTS 89
Figure 4.9: Best fit to V2 from the best single-layer models across the second CO band
head for our two epochs. The panels are as described in Fig. 4.8.
where fdata(x) is a measurement of the log-likelihood and it characterizes the discrepancy
between the image model (at a given iteration) and the observables; fprior(x) contains
the prior information included in the reconstruction, and the hyperparameter µ is used
to adjust the relative weight of the constraints between the measurements and the ones
set by the priors (Renard et al. 2011). BSMEM uses Entropy (Gull & Skilling 1984) as
regularizer in fprior(x).
As initial setup for our reconstruction, we created images of 128 × 128 pixels with a
pixel scale of 0.469 mas. The best-fit sizes of the UD models were used as starting point
in the reconstruction process. We performed the reconstructions by taking all the V2 and
CPs on each individual channel separately. Since our data are constrained by a relatively
sparse u-v coverage, this could affect the quality of the reconstructed images. Therefore,
to properly distinguish intrinsic astrophysical features from image reconstruction arti-
facts (Haubois et al. 2015; Sanchez-Bermudez et al. 2018) and to provide an estimate on
the quality of our reconstructed images, we employed a bootstrapping method (Zoubir
& Iskander 2004).
We built new data samples from the original one. The samples were created by
keeping the original number of data points unaltered, but allowing random sampling
90 CHAPTER 4. K BAND INTERFEROMETRIC ANALYSIS OF R CAR
Table 4.3: Parameters of the single-layer model for the two CO band heads and our two
epochs.
λ[µm] TL [K] RL [mas] τL Fr TL [K] RL [mas] τL Fr
January February
CO 2-0
2.2941 1316+136
−73 13.02+1.33
−1.73 1.3+0.3
−0.6 0.968+0.029
−0.030 1320+94
−60 12.61+1.09
−1.27 0.9+0.2
−0.5 0.973+0.024
−0.018
2.2944 1344 +90
−105 12.41 +1.65
−1.05 1.9 +0.5
−0.8 0.976+0.028
−0.029 1366 +120
−90 12.38 +1.35
−1.28 1.6 +0.5
−0.6 1.005 +0.033
−0.024
2.2946 1360 +106
−90 12.17 +1.22
−1.21 1.9 +0.6
−0.9 0.944 +0.027
−0.029 1413 +103
−117 12.12 +1.60
−1.09 1.7 +0.6
−0.8 0.956 +0.031
−0.032
2.2952 1358 +128
−84 12.22 +1.25
−1.34 1.6 +0.5
−0.8 0.981 +0.030
−0.026 1413 +129
−96 11.95 +1.26
−1.24 1.5 +0.6
−0.8 0.981 +0.025
−0.019
2.2954 1389 +95
−83 11.90 +1.15
−1.02 1.7 +0.5
−0.6 0.976 +0.031
−0.026 1406 +108
−83 12.07 +1.26
−1.03 1.4 +0.4
−0.6 0.989 +0.028
−0.022
2.2962 1448 +92
−111 11.28 +1.23
−0.83 1.7 +0.6
−0.9 0.962 +0.026
−0.025 1394 +195
−88 12.23 +1.39
−1.78 1.2 +0.4
−0.7 0.976 +0.029
−0.033
2.2968 1424 +148
−93 11.69 +1.17
−1.30 1.4 +0.4
−0.7 0.987 +0.023
−0.024 1504 +171
−126 10.95 +1.44
−1.18 1.3 +0.4
−0.7 0.972 +0.022
−0.023
CO 3-1
2.3229 1420 +123
−94 11.41+1.19
−1.05 1.5+0.5
−0.7 0.958 +0.026
−0.022 1467+159
−119 11.04+1.33
−1.20 1.2+0.4
−0.7 0.960 +0.021
−0.025
2.3232 1316 +131
−75 12.90 +1.43
−1.45 1.5 +0.3
−0.6 1.010 +0.027
−0.030 1406 +152
−100 11.79 +1.42
−1.40 1.2 +0.3
−0.6 0.996 +0.024
−0.025
2.3237 1309 +82
−76 12.78 +1.20
−1.07 2.1 +0.5
−0.7 0.958 +0.032
−0.033 1336 +86
−73 13.33 +1.18
−1.09 1.9 +0.5
−0.7 0.990 +0.028
−0.026
2.3240 1315 +81
−86 12.67 +1.30
−1.09 2.1 +0.6
−0.8 0.969 +0.030
−0.035 1351 +102
−79 12.80 +1.31
−1.15 1.8 +0.5
−0.7 0.971 +0.034
−0.027
2.3242 1330 +91
−81 12.18 +1.15
−1.07 1.9 +0.6
−0.8 0.960 +0.027
−0.027 1352 +116
−110 12.55 +1.88
−1.29 1.8 +0.5
−0.9 0.981 +0.034
−0.037
2.3245 1394 +112
−113 11.52 +1.51
−1.07 1.7 +0.6
−0.7 0.964 +0.026
−0.025 1401 +122
−106 12.02 +1.39
−1.17 1.5 +0.5
−0.8 0.969 +0.024
−0.024
2.3248 1393 +132
−88 11.83 +1.24
−1.26 1.4 +0.4
−0.6 0.995 +0.025
−0.024 1452 +142
−104 11.45 +1.25
−1.15 1.3 +0.4
−0.6 0.979 +0.021
−0.022
with replacement. This created data sets with points with higher weights than in the
original data set and some others with zero weights. This algorithm also produces slight
changes in the u-v plane, which allows us to trace the impact of the u-v coverage on the
reconstructed images, without loosing the statistical significance of the original number
of data points. According to Babu & Singh (1983), the statistical moments such as mean,
variance, and standard deviation of the bootstrapped samples are good approximations
to the statistical moments of the original data sets.
We employed this method by creating 50 different sampled data sets from the original
4.3. ANALYSIS AND RESULTS 91
observables, V2∗ = V2
1, ...,V
2
50 and CP∗ = CP1, ...,CP50, using the Python scikit
learn package (Pedregosa et al. 2011). We performed image reconstructions on each
sample to obtain a set of bootstrapped images: I∗ = I1..., I50. All the reconstructions
from the samples converged, and the best-fit image reproduces the trend in the sampled
observables (V2 and CPs) quite well. We calculated the mean and standard deviation
per pixel along the 50 bootstrapped images to obtain a mean image Imean and a standard
deviation map Iσ. We used Imean/Iσ to estimate the signal-to-noise ratio (S/N) of the
structures present in the mean image. We removed from the images all the features with
an S/N < 3 as they are related to reconstruction artifacts. Figs. 4.10 - 4.12 show the
Imean images for the pseudo-continuum and CO band head channels recovered. All the
observed structures in the Imean images lie within the three S/N contours. Figs. 4.20 to
4.22 show the comparison between the synthetic observables, obtained from the Imean
images, and the original data points. As can be seen, the trends are well-reproduced at
all angular scales and wavelengths.
From a visual inspection of our best reconstructed images, Imean (Figs. 4.10 - 4.12),
we note that all the pseudo-continuum images are asymmetric and show structural differ-
ences between the observed epochs. In our January images, we observe a spot toward the
south-east; while in February, a more elongated (and centered) structure is observed. Ex-
cept the change in size across the band pass, we do not see significant structural changes
with wavelength.
92 CHAPTER 4. K BAND INTERFEROMETRIC ANALYSIS OF R CAR
(a)
(b)
Figure 4.10: Best pseudo-continuum reconstructed images for the (a) January and (b)
February epochs. The white ellipse located in the right corner of the image at 2.2280
µm corresponds to the mean synthesized primary beam. The blue contours represent
10, 30, 50, 70 90, 95, 97, and 99 % of the intensity’s peak.
4.3. ANALYSIS AND RESULTS 93
(a)
(b)
Figure 4.11: Best reconstructed images per spectral channel for the (a) January and (b)
February epochs across the first CO band head. The white ellipse located at the right
corner of the image at 2.2954 µm and the blue contours are as described in Figure 4.10.
94 CHAPTER 4. K BAND INTERFEROMETRIC ANALYSIS OF R CAR
(a)
(b)
Figure 4.12: Best reconstructed images per spectral channel for the (a) January and (b)
February epochs across the second CO band head. The white ellipse located in the right
corner of the image at 2.3242 µm and the blue contours are as described in Figure 4.10.
4.3. ANALYSIS AND RESULTS 95
Principal component analysis to understand the structure of R Car
To have a more precise characterization of the asymmetries in the reconstructed images
of the CO band heads, we used the principal component analysis (PCA) described by
Medeiros et al. (2018). Those authors demonstrated that the visibilities of the principal
components are equal to those of the visibilities. This method is useful for tracing the
changes across a set of images that have the greatest effect on the visibility (amplitude
and closure phase) profile. In our case, we estimate the most significant structural
changes of the observed asymmetric structures across wavelength for each of the observed
epochs. The following procedure was applied to the ensemble of wavelength-dependent
images per epoch to extract their principal components.
For the PCA analysis, we normalized each data set by subtracting the corresponding
mean image and dividing by their standard deviation image across the wavelength range.
To perform the PCA analysis, we used the CASSINI-PCA5 package. This software allowed
us to compute the covariance matrix of the data set and to transform it into the space of
the principal components. This allowed us to determine the eigenvectors (or eigenimages
in our case) and their corresponding eigenvalues. Since we only have seven images per
data set, and the possible number of components must be smaller than or equal to
the number of images in the data set, we decided to only keep the first four principal
components. From our tests, we observe that those components explain, at least 93% of
the variance in the data sets.
The maps that correspond to the first four principal components of the CO band
heads are reported in Figs. 4.13 and 4.14. For both CO data sets, we observe that the
first two components trace the variance associated with the general changes in brightness
and size observed in the images of the CO band heads across wavelengths. However,
components 3-4 show considerably more asymmetric changes in the external structure
of the target (for contours between 50% and 10% of the peak). The observed varying
structures in components 3-4 are evenly distributed across the different position angles in
the structure of the source. We suspect that the variance explained by these components
is associated with the more asymmetric structure in the CO images. This could also
be connected with the possibility of having a clumpy CO distribution of material in
the envelope around R Car. As an example of the contribution of each component to
the asymmetric structure of the target, in Figure 4.15 we include the CPs obtained
5https://github.com/cosmosz5/CASSINI
96 CHAPTER 4. K BAND INTERFEROMETRIC ANALYSIS OF R CAR
at a wavelength of 2.2946 µm from the recovered images with cumulative components
between 1 and 4. As can be observed, the recovered images improve the recovery of
the asymmetries traced by the CPs as we add more components to them. This is an
indication that the variance observed in the eigenimages is tracing a more changing
environment in the CO distribution across the sampled wavelengths. Unfortunately, the
limited resolution of our observations is not enough to properly resolve these asymmetries
completely.
4.3.4 Pure-line CO images
Our CO reconstructed images contain a contribution associated with the underlying
stellar component. To remove such contribution, we first centered the images across the
CO band heads in a common reference position. For this purpose, we used the differential
phases (DPs) in the GRAVITY data. DPs measure the photocenter displacement of an
emission line with respect to the continuum for a given baseline length and orientation.
We only restricted the pure-line CO images to the first band head, since the DPs of the
second band head have much larger uncertainties. To compute the global 2D photocenter
displacement vector at a given spectral channel, we used the following relation (Kraus
2012):
−→p = − φi
2π
· λ−→
Bi
, (4.5)
where φi is the GRAVITY differential phase measured on baseline i,
−→
Bi is the corre-
sponding baseline vector [Bx, By], and λ is the central wavelength of the working spectral
channel. By solving the matrix represented by Eq. 4.5, we were able to find the global
[px, py] displacements of the emission line. We solved the system of equations employing
the least-squares algorithm implemented in the Python function numpy.linalg.lstsq.
This function inverts the matrix formed by each one of the observed differential phases
per baseline per channel to find the solution of the astrometric displacements. To obtain
an estimation of the precision (i.e., 1σ error-bars) of the reported displacements, we ap-
plied the same matrix inversion process to 100 different samples obtained with random
points using a Gaussian distribution with a standard deviation equivalent to the error
bars of the differential phases. The final astrometric displacements are equivalent to the
mean of the 100 different samples, while their error bars are equivalent to the standard
deviation of the displacements from the samples.We obtained the displacements for the
4.3. ANALYSIS AND RESULTS 97
different spectral channels across the CO band heads, as well as for a few continuum
channels. Fig. 4.16 shows the first CO band head offsets for our two epochs.
We can observe that the wavelengths corresponding to the continuum emission lie
close to the reference position [0,0], and the ones of the band head show displacements
within 1 mas. It is interesting to highlight that the CO centroid position of the two
epochs is different, confirming the variability of the structure of the source between
them. Once the centroid position of the images was computed, we used them to extract
the pure-line CO images using the following method. First, CPs are shift-invariant
quantities, it means that images reconstructed with this observable could be placed at
different positions in the pixel grid and, still, reproduce the CPs in the data. Therefore,
to account for this effect, we computed the continuum and CO images, and we performed
a sub-pixel shift to the corresponding centroid position obtained from the astrometric
displacements [px, py] of the differential phases. For the subsequent analysis, we took
the continuum images at λcont = 2.228 µm. Second, once all the centroids of the images
matched with the astrometric displacements, we subtracted the continuum image from
each one of the CO images across the different spectral channels. Since the total flux
in the reconstructed images is normalized to unity, before performing the subtraction
we calibrated the total flux on each image by scaling them with their corresponding
normalized flux scale from the spectra in Figure 4.2. Figure 4.17 shows the CO images
with the continuum contribution subtracted. We masked the central region to better
visualise the CO-emitting region. From a visual inspection of the images, we can identify
some features. First, there are different structures between wavelengths and epochs,
suggesting the inhomogeneous nature of CO. The images also show bright and faint
regions at different position angles. We relate them with zones of lower (faint regions)
and higher (bright regions) densities and temperatures.
98 CHAPTER 4. K BAND INTERFEROMETRIC ANALYSIS OF R CAR
(a)
(b)
Figure 4.17: Pure-line CO images of the first CO band head for the (a) January and (b)
February epochs. The region that corresponds to the subtracted continuum contribution
has been masked for better visualization of the CO region. The white ellipse located in
the lower right corner of the image at 2.2941 µm corresponds to the mean synthesized
primary beam.
4.4 Discussion
4.4.1 CO column densities from the single-layer model
The single-layer model described in Section 4.3.2 allowed us to obtain information about
the CO layer’s physical conditions (e.g., opacity, temperature, and sizes). Using the
optical thickness and temperature derived for the layer, we estimated the corresponding
column density in units of cm−2 using the following relation (Savage & Sembach 1991):
N =
τL∆v
√
π
fijλijc
πe2
mec2
, (4.6)
where τL is the optical thickness obtained from our single-layer model, ∆v is the width
of the band head in km s−1, fij is the absorption oscillator strength that depends on
4.4. DISCUSSION 99
the temperature, λij is the wavelength of the transition, c is the speed of light, e− is
the electron charge and me its mass, and the constant term π e2/me c2 = 8.8523×10−13
cm/mol. We calculated the column density at the center of the first CO band head, where
the optical thickness has its largest value, by setting λij = 2.2946 µm. Then, using Eq.
(3) from Goorvitch (1994), we calculated fij for the vibro-rotational transition V0−>2,
J20−>21 (V and J are the vibrational and rotational states, respectively). We obtained
their corresponding Einstein coefficients, A, from Chandra et al. (1996). Finally, we
calculated ∆v by fitting an inverse Lorentzian profile across the first CO band head for
the two epochs. Since ∆v is expressed in kms−1, ∆v = c∆λ/λ0.
The obtained column densities for both epochs lie in theNCO ∼ 9.19×1018 - 1.01×1019
cm−2 range. This result has two implications. On one side, these values are consistent
with the range of column densities reported by authors including Rodŕıguez-Coira et al.
(2021) and Tsuji (2006) for other evolved stars. Furthermore, these values, on average,
are high enough to shield the dust particles allowing them to grow to micron-sized scales,
where the stellar radiation could act efficiently to drive the observed winds in this kind
of object. On the other side, our current single-layer model returns large uncertainties
in the estimated parameters. Therefore, we cannot confirm if CO formation happened
between both epochs.
To obtain a direct comparison between the size of the photosphere (R∗) and the size
of the CO layer (RL), we calculated the ratio RL/R∗. We take the R∗ value reported by
Monnier et al. (2014) of 5 mas as a reference. From RL/R∗, we find a mean extension
of the layer RL ∼ 2.4 R∗. To verify the consistency of our results, we compared them
with other O-rich Miras for which the single-layer model has been applied. Table 4.4
shows the values of the layer temperature and RL/R∗ for the Mira stars R Leo, χ Cyg,
and T Cep, as well as for our R Car estimates. From the values reported, we confirm
that the parameters derived for R Car are consistent with other estimates reported in
the literature for similar type of objects.
4.4.2 Observational restrictions to the wind driving candidates for M-
type Mira stars
Given the high chemical complexity around M-type AGB stars, it is not easy to determine
the chemical composition of dust-driven winds. Even though observations tell us about
the existence of different dust species in the circumstellar environment (Molster et al.
100 CHAPTER 4. K BAND INTERFEROMETRIC ANALYSIS OF R CAR
Table 4.4: Parameters of the single-layer model for other AGBs.
Star T∗ (K) TL (K) RL/R∗
R Leo 3856±119 1598±24 2.29±0.19
χ Cyg 3211±158 1737±53 1.91±0.02
T Cep 3158±158 1685±53 2.15±0.02
R Car Jan (CO 2-0) 2800 1360 +106
−90 2.43 ± 0.20
R Car Feb (CO 2-0) 2800 1406 +108
−83 2.42 ± 0.20
R Car Jan (CO 3-1) 2800 1330 +91
−81 2.44 ± 0.20
R Car Feb (CO 3-1) 2800 1401 +122
−106 2.40 ± 0.20
Parameters of the singe-layer model extracted from Perrin
et al. (2004) for two O-rich AGB stars (R Leo and T Cep) and
one S-type (χ Cyg) compared with the parameters reported
in this work for the CO band heads and GRAVITY epochs
analyzed.
2010), not all are efficient in absorbing or scattering the star’s radiation to trigger a
stellar wind. A condition for a dust grain to be considered a wind driver is that it must
be formed closer to the stellar surface within reach of shock waves produced by the star.
The shortest distance from the star at which a dust grain can be formed depends on its
chemical and optical properties (Bladh & Höfner 2012).
As mentioned in Sect. 4.1, the observed CO vibro-rotational transitions (with ∆V =
2) originate in the deep photospheric layers where dust is forming. Therefore, with the
CO layer parameters, obtained from our single-layer model and the pure-line CO images,
we can set constraints on the types of dust grains that are considered wind-driving
candidates and that can coexist with similar physical conditions to the CO innermost
layers detected with our GRAVITY observations. It is important to remember from
Sec. 1.3.2 that for winds to be triggered, the star pulsations must levitate the gas to a
certain distance (defined as levitation distance, Rl, Eq. 1.6) where dust can condense
(known as condensation distance, Rc, Eq. 1.5). The levitation distance can be extracted
directly from the size of the layer from our single-layer model or from the pure-line
CO-reconstructed images. The single-layer model provides us an average estimate of
4.4. DISCUSSION 101
Rl = 2.43 R∗, since the model only considers a symmetric structure. However, the
reconstructed images provide us a range of Rl between 2.25 and 2.72 R∗ due to the
observed asymmetries. Assuming uesc = 44 km s−1, M = 0.87 M⊙ (Takeuti et al. 2013),
Teff = 2800 K (McDonald et al. 2012), and R∗ = 167 R⊙ (Monnier et al. 2014) for
Eq. 1.6 and that the CO band heads that we observed are formed in a region from the
stellar surface at R0 = 2 R∗, we computed an initial velocity of the CO of u0 = 13
km s−1 considering the single-layer model and u0 between 10 and 16 km s−1 using the
CO reconstructed images. These values are consistent with the ones reported in the
literature on the order of u0 ∼ 10 km s−1 (Nowotny et al. 2010).
Using Eq. 1.5, we estimated the condensation radii for several types of dust grains
assuming the condensation coefficients and temperatures reported by Bladh & Höfner
(2012) and references therein (see Sec. 1.3.2). Considering that the levitation radius has
to be larger or equal to the condensation radius, we obtained curves of possible p and Tc
values with the Rl estimates from our single-layer model and pure-line CO reconstructed
images.
Figures 4.18 and 4.19 show the Tc versus p, diagram. In both figures, the red dots
correspond to a variety of dust grains with different condensation temperatures and
values of p. The solid-black line in Fig. 4.18 displays the values of Tc and p for the
Rl value obtained from the single-layer model. The colored solid lines in Fig. 4.19
correspond to the Rl values extracted from the pure-line CO images for both observed
epochs at different orientations: 0◦, 90◦, 180◦, and 270◦. Since we have seven images
per epoch, we adopted the median value per image data set at each orientation.
Those figures show that our CO layer radius is consistent with the condensation ra-
dius for composites Mg2SiO4 and MgSiO3. However, we find layer temperatures around
1300 ± 100 K (blue and green shaded regions on Fig.4.18), which are slightly higher
than the expected condensation temperature for those two dust composites (Tc ∼ 1100
K). This temperature discrepancy could be explained by the fact that the single-layer
model considers a homogeneous and symmetric distribution of the CO layer. However,
the reconstructed images clearly show that the CO distribution is neither homogeneous
nor symmetric as shown in Fig. 4.19 where, depending on the position angle, either one
or both composites lie(s) below the estimated levitation radius. Therefore, regions with
lower temperatures and larger levitation radii could exist around R Car. This supports
the hypothesis that dust formation could occur in an anisotropic way over the CO shell
102 CHAPTER 4. K BAND INTERFEROMETRIC ANALYSIS OF R CAR
extension.
Magnesium-iron silicates (olivines and pyroxenes) with nano grain sizes (∼ 0.1 µas)
have already been detected in the circumstellar environments of M-type AGB stars
(Norris et al. 2012; Ohnaka et al. 2017), and they appear to be the main factor responsible
for the dust-driven winds observed in these stars. While our GRAVITY observations and
models suggest the possible coexistence of those dust types at the inner gaseous (CO)
shells, we require additional interferometric observations in the mid-infrared (potentially
with the beam-combiner MATISSE-VLTI) to confirm this scenario. Very recent models
of dusty winds in M-type AGBs suggest Fe enrichment of the wind as it moves far
away from the photosphere (Höfner et al. 2022). Those models suggest an inner core
of Mg2SiO4 (or of a similar composite) surrounded by more outer layers of MgFeSiO4.
New MATISSE-VLTI interferometric observations could help us to confirm this scenario
in the proximity of R Car in the near future. By confirming this scenario, mid-infrared
interferometric observations could shed new light on the problem of a lack of the 10 µm
silicate feature in the spectra of M-type AGB stars.
4.4.3 The role of CO asymmetries on mass-loss estimates
The mass-loss rate is a crucial parameter for understanding the evolution of M-type AGB
stars. Bladh et al. (2019) presented an extensive grid of DARWIN dynamic atmosphere
and wind models to establish the mass-loss rates, wind velocities, and grain sizes for M-
type AGB stars. These models assume that the principal mechanism for wind triggering
is the radiation scattering on Mg2SiO4 composites. This grid of models produce mass-
loss rates on the order of 1×10−6 M⊙ yr−1 for stellar parameters similar to R Car (e.g.,
M∗ = 0.75 M⊙; log[L∗/L⊙] = 3.55; P = 283 days).
More importantly, those models show that, for the mass (0.75 - 1 M⊙) and lumi-
nosity (log[L∗/L⊙] = 3.7) ranges of R Car, large values of the pulsation amplitude in
the models and high chemical enrichment develop stable winds based on the Mg2SiO4
composites. However, there is a region (when the piston velocity, ∆up = 2 km/s and
nd/nH = 1×1015) of low chemical enrichment and small pulsation amplitudes where
those DARWIN models do not develop stable dust-driven winds. While those models
are important to demonstrate the dynamical evolution of dust-driven winds with chemi-
cal properties based on magnesium-iron silicates, they do not incorporate the importance
of asymmetries on the development of the wind.
4.5. APPENDIX: ADDITIONAL FIGURES 103
Liljegren et al. (2018) studied the effects of complex convective structures on wind
production and mass loss in M-type AGB stars. The authors used the CO5BOLD code
(Freytag et al. 2003) to simulate the convective surface and pulsations of AGBs, using the
results as boundary conditions in their DARWIN models. Those models show that gas
is not levitated uniformly above the star’s surface (only about 70% or less of the stellar
surface is covered), causing dust not to be formed in a spherical layer around the star.
Instead, it forms clumps. This appears to have a direct implication in an overestimation
of the mass-loss rate determined with homogeneous models. In this context, our work
shows direct evidence of the asymmetric and variable gas distribution at the innermost
layers of M-type AGB stars. A specific DARWIN model of R Car is beyond the scope
of this thesis and it is left for future studies. However, by contrasting the structures ob-
served in our reconstructed images and future dedicated CO5BOLD/DARWIN models,
we could obtain better constraints on the mass-loss evolution and the role of magnesium
dust composites on the development of stable winds.
4.5 Appendix: Additional figures
In this section, we show the comparison between the synthetic and observed V2 and
CPs across the pseudo-continuum, the first and second CO band heads for our January
and February epochs (see Sect. 4.3.3). In all figures, the red dots in the main panels
correspond to the synthetic V2 and CPs extracted from the raw reconstructed images,
while the data are shown with gray dots in the main panels. The lower panels show the
corresponding residuals (in terms of the number of standard deviations) coming from
the comparison between the data and the best-fit reconstructed images.
104 CHAPTER 4. K BAND INTERFEROMETRIC ANALYSIS OF R CAR
(a)
(b)
Figure 4.13: First four principal components for the January and February epochs for the
first CO band head. The relative scale of the structures in the eigenimages is displayed
with a color bar at the bottom of each panel. Only the central 20 mas of the eigenimages
are shown in each panel. The white contours correspond to the mean image across
wavelength (per given data set) and they represent 10, 30, 50, 70, 90, 95, 97, and
99% of the intensity’s peak. The Cumulative Explained Variance (CEV) per principal
component is located in the top left side of each panel.
4.5. APPENDIX: ADDITIONAL FIGURES 105
(a)
(b)
Figure 4.14: First four principal components for the January and February epochs for
the second CO band head. The panels are as described in Fig. 4.13.
106 CHAPTER 4. K BAND INTERFEROMETRIC ANALYSIS OF R CAR
Figure 4.15: Representative example of CPs versus spatial frequencies of the data set that
corresponds to the first CO band head at 2.2946 µm (gray squares). The colored dots
show the CPs from the recovered images obtained with different numbers of principal
components (see the labels on the plot).
(a) (b)
Figure 4.16: Two-dimensional photocenter offsets in mas. The black dots correspond to
the photocenter position of the continuum emission and the labeled colored dots to the
photo-centers of the CO band head. The error bars of the reported displacements are
smaller than the used symbols.
4.5. APPENDIX: ADDITIONAL FIGURES 107
Figure 4.18: Curve of constant levitation radius (Rl = 2.43 R∗) represented with a solid-
black line. The plot shows the condensation temperature, Tc, as a function of the dust
absorption coefficient, p. The red circles indicate the positions of a variety of dust grains.
The blue and green shaded regions correspond to the CO temperature ranges estimated
from the single-layer model (see labels on the plot).
(a) (b)
Figure 4.19: The panels are similar to Fig. 4.18 but they show the the median conden-
sation radii obtained directly from the pure-line CO reconstructed images at different
position angles (see labels on the panels) for a) the January and b) February epochs.
108 CHAPTER 4. K BAND INTERFEROMETRIC ANALYSIS OF R CAR
(a)
(b)
Figure 4.20: Fit to observed V2 and CPs from the best reconstructed images across the
pseudo-continuum for the (a) January and (b) February epochs. The red dots in the
main panel correspond to the synthetic V2 and CPs extracted from the raw reconstructed
images. The data are shown with gray dots in the main panels. The lower panels show
the residuals (in terms of the number of standard deviations) coming from the comparison
between the data and the best-fit reconstructed images.
4.5. APPENDIX: ADDITIONAL FIGURES 109
(a)
(b)
Figure 4.21: Fit to observed V2 and CPs from the best reconstructed images across the
first CO band head for the (a) January and (b) February epochs. The panels are as
described in Figure 4.20.
110 CHAPTER 4. K BAND INTERFEROMETRIC ANALYSIS OF R CAR
(a)
(b)
Figure 4.22: Fit to observed V2 and CPs from the best reconstructed images across the
second CO band head for the (a) January and (b) February epochs. The panels are as
described in Figure 4.20.
5C
h
ap
te
r
Conclusions and future work
***
Although there are lots of models to explain the mechanisms involved in the mass-loss
processes during the AGB phase, the evidence to confront these models is still scarce.
One of the main reasons for this, is the lack of high angular resolution observations in
the regions were winds are formed (1–3 R∗). One way to solve the issue is trough the
analysis of high-angular resolution observations. In particular, AGB stars are bright in
the NIR making them ideal candidates to be studied with IR interferometry. In this
thesis, we employed interferometric observations to study the innermost parts of the
M-type AGB star R Car. We are confident that the results presented across these pages
will motivate novel periodic studies of other AGB stars with IR interferometry. The
main results of this work can be summarized as follows.
• The complex surface of R Car in the H−band. Current models of convec-
tion in evolved stars, although very general, predict the emergence of structures of
varying sizes on the star’s surface. Such structures can change in position and size
over intervals ranging from weeks to years. Furthermore, the relationship between
these structures and the star’s physical parameters can also be determined with
models. However, can the models accurately capture the real physics happening in
the stars? To answer this question, we have to confront those models with obser-
vations. As mentioned in Chapter 3, providing the necessary observations requires
a window where dust is transparent and the ability to resolve the object’s photo-
sphere, something achieved in this work through interferometry in the H−band.
The importance of this work relies in the fact that we were able not only to image
the star’s surface at two different pulsation phases but also to identify the surface
111
112 CHAPTER 5. CONCLUSIONS AND FUTURE WORK
structures and their changes across the time, providing new observational tests on
how convection could work in AGB stars. The main results are presented below.
⋄ To derive the size of R Car’s stellar surface, we performed a parametric fitting
to the V2 at the two epochs of observation: 2014 and 2020 (φ = 0.78 and 1.01,
respectively). This gave us an estimate of Θ2014
UD = 10.23 ± 0.05 mas (1.86
± 0.01 AU) for the 2014 epoch and Θ2020
UD = 13.77 ± 0.14 mas (2.51 ± 0.23
AU) for the 2020 epoch. According to other studies, in H−band, the object
is expected to be smaller around phase 0.8, and larger before the minimum
of the pulsation phase. The diameters estimated by us agree with this trend.
⋄ To characterize the asymmetric structure observed in the large CP variations,
we reconstructed interferometric images. Those reconstructed images reveal
the stellar surface and several bright spots lying above it. We measured
the size of these patterns using a power spectrum analysis finding that the
patterns change in size between the two pulsation phases. On average they are
smaller when the star is smaller, and larger when the star is larger. Although
this is an expected behavior, to our knowledge, this is the first time that it is
shown for an M-type Mira star using images and with a quantitative analysis.
⋄ To unveil the nature of these spots, we compared their measured size with
three different extrapolations from stellar convection models. Those models
establish convective cell sizes proportional to ∼10 times the pressure scale
height in the atmosphere of the target. Our findings suggest that the char-
acteristic sizes of these spots are consistent with the predictions. With these
results, we provide important evidence to the still scarce observational proofs
of the role of convection in the evolution of AGBs.
⋄ For a better characterization of the convective structures and their time evo-
lution, future H−band interferometric images across the pulsation cycle of R
Car are necessary. A proper comparison of those observations with state-of-
the-art hydrodynamic and radiative transfer models will be mandatory.
• The innermost gaseous layers in the K band. Understanding the mechanisms
leading to the mass loss in M-type AGB stars requires to know the types of dust
that can form and survive near the star at temperatures around the condensation
temperature limits (∼ 900-1500 K). The reason is because the first interactions
113
that lead to the subsequent mass loss occur at distances around 1–3 R∗, where the
temperatures are of this order.
State-of-the-art wind models for AGB stars allow the exploration of dust particles
with different compositions to understand how efficient they are as wind drivers
through the scattering of the star’s radiation. However, directly detecting the
dust becomes a challenge since we require the angular resolution to observe the
regions where all these interactions begin to occur, something that is not achievable
with single-dish telescopes. The problem becomes more complicated since dust is
mostly transparent in the near-infrared and, at mid-infrared wavelengths there is
a considerable degradation of the resolution compared with shorter wavelengths.
One way to address this problem is by considering one of the most stable molecules
against dissociation, which can form in the closer vicinity of the star: the CO. This
molecule has the peculiarity that some of its (vibro-rotational) transitions, mainly
those detected in the K band, form in regions where dust begins to condense.
Therefore, by directly studying these regions through high-angular resolution tech-
niques we can also provide indirect constrains on the types of dust proposed in the
wind models. The importance of this second work is that we were able not only
to image and to characterize the innermost CO layers at different pulsation phases
but also to establish physical conditions of temperature and density that the dust
candidates have to fulfill to be considered efficient wind-drivers. The main results
of the work can be summarized as follows
⋄ We obtained a first-order estimate of the size of the source by applying a
geometrical model to the V2 across the K−band for our two epochs, which
are separated by ∼ 30 days, corresponding to ∼10 % of the full pulsation
period. From this analysis, we found that the size of the target changes with
wavelength due to the presence of molecules like H2O and CO at ∼ 2µm
and ∼ 2.29 − 2.33µm, respectively. The best estimate of the stellar disk
diameter in the K band was derived at a wavelength of ∼ 2.23µm. Our
parametric model allowed us to obtain a diameter of 16.67±0.05 mas (3.03
AU) in January, 2018 and 14.84±0.06 mas (2.70 AU) in February, 2018.
⋄ We found prominent CO band heads in our GRAVITY spectra. These band
heads also appear as drops in the V2 data at the CO 2-0 (λ0 ∼ 2.2946 µm)
114 CHAPTER 5. CONCLUSIONS AND FUTURE WORK
and CO 3-1 (λ0 ∼ 2.3240 µm) vibro-rotational transitions compared with
the pseudo-continuum, which indicate that the CO line-emitting region is
more extended. The innermost CO layers are important to trace the dust
nucleation region since most of the C and O present in the atmosphere are
concentrated in the stable bound of this molecule.
⋄ We estimated the temperature, size, and optical thickness of these CO layers
by employing a parametric single-layer model. With the optical thickness, we
estimated also the column density at the minimum peak of the first CO band
head. Our estimate of the column density goes in agreement with estimates
reported for other AGB stars, and it is enough to allow the formation of dust
grains. Although our results are consistent with other works, we can provide
even better estimations with future higher angular resolution measurements
in combination with more complex models of the CO layer/s.
⋄ Using the temperature and size of the CO layer, we set constraints on the
dust candidates that could coexist with the gas. We found that the compos-
ites Mg2SiO3 and MgSiO4 can coexist at the same distances as the ones we
obtained for the CO innermost layers (∼ 2-2.4 R∗). However, the conden-
sation temperatures for these dust particles are lower than the temperatures
obtained from the CO layer model. We propose that the temperature dis-
crepancy is due to the nature of the single-layer model, which considers a
homogeneous and symmetric distribution of the CO layer, giving average
values for TL, RL, and τL. Our hypothesis is supported by the varying inter-
ferometric CPs and the posterior image reconstruction, which show that the
CO distribution is neither homogeneous nor symmetric.
⋄ To complement the information obtained from our symmetrical models, we
reconstructed interferometric images. As expected from the variation of the
CPs, we observed several asymmetries in the images across the spectral ranges
and for both epochs of observation. To discard the presence of artifacts in
our images, we employed bootstrap, PCA, and spectro-astrometric analyses.
Our findings suggest the presence of a complex, perhaps clumpy distribution
of the CO layer, which also changes within a period of 30 days or less.
5.1. FUTURE WORK 115
5.1 Future work
Throughout this thesis, we studied the convective surface and the innermost gaseous
layers of the AGB star R Car to provide observational evidence that helps us to under-
stand how mass loss and convection occurs in these types of stars, two phenomena that
are not fully understood. We also provided constraints on the types of dust that could
coexist under physical conditions similar to the surrounding gas. However, due to its
low optical thickness in the H and K bands, we were unable to detect direct signals
from the sought dust particles. As can be seen in Figs. 5.1 and 5.2, the bands where
the highest number of dust features are found are L (3.0-4.0 µm) , M (4.0 - 5.0 µm),
and N (8.0 - 13.0 µm).
Figure 5.1: ISO-SWS spectra of the prototypical oxygen-rich AGB star Mira (o Ceti),
a carbon- rich star (RY Dra) and a star with mixed chemistry (V778 Cyg): oxygen-rich
dust but carbon-rich photosphere. The most prominent spectral features are indicated.
Figure taken from Yamamura et al. (2000).
One way to study the dust responsible for these features is through the analysis
of interferometric observations. As part of the binary AGB project, of which I am
a member along with various colleagues from the interferometric community1, we are
1This work consists in the study of several M-, S- and C- type AGB stars with infrared interferometry
116 CHAPTER 5. CONCLUSIONS AND FUTURE WORK
Figure 5.2: The typical spectrum of an AGB star in the infrared from 0.5 to 25 µm. The
molecules contributing to the spectrum of an M-type AGB star are shown in red, and
the molecules contributing to a C-type star are shown in green. The black spectra shows
the molecular contributions to both type of stars. Image taken from Paladini (2011).
studying a sample of different AGB stars in the L, M , and N bands with infrared
interferometry. Those observations were obtained with the MATISSE/VLTI instrument,
using the small, medium, and large ATs configurations, resulting in angular resolutions
of λ/2Bmax = 2.25 mas and λ/2Bmax = 7.50 mas for a baseline of 140 m and for
wavelengths λ = 3µm and λ = 10µm, respectively. In particular, I am analyzing R
Car data. Although the analysis of these data is a work in progress, I can provide some
to unveil, for the first time, the asymmetric nature of the dust layers around them through image
reconstruction and atmospheric modeling.
5.1. FUTURE WORK 117
preliminary plots and describe the methodology applied.
• By applying a geometrical model to the V2 data, we will estimate the size of the
star across the L, M , and N bands, just as we did for the analyses in the H and
K bands. This will provide us with an understanding of the extent of the different
gas and dust layers surrounding the star.
• In Figs. 5.3 and 5.4 we present preliminary plots of the calibrated V2 and CPs
of R Car in the L and M bands. As we can see, the CPs vary widely between
0 - 180 degrees, indicating that the star is asymmetric in these bands. To depict
the asymmetries, image reconstruction will be necessary. Given the high spectral
resolution of our data, and the computational cost of the reconstruction algorithms,
we will select different characteristic intervals such as 3.2 µm were the contribution
of water is present, or between 4.5 - 5.0 µm where CO is present. Once we reduce,
calibrate and filter our N band data, we will repeat this process to obtain images
at this wavelength range.
• Just like with the GRAVITY instrument, MATISSE has the advantage of provid-
ing us with both, the interferometric observables and the spectrum of the star.
By employing them, we can estimate the dust properties using complex radiative
transfer models such as RADMC3D (Dullemond et al. 2012). A demonstration of
the effectiveness of this methodology can be found in Cannon et al. (2023). Here
is a brief description of the method that will be employed.
⋄ With the RADMC3D code, it is possible to obtain intensity maps from arbi-
trary 3D dust density distributions with different dust grain properties and
stellar parameters. In contrast to the work of Cannon et al. (2023) were the
authors could not perform image reconstruction of their target (α Ori), the
geometrical distribution of the dust clumps will be given by our reconstructed
images. By using those images as a starting point, we will only have to select
the chemical composition, size, and distribution of the dust grains to run the
models with these initial conditions.
⋄ Then, we can extract and compare the spectrum from the RADMC3D models
with our data to test the consistency between both of them for different dust
properties and distributions. By minimizing the differences between models
118 CHAPTER 5. CONCLUSIONS AND FUTURE WORK
and data, we will be able to estimate the dust properties that best represent
the spectra.
Figure 5.3: Calibrated V2 and CPs of R Car in the L band (3.0 - 4.0 µm) as a function
of the spatial frequency. We can see how the observables change depending of the
wavelength probed. This means that different gaseous layers might have different sizes
but also their distribution can be different. In other words, there is a chromatic behavior
in the structure of the gaseous layers around R Car both in size and in the asymmetries
traced.
Figure 5.4: Calibrated V2 and CPs of R Car in the M band (4.1 - 4.7 µm) as a function
of the spatial frequency. We can see again the chromatic effects in the observables.
5.2. PROSPECTS FOR THE FOLLOWING DECADE: THE NGVLA 119
The analysis of these L, M , and N bands data will give us new observational evi-
dence about the types of dust responsible for driving the winds in M-type AGB stars,
their physical properties and their distribution around the star. Furthermore, image
reconstruction of future interferometric observations at different epochs will allow us to
monitor how this dusty circumstellar environment is changing, something that has not
been done yet in the field with milliarcsecond resolution.
5.2 Prospects for the following decade: the ngVLA
The following decade will be crucial for our understanding of the complex phenomena
occurring at the innermost parts of AGB stars. Once the next-generation Very Large
Array (ngVLA) is finished, it will enable the imaging of astronomical sources with an
unprecedented detail by providing improvements in sensitivity and angular resolution
compared with current radio interferometers operating in the 1.2 -116 GHz range.
In the field of stellar physics, with the Karl G. Jansky Very Large Array (VLA) and
ALMA it is possible to marginally resolve some of the closest (d . 150 pc) AGBs and
RSGs whose enormous radio photospheres can span several AU and subtend up to a
few tens of mas (see e.g., Asaki et al. 2023). However, with the ngVLA, we will be able
to measure detailed properties of radio photospheres (e.g., the presence of atmospheric
temperature gradients, spots, and surface features, as well as temporal changes). With
those studies, it will be possible to complement what we know about the atmospheric
physics, including the temperature structure of the atmosphere, and to set constraints
on the mechanisms (e.g., shocks, pulsation, convection) that help to drive the observed
mass-loss rates from these stars. Together with the interferometric observations in the
infrared, it will be possible to perform more detailed comparisons with predictions of
state-of-the-art 3D dynamic atmospheric models of AGB and RSG stars. Empirically
testing the emerging new 3D models that include the effects of pulsation, convection,
shocks and dust condensation, will demand new ultra high-resolution measurements of
diverse samples of stars using instruments like the ngVLA (Akiyama & Matthews 2019).
To illustrate the capabilities of the future ngVLA, in Fig. 5.5 we present an example
of an image reconstruction of four simulated images with different characteristics. The
interferometric observables were extracted by taking into account the u − v coverage
of the ngVLA Main Array (214 18m antennas, baseline from tens of meters to ∼ 1000
km, see Fig. 5.6 ). As we can see, the reconstructed images recover the features in the
120 CHAPTER 5. CONCLUSIONS AND FUTURE WORK
simulated images with a very-high detail. The interested reader can consult Akiyama &
Matthews (2019) for more details of the reconstruction algorithms employed.
Figure 5.5: Upper row: simulated stellar models used to test the capabilities of the
ngVLA. Once the observables were simulated, an image reconstruction was performed
and the results can be seen in the lower row. The first two panels on the left were
obtained from simulations based on dynamic 3D atmospheric models and correspond
to a 1 M⊙ AGB (panel a) and a 12 M⊙ RSG (panel b). Those images were extracted
from Freytag et al. (2017) and Chiavassa et al. (2009), respectively. Panels (c) and
(d) correspond to toy models comprising a uniform circular disk with three ’spots’ of
different sizes superposed (one brighter than the underlying photosphere and two that
are cooler). These models are then located at different distances : 222 pc and 1 kpc,
respectively. Figure taken from Akiyama & Matthews (2019).
5.3. CLOSING REMARKS 121
Figure 5.6: uv-coverage of the simulated ngVLA observations. The simulated ngVLA
observations of each model in the upper panel of Fig. 5.5 were obtained using the
simobserve task in CASA. In all cases the array configuration was taken to be the
ngVLA Main Array. Figure taken from Akiyama & Matthews (2019).
5.3 Closing remarks
In this thesis, we have showed that the structures on the surface and the innermost layers
of gas/dust in AGB stars can be highly asymmetric at different spatial and temporal
scales. To characterize these changes both spatially and temporally, new interferometric
observations with baselines that allow us to map different angular scales are necessary,
and this has to be done over one or several pulsation cycles of the star.
We also showed that such detailed studies are scarce and require observations with
the highest possible quality. However, we demonstrated that they are not impossible
to achieve and that provide valuable results. Therefore, it will be worthwhile to extend
these studies to a significant sample of AGB stars in different environments. As we can
see, there are lots of new things to discover in the near future, and instruments like the
ngVLA or the VLTI and their improvements will be an essential part of these discoveries.
Eventually, we expect to build a sufficiently large sample of high-resolution observational
evidence to fine-tune the theoretical models, in order to increase their predictive power.
122 CHAPTER 5. CONCLUSIONS AND FUTURE WORK
A
A
p
p
en
d
ix
Python codes
***
Here, I show examples of the codes used to obtain the results presented in Chapter 4.
This notebook can also be found in this Github repository. The repository contains a
Jupyter notebook which is included here in a pdf format.
123
R Car GRAVITY analysis
Abel Rosales-Guzmán
September 18, 2023
1 R Car analysis tools
1.0.1 This is a Python-Jupyter notebook which contains some examples of the usage
of the different algorithms used to analyze GRAVITY-VLTI data of the Mira
star R Car.
2 Contents
Example of the MCMC model fitted to the GRAVITY-VLTI R Car V2
Example of the Molsphere model fitted to the GRAVITY-VLTI R Car V2
Example of the Bootstrapping algorithm applied to the interferometric observables of R Car
Example of the Principal Component Analysis applied to our reconstructed images of R Car
2.0.1 The following are required Python (version > 3.0) packages which need to be
installed before the execution of the current Jupyter Notebook:
Astropy
corner
emcee
lmfit
matplotlib
numpy
scipy
sklearn
oitools.py (provided in the current project/directory)
[90]: ### Author: Abel Rosales-Guzman
###The following lines load the required Python modules:
from astropy.modeling import models
import astropy.io.fits as pyfits
from astropy import constants as const
from astropy import units as u
from copy import copy
1
import corner
import cv2
import datetime
import emcee
from lmfit import Model, Parameters, minimize
import matplotlib.pyplot as plt
from matplotlib.ticker import FormatStrFormatter
from multiprocessing import Pool
import numpy as np
import oitools
import os
import scipy.special as sp
from sklearn.decomposition import PCA
from sklearn.utils import resample
#####################################
2.0.2 The following script (ln 2 and 4) provides an example of the Uniform Disk model
fitted to the GRAVITY-VLTI R Car V2. The code uses a Markov-Chan Monte-
Carlo (MCMC) method implemented through the Python emcee package. The
model uses the parametric expression of a Uniform Disk in Fourier space. This
model only contains two parameters, the diameter of the disk, and a coefficient
to account for the over-resolved flux/component (Fr) plausible present in the
data. While the model only has two parameters, we used the MCMC approach
to better estimate the statistical errors of the best-fit model. The following
lines run the Uniform Disk model for one spectral channel (2.228 µm) of R
Car’s pseudo-continuum. The same setup could be used to replicate the best-
fit models for other spectral channels in the pseudo-continuum or across the
CO bandheads.
[2]: ### MCMC auxiliary functions:
def mas2rad(x_mas):
x_rad = 4.847309743e-9 * x_mas
return x_rad
def rad2mas(x_rad):
x_mas = x_rad/4.847309743e-9
return x_mas
## The following function describes a fully illuminated uniform disk in the␣
→֒Fourier space:
def disk(u, v, diam, Fr):
"""
diam = angular diameter in mas
Fr = Scaling factor that accounts for the over-resolved flux
"""
diam = mas2rad(diam)
r = np.sqrt(u**2 + v**2)
2
v_disk = 2 * (sp.j1(np.pi * diam * r) / (np.pi * diam * r))
return Fr*(v_disk)**2
## The following function describes the log-prior probability:
def lnprior(theta):
# probability is uniform inside a given range, and zero outside
if (1. < theta[0] < 30.) and (0. < theta[1] < 2.): #and (0. < theta[2] < 3.
→֒) :
return 0.0
return -np.inf
## The following function describes the log-likelihood (i.e, the Chi2)␣
→֒probability:
def lnlike(theta, u, v, vis2, vis2err):
diam = theta[0]
Fr = theta[1]
vis2model = disk(u, v, diam, Fr)
chi2 =-0.5*np.sum((vis2-vis2model)**2/(vis2err)**2)
return chi2
## The following function describes the posterior probability:
def lnprob(theta, u, v, vis2, vis2err):
lp = lnprior(theta)
if lp == -np.inf:
return lp
return lp + lnlike(theta, u, v, vis2, vis2err)
[4]: ## The following line passes to a variable the name of the OIFITS file with the␣
→֒data:
fits_f='phasediff_18_January_a200_COMB_RCar_cont_2.228m.fits'
##############Automatic from here##########
## The data are read from the OIFITS file ####
oidata = pyfits.open(fits_f)
oi_wave = oidata['OI_WAVELENGTH'].data
oi_vis2 = oidata['OI_VIS2'].data
waves = oi_wave['EFF_WAVE']
vis2 = oi_vis2['VIS2DATA']
vis2_err = oi_vis2['VIS2ERR']
vis2_err2 = np.full(vis2.shape, 0.01)
u = oi_vis2['UCOORD']
v = oi_vis2['VCOORD']
nvis = vis2.shape[0]
3
ndim = 2
nwalkers = 250
pos1 = np.random.uniform(1., 30., size=nwalkers)
pos2 = np.random.uniform(0., 2., size=nwalkers)
fig2, ax2 = plt.subplots(ncols=1,nrows=1, figsize=(7,7))
spf_u = u/waves[0]
spf_v = v/waves[0]
bssl = np.sqrt(spf_u**2+spf_v**2)
#### The following lines run the Markov chains:
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(spf_u, spf_v,␣
→֒vis2, vis2_err2))
steps = 150
pos = np.array([pos1, pos2]).T
sampler.run_mcmc(pos,steps, progress = True)
samples = sampler.get_chain()
flat_samples = sampler.get_chain(discard=70, flat=True)
fig1, _axs1 = plt.subplots(nrows=2, ncols=1, figsize=(12,7), sharex=True)
axs1 = _axs1.flatten()
axs1[1].set_xlabel('Step number')
axs1[0].set_ylabel('Diameter')
axs1[1].set_ylabel('Fr')
inds = np.random.randint(len(flat_samples), size=50)
for i in range(ndim):
axs1[i].plot(samples[:, :, i], "k", alpha=0.3)
######## Plots ######
diam_mcmc, Fr_mcmc = map(lambda v: (v[1], v[2] - v[1], v[1] - v[0]), zip(*np.
→֒percentile(flat_samples, [16, 50, 84], axis=0)))
fig3 = corner.corner(flat_samples, labels=['D', 'Fr'], show_titles=True,␣
→֒title_fmt='0.2f',
levels=(1 - np.exp(-0.5), 1 - np.exp(-2.0)),
truths=[diam_mcmc[0], Fr_mcmc[0]], quantiles=[0.16, 0.5, 0.84],
title_kwargs={'fontsize': 12})
print('Best Values')
print('Theta_UD = {:.2f} +- {:.2f} {:.2f} mas'.format(diam_mcmc[0],␣
→֒diam_mcmc[1], diam_mcmc[2]))
print('F_r = {:.3f} +- {:.3f} {:.3f}'.format(Fr_mcmc[0], Fr_mcmc[1],␣
→֒Fr_mcmc[2]))
iind1 = np.argsort(bssl)
4
ax2.errorbar(bssl, vis2, yerr=vis2_err2, fmt='o', ms=5, label='Data', c='red')
########
ax2.plot(bssl[iind1], disk(spf_u, spf_v, diam_mcmc[0], Fr_mcmc[0])[iind1],␣
→֒label='Model', c='blue', zorder=500)
ax2.text(1.0e7, 0.8, r'$\Theta_{{UD}}={:.2f}_{{-{:.2f}}}^{{+{:.2f}}}$ mas'.
→֒format(diam_mcmc[0], diam_mcmc[1], diam_mcmc[2]), fontsize=17)
ax2.text(1.0e7, 0.65, r'$F_r={:.3f}_{{-{:.3f}}}^{{+{:.3f}}}$'.
→֒format(Fr_mcmc[0], Fr_mcmc[1], Fr_mcmc[2]), fontsize=17)
ax2.set_ylim(-0.05, 1.0)
ax2.set_xlim(0.1e7, 1.71e7)
ax2.set_xticks(np.arange(0.1e7, 1.71e7, step=0.3e7))
ax2.set_title(r'{:.3f} $\mu m$'.format(waves[0]/1.e-6))
ax2.set_xlabel(r'Spatial Frequency [$rad^{-1}$]', fontsize=18)
ax2.set_ylabel('V$^2$', fontsize=18)
plt.show()
100%|฀฀฀฀฀฀฀฀฀฀| 150/150 [00:01<00:00, 143.10it/s]
Best Values
Theta_UD = 16.67 +- 0.06 0.06 mas
F_r = 0.964 +- 0.005 0.005
5
y?
 
10 
08 
05 
04 
02 
0 
2.228 um 
  
     
   
  
 
- +0.06 Ouo =16.67*79g mal 
- +0.005 F,=0.964%0 003 
o vs 07 10 13 15 
Spatial Frequency [rad”*]
6
Di
am
et
er
 
20 
15 
Elo 
os 
0 
 
  
 
  
 
0 2 20 6 so 100 10 140 
Step number. 
 
7
Go back to top
8
2.0.3 The following part of the notebook script (ln 5, 8 and 30) provides an example
of the single layer model (the so-called MOLsphere) fitted to the GRAVITY-
VLTI R Car V2 and spectra. The code also uses a Markov-Chain Monte-Carlo
method implemented through the Python emcee package. The model uses the
parametric set of expressions (see e.g., Rodriguez-Coira et al., 2021 and/or
the main text of our manuscript) which describe a star that irradiates as a
Black-Body surrounded by an infinitesimal layer of gas, which also irradiates as
Black-Body modified by a given optical depth. The model fits simultaneously
the spectrum of the source (particularly the region of the CO bandhead in
the K-band GRAVITY data) and the squared visibilities. The best-fit model
provides an estimate of the temperature, optical depth and the size of the layer,
in combination with the size of the central star. The following lines run the
MOLsphere model for one spectral channel (2.2946 ฀m). However, the same
setup could be used to replicate the best-fit models for other spectral channels
across the CO bandheads
[5]: ###### Visibility and Intensity Functions ######
#The following function creates a pixel grid for the MOLsphere model:
def img_grid(npixels, pixelscale):
size = npixels*pixelscale
x = np.linspace(-0.5*size, 0.5*size, npixels)
y = np.linspace(-0.5*size, 0.5*size, npixels)
return x, y
# This is the function that defines the MOLsphere model:
def SLModel(Tl, Rl, tau, Fr):
"""
This function calculates the V^2 corresponding to the MOLsphere
Tl, Rl = Temperature and radius of the layer
Ts, Rs = Temperature and radius of the star
tau = optical dept of the molecular layer at a wavelength lamdbda
r = radius to evaluate the model
Fr = Scaling factor to account for the over resolved component
Wavelength is given in microns
"""
npixels = 128 #px
pixelscale = 0.46875 #[mas/px]
x, y = img_grid(npixels, pixelscale)
r=np.zeros([x.shape[0], y.shape[0]])
for j in range(128):
for k in range(128):
r[j,k]=np.sqrt(x[j]**2+y[k]**2)
I2 = np.zeros([128,128])
Rs = 5. # diameter of R Car = 10 mas, thus R = 5mas
Ts = 2800. # Kelvin, effective temperture of R Car
wavelength = slm_wave
9
bb_s = models.BlackBody(temperature=Ts*u.K)
bb_L = models.BlackBody(temperature=Tl*u.K)
B_s = bb_s(wavelength*u.micron)
B_l = bb_L(wavelength*u.micron)
# I = []
for i in range(r.shape[0]):
for j in range(r.shape[1]):
if r[i,j] <= Rs:
cos_beta1 = (1./Rl)*np.sqrt(Rl**2-r[i,j]**2)
I2[i,j] = B_s.value*np.exp(-tau/cos_beta1) +␣
→֒B_l.value*(1-np.exp(-tau/cos_beta1))
elif (r[i,j] > Rs) and (r[i,j] < Rl):
cos_beta = np.sqrt(1-(r[i,j]/Rl)**2)
I2[i,j] = B_l.value*(1-np.exp(-2*tau/cos_beta))
else:
I2[i,j] = 0.000
##Adding the Over resolved component
I3 = I2/np.sum(I2)
I3 = Fr*I3
Fr_comp = 1.-Fr
[iind0, iind1] = np.where(I3 == 0)
px_zero = iind0.shape[0]
temp = Fr_comp/px_zero
I3[iind0, iind1] = temp
###Computing the V^2
vis_synth, phi_synth, t3phi_synth = oitools.compute_obs(oi_data, I3,␣
→֒pixelscale, 6, 4)
return (np.squeeze(vis_synth))**2
def I_SLModel(Tl, Rl, tau):
"""
This function calculates the intensity (I) of the MOLsphere
Tl, Rl = Temperature and radius of the layer
Ts, Rs = Temperature and radius of the star
tau = optical dept of the molecular layer at a wavelength lamdbda
r = radius to evaluate the model
wavelength is given in microns
"""
npixels = 128 #px
pixelscale = 0.46875 #[mas/px]
x, y = img_grid(npixels, pixelscale)
r=np.zeros([x.shape[0], y.shape[0]])
for j in range(128):
for k in range(128):
r[j,k]=np.sqrt(x[j]**2+y[k]**2)
I2 = np.zeros([128,128])
10
Rs = 5. # diameter of R Car = 10 mas, thus R = 5mas
Ts = 2800. # Kelvin, effective temperature of R Car
wavelength = slm_wave
bb_s = models.BlackBody(temperature=Ts*u.K)
bb_L = models.BlackBody(temperature=Tl*u.K)
B_s = bb_s(wavelength*u.micron)
B_l = bb_L(wavelength*u.micron)
# I = []
for i in range(r.shape[0]):
for j in range(r.shape[1]):
if r[i,j] <= Rs:
cos_beta1 = (1./Rl)*np.sqrt(Rl**2-r[i,j]**2)
# print(cos_beta)
I2[i,j] = B_s.value*np.exp(-tau/cos_beta1) +␣
→֒B_l.value*(1-np.exp(-tau/cos_beta1))
# I2[i,j]=I
elif (r[i,j] > Rs) and (r[i,j] <= Rl):
cos_beta = np.sqrt(1-(r[i,j]/Rl)**2)
# print(r[i,j], cos_beta, Rl2)
I2[i,j] = B_l.value*(1-np.exp(-2*tau/cos_beta))
# I2[i,j]=I
else:
I2[i,j] = 0.000
return I2
[8]: ###Similar to the first example in this notebook, the following functions␣
→֒define the probabilities for the MCMC engine:
def resid_flux(params, Flux, fluxerr):
"""
Function used to fit the tau given Tl and Rl
"""
pars = params.valuesdict()
Tl = pars['Tl']
Rl = pars['Rl']
tau = pars['tau']
temp_I = I_SLModel(Tl, Rl, tau)
F_tau0 = flux_tau0_1[0] # Flux at zero optical depth (see Rodríguez-Coira␣
→֒et al., 2021)
F_layer = np.sum(temp_I) # Flux of the layer
F_teor = F_layer/F_tau0
return (Flux-F_teor)/fluxerr
def lnprior(theta):
#flat prior, probability is uniform inside a given range, and zero outside
11
if (1000. < theta[0] < 2800.) and (5.1 < theta[1] < 20.) and (0.01 <␣
→֒theta[2] < 2.):
# print(theta)
return 0.0
return -np.inf
def lnlike(theta, vis2, vis2err):
Tl = theta[0]
Rl =theta[1]
Fr = theta[2]
# ndof = vis2.shape[0]
# dim = len(theta)
factor = 39. #ndof
#####fitting the tau##########
params = Parameters()
params.add('tau', value=0.5, min=1.e-6, max=5.)
params.add('Tl', value=Tl, vary=False)
params.add('Rl', value=Rl, vary=False)
result = minimize(resid_flux, params, args=(norm_spec[:,1][ind2], (np.
→֒array([0.05]))**2))
taul = result.params['tau'].value
Taus_err = result.params['tau'].stderr
TChi2 = result.redchi
if TChi2 < 3: ##keeping only the smallest Chi^2
Tauss_chi2.append(TChi2)
Tauss.append(taul)
TRl.append(result.params['Rl'].value)
TTl.append(result.params['Tl'].value)
vis_synth2 = SLModel(Tl, Rl, taul, Fr)
else:
Tauss_chi2.append(TChi2)
Tauss.append(0.)
TRl.append(result.params['Rl'].value)
TTl.append(result.params['Tl'].value)
tau_no=0.
vis_synth2 = SLModel(Tl, Rl, tau_no, Fr)
return ((-0.5*np.sum((vis2 - vis_synth2)**2/vis2err**2)/(factor)))
def lnprob(theta, vis2, vis2err):
lp = lnprior(theta)
if not np.isfinite(lp):
return -np.inf
else:
return lp + lnlike(theta, vis2, vis2err)
12
[30]: ## The MCMC main routine. Notice that it has to go inside the "if" command to␣
→֒run in parallel with the Pool() module
if __name__=="__main__":
#Estimated time for the execution of this cell using 1000 iterations: 4Hrs␣
→֒40min
norm_spec = np.genfromtxt('RCar_normspec_final_Jan.txt') #January␣
→֒#norm_spec[:,0]-> wavelength, norm_spec[:,1]- > Flux
oi_data = 'phasediff_18_January_COMB_RCar_CO_2.2946m.fits' #filename of the␣
→֒OIFITS data file
specr_min = [2.2945e-6] #Spectral range defined in meters
specr_max = [2.2947e-6]
################
wave_range = [2.2945e-6, 2.2947e-6] #### Wavelength range that wants to be␣
→֒fitted
flux_tau0_1 = [0.0013643644802048] # Calculated flux if tau = 0 (See Eq. 7␣
→֒from Rodriguez-Coira et al., 2021)
lambda_co=[2.2946] #center of the band head defined in microns
################Automatic from here#############
oidata = pyfits.open(oi_data)
oi_wave = oidata['OI_WAVELENGTH'].data
oi_vis2 = oidata['OI_VIS2'].data
oi_t3 = oidata['OI_T3'].data
waves = oi_wave['EFF_WAVE']
vis2 = oi_vis2['VIS2DATA']
vis2_err = oi_vis2['VIS2ERR']
t3 = oi_t3['T3PHI']
t3_err = oi_t3['T3PHIERR']
u_coord = oi_vis2['UCOORD']
v_coord = oi_vis2['VCOORD']
u1 = oi_t3['U1COORD']
u2 = oi_t3['U2COORD']
v1 = oi_t3['V1COORD']
v2 = oi_t3['V2COORD']
nvis = vis2.shape[0]
nt3phi = t3.shape[0]
vis2_err2 = np.full(vis2.shape[0], 0.01)
##########
# To compute the UV coordinates of the closure phases:
uv_cp = np.zeros([u1.shape[0]])
u_cp = np.zeros([u1.shape[0]])
v_cp = np.zeros([u1.shape[0]])
13
u3 = -1.0 * (u1 + u2)
v3 = -1.0 * (v1 + v2)
for j in range(u1.shape[0]):
uv1 = np.sqrt((u1[j]) ** 2 + (v1[j]) ** 2)
uv2 = np.sqrt((u2[j]) ** 2 + (v2[j]) ** 2)
uv3 = np.sqrt((u3[j]) ** 2 + (v3[j]) ** 2)
if uv1 >= uv2 and uv1 >= uv3:
uv_cp[j] = uv1
u_cp[j] = u1[j]
v_cp[j] = v1[j]
elif uv2 >= uv1 and uv2 >= uv3:
uv_cp[j] = uv2
u_cp[j] = u2[j]
v_cp[j] = v2[j]
elif uv3 >= uv1 and uv3 >= uv2:
uv_cp[j] = uv3
u_cp[j] = u3[j]
v_cp[j] = v3[j]
[ind] = np.where((waves >= wave_range[0]) & (waves <= wave_range[1]))
slm_wave = waves[ind[0]]/1.e-6
print(slm_wave)
[ind2] = np.where((norm_spec[:,0] >= wave_range[0]/1.e-6) & (norm_spec[:,0]␣
→֒<= wave_range[1]/1.e-6)) ##index of the spectrum
spf_u = u_coord/waves[0]
spf_v = v_coord/waves[0]
bsl = np.sqrt(spf_u**2 + spf_v**2)
print(norm_spec[:,0][ind2][0])
Tauss = []
Tauss_chi2=[]
TRl=[]
TTl=[]
np.seterr(all='ignore')
ndim = 3
nwalkers = 20 #must be even
steps = 1000 #1000
with Pool() as pool:
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(vis2,␣
→֒vis2_err2), threads=pool)
pos0 = np.random.uniform(1000, 2800, size=nwalkers)
pos1 = np.random.uniform(5.1, 20, size = nwalkers)
pos2 = np.random.uniform(0.1, 2., size = nwalkers)
pos = np.array([pos0, pos1, pos2]).T
sampler.run_mcmc(pos, steps, progress=True)
14
samples = sampler.get_chain()
fig, axes = plt.subplots(3, figsize=(10,7), sharex = True)
fig.patch.set_facecolor('white')
labels = ['Tl', 'Rl', 'Fr']
for i in range(ndim):
ax = axes[i]
ax.plot(samples[:,:,i], 'k', alpha=0.3)
ax.set_xlim(0, len(samples))
# ax.set_ylim()
ax.set_ylabel(labels[i])
axes[-1].set_xlabel('step number');
plt.savefig('Jan_Conv_TlRlFr_{:.4f}mic.png'.format(norm_spec[:,0][ind2][0]))
plt.show()
tauss1 = np.asarray(Tauss)
[iind11] = np.where((tauss1 > 2e-6) & (tauss1 < 4.8))
Taumcmc = np.percentile(tauss1[iind11], [16, 50, 84])
q1 = np.diff(Taumcmc)
print('Tau')
print(Taumcmc[1], q1[0], q1[1])
print('Samples to discard')
discard = int(input())
flat_samples = sampler.get_chain(discard=discard, flat=True, thin=1)
prob_dist = sampler.get_log_prob(discard=discard, flat=True, thin=1)
indd = np.where(prob_dist > -6.0)
fl_samples = np.squeeze(flat_samples[indd,:])
Tl_mcmc, Rl_mcmc, Fr_mcmc= map(lambda v: (v[1], v[2] - v[1], v[1] - v[0]),␣
→֒zip(*np.percentile(fl_samples, [16, 50, 84], axis=0)))
fig = corner.corner(fl_samples, labels=['Tl', 'Rl', 'Fr'],␣
→֒show_titles=True, title_fmt='0.3e',
levels=(1 - np.exp(-0.5), 1 - np.exp(-2.0)),
truths=[Tl_mcmc[0], Rl_mcmc[0], Fr_mcmc[0]], quantiles=[0.
→֒16, 0.5, 0.84],
title_kwargs={'fontsize': 12})
fig.savefig('Jan_Corner_{:.4f}mic.png'.format(norm_spec[:,0][ind2][0]))
fig.show()
print('Best Values')
print('T_l = {:.3f} +- {:.3f} {:.3f} K'.format(Tl_mcmc[0], Tl_mcmc[1],␣
→֒Tl_mcmc[2]))
print('R_l = {:.3f} +- {:.3f} {:.3f} mas'.format(Rl_mcmc[0], Rl_mcmc[1],␣
→֒Rl_mcmc[2]))
print('Tau_l = {:.3f} +- {:.3f} {:.3f}'.format(Taumcmc[1], q1[0], q1[1]))
print('Fr = {:.3f} +- {:.3f} {:.3f}'.format(Fr_mcmc[0], Fr_mcmc[1],␣
→֒Fr_mcmc[2]))
15
print((np.sum(I_SLModel(Tl_mcmc[0], Rl_mcmc[0], Taumcmc[1]))/
→֒flux_tau0_1[0])-norm_spec[:,1][ind2])
v2_model = SLModel(Tl_mcmc[0], Rl_mcmc[0], Taumcmc[1], Fr_mcmc[0])
bsl_ord = np.argsort(bsl)
#
#######Plots########
plt.figure(figsize=(9,9))
plt.xlabel(r'Spatial Frequency $[rad^{-1}]$', fontsize=22)
plt.ylabel(r'Squared Visibility', fontsize=22)
plt.ylim(-0.05, 1.)
index = np.random.randint(len(samples), size=100)
for inds in index:
sampls=fl_samples[inds]
plt.plot(bsl[bsl_ord], SLModel(sampls[0], sampls[1],␣
→֒tauss1[iind11][inds], sampls[2])[bsl_ord],'C1', alpha=0.1)
plt.errorbar(bsl, vis2, yerr=vis2_err2, label='Data', fmt='o',␣
→֒color='limegreen')
plt.plot(bsl[bsl_ord], v2_model[bsl_ord], c='k',
label=r'Model T$_L$={:.3f} K, R$_L$={:.3f} mas, $\tau_L$={:.3f},␣
→֒F$_r$={:.3f}'.format(Tl_mcmc[0], Rl_mcmc[0], Taumcmc[1], Fr_mcmc[0]))
plt.tick_params(axis='both', labelsize=20)
plt.legend(fontsize=13)
plt.tight_layout()
plt.savefig('Jan_vis2spf_model_{:.4f}mic.png'.format(norm_spec[:
→֒,0][ind2][0]))
plt.show()
plt.figure(figsize=(9,9))
plt.xlabel(r'Spatial Frequency $[rad^{-1}]$', fontsize=22)
plt.ylabel(r'Squared Visibility', fontsize=22)
plt.ylim(-0.05, 1.)
index = np.random.randint(len(samples), size=100)
plt.errorbar(bsl, vis2, yerr=vis2_err2, label='Data', fmt='o',␣
→֒color='limegreen')
plt.scatter(bsl, v2_model, c='purple',
label=r'Model T$_L$={:.3f} K, R$_L$={:.3f} mas, $\tau_L$={:.3f},␣
→֒F$_r$={:.3f}'.format(Tl_mcmc[0], Rl_mcmc[0], Taumcmc[1],␣
→֒Fr_mcmc[0]),zorder=500)
plt.tick_params(axis='both', labelsize=20)
plt.legend(fontsize=13)
plt.tight_layout()
plt.savefig('Jan_V3_vis2spf_model_{:.4f}mic.png'.format(norm_spec[:
→֒,0][ind2][0]))
16
plt.show()
2.2946410354052205
2.294641013196623
100%|฀฀฀฀฀฀฀฀฀฀| 1000/1000 [4:40:16<00:00, 16.82s/it]
Tau
1.7735459930550008 0.4417837826432305 0.666409310755729
Samples to discard
200
<ipython-input-30-9d9a211fa08b>:121: UserWarning: Matplotlib is currently using
module://ipykernel.pylab.backend_inline, which is a non-GUI backend, so cannot
show the figure.
fig.show()
Best Values
T_l = 1334.031 +- 106.453 92.060 K
R_l = 12.489 +- 1.577 1.239 mas
Tau_l = 1.774 +- 0.442 0.666
Fr = 0.955 +- 0.033 0.030
[0.0059141]
17
Fr 
11=1334e 4031538018   
     
- +157e+00 | = 1.2498 + 0111530400   
  
  95528 0132000 
  
 
   
 
  
                  
 
181
1.0  
—— Model T,=1334.031 K, R=12.489 mas, 1,=1.774, F,=0.955 
% Data 
0.8 
  
0.4 
S
q
u
a
r
e
d
 
Vi
s 
0.2 
 
0.0     
 
0.4 0.6 0.8 10 1.2 14 
Spatial Frequency [rad”1] . 
19
Go back to top
20
2.0.4 The following example shows how the samples from the interferometric observ-
ables were created using a bootstraping algorithm. The different samples are
created to explore different best-fit images obtained with the regularized mini-
mization algorithm BSMEM to recover images. By generating different samples
from the distribution of data points in our data help us, also, to explore slightly
different u-v planes than the original distribution. This helps to ameliorate
the effect of the sparseness of the u-v plane in the mean reconstructed image.
The following part of the notebook (In 6 and 7) shows how those samples are
obtained using the random seed generator from the Python sklearn package
“resample”, which considers replacement. The example is set for one spectral
channel across the pseudo-continuum (2.228 µm) but the same code could be
used for other channels and OIFITS files. The final samples are stored into
OIFITS files, which can be used by BSMEM for imaging the source. Notice
that, for the example, only 10 different samples are created. However, this
could be expanded as required by the user.
[6]: waves_cont = np.array([2.228]) ## Wavelength used for the example (in microns)
month = ['January']
file = 'phasediff_18_'+month[0]+'_a200_COMB_RCar_cont_'+'{:.3f}'.
→֒format(waves_cont[0])+'m.fits'
path=''
directory = os.path.exists(path+'bt_method/'+'{:.3f}'.format(waves_cont[0]))
if directory == False:
print('Creating directory')
os.system('mkdir '+path+'bt_method/'+'{:.3f}'.format(waves_cont[0]))
oi_data=pyfits.open(file)
for k in range(50):
oi_data.writeto(path+'bt_method/'+'{:.3f}'.format(waves_cont[0])+'/
→֒'+'bt{}_'.format(k)+file, overwrite=True)
data= pyfits.open(path+'bt_method/'+'{:.3f}'.format(waves_cont[0])+'/
→֒'+'bt{}_'.format(k)+file, mode='update')
av2=data['OI_VIS2'].data['VIS2DATA']
at3=data['OI_T3'].data['T3PHI']
### We label each data point
numbers_v2=np.arange(av2.shape[0])
numbers_t3=np.arange(at3.shape[0])
### Using the labels we create random samples of the data points allowing␣
→֒replace and keeping the total numer of points.
### This method allows us to create as many data samples as we want. In our␣
→֒case we created 50
boot_v2 = resample(numbers_v2, replace=True, n_samples=numbers_v2.shape[0],␣
→֒random_state=k)
boot_t3 = resample(numbers_t3, replace=True, n_samples=numbers_t3.shape[0],␣
→֒random_state=k)
aa=np.array(boot_v2, dtype='int')
bb=np.array(boot_t3, dtype='int')
21
data['OI_VIS2'].data = data['OI_VIS2'].data[aa]
data['OI_T3'].data = data['OI_T3'].data[bb]
data.flush()
data.close()
[7]: fig4, _axs4 = plt.subplots(nrows=1, ncols=3, figsize=(15,6))
axs4 = _axs4.flatten()
colrs = plt.cm.gist_rainbow(np.linspace(0,1,10))
temp = np.random.randint(23)
axs4[0].set_xlabel(r'Spatial Frequency [rad$^{{-1}}$]')
axs4[0].set_ylabel(r'V$^2$')
axs4[1].set_xlabel(r'Spatial Frequency [rad$^{{-1}}$]')
axs4[1].set_ylabel(r'Closure Phase')
axs4[2].set_xlabel(r'UCOORD [m]')
axs4[2].set_ylabel(r'VCOORD [m]')
for k in range(10): ### Number of samples generated from the data
################Automatic from here#############
oidata = pyfits.open('bt_method/2.228/bt{}_'.format(k)+file)
oi_wave = oidata['OI_WAVELENGTH'].data
oi_vis2 = oidata['OI_VIS2'].data
oi_t3 = oidata['OI_T3'].data
waves = oi_wave['EFF_WAVE']
vis2 = oi_vis2['VIS2DATA']
vis2_err = oi_vis2['VIS2ERR']
t3 = oi_t3['T3PHI']
t3_err = oi_t3['T3PHIERR']
u_coord = oi_vis2['UCOORD']
v_coord = oi_vis2['VCOORD']
u1 = oi_t3['U1COORD']
u2 = oi_t3['U2COORD']
v1 = oi_t3['V1COORD']
v2 = oi_t3['V2COORD']
nvis = vis2.shape[0]
nt3phi = t3.shape[0]
vis2_err2 = np.full(vis2.shape[0], 0.01)
##########
# To compute the UV coordinates of the closure phases:
uv_cp = np.zeros([u1.shape[0]])
u_cp = np.zeros([u1.shape[0]])
22
v_cp = np.zeros([u1.shape[0]])
u3 = -1.0 * (u1 + u2)
v3 = -1.0 * (v1 + v2)
for j in range(u1.shape[0]):
uv1 = np.sqrt((u1[j]) ** 2 + (v1[j]) ** 2)
uv2 = np.sqrt((u2[j]) ** 2 + (v2[j]) ** 2)
uv3 = np.sqrt((u3[j]) ** 2 + (v3[j]) ** 2)
if uv1 >= uv2 and uv1 >= uv3:
uv_cp[j] = uv1
u_cp[j] = u1[j]
v_cp[j] = v1[j]
elif uv2 >= uv1 and uv2 >= uv3:
uv_cp[j] = uv2
u_cp[j] = u2[j]
v_cp[j] = v2[j]
elif uv3 >= uv1 and uv3 >= uv2:
uv_cp[j] = uv3
u_cp[j] = u3[j]
v_cp[j] = v3[j]
spf_u = u_coord/waves[0]
spf_v = v_coord/waves[0]
uvvis_range = np.sqrt(spf_u**2 + spf_v**2)
uvcp_range=uv_cp/waves[0]
## The plots show, with different colored dots, the 10 different sanples of␣
→֒the data generated.
## The first plot displays the samples for the squared visibilities, the␣
→֒second panel displays the closure phases
## and the third panel displays the u-v coverage explored, respectively.
axs4[0].errorbar(uvvis_range, vis2, yerr=vis2_err2, fmt='o',␣
→֒color=colrs[k], alpha=0.2)
axs4[1].errorbar(uvcp_range, t3, yerr=t3_err, fmt='o', color=colrs[k],␣
→֒alpha=0.2)
axs4[2].scatter(u_coord, v_coord, color=colrs[k], alpha=0.2)
axs4[2].scatter(-u_coord, -v_coord, color=colrs[k], alpha=0.2)
oidata.close()
plt.show()
23
2.0.5 The following part of the notebook is based on the CASSINI-PCA package.
The code (In 115, 116) shows an example of the Principal Component Analysis
(PCA) applied to our January reconstructed images across the 1st CO band
head of the R Car GRAVITY data. According to Medeiros et al., (2018), the
visibilities of the Principal Components are equal to the Principal Components
of the visibilities. This method is useful to trace the changes, across a set of
images, which have a significant effect on the visibility (amplitude and -closure-
phase) profile. In the case of the R Car data, we use the PCA to evaluate the
variability of the source across wavelengths. First, we computed the PCA of
the recovered images across the CO bandpass. Second, we evaluate how many
of those components were needed to reproduce the interferometric observables.
This allow us to filter out variable structures potentially associated with noise
and/or artifacts in the image. Hence, by selecting the number of principal
components that reproduce the observables we are able to trace “bona fide”
variablility of the target.
[115]: def observables(data, wave_do, wave_up, imrec):
"""
Function that computes the observables from the data and the synthetic␣
→֒observables
from the reconstructed image
data = The OIFITS file containing the interferometric data
[wave_do, wave_up] = Is the wavelength interval containing the spectral␣
→֒channel analysed
imrec = The reconstructed image from where we extract the synthetic␣
→֒observables
"""
scale=0.469
npix=128
######## AUTOMATIC FROM HERE ############
24
oidata = pyfits.open(data)
oi_wave = oidata['OI_WAVELENGTH'].data
oi_vis2 = oidata['OI_VIS2'].data
oi_t3 = oidata['OI_T3'].data
waves = oi_wave['EFF_WAVE']
vis2 = oi_vis2['VIS2DATA']
vis2_err = oi_vis2['VIS2ERR']
t3 = oi_t3['T3PHI']
t3_err = oi_t3['T3PHIERR']
vis2_err2 = np.full(vis2_err.shape, 0.01)
u = oi_vis2['UCOORD']
v = oi_vis2['VCOORD']
u1 = oi_t3['U1COORD']
u2 = oi_t3['U2COORD']
v1 = oi_t3['V1COORD']
v2 = oi_t3['V2COORD']
nvis = vis2.shape[0]
nt3phi = t3.shape[0]
##########
# To compute the UV coordinates of the closure phases:
uv_cp = np.zeros([u1.shape[0]])
u_cp = np.zeros([u1.shape[0]])
v_cp = np.zeros([u1.shape[0]])
u3 = -1.0 * (u1 + u2)
v3 = -1.0 * (v1 + v2)
for j in range(u1.shape[0]):
uv1 = np.sqrt((u1[j]) ** 2 + (v1[j]) ** 2)
uv2 = np.sqrt((u2[j]) ** 2 + (v2[j]) ** 2)
uv3 = np.sqrt((u3[j]) ** 2 + (v3[j]) ** 2)
if uv1 >= uv2 and uv1 >= uv3:
uv_cp[j] = uv1
u_cp[j] = u1[j]
v_cp[j] = v1[j]
elif uv2 >= uv1 and uv2 >= uv3:
uv_cp[j] = uv2
u_cp[j] = u2[j]
v_cp[j] = v2[j]
elif uv3 >= uv1 and uv3 >= uv2:
uv_cp[j] = uv3
u_cp[j] = u3[j]
v_cp[j] = v3[j]
[ind] = np.where((waves >= wave_do) & (waves <= wave_up))
uvvis_range = np.sqrt(u ** 2 + v ** 2)/waves[ind[0]]
uvcp_range=uv_cp/waves[ind[0]]
vis_synth, phi_synth, t3phi_synth = oitools.compute_obs_wave(data, imrec,␣
→֒scale, 6, 4, waves[ind[0]])
25
vis_synth=np.squeeze(vis_synth)
return vis2[:,ind[0]], vis2_err2[:,ind[0]], t3[:,ind[0]], t3_err[:,ind[0]],␣
→֒vis_synth, t3phi_synth, uvvis_range, uvcp_range
def PCA_ncomps_wave(ncomps, waves, month, interval, data_pth):
"""
Function to calculate the Principal Components of a given image data set.
ncomps = Number of desired principal components (depends on the explained␣
→֒variance that you want the components
to explain)
waves = Array with the spectral channels that you want to analyse
month = Our data was split into two months: January and February
interval = Continuum, 1st CO band head or 2nd CO band head
data_pth = Path of the OIFITS files from where we extract the␣
→֒interferometric observables (V^2 and CPs)
"""
images = np.zeros([len(waves), 128,128])
images_smooth = np.zeros([len(waves), 128,128])
for i in range(len(waves)):
img_path = 'mean_img_{}{:.4f}mic.fits'.format(month, waves[i]) # Path␣
→֒to the mean reconstructed images
img = pyfits.getdata(img_path)
images[i,:,:] = img
mn_img = np.mean(images, axis=0)
std_img = np.std(images, axis=0)
rmn_imgs = (images - mn_img)/std_img ### To remove the mean from all the␣
→֒images
images_rshp = rmn_imgs.reshape(rmn_imgs.shape[0], 128*128)
PCAs = PCA(n_comps)
PCAs.fit(images_rshp)
converted_img = PCAs.fit_transform(images_rshp)
print('Cumulative explained variance')
print(np.cumsum(PCAs.explained_variance_ratio_))
comp_data = np.zeros([len(waves), ncomps, 128, 128])
mean_PCA = PCAs.mean_.reshape(128,128)
ratio_expv_PCA = np.cumsum(PCAs.explained_variance_ratio_)
plot_ncomp(PCAs, interval, mn_img)
for i in range(len(waves)):
comp1= (PCAs.components_[0].reshape(128,128)*converted_img[i,0])
comp2 = (PCAs.components_[1].reshape(128,128)*converted_img[i,1])
comp3 = (PCAs.components_[2].reshape(128,128)*converted_img[i,2])
comp4 = (PCAs.components_[3].reshape(128,128)*converted_img[i,3])
comp_data[i,0,:,:] = (comp1 + mean_PCA)*std_img + mn_img
comp_data[i,1,:,:] = (comp1 + comp2 + mean_PCA)*std_img + mn_img
comp_data[i,2,:,:] = (comp1 + comp2 + comp3 + mean_PCA)*std_img + mn_img
26
comp_data[i,3,:,:] = (comp1 + comp2 + comp3 + comp4 +␣
→֒mean_PCA)*std_img+ mn_img
if waves[i] == 2.2946:
plot_obs_ncomp(data_pth, comp_data[i], waves[i], month, ncomps,␣
→֒interval)
def plot_ncomp(PCA_comps, month, mn_img):
"""
Function that plots the principal components with their cumulative␣
→֒explained variance
for the images corresponding to the 1st CO band head centered at 2.2946␣
→֒microns
"""
levs = np.array([0.10, 0.30, 0.50, 0.70, 0.90, 0.95, 0.97, 0.99])
fig1, _axs1 = plt.subplots(ncols=4, nrows=1, figsize=(15, 8))
axs1 = _axs1.flatten()
fig1.suptitle('Principal Components '+month, fontsize=19, y=0.75)
axs1[0].set_ylabel('Declination [mas]', fontsize=17)
axs1[0].set_xlabel('Right ascension [mas]', fontsize=17)
img_ms = np.zeros([128,128])
cv2.circle(img_ms, (64,64), 32, color=1, thickness=-1)
for i in range(4):# Number of principal components
axs1[i].set_xticks(np.arange(-30,40,10))
axs1[i].set_yticks(np.arange(-30,40,10))
axs1[i].tick_params(axis='both', labelsize=15)
temp = np.ma.masked_where(PCA_comps.components_[i].
→֒reshape(128,128)*img_ms == 0, PCA_comps.components_[i].
→֒reshape(128,128)*img_ms)
custom_cmap = copy(plt.cm.get_cmap("coolwarm"))
custom_cmap.set_bad('gray')
imm = axs1[i].imshow(temp, origin='lower', cmap=custom_cmap,␣
→֒extent=[30,-30,-30,30])#, norm=norm)
axs1[i].contour(mn_img, levels=levs*np.max(mn_img),␣
→֒extent=[30,-30,-30,30], linewidths=(1.3,), colors='white', origin='lower',␣
→֒alpha=0.9)
axs1[i].text(27, 24, 'CEV = {:.4f}'.format(np.cumsum(PCA_comps.
→֒explained_variance_ratio_)[i]), fontsize=17)
cbar=fig1.colorbar(imm, ax=axs1[i], orientation='horizontal')
cbar.formatter.set_powerlimits((0, 0))
if i != 0:
axs1[i].set_xticklabels([])
axs1[i].set_yticklabels([])
cbar.ax.tick_params(labelsize=15)
fig1.subplots_adjust(wspace=0.05, hspace=0.0)
fig1.savefig('')
plt.show()
27
def plot_obs_ncomp(data_pth, data_comp, wave, month, ncomps, interval):
"""
Function that generates a representative example of the Closure Phases␣
→֒versus
spatial frequencies of the data set that corresponds to the
1st CO band-head at 2.2946 microns (gray squares). The colored dots show
the CPs from the recovered images obtained with different numbers of␣
→֒Principal Components
data_pth = Path from where we extract the interferometric observables
data_comp = Array containing the Eigenimages (or principal components)
wave = spectral channel
month = Our data were split into to months: 'January' and 'February'
"""
fig4, _axs4 = plt.subplots(nrows=1, ncols=4, figsize=(18,5), sharey=True)
axs4 = _axs4.flatten()
axs4[0].set_xlabel(r'Sp. Freq. [$rad^{-1}$]', fontsize=18)
axs4[0].set_ylabel('Closure Phase [deg]', fontsize=18)
axs4[0].set_yticks(np.arange(-200., 201, step=50))
wave_do = (wave*1e-6)-1e-10
wave_up = (wave*1e-6)+1e-10
######
scale=0.46875
ncomps=4
colrs=plt.cm.jet(np.linspace(0,1,5))
for k in range(ncomps):
v2, v2e, t3, t3e, a, c, uvvis, uvcp = observables(data_pth, wave_do,␣
→֒wave_up, data_comp[k,:,:])
axs4[k].errorbar(uvcp, t3, yerr=t3e, color = 'gray', fmt = 'o',␣
→֒marker='s')
axs4[k].scatter(uvcp, c, zorder=50, color=colrs[k], s=10, label='N =␣
→֒{}'.format(k+1))
axs4[k].set_ylim(-201,201)
axs4[k].legend(fontsize='xx-small', title='Number of components (N)',␣
→֒title_fontsize='small')
#####
fig4.subplots_adjust(wspace = 0.05, hspace=0)
fig4.savefig('PCA_t3phi_2.2946mic.png', bbox_inches='tight')
plt.show()
[116]: ##Calling the functions to calculate and plot the Principal Components
waves_1CO = np.array([2.2941, 2.2944, 2.2946, 2.2952, 2.2954, 2.2962, 2.2968])␣
→֒#wavelengths of the 1st CO bandhead
Mths = ['January']
n_comps = 4 #Number of Eigenvectors containing >90% of the variance in the␣
→֒dataset
28
PCA_ncomps_wave(4, waves_1CO, Mths[0], 'January 1st CO', 'January_COMB_RCar.
→֒fits')
Cumulative explained variance
[0.71037513 0.84848733 0.92557209 0.96618601]
Go back to top
29
Bibliography
***
Abuter, R., Accardo, M., Amorim, A., et al. 2017, Astronomy & Astrophysics, 602, A94
Akiyama, K. & Matthews, L. D. 2019, arXiv preprint arXiv:1910.00013
Allard, F., Guillot, T., Ludwig, H.-G., et al. 2003, in Symposium-International Astro-
nomical Union, Vol. 211, Cambridge University Press, 325–332
Allard, F., Homeier, D., & Freytag, B. 2010, arXiv preprint arXiv:1011.5405
Anugu, N., Kluska, J., Gardner, T., et al. 2023, The Astrophysical Journal, 950, 149
Arroyo-Torres, B., Wittkowski, M., Chiavassa, A., et al. 2015, Astronomy & Astro-
physics, 575, A50
Asaki, Y., Maud, L. T., Francke, H., et al. 2023, arXiv preprint arXiv:2310.09664
Auer, L. & Woolf, N. 1965, The Astrophysical Journal, 142, 182
Babu, G. J. & Singh, K. 1983, The Annals of Statistics, 999
Baldwin, J., Haniff, C., Mackay, C., & Warner, P. 1986, Nature, 320, 595
Baraniuk, R., Davenport, M., DeVore, R., & Wakin, M. 2008, Constructive Approxima-
tion, 28, 253
Baron, F., Monnier, J. D., & Kloppenborg, B. 2010, in Optical and Infrared Interferom-
etry II, Vol. 7734, International Society for Optics and Photonics, 77342I
Baron, F. & Young, J. S. 2008, in Optical and Infrared Interferometry, Vol. 7013, Inter-
national Society for Optics and Photonics, 70133X
153
154 BIBLIOGRAPHY
Baud, B. & Habing, H. 1983, Astronomy and Astrophysics, 127, 73
Beckers, J., Enard, D., Faucherre, M., et al. 1990, in SPIE Conf. Ser., Vol. 1236, Inter-
national Society for Optics and Photonics, 108–124
Beltran, M. T., Cesaroni, R., Neri, R., et al. 2005, Astronomy & Astrophysics, 435, 901
Berger, J. P. & Segransan, D. 2007, New Astronomy Reviews, 51, 576
Bhatnagar, S. & Cornwell, T. 2004, Astronomy & Astrophysics, 426, 747
Bieging, J. H., Rieke, M. J., & Rieke, G. H. 2002, Astronomy & Astrophysics, 384, 965
Bik, A. & Thi, W. 2004, Astronomy & Astrophysics, 427, L13
Bladh, S. & Höfner, S. 2012, Astronomy & Astrophysics, 546, A76
Bladh, S., Liljegren, S., Höfner, S., Aringer, B., & Marigo, P. 2019, Astronomy &
Astrophysics, 626, A100
Boffin, H. M., Hussain, G., Berger, J.-P., & Schmidtobreick, L. 2016, Astronomy at High
Angular Resolution: A Compendium of Techniques in the Visible and Near-Infrared,
Vol. 439 (Springer)
Bonaca, A., Tanner, J. D., Basu, S., et al. 2012, The Astrophysical Journal Letters, 755,
L12
Brown, A. G., Vallenari, A., Prusti, T., et al. 2021, Astronomy & Astrophysics, 649, A1
Buscher, D. F. 2015, Practical optical interferometry (Cambridge University Press)
Candes, E. J., Romberg, J. K., & Tao, T. 2006, Communications on Pure and Applied
Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences,
59, 1207
Cannon, E., Montargès, M., de Koter, A., et al. 2023, arXiv preprint arXiv:2303.08892
Carrillo, R. E., McEwen, J. D., & Wiaux, Y. 2014, Monthly Notices of the Royal Astro-
nomical Society, 439, 3591
Carroll, B. W. & Ostlie, D. A. 2017, An introduction to modern astrophysics (Cambridge
University Press), 82–83
BIBLIOGRAPHY 155
Cernicharo, J., Marcelino, N., Agúndez, M., & Guélin, M. 2015, Astronomy & Astro-
physics, 575, A91
Chandra, S., Maheshwari, V., & Sharma, A. 1996, Astronomy and Astrophysics Supple-
ment Series, 117, 557
Chauvin, G., Videla, M., Beust, H., et al. 2023, Astronomy & Astrophysics, 675, A114
Chiavassa, A., Freytag, B., Masseron, T., & Plez, B. 2011, A&A, 535, A22
Chiavassa, A., Kravchenko, K., & Goldberg, J. A. 2024, arXiv preprint arXiv:2402.00187
Chiavassa, A., Kravchenko, K., Millour, F., et al. 2020, Astronomy & Astrophysics, 640,
A23
Chiavassa, A., Kravchenko, K., Montargès, M., et al. 2022, Astronomy & Astrophysics,
658, A185
Chiavassa, A., Plez, B., Josselin, E., & Freytag, B. 2009, Astronomy & Astrophysics,
506, 1351
Climent, J., Wittkowski, M., Chiavassa, A., et al. 2020, A&A, 635, A160
Colom i Bernadich, M. 2020, Measuring the Characteristic Sizes of Convection Structures
in AGB Stars with Fourier Decomposition Analyses: the Stellar Intensity Analyzer
(SIA) Pipeline.
Cowling, T. G. 1934, Monthly Notices of the Royal Astronomical Society, 94, 768
Cox, J. P. 1963, The Astrophysical Journal, 138, 487
Cox, J. P. 1968
Cox, J. P. 2017, Theory of Stellar Pulsation.(PSA-2), Volume 2, Vol. 31 (Princeton
University Press)
Creevey, O., Thévenin, F., Berio, P., et al. 2015, Astronomy & Astrophysics, 575, A26
Cruzalèbes, P., Jorissen, A., Chiavassa, A., et al. 2015, Monthly Notices of the Royal
Astronomical Society, 446, 3277
156 BIBLIOGRAPHY
Demory, B.-O., Ségransan, D., Forveille, T., et al. 2009, Astronomy & Astrophysics, 505,
205
Deutsch, A. J. 1956, The Astrophysical Journal, 123, 210
Donoho, D. L. 2006, IEEE Transactions on information theory, 52, 1289
Drevon, J., Millour, F., Cruzalèbes, P., et al. 2022, Astronomy & Astrophysics, 665, A32
Ducati, J. 2002, VizieR Online Data Catalog
Dullemond, C., Juhasz, A., Pohl, A., et al. 2012, Astrophysics Source Code Library, ascl
Duvert, G., Young, J., & Hummel, C. A. 2017, Astronomy & Astrophysics, 597, A8
Eddington, A. 1919, Monthly Notices of the Royal Astronomical Society, 79, 177
Eddington, A. S. 1918, Monthly Notices of the Royal Astronomical Society, 79, 2
Eisenhauer, F., Monnier, J. D., & Pfuhl, O. 2023, arXiv preprint arXiv:2303.00453
Elitzur, M., Brown, J. A., & Johnson, H. R. 1989, The Astrophysical Journal, 341, L95
Engels, D., Kreysa, E., Schultz, G., & Sherwood, W. 1983, Astronomy and Astrophysics,
124, 123
Eriksson, K., Nowotny, W., Höfner, S., Aringer, B., & Wachter, A. 2014, Astronomy &
Astrophysics, 566, A95
Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, Publications of the
Astronomical Society of the Pacific, 125, 306
Freytag, B. & Höfner, S. 2008, Astronomy & Astrophysics, 483, 571
Freytag, B., Holweger, H., Steffen, M., & Ludwig, H.-G. 1997, in Science with the VLT
Interferometer: Proceedings of the ESO Workshop Held at Garching, Germany, 18–21
June 1996, Springer, 316–317
Freytag, B., Liljegren, S., & Höfner, S. 2017, Astronomy & Astrophysics, 600, A137
Freytag, B., Steffen, M., & Dorch, B. 2002, Astronomische Nachrichten, 323, 213
BIBLIOGRAPHY 157
Freytag, B., Steffen, M., Ludwig, H.-G., et al. 2012, Journal of Computational Physics,
231, 919
Freytag, B., Steffen, M., Wedemeyer-Böhm, S., et al. 2003, CO5BOLD User Manual
Gaia Collaboration. 2020, VizieR Online Data Catalog, I/350
Garsden, H., Girard, J., Starck, J.-L., et al. 2015, Astronomy & astrophysics, 575, A90
Gautschy-Loidl, R., Höfner, S., Jørgensen, U. G., & Hron, J. 2004, Astronomy & Astro-
physics, 422, 289
Geballe, T., Evans, A., van Loon, J. T., et al. 2007, in The Nature of V838 Mon and its
Light Echo, Vol. 363, 110
Gehrz, R. & Woolf, N. 1971, Astrophysical Journal, vol. 165, p. 285, 165, 285
Gillett, F., Low, F., & Stein, W. 1968, The Astrophysical Journal, 154, 677
Gilman, R. C. 1969, Astrophysical Journal, vol. 155, p. L185, 155, L185
Goldreich, P. & Scoville, N. 1976, Astrophysical Journal, 205, 144
Gonzalez, C. G., Absil, O., Absil, P.-A., et al. 2016, Astronomy & Astrophysics, 589,
A54
Goorvitch, D. 1994, The Astrophysical Journal Supplement Series, 95, 535
Groenewegen, M., Baas, F., Blommaert, J., et al. 1999, Astronomy and Astrophysics-
Supplement Series, 140, 197
Gull, S. F. & Skilling, J. 1984, in IEE Proceedings F-Communications, Radar and Signal
Processing, Vol. 131, IET, 646–659
Gull, S. F. & Skilling, J. 1999, Maximum Entropy Data Consultants Ltd., Suffolk, 1
Habing, H. J. & Olofsson, H. 2013, Asymptotic giant branch stars (Springer Science &
Business Media)
Hadjara, M., Cruzalèbes, P., Nitschelm, C., et al. 2019, Monthly Notices of the Royal
Astronomical Society, 489, 2595
158 BIBLIOGRAPHY
Haubois, X., Wittkowski, M., Perrin, G., et al. 2015, Astronomy & Astrophysics, 582,
A71
Hauschildt, P. H. & Baron, E. 1999, Journal of Computational and Applied Mathematics,
109, 41
Herwig, F. 2000, A&A, 360, 952
Hinkle, K. H., Hall, D., & Ridgway, S. 1982, The Astrophysical Journal, 252, 697
Höfner, S. 2008, Astronomy & Astrophysics, 491, L1
Höfner, S., Bladh, S., Aringer, B., & Ahuja, R. 2016, Astronomy & Astrophysics, 594,
A108
Höfner, S., Bladh, S., Aringer, B., & Eriksson, K. 2022, Astronomy & Astrophysics, 657,
A109
Höfner, S. & Olofsson, H. 2018, The Astronomy and Astrophysics Review, 26, 1
Högbom, J. 1974, Astronomy and Astrophysics Supplement Series, 15, 417
Ibrahim, N., Monnier, J. D., Kraus, S., et al. 2023, The Astrophysical Journal, 947, 68
Ireland, M., Scholz, M., Tuthill, P. G., & Wood, P. R. 2004a, Monthly Notices of the
Royal Astronomical Society, 355, 444
Ireland, M., Tuthill, P., Bedding, T., Robertson, J., & Jacob, A. 2004b, Monthly Notices
of the Royal Astronomical Society, 350, 365
Ireland, M., Tuthill, P., Davis, J., & Tango, W. 2005, Monthly Notices of the Royal
Astronomical Society, 361, 337
Ireland, M. J., Monnier, J. D., & Thureau, N. 2006, in Advances in Stellar Interferometry,
Vol. 6268, SPIE, 562–569
Ireland, M. J., Scholz, M., & Wood, P. R. 2008, Monthly Notices of the Royal Astro-
nomical Society, 391, 1994
Ireland, M. J., Scholz, M., & Wood, P. R. 2011, Monthly Notices of the Royal Astro-
nomical Society, 418, 114
BIBLIOGRAPHY 159
Isbell, J. W., Pott, J.-U., Meisenheimer, K., et al. 2023, Astronomy & Astrophysics, 678,
A136
Jennison, R. 1958, Monthly Notices of the Royal Astronomical Society, 118, 276
Jeong, K. S., Winters, J. M., Le Bertre, T., & Sedlmayr, E. 2003, Astronomy & Astro-
physics, 407, 191
Jorgensen, U. G. & Johnson, H. R. 1992, Astronomy and Astrophysics, 265, 168
Joyce, M. & Chaboyer, B. 2018a, The Astrophysical Journal, 864, 99
Joyce, M. & Chaboyer, B. 2018b, The Astrophysical Journal, 856, 10
Joyce, M. & Tayar, J. 2023, Galaxies, 11, 75
Kamezaki, T., Nakagawa, A., Omodaka, T., et al. 2012, Publications of the Astronomical
Society of Japan, 64, 7
Karovicova, I., Wittkowski, M., Boboltz, D., et al. 2011, Astronomy & Astrophysics,
532, A134
Karovicova, I., Wittkowski, M., Ohnaka, K., et al. 2013, Astronomy & Astrophysics,
560, A75
Kerschbaum, F., Maercker, M., Brunner, M., et al. 2017, Astronomy & Astrophysics,
605, A116
Kippenhahn, R., Weigert, A., & Weiss, A. 1990, Stellar structure and evolution, Vol.
192 (Springer)
Klement, R., Schaefer, G. H., Gies, D. R., et al. 2022, The Astrophysical Journal, 926,
213
Kluska, J., Malbet, F., Berger, J.-P., et al. 2014, Astronomy & Astrophysics, 564, A80
Kraus, S. 2012, in Optical and Infrared Interferometry III, Vol. 8445, International
Society for Optics and Photonics, 84451H
Kraus, S., Hofmann, K.-H., Menten, K. M., et al. 2010, Nature, 466, 339
Kupka, F. 2004, Proceedings of the International Astronomical Union, 2004, 119
160 BIBLIOGRAPHY
Labdon, A., Kraus, S., Davies, C. L., et al. 2023, Astronomy & Astrophysics, 678, A6
Labeyrie, A. 1975, The Astrophysical Journal, 196, L71
Lacour, S., Thiébaut, E., Perrin, G., et al. 2009, The Astrophysical Journal, 707, 632
Lada, C. 1987, in ISU Symp. 115, Vol. 1, D. Reidel
Lapeyrere, V., Kervella, P., Lacour, S., et al. 2014, in Optical and Infrared Interferometry
IV, Vol. 9146, International Society for Optics and Photonics, 91462D
Lawson, P. 2000, Principles of long baseline stellar interferometry (JPL)
Le Bouquin, J. B., Berger, J. P., Lazareff, B., et al. 2011, A&A, 535, A67
Le Bouquin, J.-B., Berger, J.-P., Lazareff, B., et al. 2011, Astronomy & Astrophysics,
535, A67
Le Bouquin, J.-B., Lacour, S., Renard, S., et al. 2009, Astronomy & Astrophysics, 496,
L1
Lèbre, A., Aurière, M., Fabas, N., et al. 2014, A&A, 561, A85
Lebzelter, T., Hinkle, K. H., Wood, P. R., Joyce, R. R., & Fekel, F. C. 2005, Astronomy
& Astrophysics, 431, 623
Ledoux, P. & Walraven, T. 1958, Handbuch der Physik, 51, 353
Lena, P. 1979, Journal of Optics, 10, 323
Liljegren, S., Höfner, S., Freytag, B., & Bladh, S. 2018, Astronomy & Astrophysics, 619,
A47
Liljegren, S., Höfner, S., Nowotny, W., & Eriksson, K. 2016, Astronomy & Astrophysics,
589, A130
Lopez, B., Lagarde, S., Petrov, R., et al. 2022, Astronomy & Astrophysics, 659, A192
Loup, C., Forveille, T., Omont, A., & Paul, J. 1993, Astronomy and Astrophysics Sup-
plement Series, 99, 291
Ludwig, H.-G., Caffau, E., Steffen, M., et al. 2009, arXiv preprint arXiv:0908.4496
BIBLIOGRAPHY 161
Maercker, M., Mohamed, S., Vlemmings, W., et al. 2012, Nature, 490, 232
Magic, Z., Weiss, A., & Asplund, M. 2015, Astronomy & Astrophysics, 573, A89
Malbet, F., Cotton, W., Duvert, G., et al. 2010, in Optical and Infrared Interferometry
II, Vol. 7734, International Society for Optics and Photonics, 77342N
Martin, R. G., Lubow, S. H., Vallet, D., Anugu, N., & Gies, D. R. 2023, The Astrophys-
ical Journal Letters, 957, L28
Matsuura, M., Barlow, M., Zijlstra, A., et al. 2009, Monthly Notices of the Royal As-
tronomical Society, 396, 918
McDonald, I., Zijlstra, A. A., & Boyer, M. L. 2012, Monthly Notices of the Royal
Astronomical Society, 427, 343
Medeiros, L., Lauer, T. R., Psaltis, D., & Özel, F. 2018, The Astrophysical Journal, 864,
7
Molster, F., Waters, L., & Kemper, F. 2010, T. Henning (Lecture Notes in Physics, Vol.
815
Monnier, J. D. 2003, Reports on Progress in Physics, 66, 789
Monnier, J. D., Berger, J.-P., Le Bouquin, J.-B., et al. 2014, in Optical and Infrared
Interferometry IV, Vol. 9146, International Society for Optics and Photonics, 91461Q
Montargès, M., Cannon, E., Lagadec, E., et al. 2021, Nature, 594, 365
Montargès, M., Kervella, P., Perrin, G., et al. 2014, Astronomy & Astrophysics, 572,
A17
Nelder, J. A. & Mead, R. 1965, The computer journal, 7, 308
Neugebauer, G. & Leighton, R. B. 1969, Two-micron sky survey: A preliminary cat-
alog, Vol. 3047 (Scientific and Technical Information Division, National Aeronautics
and . . . )
Newville, M., Stensitzki, T., Allen, D. B., et al. 2016, Astrophysics Source Code Library,
ascl
162 BIBLIOGRAPHY
Norris, B. R., Tuthill, P. G., Ireland, M. J., et al. 2012, Nature, 484, 220
Norris, R. P., Baron, F. R., Monnier, J. D., et al. 2021, The Astrophysical Journal, 919,
124
Nowotny, W., Aringer, B., Höfner, S., & Eriksson, K. 2013, Astronomy & Astrophysics,
552, A20
Nowotny, W., Aringer, B., Höfner, S., Gautschy-Loidl, R., & Windsteig, W. 2005a,
Astronomy & Astrophysics, 437, 273
Nowotny, W., Aringer, B., Höfner, S., & Lederer, M. T. 2011, Astronomy & Astrophysics,
529, A129
Nowotny, W., Höfner, S., & Aringer, B. 2010, Astronomy & Astrophysics, 514, A35
Nowotny, W., Lebzelter, T., Hron, J., & Höfner, S. 2005b, Astronomy & Astrophysics,
437, 285
Nürnberger, D. 2003, Astronomy & Astrophysics, 404, 255
o Garatti, A. C., Stecklum, B., Weigelt, G., et al. 2016, Astronomy & Astrophysics, 589,
L4
Ohnaka, K., Hofmann, K.-H., Benisty, M., et al. 2009, Astronomy & Astrophysics, 503,
183
Ohnaka, K., Weigelt, G., & Hofmann, K.-H. 2016, Astronomy & Astrophysics, 589, A91
Ohnaka, K., Weigelt, G., & Hofmann, K. H. 2017, Astronomy & Astrophysics, 597, A20
Ohnaka, K., Weigelt, G., & Hofmann, K.-H. 2019, The Astrophysical Journal, 883, 89
Ohnaka, K., Weigelt, G., Millour, F., et al. 2011, Astronomy & Astrophysics, 529, A163
Paladini, C. 2011, Interferometry of carbon rich AGB stars (na)
Paladini, C., Baron, F., Jorissen, A., et al. 2018, Nature, 553, 310
Pashchenko, I. N. 2019, Monthly Notices of the Royal Astronomical Society, 482, 1955
Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, the Journal of machine Learning
research, 12, 2825
BIBLIOGRAPHY 163
Perrin, G., Ridgway, S., Lacour, S., et al. 2020, Astronomy & Astrophysics, 642, A82
Perrin, G., Ridgway, S., Mennesson, B., et al. 2004, Astronomy & Astrophysics, 426,
279
Perrin, G., Ridgway, S., Verhoelst, T., et al. 2005, Astronomy & Astrophysics, 436, 317
Perrin, G., Ridgway, S. T., Lacour, S., et al. 2020, A&A, 642, A82
Pluzhnik, E. A., Ragland, S., Le Coroller, H., et al. 2009, The Astrophysical Journal,
700, 114
Price, S. D. & Walker, R. G. 1976, The AFGL Four Color Infrared Sky Survey: Cat-
alog of Observations at 4.2, 11.0, 19.8, and 27.4 Micrometers (Air Force Geophysics
Laboratories, Air Force Systems Command, United States . . . )
Ragland, S., Traub, W., Berger, J.-P., et al. 2006, The Astrophysical Journal, 652, 650
Rau, G., Hron, J., Paladini, C., et al. 2017, Astronomy & astrophysics, 600, A92
Rau, G., Paladini, C., Hron, J., et al. 2015, Astronomy & astrophysics, 583, A106
Rau, U. & Cornwell, T. J. 2011, Astronomy & Astrophysics, 532, A71
Rees, N. R., Izzard, R. G., & Karakas, A. I. 2024, Monthly Notices of the Royal Astro-
nomical Society, 527, 9643
Reimers, D. 1975, Circumstellar absorption lines and mass loss from red giants, Tech.
rep., Univ., Kiel
Renard, S., Thiébaut, E., & Malbet, F. 2011, Astronomy & Astrophysics, 533, A64
Renzini, A. & Voli, M. 1981, Astronomy and Astrophysics, 94, 175
Roddier, F. 1986, Optics Communications, 60, 145
Rodŕıguez-Coira, G., Paumard, T., Perrin, G., et al. 2021, arXiv preprint
arXiv:2105.09832
Rosales-Guzmán, A., Sanchez-Bermudez, J., Paladini, C., et al. 2023, A&A, 674, A62
Rosales-Guzmán, A., Sánchez-Bermúdez, J., Paladini, C., et al. 2023, Astronomy and
Astrophysics, 674, A62
164 BIBLIOGRAPHY
Rosseland, S. 1949, Oxford
Sacuto, S., Aringer, B., Hron, J., et al. 2011, Astronomy & Astrophysics, 525, A42
Salaris, M. & Cassisi, S. 2005, Evolution of stars and stellar populations (John Wiley &
Sons)
Samus, N., Kazarovets, E., Durlevich, O., Kireeva, N., & Pastukhova, E. 2017, Astron-
omy Reports, 61, 80
Sanchez-Bermudez, J., Hummel, C., Tuthill, P., et al. 2016a, Astronomy & Astrophysics,
588, A117
Sanchez-Bermudez, J., Millour, F., Baron, F., et al. 2018, Experimental Astronomy, 46,
457
Sanchez-Bermudez, J., Millour, F., Baron, F., et al. 2018, Experimental Astronomy, 46,
457
Sanchez-Bermudez, J., Thiébaut, E., Duvert, G., et al. 2016b, Reconstruction Test Re-
port And Data Processing Cookbook, Tech. rep., MPIA, Heidelberg, Germany
Sanchez-Bermudez, J., Weigelt, G., Bestenlehner, J. M., et al. 2018, Astronomy & As-
trophysics, 618, A125
Sánchez-Monge, Á., Cesaroni, R., Beltrán, M., et al. 2013, Astronomy & Astrophysics,
552, L10
Savage, B. D. & Sembach, K. R. 1991, The Astrophysical Journal, 379, 245
Schaefer, G., Brummelaar, T. t., Gies, D., et al. 2014, Nature, 515, 234
Schöller, M. 2007, New Astronomy Reviews, 51, 628
Schwarzschild, M. 1975, The Astrophysical Journal, 195, 137
Schwarzschild, M. 1975, The Astrophysical Journal, 195, 137
Schwarzschild, M. & Härm, R. 1965, The Astrophysical Journal, 142, 855
Shore, S. N. 2003, in Encyclopedia of Physical Science and Technology (Third Edition),
third edition edn., ed. R. A. Meyers (New York: Academic Press), 737–747
BIBLIOGRAPHY 165
Skilling, J. & Bryan, R. 1984, Monthly Notices of the Royal Astronomical Society, Vol.
211, NO. 1, P. 111, 1984, 211, 111
Solomon, P., Jefferts, K. B., Penzias, A. A., & Wilson, R. W. 1971, The Astrophysical
Journal, 163, L53
Sonoi, T., Ludwig, H.-G., Dupret, M.-A., et al. 2018, in PHysics of Oscillating STars.
Proceedings from the PHOST (PHysics of Oscillating STars) symposium hosted by
the Oceanographic Observatory in Banyuls-sur-mer (France) from 2-7 September 2018.
This conference honours the life work of Professor Hiromoto Shibahashi, 27
Strong, D. & Chan, T. 2003, Inverse problems, 19, S165
Sutton, E. C. & Wandelt, B. D. 2006, The Astrophysical Journal Supplement Series,
162, 401
Takeuti, M., Nakagawa, A., Kurayama, T., & Honma, M. 2013, Publications of the
Astronomical Society of Japan, 65, 60
Tallon-Bosc, I., Tallon, M., Thiébaut, E., et al. 2008, in Society of Photo-Optical Instru-
mentation Engineers (SPIE) Conference Series, Ser
Tayar, J., Somers, G., Pinsonneault, M. H., et al. 2017, The Astrophysical Journal, 840,
17
Thiébaut, É. 2009, New Astronomy Reviews, 53, 312
Thiébaut, É. 2013, European Astronomical Society Publications Series, 59, 157
Thompson, A., Moran, J. M., & Swenson Jr, G. W. 2017, Interferometry and synthesis
in radio astronomy (Springer Nature)
Trampedach, R., Asplund, M., Collet, R., Nordlund, Å., & Stein, R. F. 2013, The
Astrophysical Journal, 769, 18
Trampedach, R., Stein, R. F., Christensen-Dalsgaard, J., Nordlund, Å., & Asplund, M.
2014, Monthly Notices of the Royal Astronomical Society, 445, 4366
Tremblay, P.-E., Ludwig, H.-G., Freytag, B., Steffen, M., & Caffau, E. 2013, A&A, 557,
A7
166 BIBLIOGRAPHY
Tsuji, T. 2006, The Astrophysical Journal, 645, 1448
Ulrich, R. K. 1970, Astrophysics and Space Science, 7, 71
Vassiliadis, E. & Wood, P. 1993, Astrophysical Journal, Part 1 (ISSN 0004-637X), vol.
413, no. 2, p. 641-657., 413, 641
Vehoff, S., Hummel, C., Monnier, J., et al. 2010, Astronomy & Astrophysics, 520, A78
Viani, L. S., Basu, S., JM, J. O., Bonaca, A., & Chaplin, W. J. 2018, The Astrophysical
Journal, 858, 28
Whitelock, P. A., Feast, M. W., & Van Leeuwen, F. 2008, Monthly Notices of the Royal
Astronomical Society, 386, 313
Wickramasinghe, N., Donn, B., & Stecher, T. 1966, The Astrophysical Journal, 146, 590
Wilson, W. J. & Barrett, A. H. 1968, Science, 161, 778
Wittkowski, M., Boboltz, D. A., Driebe, T., et al. 2008, Astronomy & Astrophysics, 479,
L21
Wittkowski, M., Boboltz, D. A., Ireland, M., et al. 2011, Astronomy & Astrophysics,
532, L7
Wittkowski, M., Chiavassa, A., Freytag, B., et al. 2016, Astronomy & Astrophysics, 587,
A12
Wittkowski, M., Hauschildt, P., Arroyo-Torres, B., & Marcaide, J. 2012, Astronomy &
Astrophysics, 540, L12
Wittkowski, M., Hofmann, K.-H., Höfner, S., et al. 2017, Astronomy & astrophysics,
601, A3
Wittkowski, M. & Paladini, C. 2014, European Astronomical Society Publications Series,
70, 179
Wittkowski, M., Rau, G., Chiavassa, A., et al. 2018, Astronomy & Astrophysics, 613,
L7
Woolf, N. & Ney, E. 1969, Astrophysical Journal, vol. 155, p. L181, 155, L181
BIBLIOGRAPHY 167
Xiong, D. & Deng, L. 2007, Monthly Notices of the Royal Astronomical Society, 378,
1270
Yamamura, I., Dominik, C., De Jong, T., Waters, L., & Molster, F. 2000, Astronomy
and Astrophysics, 363, 629
Zoubir, A. M. & Iskander, D. R. 2004, Bootstrap techniques for signal processing (Cam-
bridge University Press)