Challenges and Techniques for Simulating Line Emission

Modeling emission lines from the millimeter to the UV and producing synthetic spectra is crucial for a good understanding of observations, yet it is an art filled with hazards. This is the proceedings of"Walking the Line", a 3-day conference held in 2018 that brought together scientists working on different aspects of emission line simulations, in order to share knowledge and discuss the methodology. Emission lines across the spectrum from the millimeter to the UV were discussed, with most of the focus on the interstellar medium, but also some topics on the circumgalactic medium. The most important quality of a useful model is a good synergy with observations and experiments. Challenges in simulating line emission are identified, some of which are already being worked upon, and others that must be addressed in the future for models to agree with observations. Recent advances in several areas aiming at achieving that synergy are summarized here, from micro-physical to galactic and circum-galactic scale.


Introduction
Line emission from the Interstellar Medium (ISM) of galaxies carries information that is crucial in understanding galaxy evolution. Observing line emission across the electromagnetic spectrum allows us to characterize the mass, composition, and chemical state of the ISM, as well as to trace galaxy properties such as star formation rate (SFR), metallicity and dynamics. For example, the emission from major cooling lines, such as Hα or [C II], is sensitive to the physical conditions (densities, radiation field) and dynamics of the ISM. In addition, emission lines work on all physical scales, from galaxy dynamics and inflows to turbulent and collapse motions in star-forming clouds and cores. By systematically comparing spectral-line signatures of different physical models, one can correctly identify the physical processes occurring in these regions. Furthermore, the emission from ionized interstellar gas contains particularly valuable information about the nature of the ionizing radiation sources in a galaxy. In fact, prominent optical emission lines are routinely used to estimate whether ionization is dominated by young massive stars (tracing SFR), an AGN or evolved, post-asymptotic giant branch (post-AGB) stars. Three of the most widely used line-ratio diagnostic "BPT" diagrams 1 , relate the [OIII]/Hβ ratio to the [NII]/Hα, [SII]/Hα and [OI]/Hα ratios. These diagrams have proven useful in identifying the nature of the ionizing radiation in large samples of galaxies in the local Universe [2,3]. Complementary to line emission are the observations of absorption lines of the circumgalactic medium (CGM), which can give key information on the history of the feedback, in terms of chemical, ionization, and thermodynamical state of the outflowing/inflowing gas, that regulates the star formation process. Gas kinematics, from both emission and absorption, give information about large scale gas flows. Thus galactic outflows, from active galactic nuclei (AGN) and starbursts, can be combined with CGM absorption line observations, to study the star formation history, AGN activity history, and feedback processes that regulate both the evolution of the galaxy and its environment.
Looking back on the past three decades, the approach to creating synthetic observations of line emission has gone from simplified analytical modeling to complex simulations, increasing the reliability of the results. Since many important cooling lines emerge from the photodissociation regions (PDRs) of the ISM, modeling line emission from the PDRs has been an active field of research since the basic 1D PDR models came into place in the 80's (e.g., [4][5][6][7][8]). This early modeling work included the reprocessing of starlight in the UV to infrared continuum by dust grains and polycyclic aromatic hydrocarbons (PAHs), and the stratified layers of increasingly photo-dissociated species as one moves through the neutral gas of a cloud towards the HII region. This picture is illustrated in the top panel of Figure 1. In the late 90s, the models presented in the 80's [4,9] were used as basis for a modeling effort to create line ratio diagnostic plots of the [O I], [C II], C I and CO FIR emission lines [10]. Later the modeling was improved with the online PDR Toolbox 2 as a result [11,12].  [13] for model V1 of [14], adapted to extend to A V = 100. Still, one of the main problems when comparing these models to actual observations is that the models are based on patches of gas with constant density, metallicity, and radiation field strength whereas any given galaxy (or region in a galaxy) will be a superposition of different gas states and radiation conditions as illustrated in Figure 2 (Ideally, each cloud should also have a radial gradient, but this feature has been omitted from the Figure for the sake of clarity). With improved simulations of galaxy (and cloud) formation, it became possible in the 21st century to use galaxy/cloud simulations as direct input to photoionization codes, that also calculate the line emission (e.g., [15][16][17][18][19][20][21][22][23]). This new approach gave way for a more realistic picture in which each region of the ISM is approximated by a combination of different gas conditions. Some of these simulations include radiative transfer and the non-equilibrium impact on the gas temperature of the local radiation field generated by nearby star formation (e.g., [24]), allowing for more accurate estimation of the nebular line emission. This method is being extended to simulations of the CGM, with increased resolution to capture the impact of outflows on accreting structure (see Section 4.4). Improved numerical techniques such as those mentioned above have also meant a slight division in the field of line modeling, between groups would who specialize in creating the photoionization codes and groups that apply these codes to hydrodynamical simulations. The photoiononization codes have evolved to also calculate line emission and are kept updated as experimental values for collision strengths, chemical rate coefficients and photoelectric heating rates are being revised and databases of atomic and molecular line data increase in size (e.g., [25][26][27][28][29][30][31][32]).
Similar to the case of PDRs, significant progress has been achieved over the last few years in modelling nebular emission from ionized regions around young star clusters, AGN and post-AGN stellar populations in a full cosmological context. Yet, fully self-consistent models of this kind are currently limited by the performance of cosmological radiation-hydrodynamic simulations and insufficient spatial resolution on the scales of individual ionized regions around stars and active nuclei. To circumvent these limitations, some pioneer studies proposed the post-processing of cosmological hydrodynamic simulations and semi-analytic models with photoionization models to compute the cosmic evolution of nebular emission [33][34][35]. Some further improvement has been achieved recently by [36], not only accounting for the integrated nebular emission from young stars (as in [34,35]), but also from AGN and post-AGB stars, based on the star formation, AGN luminosity, gas density, and chemical enrichment histories of the simulated galaxies.
Across different aspects of simulating line emission, key questions remain to be answered: • What are the best emission lines to trace various ISM properties and ionizing sources in galaxies? • How can we use emission lines to trace feedback and ISM evolution with redshift? • How should absorption features be correctly interpreted? • Where do we stand in deriving sub-grid physics and comparing codes? • How do we coordinate our efforts?
To address such questions, the conference "Walking the Line" was held in Phoenix, Arizona on 14-16 March 2018. The main topic was Simulating Line Emission from Galaxies 3 and the program 4 included 35 min presentations from three invited speakers and 15 min presentations from 24 of the remaining 27 participants. In total, nine participants were PhD students and 1/3 of all participants were women. Most talks are now publicly available with video recordings online 5 . In addition, the conference had daily discussion sessions stimulating in-depth exploration of new topics and sparking new connections. This way, the meeting was an opportunity to discuss the challenges that need to be solved in order to make more realistic simulations and a fairer comparison with observations.
To take a specific case study where line emission has come to play an important role, one can consider the problem of measuring gas mass, which became a common thread for a large part of the workshop. Measuring the gas mass of galaxies, the fuel for star formation, as a function of cosmic time is a crucial component in understanding the formation and evolution of galaxies. Different groups have observed the gas content of main-sequence galaxies through cosmic time (the bulk of the galaxy population that is responsible for forming new stars at any epoch in cosmic history) either through their 12 CO (hereafter CO) emission or the dust continuum (e.g., [37][38][39][40][41][42][43][44][45][46][47]). Both approaches (CO and dust) rely on uncertain conversion factors between the observed luminosity and an estimated gas mass. ISM physical processes can further complicate the use of CO and dust continuum to reliably measure the gas mass of galaxies. For example, the contrast of the CO emission line and dust continuum against the cosmic microwave background becomes lower at higher redshifts (e.g., [23,48]) and under the influence of cosmic rays the CO molecule can be destroyed [49,50]. Reliable alternative measures of the gas mass in a well defined sample of galaxies are therefore essential to overcome the systematic uncertainties in the CO and dust conversion factors and to overcome the CMB contrast. Alternative options seen in the literature include for instance the emission from [CI] [51,52], [CII] [53], from H 2 O (e.g., [54]) and PAHs features (e.g., Cortzen et al. submitted 2018), and from optically thin isotopologues (e.g., [55]).The synergy between galaxy formation simulations, ISM chemistry simulations, and radiative transfer codes has the power to pave the way towards reliably measuring gas masses. This synergy allows for a controlled setting in which the conversion between sub-mm line/continuum strength and H 2 mass can be explored. This is important to (1) quantify under which physical conditions classical approaches to measure gas masses such as CO and dust continuum emission break down and (2)  In this paper we summarize the topics covered in the conference, with particular attention to the break-out sessions of this workshop. We report our conclusions on the state of the art of modelling emission lines, by going to progressively larger physical scales (Sections 2-4), concluding by discussing the possible ways to move forward as a community (Section 5).
Note that the number of emission (and absorption) lines observed from the ISM of galaxies is ever increasing, and we did not attempt to cover them all with this workshop. However, to give the reader a quick overview of the lines that were discussed at the workshop and are commonly used to diagnose the ISM and/or CGM, we present in Table 1 a list of lines that will be mentioned in this paper.  (4) (1) We give air wavelengths above 2000 Å and below 20,000 Å, and vacuum wavelengths otherwise; (2) Collisionally excited line; (3) https://physics.nist.gov/PhysRefData/ASD/lines_form.html; (4) http://home. strw.leidenuniv.nl/~moldata/.

