We use the gas-grain chemistry code UCLCHEM to explore the impact of cosmic-ray feedback on the chemistry of circumstellar disks. We model the attenuation and energy losses of the cosmic-rays as they propagate outwards from the star and also consider ionization due to stellar radiation and radionuclides. For accretion rates typical of young stars, $\dot M_* \sim 10^{-9}-10^{-6}$ M_\odot yr$^{-1}$, we show that cosmic rays accelerated by the stellar accretion shock produce a cosmic-ray ionization rate at the disk surface $\zeta \gtrsim 10^{-15}$ s$^{-1}$, at least an order of magnitude higher than the ionization rate associated with the Galactic cosmic-ray background. The incident cosmic-ray flux enhances the disk ionization at intermediate to high surface densities ($\Sigma > 10$ g cm$^{-2}$) particularly within 10 au of the star. We find the dominant ions are C$^+$, S$^+$ and Mg$^+$ in the disk surface layers, while the H$_3^+$ ion dominates at surface densities above 1.0 g cm$^{-2}$. We predict the radii and column densities at which the magneto-rotational instability (MRI) is active in T Tauri disks and show that ionization by cosmic-ray feedback extends the MRI-active region towards the disk mid-plane. However, the MRI is only active at the mid-plane of a minimum mass solar nebula disk if cosmic-rays propagate diffusively ($\zeta \propto r^{-1}$) away from the star. The relationship between accretion, which accelerates cosmic rays, the dense accretion columns, which attenuate cosmic rays, and the MRI, which facilitates accretion, create a cosmic-ray feedback loop that mediates accretion and may produce luminosity variability.

We utilize a modified astrochemistry code which includes cosmic ray attenuation in-situ to quantify the impact of different cosmic-ray models on the CO-to-H$_2$ and CI-to-H$_2$ conversion factors, $X_{\rm CO}$ and $X_{\rm CI}$, respectively. We consider the impact of cosmic rays accelerated by accretion shocks, and show that clouds with star formation efficiencies greater than 2\% have $X_{\rm CO} = (2.5 \pm 1)\times10^{20}$ cm$^{-2}$(K km s$^{-1}$)$^{-1}$, consistent with Milky Way observations. We find that changing the cosmic ray ionization rate from external sources from the canonical $\zeta \approx 10^{-17}$ to $\zeta \approx 10^{-16}$ s$^{-1}$, which better represents observations in diffuse gas, reduces $X_{\rm CO}$ by 0.2 dex for clusters with surface densities below 3 g cm$^{-2}$. We show that embedded sources regulate $X_{\rm CO}$ and decrease its variance across a wide range of surface densities and star formation efficiencies. Our models reproduce the trends of a decreased $X_{\rm CO}$ in extreme cosmic ray environments. $X_{\rm CI}$ has been proposed as an alternative to $X_{\rm CO}$ due to its brightness at high redshifts. The inclusion of internal cosmic ray sources leads to 1.2 dex dispersion in $X_{\rm CI}$ ranging from $2\times10^{20} < X_{\rm CI} < 4\times10^{21}$ cm$^{-2}$ (K km s$^{-1}$)$^{-1}$. We show that $X_{\rm CI}$ is highly sensitive to the underlying cosmic ray model.

