Anisotropic diffusion cannot explain TeV halo observations

TeV halos are regions of enhanced photon emissivity surrounding pulsars. While multiple sources have been discovered, a self-consistent explanation of their radial profile and spherically-symmetric morphology remains elusive due to the difficulty in confining high-energy electrons and positrons within ~20 pc regions of the interstellar medium. One proposed solution utilizes anisotropic diffusion to confine the electron population within a"tube"that is auspiciously oriented along the line of sight. In this work, we show that while such models may explain a unique source such as Geminga, the phase space of such solutions is very small and they are unable to simultaneously explain the size and approximate radial symmetry of the TeV halo population.

TeV halos are regions of enhanced photon emissivity surrounding pulsars. While multiple sources have been discovered, a self-consistent explanation of their radial profile and spherically-symmetric morphology remains elusive due to the difficulty in confining high-energy electrons and positrons within ∼20 pc regions of the interstellar medium. One proposed solution utilizes anisotropic diffusion to confine the electron population within a "tube" that is auspiciously oriented along the line of sight. In this work, we show that while such models may explain a unique source such as Geminga, the phase space of such solutions is very small and they are unable to simultaneously explain the size and approximate radial symmetry of the TeV halo population.
Our understanding of TeV halos hinges on one key question: are TeV halos produced in peculiar regions of the ISM that have pre-existing conditions ripe for halo formation? Or, conversely, are TeV halos produced throughout the bulk of the ISM and powered by the natal SNR, PWN, or potentially supergiant star which produces a local environment necessary for halo formation? In the former case, only a small fraction of pulsars will produce observable TeV halos. In the latter, TeV halos are expected to surround most energetic pulsars.
Observations support the latter case. Ref. [1] ranked all ATNF catalog [13] pulsars in terms of their "spin-down flux" (spin-down power divided by the pulsar distance squared). Assuming that all pulsars convert an equivalent fraction of their spin-down power into γ-ray emission, the TeV halo flux should be proportional to spin-down flux. Indeed, Ref. [1] found that five of the seven pulsars with the highest spindown-flux were detected as TeV sources, while none of the 48 dimmer objects in the HAWC field of view were detected. Subsequent observations detected TeV emission from another relatively high (the 11th brightest) pulsar [14]. Additional TeV halo observations by the HAWC, H.E.S.S. and LHAASO collaborations have provided additional support for the conclusion that TeV halos can be found throughout the Milky * pedro.delatorreluque@fysik.su.se, ORCID: 0000-0002-4150-2539 † ottavio.fornieri@gssi.it, ORCID: 0000-0002-6095-9876 ‡ linden@fysik.su.se, ORCID: 0000-0001-9888-0971 Way [8,10,15,16]. The theoretical arguments for "innate" or "sourceproduced" turbulence are less certain. Early studies focused on the potential that Geminga and Monogem exist in a "lowdiffusion pocket", that would potentially extend all the way to the solar position [17]. However, such a model is incompatible with both local observations, which indicate that the diffusion coefficient near Earth is not abnormally low [18], and global observations, because the existence of many large lowdiffusion regions surrounding pulsars would be incompatible with observed CR secondary-to-primary [4,6,19].
Several classes of models have been proposed to explain TeV halos. One popular model focuses on the potential for CRs accelerated by the pulsar or associated SNR [20] to excite a resonant streaming instability that self-confines the CRs near the source [21,22]. These models potentially explain the evolution of halos, but many complexities of CR turbulence must be solved to make precise predictions. Rectilinear propagation models argue that diffusion is not inhibited, but particle propagation is instead ballistic on small scales, which produces an effective suppression of high-angle emission [23]. However, such models may require an unphysically high efficiency for the pulsar e + e − production [24].
Finally, another kind of model argued that the apparent angular size of Geminga and Monogem are consequences of anisotropic diffusion with a maximal diffusion constant similar to the galactic average. In this scenario, the direction of efficient diffusion is oriented along the line-of-sight (LoS) towards Earth, while diffusion is strongly inhibited in the two visible directions perpendicular to the LoS [25] (see also Ref. [26]). This model is theoretically motivated by synchrotron polarization measurements which indicate that local diffusion is dominated by flux tubes on scales between 1-100 pc [27,28]. However, such a model does not predict that many TeV halos would be seen, as observable halos would only be expected from sources that have flux tubes that are fortuitously aligned towards Earth.
In this paper, we systematically re-examine the class of anisotropic diffusion models. We show that they cannot simultaneously account for the radial size and approximate spherical symmetry of the observed TeV halo population. We note that this conclusion holds for any CR-powered source (hadronic or leptonic), implying more generally that anisotropic diffusion does not dominate the propagation of particles near energetic sources.

