We propose a new method to search for axion-like particles (ALPs) based on the gamma-rays produced concomitant with high-energy astrophysical neutrinos. The existence of high-energy neutrinos implies production of gamma-rays in the same sources. Photons can convert into ALPs in the sources' magnetic fields, and will travel as ALPs through extragalactic space. Back-conversion in the Milky Way's magnetic field leads to a diffuse anisotropic high-energy photon flux that existing and upcoming gamma-ray detectors, like HAWC, CTA, and LHAASO can detect. This method probes unexplored ALP parameter space, with LHAASO being realistically sensitive to couplings above $10^{-11}\, \rm{GeV^{-1}}$ and masses up to $3\times 10^{-6} \, \rm{eV}$ in ten years. Our technique also explores viable ALP dark matter parameter space.

Magnetic fields in extragalactic space between galaxy clusters may induce conversions between photons and axion-like particles (ALPs), thereby shielding the photons from absorption on the extragalactic background light. For TeV gamma rays, the oscillation length ($l_{\rm osc}$) of the photon-ALP system becomes inevitably of the same order as the coherence length of the magnetic field ($l$) and the length over which the field changes significantly (transition length $l_{\rm t}$) due to refraction on background photons. We derive exact statistical evolution equations for the mean and variance of the photon and ALP transfer functions in the non-adiabatic regime ($l_{\rm osc} \sim l \gg l_{\rm t}$). We also make analytical predictions for the transfer functions in the quasi-adiabatic regime ($l_{\rm osc} \ll l, l_{\rm t}$). Our results are important in light of the upcoming Cherenkov Telescope Array (CTA), and may also be applied to models with non-zero ALP masses.

We study the impact on the Sun of an exotic energy-loss channel caused by plasmon decay into fermionic minicharged particles with charge $\epsilon e$ and mass $m_f$. We compare solar models with this extra emission to helioseismological and neutrino data, obtaining a bound $\epsilon <2.2\times 10^{-14}\ ({\rm 95\%\, CL})$ for $m_f\lesssim 25 \, \text{eV}$. Our result is comparable to previous limits from the cooling of globular cluster stars, while at the same time it is better understood and takes theoretical and observational errors into account.