{ "cells": [ { "cell_type": "markdown", "id": "02b47774", "metadata": {}, "source": [ "# Fractal Aggregate Spectra in PICASO\n", "\n", "In this tutorial, we will learn:\n", "\n", "1. How to run virga for non-spherical cloud particles. \n", "2. How to generate and compare spectra for fractal aggregates vs. spheres in PICASO\n", "3. How to compare the particle sizes and optical properties for aggregates vs. spheres\n", "\n", " This tutorial builds directly on the methods in tutorial \"Fractal Aggregate Aerosols in a Hot Jupiter\", using the same planet and conditions as an example. However in this tutorial, the atmosphere is set-up within PICASO and virga is called using PICASO's built-in function. The varied spectra from each choice of particle shape are then compared. \n", "\n", "Aerosols can come in a variety of shapes and sizes, which we typically characterise by their fractal dimension (Df). Under this framework, aerosols are made of collections of small spheres called monomers, and Df=3 would represent the most compact possible arrangement, whereas Df=1 would be something perfectly linear. Real aerosols are usually somewhere in-between 1 and 3, and the ultimate shape of an aerosol is usually a result of their formation conditions. VIRGA v2.0 allows us to silmultaneously model the effects of a range of different particle shapes on spectra so we can find the best fit to observational data. For further details and example of typical shapes, see Lodge & Moran et al. (2025).\n", "\n", "Striving to maintain the simplicity and ethod of virga, aerosols can be added by providing two simple parameters:\n", "\n", " a) The fractal dimension (general shape) of the aerosols (which can be an array).\n", " b) EITHER the number of monomers (N_mon) OR the radius of the monomer (r_mon).\n", "\n", "Figure 2 of Moran & Lodge et al. (2025) demonstrates both of the choices for (b) and explains them in detail -- see this for more guidance. in this tutorial we will choose an array of fractal dimensions, and we choose to fix the number of monomers ($N_{mon}$=1000).\n", "\n", "Before beginning, make sure you have downloaded the new optical database for aggregates from here: https://zenodo.org/records/16581692 (which are all calculated with the assumption that $N_{mon}$=1000) and update the filepaths to the VIRGA/PICASO reference data in the first cell:" ] }, { "cell_type": "code", "execution_count": null, "id": "9fb3885c", "metadata": {}, "outputs": [], "source": [ "# --------------------------------------------- UPDATE THESE FILEPATHS FOR YOUR OWN SETUP ------------------------------------------------------------------\n", "\n", "# first, set the locations of databases needed for picaso if you have not done this already via the code installation\n", "#import os\n", "#os.environ['picaso_refdata'] = '/PICASO/reference' # location of reference data\n", "#os.environ['PYSYN_CDBS'] = '/PICASO/PYSYN_CDBS' # location of stellar data\n", "mieff_dir = '/Users/nbatalh1/Documents/data/VIRGA_2_mieff_files' # set file paths for VIRGA v2 optical database (.mieff files) -- make sure you have downloaded the new database that includes aggregates (v2!)\n", "\n", "# --------------------------------------------------------------------------------------------------------------------------------------------" ] }, { "cell_type": "code", "execution_count": null, "id": "944f9007", "metadata": {}, "outputs": [], "source": [ "#standard import statements\n", "import numpy as np\n", "import pandas as pd\n", "import astropy.units as u\n", "\n", "#radiative transfer and atmosphere code\n", "from picaso import justdoit as pdi\n", "\n", "#plotting tools\n", "from picaso import justplotit as ppi\n", "from virga import justplotit as vpi\n", "ppi.output_notebook()\n", "from bokeh.plotting import show, figure\n", "from bokeh.io import output_notebook\n", "import matplotlib.pyplot as plt " ] }, { "cell_type": "markdown", "id": "aa4a7158", "metadata": {}, "source": [ "## Run PICASO to create the clear (without condensate) transmission spectrum\n", "\n", "To start with, let's create a clear (cloud-free) transmission spectrum using the built-in Hot Jupiter profile within Picaso." ] }, { "cell_type": "code", "execution_count": null, "id": "33c11958", "metadata": {}, "outputs": [], "source": [ "# First, we import the opacity values for the models of the stars. We do this using a database stored locally on our computer, and it just imports the values \n", "# for the specified range of wavelengths that we are interested in (rather than all possible opacities for all star types at all wavelengths)\n", "\n", "opacity = pdi.opannection(wave_range=[0.3,15]) # import opacities for stellar models between the range: 0.3 um -> 15 um\n", "jupiter = pdi.inputs() # initialise planet atmosphere in PICASO\n", "\n", "jupiter.phase_angle(0) # set phase angle to 0\n", "\n", "# SET PLANET PARAMETERS: \n", "\n", "planet_mass = 1 # M_Jup\n", "planet_radius = 1 # R_Jup\n", "\n", "jupiter.gravity(mass=planet_mass, mass_unit=pdi.u.Unit('M_jup'), radius=planet_radius, radius_unit=pdi.u.Unit('R_jup')) # order: planet mass, unit, planet radius, unit\n", "\n", "# INPUT STELLAR PARAMETERS: \n", "\n", "stellar_temperature = 5700 # K\n", "stellar_metallicity = 0 # (Fe/H) dex\n", "stellar_surface_gravity = 4.5 # (log g) cgs units\n", "stellar_radius = 1 # R_Sun\n", "\n", "jupiter.star(opacity, stellar_temperature, stellar_metallicity, stellar_surface_gravity, radius=stellar_radius, radius_unit = pdi.u.Unit('R_sun') ) # order: Opacity database, Temperature, Metallicity, Surface gravity (log g), Radius, Radius unit\n", "\n", "# ADD THE ATMOSPHERIC CHEMISTRY \n", "\n", "jupiter.atmosphere(filename=pdi.HJ_pt(), sep=r'\\s+')" ] }, { "cell_type": "markdown", "id": "8898ace5", "metadata": {}, "source": [ "Our planet is complete! Let's compute and plot our cloud-free transmission spectrum:" ] }, { "cell_type": "code", "execution_count": null, "id": "72faf532", "metadata": {}, "outputs": [], "source": [ "# Ask PICASO to do all of the radiative transfer calculations and compute the transmission spectrum for the clear model\n", "clear_spectrum= jupiter.spectrum(opacity, full_output=True,calculation='transmission')\n", "\n", "# create list to store all spectra\n", "all_spectra=[]\n", "\n", "# regrid to resolution R=300\n", "wno, rprs2 = pdi.mean_regrid(clear_spectrum['wavenumber'] , clear_spectrum['transit_depth'], R=300)\n", "full_output = clear_spectrum['full_output']\n", "\n", "# store cloudy spectrum for clear case\n", "all_spectra.append(rprs2*1e6)\n", "\n", "# set custom colors for cloud species\n", "cloud_colors = ['#ce3fce', '#7cf2fb', '#33b9cc', '#1c7bb8', '#295091', '#212255'] # order: spheres, 1.2, 1.6, 2.0, 2.4, 2.8\n", "spectra_colors = ['#bcbcbd'] + cloud_colors # add an extra element for the clear model for when plotting spectra. New order: clear, spheres, 1.2, 1.6, 2.0, 2.4, 2.8\n", "\n", "# plot spectrum\n", "ppi.show(ppi.spectrum(wno,rprs2*1e6,plot_width=900,plot_height=300, palette=spectra_colors[0:1]))" ] }, { "cell_type": "markdown", "id": "070aa3e5", "metadata": {}, "source": [ "## Adding Mg2SiO4 Clouds with Spherical Particles\n", "\n", "Mg2SiO4 is a typical cloud species in Hot Jupiter environments -- let's add some spherical cloud particles to start with. This time, we pass all important cloud variables directly to the .virga function from within PICASO." ] }, { "cell_type": "code", "execution_count": null, "id": "d9b200e0", "metadata": {}, "outputs": [], "source": [ "# Set Kzz, metallicity and mean molecular weight \n", "metallicity = 1 #atmospheric metallicity relative to Solar\n", "mean_molecular_weight = 2.3 # atmospheric mean molecular weight\n", "constant_kzz = 10**10 # constant Kzz\n", "f_sed = 0.5 # sedimentation efficiency\n", "sigma_width = 2 # width of lognormal distribution\n", "aggregate_species = ['Mg2SiO4'] # Aggregate species\n", "\n", "kz = np.zeros(len(jupiter.inputs['atmosphere']['profile']))+ constant_kzz # make a constant Kzz profile as an array\n", "jupiter.inputs['atmosphere']['profile']['kz'] = kz # add this Kzz profile to the PICASO inputs\n", "\n", "# make a copy of the clear spectrum\n", "cloudy_jupiter_spheres = jupiter\n", "\n", "# run virga: this will determine where the clouds are, what particles sizes they are made from, and how opaque they are at each pressure level\n", "clouds_from_virga_spheres = cloudy_jupiter_spheres.virga(aggregate_species, directory=mieff_dir, \n", " fsed=f_sed, mh=metallicity, mmw = mean_molecular_weight, \n", " sig=sigma_width, aggregates=False) # Add the Mg2SiO4 clouds! The function '.virga' runs VIRGA from within PICASO, using all of the same atmosphere values etc, so that we don't have to set them again\n", "\n", "# calculate a transmission spectrum for this new (cloudy) model (with added Mg2SiO4) and regrid to a resolution of R=300\n", "cloudy_spectrum = cloudy_jupiter_spheres.spectrum(opacity, full_output=True, calculation='transmission')\n", "cloud_wno,cloud_rprs2 = pdi.mean_regrid(cloudy_spectrum['wavenumber'] , cloudy_spectrum['transit_depth'], R=300)\n", "\n", "# store cloudy spectrum for spheres\n", "all_spectra.append(cloud_rprs2*1e6)\n", "\n", "# show transmission spectrum\n", "print(f'For Spheres:')\n", "ppi.show(ppi.spectrum([wno,cloud_wno],[(rprs2*1e6),(cloud_rprs2*1e6)],\n", " legend=['No Mg2SiO4','With Mg2SiO4'],plot_width=900,plot_height=300, palette=spectra_colors[0:2])) # palette argument simply assigns custom colors to match plots below" ] }, { "cell_type": "markdown", "id": "1ffb7f72", "metadata": {}, "source": [ "## p-T profile\n", "\n", "To view the p-T profile and get an idea of where in the atmosphere our Mg2SiO4 condensates are forming into clouds, we can use virga's usual pt() function." ] }, { "cell_type": "code", "execution_count": null, "id": "40b7bc81", "metadata": {}, "outputs": [], "source": [ "# Plot pressure-Temperature profile\n", "show(vpi.pt(clouds_from_virga_spheres,plot_height=450))" ] }, { "cell_type": "markdown", "id": "a624c429", "metadata": {}, "source": [ "## ADDING FRACTAL AGGREGATE CLOUDS\n", "\n", "Now let's add fractal aggregate clouds to our clear model instead of spheres." ] }, { "cell_type": "code", "execution_count": null, "id": "14d2195b", "metadata": {}, "outputs": [], "source": [ "# Set a range of fractal dimensions to explore (a different model will be run each time). This can be a single value or an array. \n", "# In this example, we will choose the values [1.2, 1.6, 2.0, 2.4, 2.8]:\n", "d_f_values = np.linspace(1.2,2.8,5)\n", "\n", "# make a copy of the clear spectrum for each d_f value\n", "cloudy_jupiter_list=[]\n", "for i in range(len(d_f_values)):\n", " cloudy_jupiter_list.append(jupiter) # copy the clear spectrum\n", "\n", "# create a list to store all of the different cloud models for each fractal dimension\n", "clouds_from_virga_fractals=[]\n", "\n", "# loop over all particle shapes\n", "for i in range(len(d_f_values)):\n", "\n", " d_f = d_f_values[i]\n", " N_monomers = 1000 # Here we provide the number of monomers. If you prefer to provide r_mon instead, set it here (but note that it needs to be in cm, not um), and submit r_mon instead of N_mon in the .virga() function below\n", "\n", " # That's it -- Df and N_mon are all you need to prescribe for aggregates! More details on other typical aggregate parameters below for experts:\n", " # k0 will usually be calculated using Eq. 14+15 of Moran & Lodge et al. (2025), but this can be prescribed here if you wish to use a fixed value.\n", " # Radius of gyration will be calculated within the code: virga can calculate the mass of the each aggregate, meaning that r_mon can be calculated\n", " # from N_mon (or vice versa). k0 is calculated as stated above, so including Df we have everything we need to find R_gyr using Eq. 2 of Moran & Lodge \n", " # et al. (2025)\n", "\n", " # take the clear atmosphere (cloudy_jupiter_list[i]) and run virga to find the clouds for each fractal dimension (using aggregates=True and our aggregate parameters). \n", " # this method saves each cloud model as a new model in a list (clouds_from_virga_fractals).\n", " clouds_from_virga_fractals.append(cloudy_jupiter_list[i].virga(aggregate_species, # Add the Mg2SiO4 clouds! The function '.virga' runs VIRGA from within PICASO, using all of the same atmosphere values etc, so that we don't have to set them again\n", " directory=mieff_dir, # location of .mieff files\n", " fsed=f_sed, # f_sed\n", " mh=metallicity, # metallicity\n", " mmw = mean_molecular_weight, # mean molecular weight of atmosphere\n", " sig=sigma_width, # width of lognormal distribution\n", " aggregates=True, # model cloud particles as aggregates\n", " Df=d_f, # fractal dimension of aggregates\n", " N_mon=N_monomers)) # number of monomers in an aggregate\n", "\n", " # calculate a transmission spectrum for this new (cloudy) model (with added Mg2SiO4) and regrid to a resolution of R=300\n", " cloudy_spectrum = cloudy_jupiter_list[i].spectrum(opacity, full_output=True, calculation='transmission')\n", " cloud_wno,cloud_rprs2 = pdi.mean_regrid(cloudy_spectrum['wavenumber'] , cloudy_spectrum['transit_depth'], R=300)\n", " \n", " # store spectrum for each fractal\n", " all_spectra.append(cloud_rprs2*1e6)\n", " \n", " # show transmission spectrum for each fractal dimension\n", " print(f'For d_f = {d_f}:')\n", " ppi.show(ppi.spectrum([wno,cloud_wno],[(rprs2*1e6),(cloud_rprs2*1e6)],\n", " legend=['No Mg2SiO4','With Mg2SiO4'],plot_width=900,plot_height=300, palette=list(spectra_colors[i] for i in [0,(i+2)]))) # the palette argument iterates through the list of custom colors so that fractal has a new color" ] }, { "cell_type": "markdown", "id": "f3ff32d1", "metadata": {}, "source": [ "# Aggregate Analysis Tools\n", "\n", "We've made our clouds! Now let's see how the aggregates compared to spheres. There are a number of tools that have been added to compare between the different models, and each one can give you interesting diagnostic information. Let's explore them one-by-one. First, set-up the figure size, location of the directories that hold the data, and the min/max wavelength to plot." ] }, { "cell_type": "code", "execution_count": null, "id": "b0597dad", "metadata": {}, "outputs": [], "source": [ "# set figure size in notebook\n", "plt.rcParams['figure.figsize'] = [15, 10]\n", "\n", "# set min and max wavelength to show on plots\n", "min_wavelength = 0.3 # in um\n", "max_wavelength = 15 # in um" ] }, { "cell_type": "markdown", "id": "b04b4b71", "metadata": {}, "source": [ "## Comparing the Spectra\n", "\n", "This tool plots the transmission spectrum for each of the different fractal dimensions, on log-axes this time (so they look a little different to the plots above)." ] }, { "cell_type": "code", "execution_count": null, "id": "c2cc29a5", "metadata": {}, "outputs": [], "source": [ "# plot spectra - compare the transmission spectrum for each of the different fractal dimensions, on log x-axis this time\n", "\n", "ppi.show(ppi.spectrum([wno]*7,[all_spectra],\n", " legend=['clear', 'spheres', '1.2', '1.6', '2.0', '2.4', '2.8'],\n", " plot_width=1100,plot_height=600, \n", " x_axis_type='log', \n", " palette=spectra_colors))" ] }, { "cell_type": "markdown", "id": "14ee5289", "metadata": {}, "source": [ "Interpretation: In this example, the Df=1.2 shapes are the least opaque and have the smallest transit depth of all of the different shapes modeled. Note that this isn't always the case (see Lodge & Moran 2025 for examples where the opposite is true!)." ] }, { "cell_type": "markdown", "id": "80550169", "metadata": {}, "source": [ "## Visualising Particle Sizes in each Pressure Layer\n", "\n", "Virga will typically find different particle sizes for different aggregate shapes because of how they behave dynamically (see Moran & Lodge et al. 2025). Lower fractal dimensions tend to form into larger particles (because they have slower fall speed vs compact spheres. However, the upwards forces are the same in both versions of the model. Therefore, the aggregates need to form into more massive particles to stay stable against the same upwards velocity from Kzz/convection). The mean particle sizes that virga found for each shape can be seen and compared using the tool below. " ] }, { "cell_type": "code", "execution_count": null, "id": "48f1c63d", "metadata": {}, "outputs": [], "source": [ "vpi.aggregates_pressure_vs_radius(clouds_from_virga_spheres, clouds_from_virga_fractals, d_f_list=d_f_values, colors=cloud_colors)" ] }, { "cell_type": "markdown", "id": "1426452b", "metadata": {}, "source": [ "\n", "Interpretation: The Df=1.2 model has the largest particles, which are at the bottom of the cloud-deck. In all models, cloud particles get smaller as altitude increases." ] }, { "cell_type": "markdown", "id": "ff27cb13", "metadata": {}, "source": [ "## Probing the relationship between Pressure, Radius and Number Density\n", "\n", "Because the same amount of condensate is condensing out in each model in virga (regardless of particle shape), if the low fractal dimensions form into larger particles, there must be fewer of them. Plotting the number density at each pressure layer shows demonstrates this balance between particle size and number density visually with altitude. The marker size indicates particle size. " ] }, { "cell_type": "code", "execution_count": null, "id": "52ca3460", "metadata": {}, "outputs": [], "source": [ "vpi.aggregates_pressure_vs_number_density(clouds_from_virga_spheres, clouds_from_virga_fractals, d_f_list=d_f_values, colors=cloud_colors)" ] }, { "cell_type": "markdown", "id": "c8bf4e5f", "metadata": {}, "source": [ "Interpretation: The Df=1.2 clouds are shown to be made of fewer (but larger) particles than in the other models." ] }, { "cell_type": "markdown", "id": "5aed211d", "metadata": {}, "source": [ "## Visualising .mieff files and Optical Properties\n", "\n", "This tool allows you to see how the optical properties of the cloud particles vary with shape. $Q_{ext}$ represents how efficiently they will absorb and scatter radiation, and it is a function of particle size, wavelength and shape. This reads all of the data from the .mieff database and plots $Q_{ext}$ for all particle sizes (blue = smallest particles in the grid, green = largest particles in the grid).\n" ] }, { "cell_type": "code", "execution_count": null, "id": "34549e44", "metadata": {}, "outputs": [], "source": [ "# plot Q_ext as a function of wavelength for each fractal dimension\n", "vpi.aggregates_optical_properties(aggregate=aggregate_species[0], mieff_dir=mieff_dir, d_f_list=d_f_values, min_wavelength=min_wavelength, max_wavelength=max_wavelength)" ] }, { "cell_type": "markdown", "id": "2d9145b9", "metadata": {}, "source": [ "Interpretation: This is usually a busy plot, but gives some intuition into how linear and compact aggregates compare to their spherical counterparts. Here the linear fractals with Df = 1.2 are the most opaque at short wavelengths. We saw that these particles are the least opaque in the transmission spectrum, which is because their number density is so much lower as shown in the previous figure.\n", "\n", "These plots can highlight resonance features and aid in checking that your optics are behaving as expected for a particular chemical species." ] }, { "cell_type": "markdown", "id": "fc8990e6", "metadata": {}, "source": [ "## Checking that Particle Radii are within the .mieff Grid\n", "\n", "It's important that virga finds particles that are within the range of optical properties in the .mieff grid, otherwise virga may underpredict or overpredict the opacity of the particles. Virga will give warnings if it finds particles off-grid, but this plot is a quick visual check to see if and where that happens:" ] }, { "cell_type": "code", "execution_count": null, "id": "37636d6b", "metadata": {}, "outputs": [], "source": [ "vpi.aggregates_wavelength_radius_grid(clouds_from_virga_spheres, clouds_from_virga_fractals, aggregate=aggregate_species[0], d_f_list=d_f_values, mieff_dir=mieff_dir, colors=cloud_colors)" ] }, { "cell_type": "markdown", "id": "ec2239e0", "metadata": {}, "source": [ "Interpretation: \n", "- The red outline box marks the extremes of the .mieff grid, with faint red vertical lines representing the grid spacing in radius.\n", "- The vertical extent of the box is just to show the wavelength range in the .mieff files. The colored markers show the column density so that we have an indication of how many particles were in each radius bin for each model.\n", "- The \"radius\" of an aggregate is defined as the compact radius (the radius of a sphere made from the same volume of material) and this is how they are stored in the .mieff files. The particle sizes obtained by virga at each pressure layer and for each fractal dimension are plotted on this grid as markers, and these should all ideally be inside the horizontal bounds of the red box. If not, you may need to adjust your Kzz/fsed values (both of these can effect particle size), or make a new .mieff grid (see next tutorial!). It's ok for the particles to be near the vertical bounds of the red box -- the marker y-positions just scale with column density, so they will always sit on the top/bottom edge of the box." ] } ], "metadata": { "kernelspec": { "display_name": "pic312", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.2" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": true, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 5 }