A. Theory
To study the lepton distribution, u(r, t, E e ), around pulsars, we make use of the standard transport equation [29]: where D ij is the diffusion tensor, the energy derivative accounts for the energy losses and S(r, t, E e ) = Q(E)L(t)δ(r) is the source term, representing the flux injected by a point source located at r = 0 as a function of time.
The diffusion of charged particles depends on the local magnetic field, B tot . Magnetic fields in astrophysical plasmas can be described as B tot = B 0 + δB, namely the sum of a large-scale background field (B 0 ), with a coherence scale between ∼ 1 − 100 pc [30], and a small-scale field (δB) that depends on the size of the source that is powering turbulence. For TeV halos, typical turbulence scales are L inj ∼ O(10) pc, while SNRs may inject turbulence up to L inj ∼ O(100) pc.
To account for the effect of the magnetic field structure in particle transport, it is useful to decompose the diffusion tensor into directions parallel and perpendicular to the large-scale magnetic field lines as y, z, in a Cartesian reference frame. Placing the background magnetic field along the z-axis (B 0 = (0, 0, B 0 )) we exploit the axisymmetric nature of the problem and write D xx = D yy = D ⊥ , z |B 0 | 2 = D , and all D ij = 0. This allows us to solve Equation (1) in cylindrical coordinates (r, z, φ), for a cylinder oriented along the z axis, such that D zz = D and D rr = D ⊥ , where r is the polar coordinate x 2 + y 2 = r. The diffusion equation becomes: where, because of the cylindrical symmetry, the gradients involving the azimuthal coordinate φ vanish. The physics behind parallel and perpendicular diffusion on scales similar to the Larmor radius is different [31]. Parallel diffusion is the result of the scattering of particles against δB, while (mainly) the random walk of the lines themselves (fieldline random walk) is responsible for perpendicular diffusion. Therefore, if the injected turbulence is strong enough to considerably affect the preferential direction of the background field B 0 on small scales, particle motion tends not to have a privileged direction, and is instead isotropic. Conversely, a weak turbulence does not alter the direction of B 0 . The intensity of the injected turbulence is represented by the so-called Alfvénic Mach number, defined as M A ≈ δB B 0 Linj at the turbulence injection length-scale, L inj .
Anisotropic diffusion is an enticing explanation for the TeV halo morphology because it can explain the spatial extension of halos without invoking a diffusion coefficient that is orders of magnitude below standard ISM diffusion [2]. Whenever the background magnetic field direction is oriented with our LoS, we observe the TeV halo in only the directions where diffusion is inhibited, and the low-diffusion coefficient becomes a projection effect [25]. However, the overall diffusion coefficient, which is uninhibited along the LoS, remains consistent with global cosmic-ray measurements.
Quantitatively, we set D to match cosmic-ray measurements (e.g. the boron-over-carbon ratio) and set the perpendicular diffusion coefficient using the model D ⊥ = D M 4 Aderived by Ref. [32] for particle energies below ∼ 10 TeV, corresponding to Larmor radii smaller than the injection scale L inj ∼ O(10 − 100) pc. TeV halo observations constrain our models to 0 < M A ≤ 1, which spans from the anisotropic case (δB B 0 , or M A 0.1) to the isotropic one (δB = B 0 , or M A = 1). This implies that particle diffusion perpendicular to the local field can be strongly inhibited, depending on the turbulence strength and injection scale.
We parameterize the energy scaling of parallel diffusion as D = D 0 E E 0 δ , where D 0 is set at a chosen normalization energy E 0 and δ is derived from the spectral index of the turbulent power spectrum. We fix D 0 = 3.8×10 28 cm 2 s −1 at E 0 = 1 GeV and consider a Kolmogorov spectrum for which δ = 0.33. We note that these values are standard [33,34] and compatible with the first-principle calculations in Ref. [35]. Leptons in the halo interact with their environment to produce bright γ-ray emission, predominantly through inverse-Compton scattering (IC) of the surrounding Interstellar Radiation Field (ISRF) [17]. The IC emissivity results from the convolution of the photon field and the CR spectrum [36]: where dσ IC dE γ is the differential production cross-section of γ-rays with energy E γ resulting from the collision of an electron with energy E e and a background photon with energy E ph , dn ph dE ph is the spectral density of ISRF photons and Φ e is the differential electron flux, which, for isotropic emission, is Φ e (E e , r) = c 4π × u(E e , r). The differential cross section for IC is given by [36]: where σ T is the Thomson cross section, m e is the electron mass. Here p = 4E ph E e m 2 e and where the cross-section vanishes outside m 2 e 4E 2 e ≤ 1.