Tools for Solving Level Populations and Line Excitation
Several software tools are available online for the calculation of line emission. While not covering all, the tools presented and discussed at this workshop are listed alphabetically in Table 2 for quick reference and comparison. Before comparing the tools, it is worth noting that they are built to achieve different levels of precision at the cost of computation time. As such, codes like LIME, MOLLIE, and RADMC-3D are very accurate in 3D, but slow compared to CLOUDY and DESPOTIC which solve for chemistry and temperature in 1D and are intended for larger parameter studies that can be analyzed fast.
In terms of chemical networks, CLOUDY is the most comprehensive tool with 625 species and a large range of permitted densities that allows for application to models on scales from clouds and HII regions to galaxies (e.g., [58][59][60]). DESPOTIC has been applied in post-processing to galaxy simulations mainly (e.g., [21,22,61]) but also to smaller regions such as the Galactic central molecular zone [62]. As mentioned above, CLOUDY and DESPOTIC are restricted to problems in 1D whereas the remaining tools in Table 2 work in 3D. LIME, MOLLIE, and RADMC-3D do exact radiative transfer in 3D and are typically used on smaller scales where non-symmetric features such as filaments and turbulence become important (e.g., [63,64]). ART 2 also does radiative transfer in 3D including the continuum from far-UV to radio wavelengths and incorporating the resonant line Lyα (e.g., [29,65,66]) (an updated version of the code including CO and some prominent far-infrared fine-structure lines such as CII, OI, OIII and NII will be out later this year). Finally, MAIHEM is unique in taking into account the non-equilibrium effects of turbulence on line emission (c.f. Section 3.3).  (8) http://www.ita.uni-heidelberg.de/~dullemond/software/radmc-3d/; (9) Escape probability formalism.