We present astrochemical photo-dissociation region models in which cosmic ray attenuation has been fully coupled to the chemical evolution of the gas. We model the astrochemical impact of cosmic rays, including those accelerated by protostellar accretion shocks, on molecular clouds hosting protoclusters. Our models with embedded protostars reproduce observed ionization rates. We study the imprint of cosmic ray attenuation on ions for models with different surface cosmic ray spectra and different star formation efficiencies. We find that abundances, particularly ions, are sensitive to the treatment of cosmic rays. We show the column densities of ions are under predicted by the `classic' treatment of cosmic rays by an order of magnitude. We also test two common chemistry approximations used to infer ionization rates. We conclude that the approximation based on the H$_3^+$ abundance under predicts the ionization rate except in regions where the cosmic rays dominate the chemistry. Our models suggest the chemistry in dense gas will be significantly impacted by the increased ionization rates, leading to a reduction in molecules such as NH$_3$ and causing H$_2$-rich gas to become [C II] bright.

We construct a model for cosmic ray acceleration from protostellar accretion shocks and calculate the resulting cosmic ray ionization rate within star-forming molecular clouds. We couple a protostar cluster model with an analytic accretion shock model to calculate the cosmic ray acceleration from protostellar surfaces. We present the cosmic ray flux spectrum from keV to GeV energies for a typical low-mass protostar. We find that at the shock surface the spectrum follows a power-law trend across 6 orders of magnitude in energy. After attenuation, the spectrum at high energies steepens, while at low energies it is relatively flat. We calculate the cosmic ray pressure and cosmic ray ionization rate from relativistic protons at the protostellar surface and at the edge of the core. We present the cosmic ray ionization rate for individual protostars as a function of their instantaneous mass and final mass. The protostellar cosmic ray ionization rate is $\zeta \approx 0.01 - 1$ s$^{-1}$ at the accretion shock surface. However, at the edge of the core, the cosmic ray ionization rate drops substantially to between $\zeta \approx 10^{-20}$ to $10^{-17}$ s$^{-1}$. There is a large spatial gradient in the cosmic ray ionization rate, such that inner regions may experience cosmic ray ionization rates larger than the often assumed fiducial rate, $\zeta = 3\times10^{-17}$ s$^{-1}$. Finally, we calculate the cosmic ray ionization rate for protostellar clusters over 5 orders of magnitude of cluster size. We find that clusters with more than approximately 200 protostars produce a higher cosmic ray ionization rate within their natal cloud than the fiducial galactic value.

We construct a semi-analytic model to study the effect of far-ultraviolet (FUV) radiation on gas chemistry from embedded protostars. We use the Protostellar Luminosity Function (PLF) formalism of Offner & McKee (2011) to calculate the total, FUV, and ionizing cluster luminosity for various protostellar accretion histories and cluster sizes. We compare the model predictions with surveys of Gould Belt star-forming regions and find the Tapered Turbulent Core model matches best the mean luminosities and the spread in the data. We combine the cluster model with the photo-dissociation region astrochemistry code, {\sc 3d-pdr}, to compute the impact of the FUV luminosity from embedded protostars on the CO to H$_2$ conversion factor, $X_{\rm CO}$, as a function of cluster size, gas mass and star formation efficiency. We find that $X_{\rm CO}$ has a weak dependence on the FUV radiation from embedded sources for large clusters due to high cloud optical depths. In smaller and more efficient clusters the embedded FUV increases X$_{\rm CO}$ to levels consistent with the average Milky Way values. The internal physical and chemical structure of the cloud are significantly altered, and $X_{\rm CO}$ depends strongly on the protostellar cluster mass for small efficient clouds.

We investigate the spectral correlations between different species used to observe molecular clouds. We use hydrodynamic simulations and a full chemical network to study the abundances of over 150 species in typical Milky Way molecular clouds. We perform synthetic observations in order to produce emission maps of a subset of these tracers. We study the effects of different lines of sight and spatial resolution on the emission distribution and perform a robust quantitative comparison of the species to each other. We use the Spectral Correlation Function (SCF), which quantifies the root mean squared difference between spectra separated by some length scale, to characterize the structure of the simulated cloud in position-position-velocity (PPV) space. We predict the observed SCF for a broad range of observational tracers, and thus, identify homologous species. In particular, we show that the pairs C and CO, C$^{+}$ and CN, NH$_3$ and H$_2$CS have very similar SCFs. We measure the SCF slope variation as a function of beam size for all species and demonstrate that the beam size has a distinct effect on different species emission. However, for beams of up to 10", placing the cloud at 1 kpc, the change is not large enough to move the SCF slopes into different regions of parameter space. The results from this study provide observational guidance for choosing the best tracer to probe various cloud length scales.