B. Numerical setup
Equation (2) cannot be solved analytically. In this paper, we numerically evolve our system using a Crank-Nicolson scheme [37] as detailed in Appendix B. We examine values of M A spanning from 0.1 to 1 and align the large-scale magnetic field with the z-axis. The diffusion equation is evolved on a 2D grid of radius 60 pc and height [−60 pc, +60 pc] with (N r , N z ) = (200, 200) points. Adopting a reference distance to Geminga of d Gem = 250 pc, this corresponds to a window ∆φ [−13 • , +13 • ] and angular resolution φ res 0.1 • .
We assume that the leptons are injected from a point-like pulsar source, noting that the assumed source size does not affect our results. Drawing on pulsar observations, we set the luminosity function to where L 0 is the luminosity of the source at t = 0, n is the braking index and τ 0 is the pulsar spin-down timescale. In our simulations we set L 0 = 2.8 × 10 37 erg s −1 , τ 0 = 12 kyr and n = 3, and normalize our results by imposing that the total energy released by the pulsar since its birth is W e = 1.1 × 10 49 erg, consistent with previous studies [21,25,38,39].
We convert this spin-down power into electron and positron pairs with an efficiency η, which is fit to data, but cannot exceed 1, yielding an injection spectrum Q(E e ) = η Q 0 E e GeV −α × e −Ee/Ecut , where α = 1.6 [19] and E cut = 200 TeV, and Q 0 is a normalization constant. We compute the electron flux from 0.1 − 300 TeV and the IC-produced γ-ray flux from 0.1 − 200 TeV. In our setup, we stop the simulation at the age of Geminga, t ch ∼ 342 kyr.
We point out that, by comparing the different relevant timescales, we can conclude that our final result is robust with respect to the age of the pulsar. There are three relevant timescales: (1) the age of the pulsar, (2) the spin-down timescale of the pulsar, which we compute by determining the time period over which the pulsar luminosity changes by a factor of e, and (3), the diffusion timescale, which determines the rate at which the leptons produced by the pulsar leave the simulation volume. The age of the pulsar in our simulation is 342 kyr. We note that the pulsar spindown timescale is a power-law, and not an expontential, so it changes significantly as a function of pulsar age. At 342 kyr, the effective pulsar spindown timescale is ∼180 kyr. The diffusion timescale, on the other hand, is only ∼2 kyr (for D 0 = 3.8 × 10 28 cm 2 /s). Over this timescale, the spin-down power only changes by ∼ 1%. This means that, while the spindown timescale affects the normalization of the particle density, it will have no effect on the morphology of the TeV halo in our simulation.
In fact, the axial ratio of diffusion parallel and perpendicular to the magnetic field lines is even more robust to changes in the pulsar age or spin-down timescales, because it is based on the ratio for particles to diffuse in each direction, which is independent of the instantaneous pulsar power. Thus, changes in the age of the modeled pulsar will produce negligible changes in the results of our study. On top of this, we also stress that the injection energy dependence assumed does not affect any of our conclusions on the morphology and radial profile of the γ-ray emission.
Synchrotron and IC energy losses are calculated as: where σ T 6.65 × 10 −25 cm 2 is the Thomson crosssection, and (U B , U i ) are the magnetic field and ISRF energy densities. We set U B to be equal to the dominant ordered field,  [41], an approximation that is valid for the energy-range and accuracy we require, but may fail for high-precision GeV measurements [42].