Correcting Heavy-Element Energy Levels
Predictions on line emission modeling depends on the reliability of fundamental atomic data. Software tools such as those discussed in the previous section use theoretical and experimental atomic and molecular rates and energies to calculate emissivities and line wavelengths respectively, mostly from external databases as CHIANTI [69,70] or NIST. Nowadays, theoretical energy levels of heavy many-electron ions disagree with experimental values, leading to incorrectly predicted wavelengths. This happens in all the energy and wavelength ranges. For instance, there are differences between the energies of Fe II reported by the NIST and CHIANTI databases, where several terms reported by NIST are missing in CHIANTI. Similarly, efforts to fit the spectra of X-ray satellite observations have accentuated the necessity of including transitions between levels with high main quantum number n in the theoretical calculations (e.g., [71][72][73][74]). A displacement of about 6 Å has been found in some X-ray lines due to theoretical and experimental differences for energy levels [74], and observation-theory intensity ratios as large as 10 have been reported for several lines in Fe XVI. It is believed that some of the differences could be explained by the blending of the lines with emissions from other ions [75]. Moreover, experimental measurements using ion traps (EBIT) have revealed significant differences with theoretical energies (see e.g., [74,76,77]). This situation has improved over the last decade. The high computational power reachable nowadays allows calculating energies and collision strengths for groups of ions with the same iso-sequence (e.g., [78][79][80][81]), or the same type of transitions for all ions [82]. As different codes calculate new sets of data as well as fit spectra (e.g., FAC [83], GRASP [84], AUTOSTRUCTURE [85]) and new data are stored in databases such as CHIANTI [69,70]) and HULLAC [73], there is better general agreement between theory and experiments (see for example [86]). However, there are still lines which remain unidentified, while others are poorly described because they are blended or have wrongly assigned transitions [87].
To illustrate this, two possible cases of disagreement between theory and experiment, both taken from CHIANTI (v7, [88]), are plotted in Figure 3. In the upper panel, the theoretical energies for Fe 16+ are, at most, roughly 1% off the experimental values. However, theory predicts many more levels that have not yet been observed. In the lower panel, the few levels predicted for Mg + are systematically lower than experimental values. Although the error for Mg + levels is around 3%, in the worst case, a line to the ground level can be shifted by as much as 25 Å.  [89]; observed energy levels (n = 6,7) from the wavelength measurements of [90]. (Bottom) Mg + : theoretical energies: [78]; observed energies from NIST: [91].

Radiative Transfer of UV to Infrared
The reprocessing of UV light into infrared via dust grains is crucial for the energy balance in the ISM and thereby the line emission. In addition, it allows for a determination of the dust continuum emission which is as important an observable as line emission. However, radiative transport (RT) is a highly complex, non-local, and non-linear problem. It is considered one of the four grand challenges in computational astrophysics. 6 Considerable progress has been achieved over the last 20 years in terms of carrying out exact 3D radiative transfer simulations of dust, as outlined in [92]. Yet, modeling dust radiative transfer in large hydrodynamical simulations is very computationally expensive, and it is usually ignored. A few pioneering studies that account "on the fly" for radiative transfer and the effects of dust are (e.g., [24,[93][94][95][96][97][98]). RAMSES-RT [99,100] implements radiative transfer in a cosmological code, albeit without a full treatment for dust emission. Examples of other codes that transfer radiation on the fly include EMMA [101] in a cosmological setting, as well as ORION [102] and Athena [103] in an ISM setting. Generally speaking, the impact of dust on the energy balance of simulated galaxies is typically not accounted for when computed on the fly. Alternative approaches to modeling dust are therefore worth exploring.
Ideally, a full simulation of the hydrodynamical evolution of astrophysical clouds, including galaxies, needs to account for the propagation of radiation through the simulated volume, and for its interaction with matter. The latter requires detailed knowledge of the quantum structure of all ions, molecules, and dust present in the simulated gas, as well as keeping track of the chemical reactions that couple the various species. For a simulation with a few million cells or particles, the problem becomes extremely computationally expensive, and the non-local couplings induced by radiation render it practically intractable with the computational means available at present.
An alternate solution is to solve for the hydrodynamics and the microphysics separately, in post-processing. In this approach, a code such as those listed in Table 2 is required to simulate the microphysics. CLOUDY [13] performs 1D spectral synthesis simulations of astrophysical clouds, accounting for all ions of elements up to Zn, a large number of molecular species, and all important microphysical processes. The propagation of radiation is done with the escape probability formalism, which is an expedient, and generally accurate approximation [104]. However, a few pathological cases exist for which escape probabilities are known to not give correct results [105]. Exact radiative transfer methods are required if a code is to be used in conjunction with hydrodynamical codes. The CLOUDY team will seek to implement exact radiative transfer, eventually with full treatment of the interactions of dust and radiation. Even with these improvements, however, due to their 1D nature these simulations will not capture the time-varying three-dimensional radiative fields that exist in nature. An exact 3D treatment of radiative transfer will be needed to connect the microphysics and 3D hydrodynamical models.

The Internal Density and Velocity Profile of Molecular Clouds
In calculating emission line strengths from molecular clouds, the adopted radial profiles of densities and radial velocity profile of the gas can drastically change the result. For example, at the scale of molecular clouds and their substructure, the so-called "infall profiles" are routinely used in identifying collapsing dense cores (see, e.g., [106]). These profiles consist of self-absorbed lines with a blue excess in moderately optically thick lines, in conjunction with single-peaked profiles of similar total width in optically thin lines. These profiles have been interpreted in terms of a dense molecular cloud core undergoing "inside-out collapse" [107], where the inner regions of the core are undergoing collapse, while the external envelope is in hydrostatic equilibrium, the two regions mediated an expanding rarefaction front. The assumed underlying radial velocity profile is essential for the generation of the lines. Recent studies (e.g., [108] Loughnane et al. 2018, submitted) have focused on the effect different assumed radial density and velocity profiles have on the emitted lines, and favored the collapse of marginally-unstable structures, characterized by "outside-in" velocity profiles rather than the canonical inside-out one. Thus, traditional line modeling has to assume a velocity profile, but there may be a degeneracy such that different combinations of density and velocity profiles may produce similar spectral lines. Synthetic observations of numerical simulations are thus essential in order to break the degeneracy, since they produce dynamically self-consistent density and velocity profiles, restricting the range of physically acceptable choices and thus facilitating the identification of the actual physical conditions that produce the signature spectral lines.

Simulating the Ionizing UV Field that Clouds Are Embedded in
The standard way to model the ionizing UV radiation, which ultimately determines the chemistry of the ISM, is to use stellar population synthesis (SPS) codes. In the SPS methodology, the integrated spectrum from populations of stars with given characteristics (IMF, star formation history, metallicity, etc.) is derived from stellar evolutionary tracks, which give the H-R diagram positions of stars of given masses with time, and from individual stellar spectra, which can be empirical or theoretical. If theoretical, they are computed with stellar atmosphere codes. Figure 4 shows examples of model predictions for populations of stars resulting from an instantaneous burst of star formation, which were computed with starburst99 7 .
However, using SPS codes comes with a few cautions worth knowing. For example, SPS codes from independent research groups account for different stellar astrophysics (rotation, close-binarity, etc., see e.g., [109]) and predict different ionizing spectra (e.g., Figure 2 of [110]). In addition, the ionizing continuum from individual ionizing stars cannot be observed directly, because it is strongly absorbed in the ISM. This makes constraining model predictions of the ionizing spectrum of populations of stars difficult. The escape fractions of ionizing photons from a few distant galaxies have been estimated from direct detections of ionizing photons (e.g., [111]). However, these are insufficient to constrain the ionizing spectrum from massive star populations, because they are just escape fractions. Thus, we must rely on theory. However, there remain large uncertainties in massive star evolution (trajectories in H-R diagram and stellar atmospheres), particularly for post-MS phases and stars more massive than 20 M .

Direct Observations of Low-Metallicity Massive Stars
Massive (M > 8 M ) extremely metal-poor stars (Z/Z ≤ 0.1) are expected to be present in extremely metal-poor galaxies, i.e., galaxies with ionized-gas oxygen abundances of ≤solar/10. Since such galaxies are farther than ∼10 Mpc away, they are too distant to enable observations of their individual massive stars. This is a problem for understanding a variety of astrophysical objects which range from the first stars and galaxies which re-ionized the universe to the metal-poor stars that give rise to the formation of heavy binary black holes. Direct observations of individual, massive, extremely metal-poor stars will only be possible if telescopes such as LUVOIR 8 and HabEx 9 , which are being proposed to NASA, ever see the light.  [113] were used. The shaded area indicates the region in wavelength used to integrate the SED for obtaining the FUV flux in e.g., Habing units [114], setting the chemical state of the gas. We note that comparing these models of different metallicities at the same age can be misleading, because stars of different metallicities have different lifetimes, i.e., a star of 3 Myr with solar metallicity will not be in the same evolutionary stage as a star of 3 Myr with 0.1 times solar metallicity.

Implementing Turbulence and Shocks in Simulations of the ISM
In addition to the importance of modeling stars as described in the previous two subsections, another important source of heating and ionization is turbulence and shocks. The introduction of density gradients and additional sources of ionization through these processes can have a significant impact on the observed spectrum, and hence the emission line diagnostics. Examples of how different levels of turbulence affect the density structure of ISM gas is shown in Figure 5. One solution to understanding the effects of turbulence is through the publicly available MAIHEM 10 code. This FLASH (version 4.4) based code tracks ionizations, recombinations, and species by species radiative cooling for hydrogen, helium, carbon, nitrogen, oxygen, neon, sulfur, calcium, and iron [32,67,68]. By including a solenoidal stirring mechanism, turbulence is introduced into the gas that then propagates through all spatial scales.
At the meeting, there was an interest in applying MAIHEM to Baldwin-Phillips-Terlevich (BPT) diagrams [1,115] and other UV diagnostics, as shown in Figure 6, and other UV diagnostics. Specifically, turbulence in the gas might help explain the observed dispersion in some spatially-resolved BPT diagrams of nearby AGN [116],  In addition, the treatment of shocks is particularly complex as it requires understanding the physical conditions in the emission line gas before, during, and after the shock event. These events can cause collisional ionization and the resulting recombinations may release additional ionizing photons that interact with surrounding material, depending on the shock velocity. A new database of shock models calculated using the MAPPINGS V 11 code [117] was presented at the workshop, and details can be found by consulting the documentation available on the Mexican Million Models database (3MdB) website 12 .

From Cloud to Galaxy Scale Simulations
Cloud scale simulations starting at masses of 10 4 -10 6 M have shown us how gas structures emerge and evolve for different inflow rates from converging streams of diffuse gas, corresponding to large-scale velocity dispersions in galaxy simulations (e.g., [118][119][120][121][122][123][124]). However, these simulations are mostly made at CNM initial conditions, i.e., in a limited temperature range. Projects dealing with line emission on galaxy scales would benefit greatly from seeing the impact on small-scale simulations from different temperatures reflecting different radiation backgrounds and metallicities. The community could also benefit from more ISM simulation work on scales in between cloud and galaxy sizes (e.g., [125]) and also observations of individual molecular clouds in external galaxies (e.g., PHANGS survey 13 and the LMT effort 14 ).
When simulating line emission from objects extracted from analytical/cosmological simulations, it is often a challenge to identify the minimum amount of information you need in order to get a decent estimate of the emission. For example, some simulations calculate abundances of specific elements on the fly, whereas others only contain one number for the metallicity and an abundance table must be adopted. Including several elements separately instead of a fixed abundance table, has been shown to make a significant difference for modeling of the fine-structure [C II] line [126].
Additionally, other possible uncertainties are present when considering galaxy simulations that include in their modeling both chemical networks (e.g., [127,128]) and radiative transfer (e.g., [129,130]); while such simulations can capture dynamical and thermodynamical effects due to photochemistry, they are not refined enough to calculate line emission, which is to be done in post-processing. However, in post-processing line calculations typically consider photoionization equilibrium and it is difficult to account for the past thermal and ionizational history of the simulated gas.
One of the more valuable conclusions from the discussions on galaxy-scales simulations, was the importance of simulating more than one emission line simultaneously. By simulating different lines, arising in different ISM phases, and comparing with observations, one ensures that the post-process recipes not only satisfy what is seen in one ISM phase, but is consistent across the entire galaxy.

Mapping Simulated Galaxy Samples to Observed Samples
An important step in simulating line emission from galaxies and from interstellar clouds is to compare with observed analogues of the model systems. Here, a number of issues continue to stand out, of which we discussed the following: 1 Observed SFRs come from tracers such as Hα, UV, IR and radio sensitive to a time averaged SFR of 10 to a few hundred Myr. Model SFRs are often the instantaneous SFRs. The question here is whether models should use time averaged SFRs or produce continuum and line emission to measure SFRs like observers do? One of the better solutions is to do both and identify any possible biases that are introduced by using observational methods. 2 Model metallicities are often M metal /M gas|star , whereas observed metallicities are determined using a number of ionized emission lines and presented as e.g., 12 + log(O/H). As with SFR, the best option at the moment might be to directly calculate 12 + log(O/H) if the oxygen abundance is tracked by the simulation used, although even this approach has problems as the observed 12 + log(O/H) will ultimately be luminosity-weighted which is hard to replicate in a simulation. 3 A realistic comparison between modeled and observed UV/optical emission lines (and continuum) requires the correct treatment of dust absorption within simulations. The amount of intervening dust in galaxy scale simulations between the photon source and the observer is ideally calculated self-consistently using dust chemical networks (e.g., [131][132][133][134]). However, often the amount of intervening dust is scaled as a function of the gas-phase metallicity by assuming a fixed grain size distribution. Finally, the geometry of the galaxy and the relative star-to-dust location is critical in correctly determining the observed properties (e.g., [135,136]), including sub-grid modeling of the dust properties of the birth-clouds of young stars.

4
Atomic hydrogen emission at 21cm is out of reach of current instrumentation beyond at z > 0.3 (SKA and SKA pathfinders such as APERTIF, MEERKAT, and ASKAP will be a remedy for this). This hampers the validation of the models that predict an evolution of atomic to molecular gas fraction with molecular gas dominating the gas budget within the optical scale of the galaxies at higher redshifts. 5 Analyzing a large sample of strongly star-forming SDSS galaxies with nebular He II λ4686 emission has shown that the luminosity of this line can only be reproduced with single bursts of star formation of 20 percent solar metallicity or higher, and ages of 4-5 Myr, when the extreme UV continuum is dominated by Wolf-Rayet stars [137]. He II λ1640 in the UV is 10 times stronger than the optical He II λ4686 line, and has been observed in broad (FWHM 1000 km s −1 ) and narrow emission for tens of dwarf star-forming galaxies also selected from SDSS. Reproducing the luminosity of the narrow nebular He II λ1640 emission from these galaxies has been challenging, even with state-of-the-art spectral synthesis models which combine the newest Charlot & Bruzual population synthesis models, which include very massive (300 M ) stars, with photoionization models, as described in [138]. This failure is reported for instance in [139]. In this workshop, Aida Wofford presented the case of one of the most metal poor nearby galaxies known, SBS 0335-052E, which has a metallicity of solar/20. None of the models which they tried were able to reproduce the luminosity of the He II λ1640 line. This is a problem because this line will be one of the only diagnostics of massive stars in future observations with large telescopes such as JWST (see a recent technical and scientific description in [140]), TMT [141], and e-ELT 15 . These telescopes will obtain rest-frame UV spectra for thousands of galaxies, at redshifts between 10 and 15, in the era of re-ionization when the first stars and galaxies formed. 6 How do AGN affect the comparison between observed and simulated galaxies? Comparing galaxies of similar SFR can become problematic, as the ionizing radiation and presence of a radio jet will enhance Hα and radio emission that is typically used as a SFR diagnostic. In addition, the radiation will heat dust that is often used in mass estimations. See the following section for more on the effect of AGN and mass outflows. Although not directly discussed at this workshop, X-ray dominated regions (XDRs) are another important result of having an AGN present, and they must be modeled in order to reproduce certain emission lines. For example, this is extremely relevant for high-J CO lines, and it is in fact still not clear if high-J CO lines are influenced more by the presence of X-ray Photons or shocks [142]. Important theoretical work on modeling XDRs was published in 2005 [143] and CLOUDY has also been used to model XDRs (e.g., [144]).

The Effects of AGN and Mass Outflows on Line Emission
The presence of an AGN can profoundly affect every aspect of observed emission line profiles including their velocities, widths, symmetry, and relative flux ratios. The AGN emits radiation anisotropically, and under the proper physical conditions can effectively couple with the ISM. This gives rise to outflows [145] of molecular and ionized gas that are likely radiatively driven. Understanding these outflows is critical as they may deliver feedback to the galaxy by clearing the inner regions of star forming gas, and assist in establishing the observed scaling relationships between galaxies and their supermassive black holes (SMBHs) [146][147][148][149].
At the workshop a subset of topics were covered, highlighting the physical mechanisms that couple the ISM to the AGN radiation field and affect the observed emission line properties. As an example, recent results from post-processing high-resolution, cosmological zoom-in simulations of massive galaxies with nebular emission models of galaxies and AGN were presented at the workshop [36,150]. One of the main results is that at least for massive galaxies, AGN-driven outflows seem to be a necessary component in order to reproduce the observed cosmic evolution of [OIII]/Hβ at fixed [NII]/Hα via strongly regulating the star formation history, which controls nebular emission from young stars via the ionization parameter (see Figure 7). There was also a particular focus on what can be gleaned from spatially resolved observations and models. In the local universe, spatially resolved BPT diagrams have proven their value in understanding how the physical conditions of the emitting gas change with distance from the SMBH, and can be used to constrain CLOUDY photoionization models of ionized mass outflows [116]. A first theoretical modeling of such spatially resolved BPT diagrams, using cosmological zoom-in simulations of massive galaxies shows that central low ionization emission in massive galaxies is not necessarily caused by an AGN, but can be also due to the ionizing radiation coming from post-AGB stars (Hirschmann et al., in prep.). For some galaxies, direct signatures for AGN-driven outflows are also visible in nebular emission line maps. These studies have primarily focused on the optical nebular emission line gas; however, observations reveal that powerful outflows of molecular gas may equal or even dominate the feedback in some galaxies. A conundrum in this area has been understanding how molecular outflows remain stable over time, as the physical conditions in the gas suggest many molecules should be destroyed. New modeling with advanced chemistry networks has found that molecules can form in-situ within the outflow, providing a significant opacity source that helps to drive the molecular outflows to the observed velocities [151]. Additional work in this area is needed to fully understand how the presence of an AGN modifies the standard emission line diagnostics for determining properties such as dust temperature, metallicity, and the inferred star formation rates.
A special case combining the effect of AGN, going from the host galaxy to the surrounding CGM, is the high-redshift radio galaxy MRC 0943-242. S. Kolwa presented nuclear emission and CGM absorption observations of MRC 0943-242 (see Figure 8), whose emission is spatially integrated over a 0.6 arcsec (4.75 kpc) radius aperture centered at what is assumed to be an AGN. Hence, the lines reveal gas in the ISM in close proximity to the nucleus of the galaxy and ionized by the AGN radiation. The emitting gas is turbulent with a velocity width larger than 1000 km s −1 . Models and simulations are important in understanding the physical mechanisms behind the turbulence in the ionized gas caused by processes surrounding the AGN. Further, Figure 9 shows absorption from the gas of various column densities in the CGM. The kinematics of this gas will be studied soon with IFU observations. What impact the AGN has on heating or removing the galaxy gas may be uncovered by combining observations of the ISM emission excited by the AGN with emission from CGM absorbers possibly resulting from AGN outflows. Further, these observations may discriminate between different models of the origin and fate of CGM gas.

Simulating Line Absorption from the CGM
While the workshop focus was not specifically on absorption line profiles, or the CGM, several talks and break-out discussions emphasized the connection between modeling emission lines and absorption lines, and the value in collaborations between these two fields. In particular, the CGM can give further constraints on the nature of the galaxy that could be important for converting emission line observations to other characteristics of the galaxy, such as gas mass.
Simulations of the CGM and its corresponding absorption line profiles are progressing, even given the difficulty of resolving this medium. Observations of distant quasars with absorption from the CGM of foreground galaxies provide useful metrics to help improve simulations. The characteristics of these observations are described well in [152] and references therein. In particular, observations find an abundance of OVI in the CGM of galaxies of a range of masses, only well matched in simulations at high galaxy masses above 10 10.5 M [153]. Other simulations post-process their outputs with Trident [154] to produce mock absorption line profiles and compare them with quasar absorption line observations. The predicted CGM column densities and covering fractions of ions such as OVI depend on the models of star formation and stellar feedback (see Figure 10; Chuniaud et al. in prep). Star formation models influence the amount of metals produced in the simulations, while feedback models impact their ability to leave the galaxy and persist in the CGM. Figure 10. CGM OVI covering fractions in a sub-L * simulated galaxy at z = 3 for different subgrid models. SF#FB# corresponds to whether the simulations uses star formation model 1 or 2 (SF1 and SF2), and stellar feedback model 1 or 2 (FB1 and FB2). SF1 uses a density threshold for triggering star formation. FB1 corresponds to feedback from a stellar population with a Salpeter IMF and the massive stars exploding with 10 51 erg per 10 M and the energy injected in the simulation matching the Sedov solution [155]. Both of these subgrid models depend on the resolution of the simulations (here 10 pc). SF2 uses a variable star formation efficiency for gas that is not supported against collapse by turbulence, following [156,157]. FB2 is consistent with a Chabrier IMF, and either the energy or momentum solution is injected depending on which is resolved by the simulations [158].
Observations also find that these OVI lines typically have suprathermal line widths [159]. However, this could be due to blending of multiple components along the line of sight. Simulations could reveal the nature of these broad components, but this may depend on resolution. The dependence of the simulated CGM on the numerical resolution is a debated topic. It is unclear what are the scales of the cool gas in the CGM, which could range from tens of pc down to sub-pc scale (e.g., [160]). Increasing the resolution in the CGM leads to better resolving the central density peaks of outflowing and inflowing material in the CGM. The initially slightly higher density has faster cooling rates, which is amplified over time. Thus, simulations with better resolution in the CGM, such as in NEPHTHYS presented by M. L. A. Richardson, form more gas clumps, mostly in inflowing gas. Inflows are enriched from outflows, leading to enhanced cooling. However even the high resolution CGM simulations shown by Richardson et al. (NEPHTHYS) which have 100 physical pc resolution, may by numerically supported. Observed clumps in their simulations are roughly 10-20 cells across.
Many models of the CGM assume the ionizing radiation source to be the background UV reionization field, but the central galaxy can dominate the background over periods of time, including during AGN flaring or starburst episodes (e.g., [153,161]). Non-equilibrium treatment of the ionization state could make this effect be longer-lasting. Future work, including the production of look-up tables to convert simulation data to predicted ionization states, will need to consider this.
Finally, the CGM will be better constrained in the future with the ELT and other next-generation observatories that will be able to use Lyman-break galaxies, far more numerous than quasars, to have multiple line-of-sight absorption profiles per foreground galaxy. This will allow for spatially resolved study of the CGM, and give more instruction for combining different foreground galaxy column density to make general CGM column density profiles.

Discussion: How Can We as a Community Move Forward?
This workshop was the first of its kind (to our knowledge) to attempt to create a community of astronomers simulating line emission and working on the comparison with observations. In order to continue the sharing of knowledge within the community, we discussed ways in which to make the material of the workshop available to a larger group and encourage collaborations in the future. For this workshop, it was decided to use zenodo.org where all speakers had the option to upload their slides together with a recording of their talk (where available) to a community space 16 . This approach automatically provides a DOI for each upload for easy reference and we hope that future workshops of this kind will likewise attempt to make the presentations publicly available.
It should also be mentioned that the specific problems of simulating line emission from plasma are not only faced by the astronomical community. Within the field of plasma physics, collisional-radiative models are used to analyze plasma kinetics and spectra as well as to design plasma experiments in the laboratory. Yet, the two fields, i.e., that of astronomy and that of plasma physics, are currently far from fully benefiting from each other. We therefore wish to highlight corresponding established workshop series for non-LTE modeling, namely the "NLTE Code Comparison" workshops 17 and the "Spectral Line Shapes in Plasmas" workshops 18 .

Getting Involved With the Community of Developers
As Table 2 shows, there are several software tools available for simulating line emission, each with their own group of developers behind. However, in order to achieve the best usage of these tools and adapting them to astrophysically relevant scenarios, it is crucial that users are able to connect with the developers to request support and updates, or even to adjust the underlying code themselves. Here we detail how this is done for two examples listed in Table 2.

The CLOUDY Community
Small tutorials on using CLOUDY have been uploaded to youtube by individuals not part of the CLOUDY developer team and such video uploads are much appreciated by many and strongly encouraged by the CLOUDY community. It was suggested to establish a CLOUDY-branded youtube channel or youtube playlist for these videos. Also, many people are not aware that you can edit the databases used by CLOUDY to for example include new species. The CLOUDY quickstart guide has been noted to be an excellent tool for introducing new researchers and students to making predictions with the software. Maintaining it and keeping it up-to-date will be a priority. The CLOUDY team is always looking for more participation on the Yahoo Groups forum. A goal is that questions and discussions will take place there among users, rather than solely between a user and cloudy developers. An interactive wiki has also been suggested as a place for users to post about specific cases or issues in their CLOUDY research.

LIME
LIME resides on the code hosting platform GitHub (https://github.com/lime-rt/lime) from where it can be downloaded. A small team of active developers maintain and develop the package and minor releases appear regularly, approximately once per year. Users that experience problems when using LIME or have ideas to improve on the code, are welcome to create issues on GitHub where they will be answered by one or more developers. It is also possible to fork the code and develop it for specific needs. Documentation exists in the form of a written PDF manual as well as an issue tracker. Both are available on GitHub. A bi-annually summer school on Monte Carlo techniques in radiative transfer problems, aimed at PhD-students, is held at University of St Andrews School of Physics & Astronomy and this school offers both a lecture and tutorials on how to use LIME.

Conclusions
This paper reviews recent progress and identifies current problems in the field of simulating line emission from the ISM and CGM. Conclusions presented here are a result of discussion sessions at the workshop "Walking the Line" held in Phoenix Arizona on 14-16 March 2018. This international meeting was the first of its kind to gather astronomers specializing in the modeling of line emission and brought together 30 scientists from grad student to professor level. Most talks with follow-up questioning sessions have been made available online, but as a more digestible release of knowledge, we have condensed the outcome on the most relevant topics in this paper. Returning to the questions posed in the introduction, below are the main conclusions, focusing on the methods used to give the most reliable answers to each question.
What are the best emission lines to trace various ISM properties and ionized sources in galaxies? The best method to answer this question of course depends on the type and scale of the emitting source in the ISM considered. In order to provide the reader with an overview, though not complete, of tools used to calculate line emission either on-the-fly in a simulation or in post process, we listed the software tools that were discussed at the workshop. One outstanding issue that can affect these tools, mostly when considering line emission in the X-ray but also at longer wavelengths [162], is to correct theoretical heavy-element energy levels that do not agree with experimental values.
Key to deriving the line emission from any cloud is the UV field irradiating that cloud. However, our knowledge on the ionizing continuum of stars is severely hampered by the observational difficulties and thus currently restricted to theoretical predictions. Direct observations of individual, massive, low-metallicity stars will only be possible with telescopes such as LUVOIR and HabEx. Once the UV field is known, radiative transfer of UV into infrared light is a crucial part of any method used, especially when studying neutral or molecular regions. While still too computationally expensive to run alongside galaxy-scale simulations, work is ongoing to include it more regularly in post-process as part of any line emission simulation.
In order to compare simulated galaxies with observed ones, we need to put more thought into the choice of physical properties being compared. For example, a SFR that is modeled as instantaneous can hardly be compared to an observed SFR that is essentially derived as an average over 10 to a few-hundred Myr. Many other parameters such as the effects of an AGN or the inclination of a galaxy can distort an attempt to compare model galaxies with observed ones. Finally, it was concluded that the aim should always be to simulate more than one emission line simultaneously, in order to assure consistency across an the entire galaxy.
How can we use emission lines to trace feedback and ISM evolution with redshift? Feedback from AGN in the form of radiation and outflows has a profound effect on the nebular line emission. As an example, simulations of massive galaxies with and without AGN were presented in which the presence of an AGN could clearly be identified from the [OIII]/Hβ and [NII]/Hα emission-line ratios. AGN-driven outflows thus strongly regulates the star formation history, at least for massive galaxies, which controls nebular emission from young stars via the ionization parameter. As an observational example, the high-redshift galaxy MRC 0943-242 is believed to contain an AGN, with line emission (and absorption) revealing ionized ISM in close proximity to the nucleus as well as outflows and high levels of turbulence.
In close relation to stellar and AGN feedback, turbulence and shocks are another important source of heating and ionization. Recently developed tools for the treatment of turbulence and shocks were presented to aid the line emission simulations in cases where necessitated.
On cloud-scale level, simulating the line emission from collapsing dense cores can reveal whether the molecular cloud is undergoing an "inside-out" or "outside-in" collapse. Work is undergoing to break this degeneracy by making synthetic observations of signature spectral lines for comparison with observations.
How should absorption features be correctly interpreted? Simulations of line absorption through the CGM are faced by its own set of problem, where for example the treatment of the ionizing field from the galaxy source can play a crucial role but has yet to be fully treated in a non-equilibrium manner. Another choice to be made when simulating the CGM is the level of spatial resolution needed, as better resolution leads generally to more gas clumps and a point of convergence must be found.
Where do we stand in deriving sub-grid physics and comparing codes? When studying ISM evolution with redshift, line emission must be derived from cosmological simulations which capture the evolution but comes at low spatial resolution compared to semi-analytical models. In addition, full chemistry and radiative transfer must typically be applied in post process. The uncertainties in the approach could be alleviated by more work on scales between galaxies and clouds, because simulations of the latter come at higher resolution and could hence be used as a benchmark on larger scales.
How do we coordinate our efforts? We hope to have set an example for more theoretical workshops of this kind to come, and look forward to seeing the results of the ongoing projects presented here. The developers of the software tools used to derive line emission, increasingly encourage their users to report problems or suggest improvements to their code via online platforms. With the list of caveats identified in this workshop, albeit long, we continue to believe that the future of simulating line emission is bright, especially when research groups in the field come together to find and solve common problems.