III. RESULTS
Using the modelling described above, we produce mock observations for Geminga-like TeV halos for various quantities of the parameters M A and the inclination angle of the simulation with respect to the LoS, ψ incl . We note that these models are produced in a 2D simulation with dynamics that are dependent on the parameter M A . The parameter ψ incl , on the other hand, corresponds to mock observations taken from different angles with respect to the axes of the simulation. We utilize the lower-case letters r and z when we discuss the physical coordinates of the simulation, and the upper-case letters R and Z to denote the projected coordinate system along our line of sight (see Figure A). We note that R = Z when ψ incl = 0 • and that R = r, Z = z when ψ incl = 90 • . Figure 1 (top) shows the morphology of the γ-ray emissivity as a function of r and z at 20 TeV, as computed in Equation 3. In Figure 1 (bottom), we show the observed extension of the simulated halo as a function of the angle ψ incl , which corresponds to rotations of our simulated cylinder with respect to our LoS along the r-axis, and compare our results to the 68% and 82% of the flux contained in the Geminga TeV halo as reported by Ref. [2]. We note that rotations around the zaxis do not change the morphology of the halo with respect to our LoS due to the cylindrical symmetry of the system, while rotations around the r-axis change the morphology that is projected on the plane-of-the-sky (c.f. Figure 2 of Ref. [25]). Figure 1 demonstrates that if anisotropic diffusion produces TeV halos, we should detect a variety of both highly extended and asymmetric objects (as seen at different inclination angles, ψ incl ). This is in tension with the fact that observed TeV 5 the halo appears roughly spherically symmetric, but the lack of inhibited diffusion makes the halo too large to explain observed systems. Notably, we see that for M A ≥ 0.5 there is no value of ψ incl for which the containment angle along the z-axis is consistent with HAWC observations of Geminga. We can formalize the excluded ψ incl angles based on the morphology and symmetry of simulated TeV halos by imposing two conditions: (i) that the emission should not be very asymmetric (i.e. the extension of the halo in any direction should not be much larger than the extension in the perpendicular one). (ii) the size of the emission should not be much larger than 5.5 • (i.e. 24 pc, given the distance from Geminga), to be consistent with the size of Geminga reported in Ref. [2], which corresponds to ∼ 82% of the flux contained around the source. We additionally calculate the size the halo at ∼ 68% (∼ 1σ) containment, which is 4.3 • (∼ 19 pc) for Geminga.
To quantify the first condition (hereafter, the symmetry condition) we impose that the projected extension of the halo in one direction must not be more than 100% larger than in the other direction (Z/R < 2), which is a very conservative choice. Figure 6 of the supplementary material (SM) shows the Z/R ratio for different inclination angles. The sec-ond condition (hereafter, the size condition, see bottom panels of Figure 1) imposes that the extension of the halo projected on the plane-of-sky along Z is within the size uncertainty reported by HAWC, which is 5.5±0.7 • (∼ 24±3 pc). While the first condition only depends on the ratio D ⊥ D = M 4 A , the second depends on both such ratio and the normalization of D , which is fixed to the diffusion coefficient obtained from analyses of CR secondary-to-primary ratios. Since the normalization of the D in the Galaxy is uncertain by at least ∼ 30%, mainly due to cross sections uncertainties [43][44][45], we have also tested other values of the normalization of D around D 0 = 3.8 × 10 28 cm 2 s −1 , as we discuss below.
In Figure 2, we show the constraint on the TeV halo population in the parameter space of M A and ψ incl . Only a very reduced space of inclination angles (ψ incl < 5 • ) is able to simultaneously account for the radial size and measured symmetry of a typical TeV halo. This means that, unless there is a reason to believe that all existing TeV halos are aligned with our LoS, the anisotropic model is not able to explain the observation of multiple symmetric TeV halos and the lack of observed asymmetric ones. Quantitatively, we note that if all inclination angles are equally probable, then the probability of observing a TeV halo with an inclination angle less than The parameter space of inclination angles (ψ incl ) and asymmetric diffusion parameters which produce TeV halos that fulfill the symmetry (Z/R < 2) and size (θc ∼ θGeminga ± 1σ) conditions at both 68% and 82% containment. In Figure 7 of the appendix, we show the resulting TeV halo parameter space for values of D , D 0 , covering a one order of magnitude (from 10 28 to 10 29 cm 2 s −1 ). As expected, larger values of D 0 further reduce the allowed parameter space, because the halo becomes more extended, while lower D 0 values only slightly increase the fraction of inclination angles that satisfy both conditions, even for our extreme values that are in significant tension with galactic secondary-to-primary ratios if these values are standard for the Galaxy.
The relatively simple symmetry and size conditions already rule out the vast majority of the M A ψ incl parameter space. Additionally, the integrated profile is expected to show clear signatures of anisotropic diffusion when performed following independently perpendicular axes (i.e. the integrated γ-ray emission along the different axes is expected to be different). These signatures would be detectable, although no collaboration reported this kind of differences yet. In fact, both the results from Fig. 2 and the integrated profile are complementary information that must be used to prove or discard the observation of asymmetric TeV halos.
In Figure 3, we show the emission profiles in both the r (perpendicular to the magnetic field) and z (parallel to the field) directions for values of M A = 0.1, M A = 0.3 and M A = 1. These are computed integrating the emissivity obtained from Eq. 3 along our LoS at each projected axis. Given the anisotropic structure of the predicted halos, the profile is computed along the Z and R axes separately. A crucial point here is that the computation of the profile by the HAWC collaboration was done taking the average emission at circular rings around the source, which does not allow the observation of any feature of anisotropy or asymmetry from the halo. Therefore, this kind of profile is not appropriate for asymmetric (anisotropic) objects. We point out that the recent works studying TeV halos compute their profile assuming circular symmetry as well [9,16,25,38,46]. However, we go beyond previous works and calculate the independent profiles in the r and z direction, which allows us to gauge the asymmetry of our model. To have a qualitative comparison, HAWC's surface brightness [17] is also shown in Fig. 3. As discussed, in the case of an asymmetric halo in the ψ incl = 90 • case, observations would detect a profile that is starkly different in each direction (at least for M A ≤ 0.5). This remains valid for angles ψ incl > 0 • .
We stress that our results are not in contradiction with the results found by the authors of Ref. [25], since, in fact, the values of M A and ψ Incl that are compatible with our conditions are M A 0.3 and ψ incl 5 • Instead, our results indicate that the phase space for this solution is small, and the probability of having multiple systems in such a configuration is extremely low.

IV. DISCUSSION AND CONCLUSIONS
TeV halos constitute a new class of astrophysical objects which have the capability to significantly advance our understanding of galactic diffusion [47]. In this work, we have demonstrated that one of the more popular models, where anisotropies in local diffusion explain the TeV halo morphology, is inconsistent with TeV halo observations. Specifically, we have explored and analyzed different morphological signatures of anisotropic diffusion that are predicted by this model but are not observed in detected TeV halos.
We have analyzed the morphology of anisotropic TeV halos as a function of two key parameters: M A , which controls the ratio of the diffusion coefficients perpendicular to and along the background magnetic field, and ψ incl , which controls the angle between the magnetic field and the observer's LoS. Our results constrain anisotropic TeV halo models in three ways: (1) we constrain M A to be smaller than ∼ 0.5 to prevent the TeV Halos from becoming too large compared to current measurements, (2) we constrain ψ incl to be less than ∼ 5 • in order to prevent observed TeV halos from having a significant visual asymmetry that would appear oblong or "spaghetti shaped" on the sky, (3) we show that the expected surface brightness along different axes is significantly different for asymmetric objects, which could lead to easily discard values of M A smaller than ∼ 0.3. In this context, we stress that it would be crucial to experimentally measure the integrated emission projected along the different axes, since such a test would be able to unequivocally detect signatures of anisotropic diffusion.
To be precise, our analysis specifically indicates that anisotropic diffusion cannot explain the observation of several TeV halos in any scenario where the diffusion coefficient in the uninhibited direction is compatible with best-fit values from galactic secondary-to-primary ratios. Our models leave open the possibility that the diffusion coefficient surrounding TeV halos is mildly anisotropic. However the diffusion coefficient in every direction must be significantly inhibited compared to the average diffusion coefficient of the Milky Way. As a result, TeV halos must occupy or generate regions with unique diffusion characteristics compared to the Milky Way average. We note that these conclusions are generically true for any CR source in the galaxy, and indicate that particle diffusion near sources with observable γ-ray emission cannot simultaneously be strongly anisotropic and have uninhibited diffusion along the preferential direction. Although our analysis is based on the anisotropic-diffusion model put forward in Yan and Lazarian [32], our conclusions remain valid for any anisotropic model where the scalings of the perpendicular and parallel diffusion coefficients are similar -namely δ δ ⊥ , given the typical parameterization D ,⊥ ∝ E δ ,δ ⊥ -which is supported by numerical simulations [48].
In conclusion, observations of suppressed and spherically symmetric diffusion provide further credence in favor of models where the diffusivity is reduced not due to a geometrical effect, but rather intrinsically inhibited/subdominant due to subtle mechanisms, either generated by the compact object or pre-existing in the region. These models include, for instance, (i) models with self-generated turbulence that efficiently con-fine CRs more than typical for the ISM [20][21][22], (ii) models with rectilinear propagation in the first stage of the particle injection [38], or (iii) models where the correlation length for the magnetic field is extremely small (∼ 1 pc), such that particles are trapped within the magnetic field structure of the halo on timescales equivalent to HAWC-observations [49].
This appendix aims at detailing the numerical scheme implemented to solve the transport equation described in the text. We use the Crank-Nicolson (CN) expansion, which is second order in energy, space and time. A detailed description of our numerical discretization, an example script with this numerical prescription and a document clarifying some of the relevant numbers for the Geminga TeV halo based on Ref. [2] are publicly available at https://github.com/tospines/ Analyses-and-plotting-codes/tree/main/Anisotropic_TeV_Halos. The numerical algorithm for our scheme is given by: Here u is the density of electrons, E represents their energy and the subscripts i, k and are the indexes of the spatial step in r and z direction and energy bin, respectively, and τ is the time-step index. The diffusion coefficients in the perpendicular and parallel direction to the magnetic field (defined as the z-axis) are denoted as D and D ⊥ , respectively. σ T = 6.652 × 10 −29 m 2 is the Thomson cross section, the energy density of the magnetic field is denoted as U B = B0 4π , the energy density of the different photon fields included in this study is U ph = U CMB + U IR + U opt + U UV and the term b = dE dt = − where T i is the temperature of each ISRF component. Recently, an analysis by Ref. [42] showed that this approximation of the true Klein-Nishina cross-section is now insufficient to describe highly-precise AMS-02 data at the percent level. However, the accuracy of this fit is more than sufficient to describe the electron cooling parameters which are shown here.
To speed up the computations, we differentiate at first order over energy (O(dE)), which provides sufficient accuracy so long as the spacing of our energy steps (δE) is larger than the energy-loss in a given time-step (δt). Our results are accurate to second order in time and space (O(dx 2 dt 2 )). Equation (B1) is not directly solvable in its form, so we solve it applying the alternatingdirection implicit (ADI) method. This requires converting Equation (B1) into two equations (concretely, we follow the Peaceman and Rachford scheme [50]), which differentiate in steps of ∆t/2 (implying that in the injection ∆t → ∆t/2). Eq. B3 solves the equation implicitly in the z-direction while explicitly in the r-direction, and Eq. B4 evolves solving the equation explicitly in the z-direction and implicitly for the r-direction. The coefficients involved, from A to S (A' to S'), result from rearranging the new equations and keeping all the discretized terms at time τ + dt/2 (τ + dt) on the left-hand side and all the terms at τ (τ + dt/2) on the right-hand side. These coefficients are commonly referred to as Crank-Nicolson coefficients. Finally, each of these equations is easily solvable using (tri-diagonal) matrix operations.
These two equations can be expressed as:  In each case the pulsar efficiency is set to 1, which means that models which exceed the HAWC data are potentially compatible with HAWC observations, while models that fall below the HAWC data are in tension with HAWC observations due to energetic considerations.
where the terms from A to G are the Crank-Nicolson coefficients associated to each of the density bin indexes. The left panel of Figure 5 shows the propagated spectrum of electrons.This emission is proportional to the flux of gamma rays generated from IC process, as shown in Figure 5 (right), where we also include the Geminga gamma-ray flux measured by HAWC [17]. Here we observe that at about 50 TeV, the emission starts to be suppressed by the Klein-Nishina effect and the injection cut-off. In this figure the efficiency is set to 1. In Figure 6 we show the ratio of the extension of the simulated TeV halo in the projected Z-axis divided by the extension of the halo in the projected R-axis for different inclination angles (ψ incl ) with respect to our LoS. As this ratio becomes bigger, bigger is the asymmetry of the object. Here, we see that models with M A smaller than ∼ 0.6 become increasingly incompatible with the isotropy of TeV halo observations as seen from most inclination angles. In Figure 7 we show the constraint on the TeV halo population in the parameter space of M A and ψ incl . Only a very reduced space of inclination angles (ψ incl < 5 • ) is able to simultaneously account for the radial size and measured symmetry of a typical TeV halo. The left panel shows the allowed M A and ψ incl parameter space assuming a normalization of the diffusion coefficient D 0 = 10 28 cm 2 s −1 while the right panel assumes a normalization D 0 = 10 29 cm 2 s −1 . As discussed in the main text, we observe that even considering slightly different diffusion coefficient normalizations the allowed parameter space able to explain the size and symmetry conditions imposed is very limited. We also observe here that as we go to lower normalizations, the parameter space increases, meaning that inhibited diffusion is preferred.  Figure 2, but for normalization of the diffusion coefficient set to D0 = 10 28 cm 2 s −1 (left panel) and D0 = 10 29 cm 2 s −1 (right panel). These models are both extreme, already residing in significant tension with cosmic-ray secondary-to-primary ratios.

Appendix C: Line-of-sight integrated emission -2D Projected halos
In this section, we report the 2D projected surface brightness of the simulated halos for the M A = 0.3 case, both for ψ incl = 0 • and ψ incl = 90 • , since these images can be directly compared to the ones reported by experiments (Fig. 8).