Inferring the Interstellar Magnetic Field Direction from Energetic Neutral Atom Observations of the Heliotail

Determining the magnitude and direction of the interstellar magnetic field (B ISM) is a long-standing problem. To date, some methods to infer the direction and magnitude have utilized best-fit models to the positions of the termination shock and heliopause measured by Voyager 1 and 2. Other models use the circularity of the Interstellar Boundary Explorer (IBEX) ribbon assuming a secondary energetic neutral atom (ENA) mechanism. Previous studies have revealed that the B ISM organizes the orientation of the heliotail with respect to the solar meridian. Here we propose a new way to infer the direction of the B ISM based on ENA observations of the heliotail. IBEX observations of the heliotail have revealed high-latitude lobes of enhanced ENA flux at energies >2 keV. Analyses showed that the high-latitude lobes are nearly aligned with the solar meridian, while also exhibiting a rotation with solar cycle. We show, using steady-state solar wind conditions, that the inclination of the lobes reproduced with commonly used values for the angle (α BV ) between B ISM and the interstellar flow in the hydrogen deflection plane (40° < α BV < 60°) is inconsistent with the IBEX ENA observations. We report that 0° < α BV < 20° best replicates the heliotail lobe inclinations observed by IBEX. Additionally, our model results indicate that the variation of the solar magnetic field magnitude with solar cycle causes the longitudinal rotation of the lobes observed by IBEX by affecting the inclination of the lobes.


Introduction
The heliosphere is created through the interaction of the solar wind and interstellar medium (ISM; Parker 1961; Baranov et al. 1971;Holzer 1989;Baranov & Malama 1993;Dialynas et al. 2022;Galli et al. 2022;Kleimann et al. 2022, and references therein).Embedded in the ISM is the B ISM , which is characteristic of the local galactic environment through which the Sun is traveling.With Voyager 1 and 2 measurements at only 38 and 14 au beyond the heliopause (Gurnett et al. 2013;Krimigis et al. 2013Krimigis et al. , 2019;;Stone et al. 2013;Burlaga et al. 2019;Gurnett & Kurth 2019;Richardson et al. 2019), the B ISM is still mediated by the heliosphere (Burlaga et al. 2022).With no direct observations of the B ISM in the unperturbed ISM, indirect methods are required to decipher its orientation and magnitude.Aside from Voyager observations, there are several ways to infer the magnitude of the magnetic field of the ISM, one such being through a pressure balance at the heliopause and/or between the local ISM and the local interstellar cloud (e.g., Dialynas et al. 2019;Dialynas et al. 2020;Linsky et al. 2022), but obtaining the direction of the B ISM has been more challenging.
While neutral hydrogen from the ISM is affected by charge exchange with the interstellar plasma piled up in front of the heliosphere, neutral helium can pass through the heliosphere relatively unimpeded.Lallement et al. (2005Lallement et al. ( , 2010) ) found that the inflow direction of hydrogen (λ = 252°.5 ± 0°.7, β = 8°.9± 0°.5, where λ and β are the ecliptic longitude and latitude, respectively) and the direction of the pristine interstellar helium (λ = 255°.4± 0°.5, β = 5°.2± 0°.2; Witte 2004) define the hydrogen deflection plane (HDP), which constrains the direction of the B ISM .Using measurements from the Solar Wind ANisotropy experiment aboard the Solar and Heliospheric Observatory spacecraft, they found that there is a deflection of the inflowing hydrogen from that of the pristine local interstellar gas flow by 4°± 1°.Lallement et al. (2005) suggested that the HDP reflects an asymmetry between the pristine ISM and the B ISM .This was supported by Izmodenov et al. (2005), who confirmed the deviation of the interstellar neutral hydrogen flow from the pristine local interstellar gas flow when the B ISM is offset from the interstellar flow.Pogorelov et al. (2008) also demonstrated via modeling that the deflection of hydrogen atoms affected by the presence of the heliosphere occurs in the plane defined by the vectors of the B ISM and interstellar flow velocity.Additionally, the direction of the secondary helium (the so-called "Warm Breeze") co-aligns with the deflection plane (Kubiak et al. 2016).
While Lallement et al. (2005) were able to constrain the asymmetry of the pristine ISM flow and B ISM direction to both lie in the HDP, the angle between the B ISM and ISM flow velocity vectors (α BV ) at this point remained an unknown quantity along with the magnitude of the B ISM .When Izmodenov et al. (2005) first modeled the deflection of neutral hydrogen from the ISM flow, they found that a B ISM magnitude of B ISM ∼ 2.5 μG with an α BV = 45°could reproduce the deviation measured by Lallement et al. (2005).Opher et al. (2006) were the first to show quantitatively that a B ISM in the HDP reproduces the termination shock asymmetries observed by Voyager 1 and 2. Opher et al. (2007) found that α BV = 30°with B ISM = 1.8 μG produced MHD model results most consistent with Voyager radio observations that span across the sky as a ribbon (Gurnett et al. 2006).Pogorelov et al. (2007) found that when including the solar magnetic field the HDP does not exactly coincide with the plane defined by the interstellar neutral helium and hydrogen velocity vectors in the supersonic solar wind region, but a B ISM magnitude of B ISM > 3 μG can account for the observed asymmetry.Pogorelov et al. (2008) found that to reproduce the deflection of neutral hydrogen from the pristine ISM flow α BV = 30°and B ISM = 3 μG was needed.Izmodenov (2009) used kinetic-MHD modeling to find that α BV = 15°-30°a nd B ISM = 2.5-3.5 μG matched best with Voyager termination shock crossings but neglected the solar magnetic field in their study.Opher et al. (2009) constrained the values of α BV to lie between 20°and 30°with B ISM = 3.7-5.5 μG based on multifluid MHD modeling of the Voyager termination shock crossings.
When the Interstellar Boundary Explorer (IBEX; McComas et al. 2009b) started observing energetic neutral atoms (ENAs) in 2009, a ribbon of enhanced ENA flux was discovered (McComas et al. 2009a) that formed a circular arc centered on (λ, β) ∼ (221°, 39°), which was thought to be near the direction of the B ISM (Funsten et al. 2009).The source of this ribbon, referred to as the IBEX ribbon, had been an open question in the community for some time (see references within McComas et al. 2014).The source for the IBEX ribbon is a secondary ENA mechanism, where ENAs created in the solar wind are able to leave the heliosphere, where they charge exchange with the ISM plasma and become pickup ions (PUIs); however, the microphysics of PUIs as they gyrate along the B ISM is still under debate.These PUIs then gyrate around the B ISM until they charge exchange again, becoming ENAs (so-called secondary ENAs), which are observed by IBEX. McComas et al. (2017) showed that the IBEX ribbon and globally distributed flux (GDF) from the heliosheath evolve differently with solar cycle, thereby supporting the secondary ENA mechanism.The ribbon was interpreted as the location of the B ISM • r = 0 surface (Schwadron et al. 2009(Schwadron et al. , 2011)).Using the circularity of the IBEX ribbon and this assumed secondary ENA mechanism, Zirnstein et al. (2016b) found that the B ISM orientation and magnitude that best replicated the IBEX ribbon shape and structure was B ISM = 2.9 μG and α BV = 40°, with an ecliptic longitude and latitude of (λ, β) = (227°, 34°).
At the same time, higher-energy ENAs were measured by the Ion and Neutral Camera (INCA; Krimigis et al. 2009) on the Cassini mission around Saturn.These observations showed the existence of a high-intensity, relatively wide, latitudedependent ENA region that loops through a rough circle in the sky in ecliptic coordinates.Dialynas et al. (2013) considered the interpretation of the INCA ENA measurements with respect to the B ISM .They showed that the deviation of the INCA ENA emissions from the ecliptic equator is effectively minimized in a rotated frame where the B ISM is normal to the radial direction outward from Saturn, where its north pole points toward (λ, β) = (190°, 15°).
More recent works have suggested alternative B ISM configurations.Opher et al. (2020) used a multi-ion (separating the cold and hot components of the solar wind) multifluid MHD model of the heliosphere and found that using α BV = 40°from the derived IBEX ribbon direction of Zirnstein et al. (2016b) and B ISM = 3.2 μG produced reasonable agreement with the Voyager 1 and 2 magnetic field observations outside of the heliopause.Izmodenov & Alexashov (2020) used time-dependent kinetic-MHD modeling and showed that α BV = 60°and B ISM = 3.75 μG best replicated the Voyager 1 and 2 termination shock and heliopause asymmetries.Frisch et al. (2022) used dust grain polarization and found an B ISM orientation close to that of Zirnstein et al. (2016b); however, the dust polarization is sensitive to much larger scales (∼parsecs) and not to the field closer to the heliosphere on astronomical unit scales where a majority of the ENAs observed by IBEX are produced.Rankin et al. (2023) suggested that the interstellar magnetic field derived from the IBEX ribbon from Zirnstein et al. (2016b) matched well with the overall B T and BN (i.e., tangential and normal, respectively) components of the interstellar magnetic field as observed in the Voyager 1 and 2 spacecraft frames, though the ribbon-derived direction was unable to match the B R (i.e., radial) component of the field.While the conclusions of Rankin et al. (2023) are currently being debated within the community, one important result is the need to also compare modeled results to the observed magnetic field components in addition to the observed azimuthal and elevation angles.McComas et al. (2013) presented the first tail-centered ENA maps of the heliosphere from IBEX and discovered that the tilt of the port/starboard lobes is caused by the magnitude and direction of B ISM .Schwadron et al. (2014) isolated the ENA contribution from the heliosheath in the form of the GDF and found the presence of two high-latitude lobes at energies >2 keV in the tail of enhanced ENA flux.Opher et al. (2015) used global MHD modeling to first suggest that these highlatitude lobes could result from the confinement of the solar magnetic field into two high-latitude lobes, in accordance with the simplified analytic model by Yu (1974), which suggested the importance of the solar magnetic field in shaping the heliotail.Zirnstein et al. (2016aZirnstein et al. ( , 2017) ) used a time-dependent simulation of the heliosphere with solar cycle variations and concluded that the high-latitude heliotail lobes reflected the fast/slow solar wind profile, as they were able to qualitatively reproduce the heliotail ENA maps shown by IBEX.One limitation of their modeling was the combined effects of the solar magnetic field and solar wind profile, which could not be separated.Kornbleuth et al. (2018Kornbleuth et al. ( , 2020) ) investigated the effect of the solar magnetic field versus the solar wind latitudinal profile on ENA maps using steady-state MHD solutions and concluded that the solar magnetic field also plays an important role in the emergence of the high-latitude heliotail lobes.This indicated that the appearance of the high-latitude lobes in the IBEX observations is a function of both the solar cycle variations of the solar wind and the solar magnetic field.
In addition to the previously mentioned theoretical and modeling efforts studying the high-latitude lobes, Dayeh et al. (2022) tracked the size and location of the high-latitude heliotail lobes (herein referred to as heliotail lobes) over a solar cycle of observations and found a cyclic behavior of the lobes, with a change in their inclination angle, defined here as the angle of each lobe center from the bisecting solar meridian.The important conclusion from Dayeh et al. (2022) is that the lobe inclinations were nearly aligned with the solar meridian, while previous ENA modeling efforts (Zirnstein et al. 2017;Baliukin et al. 2020;Kornbleuth et al. 2020Kornbleuth et al. , 2023a) ) have shown a much more inclined lobe profile.Fahr et al. (1986Fahr et al. ( , 1988) ) were the first to study the effect of the B ISM on the exterior heliosphere using a simplified model that utilized a single surface to reflect the termination shock, heliopause, and bow shock.Following this study, Ratkiewicz et al. (1998) used MHD modeling to investigate this effect and found asymmetries arising from various α BV angles.Previously, Opher et al. (2007Opher et al. ( , 2009) ) showed that the heliosphere orients itself along the direction of the B-V plane between the B ISM and interstellar flow.Prested et al. (2010) used ENA modeling to demonstrate that different B ISM directions will manifest in ENA maps by changing the orientation of observable profiles.Both Heerikhuisen et al. (2014) and Kleimann et al. (2016) performed a parameter study involving different magnitudes of B ISM and found that stronger intensities led to more highly inclined, ellipsoidal heliotail configurations.Opher et al. (2017) demonstrated through MHD modeling that the B ISM can reconnect with the solar magnetic field and drape along the heliopause, yielding an additional magnetic pressure as also demonstrated in Kornbleuth et al. (2021).Therefore, the interstellar magnetic field has the ability to orient the heliosphere, and by extension the heliotail, which may be observable through ENAs.
In this work, we attempt to address the influence of the B ISM on the heliosphere by investigating its effects on the highlatitude heliotail lobe profile seen in IBEX ENA observations.Using steady-state MHD modeling (i.e., no time variations), we compare the effects of the solar wind profile, solar magnetic field, and interstellar magnetic field on the heliotail lobe inclinations inferred from ENA flux.In Section 2, we present the MHD model (Section 2.1) and the ENA model (Section 2.2) used in this work.In Section 3, we present the effects of the solar wind and solar magnetic field on the heliotail lobes.In Section 4, we present the effects of the B ISM on the heliotail lobes.Lastly, in Section 5, we discuss the results.

Models
For the MHD models, we use a kinetic-MHD model, which we describe in Section 2.1.Additionally, we use the ENA model from Kornbleuth et al. (2023b), which is based on the model from Kornbleuth et al. (2018Kornbleuth et al. ( , 2020) ) and described briefly in Section 2.2.

Global Model of the Heliosphere
We use the BU kinetic-MHD model of the heliosphere that is based on Opher et al. (2009Opher et al. ( , 2015) ) and Chen et al. (2024), developed as the outer heliosphere component of the Space Weather Modeling Framework (Tóth et al. 2005).The model is a 3D kinetic-MHD model that solves for a single-ion plasma that combines the cold solar wind plasma and hot PUIs as one fluid.The BU model allows for numerical reconnection at the heliopause between the solar magnetic field and B ISM , although reconnection is suppressed in the nose while allowing it in the flanks (see Opher et al. 2017 for more details).
To support efficient time-accurate kinetic-MHD simulations of the heliosphere, we have implemented a new kinetic model as the Particle Tracker (PT) component (Chen et al. 2024).The new model is based on the FLexible Exascale Kinetic Simulator (FLEKS; Chen et al. 2023), which was originally developed as a particle-in-cell code.We have extended the capacity of FLEKS to support solving the Boltzmann equation for neutrals, i.e., tracking the motion of neutral particles and simulating the charge exchange between the kinetic neutrals and the plasma fluid.Compared to the previous implementation based on AMPS (Michael et al. 2022), we have developed two advanced features into FLEKS to enable time-accurate simulations: adaptive mesh refinement (AMR) and particle resampling algorithms.AMR mesh reduces the total number of cells and particles while maintaining resolution in areas of interest.Since the charge-exchange process generates new neutrals, which requires launching new macroparticles in simulations, the number of macroparticles in the simulation domain will increase rapidly.To solve this problem, we have designed particle resampling algorithms, including a splitting algorithm and a merging algorithm, to regulate the macroparticle number of each cell.
When running the BU model, the Outer Heliosphere (OH) stand-alone component is used first with the multifluid neutral approximation to achieve a steady state.Then, the OH component drops the neutral fluids and only simulates the plasma fluid, and the PT component launches macroparticles based on the steady-state neutral profiles as its initial conditions.The OH and PT components are coupled and the source terms for the plasma fluid of the OH component are obtained from the PT component by calculating the charge exchange between each macroparticle and the plasma fluid.We note that for the simulations included in this work, we first produce steady-state plasma solutions with the multifluid neutral approximation.After obtaining a steady-state plasma solution with multifluid neutrals, we generate a kinetic neutral solution using FLEKS corresponding to the steady-state plasma solution by cycling FLEKS for 100 yr to obtain sufficient statistics (similar to the process used in Michael et al. 2022).We do this to save computational time owing to the number of cases used in this work.As shown in Michael et al. (2021), a steady-state kinetic-MHD solution leads to a shrinking of the heliosphere, notably in the heliotail; however, the inclination of the heliosphere traveling in the ISM does not change.Therefore, while allowing the plasma solution to converge to steady state with a kinetic neutral solution could affect the latitude of the heliotail lobes, it would likely not affect the longitude of the lobes.
We consider 11 different simulated steady-state MHD solutions with different solar wind and interstellar conditions.The goal is to identify physical processes within the heliosphere that can affect the inclination of the high-latitude heliotail lobes observed by IBEX.Table 1 summarizes all 11 cases studied.
In Cases 1-4 we apply the same B ISM orientation as derived from the secondary ENA mechanism for the creation of the IBEX ribbon, and it is α BV = 40° (Zirnstein et al. 2016b)  For the solar wind conditions in Cases 1 and 2, we use yearly averaged data corresponding to times of solar minimum and solar maximum.Since the heliotail contains a layering of wind from different years, there is not one particular time that would precisely correlate with the observed ENA fluxes in our steadystate modeling.Therefore, we use solar wind conditions from the years 2008 and 2014 for Cases 1 and 2, respectively, as they reflect solar minimum and solar maximum conditions during the IBEX observational window.While this method does not truly replicate time-dependent modeling owing to the absence of fast and slow solar wind layering over time, our alternative method is likely to yield similar results.We do note that the traceback times of solar minimum and solar maximum for the IBEX ENAs noted in this paper come from Dayeh et al. (2022), who obtain their simulated time-delay estimates from Zirnstein et al. (2021).
Case 1 has a plasma density and speed of 5.1 cm −3 and 437.7 km s −1 , respectively, at 1 au in the ecliptic plane and a faster, less dense wind at the poles (see Figure 1 of Kornbleuth et al. 2020) with the assumption of heliolongitudinal symmetry (Sokół et al. 2015).For the solar magnetic field, we assume a Parker spiral and assume a unipolar configuration as in Opher et al. (2015) to minimize spurious numerical reconnection.We use a B r,SW = 29.4μG for Case 1, which is the averaged radial magnetic field magnitude during the year 2008 based on OMNI magnetic field data at 1 au (King & Papitashvili 2005).Case 2 has a plasma density and speed of 6.05 cm −3 and 407.9 km s −1 , respectively, at 1 au in the ecliptic plane, as well as similar densities and speeds at higher latitudes (Sokół et al. 2015).Like in Case 1, the solar magnetic field magnitude for Case 2 is derived from OMNI magnetic field data corresponding to the year 2014, with B r,SW = 44.4μG at 1 au.For Case 3, we use the same solar magnetic field magnitude as Case 1, but with a uniform solar wind profile akin to solar maximum conditions.Case 3 has a plasma density and speed of 7.9 cm −3 and 417.1 km s −1 , respectively, at all latitudes and longitudes at 1 au, corresponding to OMNI data from 1984.5 to 2006.5 from Izmodenov et al. (2008) and Opher et al. (2015).For Case 4, we apply the same uniform solar wind plasma profile as in Case 3, but with a stronger solar magnetic field magnitude, with B r,SW = 64.5 μG at 1 au.
We implement the above solar wind conditions at 10 au within our MHD model.In order to do so, we extrapolate the solar wind conditions from 1 au to the inner boundary at 10 au by assuming adiabatic expansion.We also extrapolate the solar magnetic field from 1 to 10 au by treating the solar magnetic field as a Parker spiral (Parker 1958).
For Cases 5-8, we only change the B ISM orientation.For the solar wind conditions of these cases, we use 22 yr averaged solar cycle conditions (1995-2017) from Izmodenov & Alexashov (2020) that were first applied in the BU MHD model in Kornbleuth et al. (2021).For the solar magnetic field a Parker solution is assumed, with the radial component of the magnetic field strength being B r,SW = 37.5 μG at 1 au.For the B ISM , we change α BV while keeping the B-V plane as the HDP for each case, with Cases 5-8 reflecting an α BV of 60°, 40°, 20°, and 10°, respectively.We note that the direction of the primary interstellar flow remains unchanged in our modeling, and only the orientation of B ISM is changed.
For Cases 1-8 we assume a B ISM magnitude of B ISM = 3.2 μG.For Cases 9 and 10 we vary the B ISM magnitude while keeping α BV = 40°fixed to investigate the effect of magnitude.For Case 9 we use a lower B ISM magnitude of B ISM = 2.9 μG based on the best fit to the IBEX ribbon from Zirnstein et al. (2016b).For Case 10 we use a higher B ISM magnitude of B ISM = 4.4 μG, which is a value originally from Izmodenov (2009) that best fit the Voyager termination shock crossings, but with a different value of α BV .
Cases 1-10 use the same ISM conditions with the exception of the B ISM .All the neutral and ionized populations in the ISM are assumed to have the same bulk velocity v ISM = 26.4km s −1 (longitude = 75°.4,latitude = −5°.2 in ecliptic J2000 coordinate system) and temperature T ISM = 6519 K at the outer boundary 1500 au from the Sun, where the pristine ISM is not mediated by the heliosphere.In the ISM, the proton density is assumed to be n p,ISM = 0.06 cm −3 and the neutral H atom density is n H,ISM = 0.18 cm −3 , which is in agreement with the pristine 1°. 2 ± 2°. 5 Data (14-15) 1°. 0 ± 1°. 9 Data (16-17) 0°. 1 ± 0°. 4 Note.The inclination angle (ω) is determined by the angle of the lobe centers from the bisecting solar meridian (see Figure 1), while α BV is the angle between the B ISM and ISM flow directions in the HDP.Data from For Case 11, we increase the ISM proton density to be n p,ISM = 0.09 cm −3 , while keeping other ISM parameters the same.We consider this case because helium ions are not included in our modeling, so n p,ISM = 0.09 cm −3 accounts for the remainder of the effective mass density (Bzowski et al. 2019).While this case overpredicts the ISM pressure on the heliosphere because protons charge exchange at a significantly higher rate than helium ions, it offers insight into the effect of the ISM ram and thermal pressures on the heliotail lobes.For Case 11 we also use B ISM = 3.2 μG and α BV = 40°.

ENA Model
Here we present a brief overview of our ENA model, which is described in Kornbleuth et al. (2023aKornbleuth et al. ( , 2023b)).We model the ion flux at the termination shock as initially presented in Gkioulidou et al. (2022) and find the best-fit ion populations and distributions to match the hybrid simulation results at the termination shock obtained by Giacalone & Decker (2010) and Giacalone et al. (2021).These ion populations and distributions then propagate along plasma velocity streamlines in the heliosheath, where they undergo charge exchange to become ENAs and therefore also experience preferential losses at higher energies as they are replaced by lower-energy ions (i.e., extinction).
We capture the ENA flux along a particular line of sight (LOS) by integrating where ¢ r is the vector along a particular LOS as a function of θ and f, s is the length of the ¢ r vector, m p is the mass of a proton, and f p is the phase-space velocity distribution, which has a different form depending on the ion population being considered (see Kornbleuth et al. 2023b for further details).The velocity of the parent ion in the frame of the plasma is given by v plasma = |v p − v i |, with v p and v i being the velocities of the bulk plasma and the parent proton, respectively.The ion density, n i , is derived from the ion density fraction for a given population, while the ion temperature, T i , is derived from the ion density and energy fractions as T i = E i /n i .T i is a fraction of the local MHD temperature based on the thermal pressure fraction of the ion species.We assume quasi-neutrality and use the approximation that the electrons have the same temperature as the thermal solar wind ions.Additionally, n H is the neutral H density, and σ ex is the charge-exchange cross section (Lindsay & Stebbings 2005).The survival probability S(E) of ENAs in our model represents the likelihood that an ENA created in the heliosheath will charge exchange prior to being observed by IBEX and is calculated based on the work of Bzowski (2008).We calculate the survival probability of an ENA along a radial LOS because the trajectories of H atoms with the energies modeled here are almost straight.We place the observer at the termination shock such that the survival probability correction is out to 100 au.Considering the distance from the Sun, we only calculate the survival probability due to charge exchange.To compare with IBEX-Hi observations, we use the IBEX-Hi energy transmission functions (http://ibex.swri.edu/ibexpublicdata/CalData/Hi/)for each IBEX-Hi energy step in our modeling.
One key aspect of our ENA modeling for this study is the ability to identify the centers of the high-latitude heliotail lobes.Identifying the centers of the lobes is critical for determining the inclination of the lobes.We employ a similar method to Zirnstein et al. (2016a) and Dayeh et al. (2022), where, after defining the boundaries of the lobes by isolating the region of enhanced flux relative to the GDF, we determine the lobe centers that are weighted by the fluxes.We do this by converting the location of each pixel into a latitudinal/ longitudinal unit vector weighted by its flux, and then we find the lobe centers by combining these vectors (similar to a centerof-mass calculation).We then find the inclination of the lobes in each modeled case by evaluating the angle of the north/ south lobe centers relative to the solar meridian (see Figure 1).For each case, we take the arctangent of the difference in the ecliptic longitudes and the difference in the ecliptic latitudes of the north/south lobe centers to obtain the inclination angle.

Effect of Solar Wind Conditions on Heliotail Lobe Inclinations
Here we consider how varying solar wind conditions can affect the high-latitude heliotail lobes as seen in 4.29 keV IBEX ENA observations and ENA modeling at the corresponding energy band.We choose this energy band because it is least affected by the IBEX ribbon.We separate the two effects to determine whether the solar wind profile or solar magnetic field has the ability to affect the lobe inclination.While the ENA model in Kornbleuth et al. (2020) does not reflect the updated PUI distributions presented in Kornbleuth et al. (2021Kornbleuth et al. ( , 2023a)), the results of Kornbleuth et al. (2023b) demonstrated that the improved PUI modeling only affects the quantitative flux intensity and not the qualitative profile of the modeled ENA maps.Therefore, using the ENA results from Kornbleuth et al. (2020) for Cases 1-4 is sufficient to understand the effect of the solar wind effects on the high-latitude heliotail lobes.
For Cases 1 and 2, which reflect steady-state solutions using realistic, time-averaged solar wind profiles corresponding to solar minimum and maximum, respectively, we note a rotation of the high-latitude heliotail lobe centers as seen in Table 1 and Figure 2. The inclination of the lobes shifts by approximately 9°from 2008 to 2014.While this rotation is notably larger than is seen in IBEX observations noted in Dayeh et al. (2022) of approximately 1°, the decrease in inclination angle from solar minimum to solar maximum in our simulations is in qualitative agreement with the IBEX observations when accounting for transit time delays of ENAs corresponding to times of solar minimum and maximum.However, considering that the two cases differ in both solar wind profile and solar magnetic field magnitude, the effects of each need to be isolated.
For Cases 1 and 3 in Table 1, we present a comparison of the lobe inclination angles when changing the solar wind profile, 2. Evolution of the north lobe (NL) and south lobe (SL) centers with changing solar wind profile and solar magnetic field.ENA sky maps show flux centered on the downwind (tail) direction that is normalized to the maximum flux in each map corresponding to Cases 1-4 listed in Table 1.For each map, the lobe center is displayed as the zeroed (blue) pixel within each north/south heliotail lobe.Top: Case 1 (left) corresponds to 2008-averaged solar minimum conditions with a weak solar magnetic field (B r,SW = 2.94 μG), while Case 2 (right) corresponds to 2014-averaged solar maximum conditions with a weak solar magnetic field (B r, SW = 44.4μG).Bottom: Case 3 (left) corresponds to a uniform slow solar wind (e.g., solar-maximum-like conditions) with a weak solar magnetic field, and Case 4 (right) also has a uniform slow solar wind, but with a strong solar magnetic field (B r,SW = 64.5 μG).
but keeping all other conditions (i.e., solar magnetic field, B ISM , etc.) identical.Here, as seen in Figure 2, we find that, despite very different solar wind conditions, with Case 1 corresponding to latitudinally varying solar minimum conditions and Case 3 corresponding to uniform (i.e., solar-maximum-like) conditions, the lobe center locations and inclinations are identical for the two cases.In comparing Cases 2 and 4, which both have strong solar magnetic fields, we see that the two cases display similar longitudinal centers for the north and south lobes with different latitudinal centers.While heliolongitudinal symmetry is assumed for the solar wind profile of all cases, Cases 1 and 2 have distinctly different latitudinal solar wind profiles than the uniform solar wind profiles of Cases 3 and 4. The latitudinal profile of the wind therefore has the potential to affect the latitudinal profile of ENAs and also the intensity in the lobes.As shown by Zirnstein et al. (2020), the layering of the fast and slow solar wind in the high-latitude heliotail affects the latitudinal center of the high-latitude lobes seen in ENA observations, which is an effect absent in our steady-state solutions.However, this will not significantly affect the longitudinal offset of the lobes.This is because there is heliolongitudinal symmetry of the solar wind profile in the outer heliosphere due to the rotation of the Sun over a 27-day period.Therefore, longitudinal shifts of the high-latitude lobe centers are less affected by the solar wind latitudinal profile as noted between Cases 1 and 3.Although our conclusions suggest the independence of the lobe centers on the solar wind profile, it is worth noting that the layering of slow and fast solar wind over long times (years) would result in a nonsymmetric SW profile at high latitudes.This could have potential effects on the lobe centers that need to be further investigated.
For Cases 3 and 4 in Table 1, we compare the effect of changing solar magnetic field strength when all other factors (such as solar wind profile) are identical.Here, as seen in Figure 2, we find that an increase in the solar magnetic field strength leads to a decrease in the heliotail lobe inclination angle.Relative to Case 3, in Case 4 the lobe centers move in toward one another owing to the changing solar magnetic field, leading to a lobe inclination angle of 5°.44 ± 2°.72, in contrast to 14°.04 ± 3°.48 for Case 3.This result can be seen in the results of Kornbleuth et al. (2020) as well, where it was shown that the solar magnetic field resists the influence of the B ISM , which results in an increased alignment with the solar meridian.We note that in the results of Dayeh et al. (2022) the effect of solar cycle appears to correlate with the changing heliotail lobe inclination.As the IBEX ENA observations reflect a change from solar minimum to solar maximum, the inclination angle decreases with a trend toward an inclination angle of 0°.This is noted in Table 1, where the 2012-2013 time period observed by IBEX, which corresponds to solar minimum after accounting for transit time delays, shows an inclination angle of 1°.19 ± 2°.47.In contrast, the 2016-2017 time period, which corresponds to solar maximum, shows an inclination angle 0°.07 ± 0°.40.With these results, Dayeh et al. (2022) identified a rotation of the high-latitude heliotail lobes in IBEX-Hi observations, which appears to be consistent with our conclusions that since the solar magnetic field magnitude will increase as the solar cycle progresses from solar minimum to solar maximum, the lobes will experience a rotation.However, due to uncertainties, we cannot make any firm conclusions from the data.In the top panels of Figure 3, we overlay the results of Table 1 for Cases 1-4.

Effect of the B ISM Orientation on Heliotail Lobe Inclinations
While in Section 3 we evaluated the effect of the changing solar wind profile and solar magnetic field strength on the inclination of the high-latitude heliotail lobes, here we isolate the effect of the B ISM and interstellar ram pressure.We consider Cases 5-11 in Table 1, which correspond to the maps shown in Figure 4. We note here that the goal of this study is to determine a general range for the B ISM orientation and that we cannot make any firm conclusions on the exact B ISM conditions, as the effects of the magnitude and orientation are not separable.Here we present initial constraints, with a larger parameter study including more B ISM magnitudes to be conducted in the future.
To investigate the effect of the B ISM orientation, we test four cases: α BV = 60°(Case 5), which corresponds to the best-fit orientation for matching the Voyager crossings from the kinetic-MHD modeling of Izmodenov & Alexashov (2020); α BV = 40°( Case 6), which corresponds to the best-fit orientation for the IBEX ribbon secondary ENA mechanism; α BV = 20°(Case 7), which corresponds to the best-fit orientation from Opher et al. (2009); and α BV = 10°(Case 8), which serves as a lower limit in this phase-space exploration of α BV .
The solar magnetic field strength we use in Cases 5-8 corresponds to a higher solar magnetic field strength than is used in Cases 1 and 3.This value of B r,SW = 37.5 μG for Cases 5-8 more closely reflects the transition between solar minimum and solar maximum.While the 22 yr averaged solar wind conditions used in these cases do not entirely align with the times observed by IBEX-Hi, as noted in Section 3, we believe that the solar magnetic field intensity is more important than the solar wind profile in dictating the longitude of the lobe centers, which influences the inclination of the lobes.Comparing to the heliotail lobe center locations from Dayeh et al. (2022), we focus on the IBEX ENA observations from 2014 to 2015, as this best correlates with the solar cycle transition period.During this time, the observed inclination angle for the lobes is 0°.96 ± 1°.94.Accounting for model uncertainties, this inclination angle is reproduced for α BV = 10°.However, if we considered the period of IBEX ENA observations that best correlated with solar minimum conditions when accounting for approximately 5 yr time delay (2012-2013) from Dayeh et al. (2022), the inclination angle of the lobes is 1°.19 ± 2°.47, which is best matched by α BV = 10°-20°when comparing the inclination angles and uncertainties.We use the B r,SW = 37.5 μG and compare with the transition period owing to the fact that this allows for the best method of comparison between our steady-state modeling and observations.However, the influence of the changing solar cycle, and by extension the solar magnetic field, will lead to different amounts of resistance from the solar magnetic field to the influence of the B ISM .
An interesting result of this study is the similar inclination angle of Case 5 as compared to Cases 1 and 3.As mentioned in Section 3, Cases 1 and 3 have the same inclination angle in our study owing to having the same solar magnetic field magnitude at 1 au and the same interstellar magnetic field magnitude and orientation.For Case 5, which has α BV = 60°, it manages to have the same inclination angle as Cases 1 and 3, which have α BV = 40°.This result follows from our conclusions in Section 3.While α BV = 60°should produce a more inclined profile of the heliotail lobes relative to α BV = 40°, Case 5 also has a stronger solar magnetic field magnitude of B r,SW = 37.5 μG relative to B r,SW = 29.4μG in Cases 1 and 3. Therefore, the result reflects the interplay of both the interstellar and solar magnetic field in orienting the lobes.
In the middle panels of Figure 3, we overlay the results of Table 1 for Cases 5-8 to show the changing inclination angles of the lobes with α BV .We find that α BV = 10°provides the best-fit replication for the observed ENA lobe inclinations.One case we neglected to model here was α BV = 0°.In this case, the B ISM and interstellar flow are aligned, which would yield a heliotail orientation aligned with the solar meridian.Decreasing α BV leads to a decrease in the inclination of the heliotail.Considering that the IBEX ENA observations indicate a rotation of the lobe inclination over time and on heliospheric timescales we do not anticipate any changes in the B ISM orientation or magnitude, our results suggest that α BV > 0°in the HDP, which is in agreement with Voyager observations (e.g., Burlaga et al. 2013Burlaga et al. , 2019)).Considering that the uncertainty of the IBEX ENA observations from 2012 to 2013 falls within the limits of α BV = 20°, we conclude that the best-fit interstellar magnetic field orientation to the IBEX ENA observations has 0°< α BV < 20°in the HDP.
While 0°< α BV < 20°in the HDP yields the best comparison with respect to the inclination of the high-latitude lobes as seen in IBEX-Hi observations, we do note a discrepancy in the nose direction.As seen in Figure 4 for Cases 7 and 8, with a decrease α BV we note a decrease in the ENA flux for the noseward direction at 4.29 keV.This likely arises from the decreasing tangential component of the interstellar magnetic field with decreasing α BV , which affects the magnetic pressure on the noseward heliosphere.Previous studies, such as Opher et al. (2009) and Zieger et al. (2013), which employed an α BV = 20°, suggested a stronger interstellar magnetic field intensity (4.4 μG) than in Cases 1-8, which would increase the magnetic pressure and potentially increase the noseward ENA flux relative to Case 7 shown here.Additionally, inclusion of a time-dependent solar wind is critical to further analyze this phenomenon.
To investigate the effect of the B ISM magnitude, we probe two cases using the same solar wind conditions and interstellar plasma and neutral conditions as in Cases 5-8 and α BV = 40°.We use this particular configuration for α BV because it provides the best fit to the IBEX ribbon (Zirnstein et al. 2016b) and was also used to replicate Voyager 2 observations in Opher et al. (2020).While in Case 6 we used the B ISM magnitude from Opher et al. (2020) with α BV = 40°, in Case 9 we present the IBEX-ribbon-derived magnitude of 2.9 μG from Zirnstein et al. (2016b).In Case 10, we use B ISM = 4.4 μG, which reflects an upper limit on B ISM based on previous modeling efforts (Izmodenov 2009;Opher et al. 2009Opher et al. , 2015;;Zieger et al. 2013), though in those models a value of α BV < 20°was used.2, but with non-normalized ENA flux in units of (cm 2 s sr keV) −1 .Cases 5-11 all use the same solar wind profile (1995-2017 averaged) and solar magnetic field strength (B r,SW = 37.5 μG).For Cases 5-8, we vary α BV from 60°to 10°.For Cases 9 and 10 we replicate Case 6, but with a lower (B ISM = 2.9 μG) and higher (B ISM = 4.4 μG) interstellar magnetic field intensity, respectively.For Case 11, we also replicate Case 6, but with a higher interstellar plasma density (n p,ISM ) to investigate the effect of higher ISM ram pressure on the lobes.
The modeled ENA maps for Cases 9 and 10 are presented in Figure 4, and the lobe centers derived from these maps are presented in Figure 3.In comparing the inclination angles included in Table 1, it is shown that the inclination angle of the lobe centers increases with increasing B ISM magnitude.With increasing B ISM magnitude, the magnetic pressure in the ISM increases, which further counteracts the effect of the solar magnetic field.For the cases of α BV = 40°ranging from 2.9 to 4.4 μG, we find that the inclination angle increases from 7°.1 ± 3°.6 for 2.9 μG (Case 9) to 10°.6 ± 3°. 5 for 3.2 μG (Case 6) to 12°.8 ± 2°. 5 for 4.4 μG (Case 10).An interesting aspect of this trend is the migration of the northern lobe, as seen in the bottom left panel of Figure 3.For Cases 6, 9, and 10, the lobe center location in the south maintains a relatively constant longitude, whereas in the north the lobe center migrates closer to the observed IBEX data.Coupled with the model results of changing α BV from the middle left panel of Figure 3, these results are consistent with the findings from Powell et al. (2024), who compared modeled Lyα profiles to observations and found that a lower modeled B ISM magnitude and a smaller α BV provided the best trends to match observations.Considering that we do not include helium in our MHD model, our solutions for Cases 1-10 do not take into account the effect of a higher ISM ram pressure resulting from the added mass density of helium.To investigate the effect of a higher ISM ram pressure due to the presence of helium, which is not accounted for in our modeling, in Case 11 we present the same conditions as in Case 6 (α BV = 40°, B ISM = 3.2 μG), but with n p = 0.09 cm −3 instead of n p = 0.06 cm −3 as an approximation.Besides helium ions in the ISM, there are also alpha particles in the solar wind, which (to some degree) balance the additional pressure from ISM helium ions.The effect of the solar wind and ISM ionized helium was studied in Izmodenov et al. (2003).We present our results in Figures 4 and 3, as well as in Table 1.As shown in the table, increasing the ISM ram pressure by 50% only has a slight effect on the inclination of the lobes.The inclination angle with a higher ISM ram pressure changes from 10°.6 ± 3°.5 for Case 6 to 11°.3 ± 3°. 7 for Case 11.This is also well represented in the bottom panels of Figure 3, where Cases 6 and 11 have nearly identical lobe centers, with the only difference being a slight shift in the southern lobe center northward.Therefore, the ISM ram pressure does not appear to have a strong effect on the inclination of the heliotail lobes.

Discussion
We present an alternative method to determining the orientation of the B ISM based on the IBEX ENA observations using the BU MHD model and the ENA model of Kornbleuth et al. (2023a).We find that our model results indicate that the high-latitude heliotail lobe rotation in longitude with respect to the solar meridian seen in IBEX observations and identified in Dayeh et al. (2022) is due primarily to changes in the solar magnetic field and not the solar wind fast/slow profile.Additionally, the inclination of the lobes is determined by the B ISM orientation, which, based on our model uncertainties, gives a range of 0°< α BV < 20°in the HDP.As was shown previously by Opher et al. (2007Opher et al. ( , 2009)), the heliosphere orients itself along the B ISM direction, which can be seen through the orientation of the high-latitude heliotail lobes of enhanced ENA flux at high energies (>2 keV) resulting from solar wind being deflected to the heliotail.The results here also demonstrate closer agreement to ENA observations of the heliotail from INCA, which indicated a vertical alignment of the heliotail structure (see, e.g., Kornbleuth et al. 2023a).In this work, we also considered B ISM strengths of 2.9 and 4.4 μG.Using a stronger B ISM strength (∼4 μG), the resistance by the solar magnetic field decreases, leading to a more inclined heliotail.As such, our results indicate that a low B ISM (∼3 μG) is required to match the IBEX ENA observations of the highlatitude heliotail lobe inclination and rotation (Dayeh et al. 2022).However, a low B ISM may contradict the magnetic field intensity observed in situ by Voyager 1 and 2 in the ISM.Further constraining of the B ISM orientation and magnitude is left to a future study.
One key limitation of this work is the use of steady-state MHD modeling instead of time-dependent MHD modeling in comparing to ENA observations that reflect changing solar wind conditions over time.In Section 3, we run cases both for solar minimum and maximum and for solar maximum with different solar magnetic field conditions.By testing these limits, we are able to place approximate constraints on the variation of the heliotail lobe inclinations over the course of the solar cycle.In particular, for Case 1 (i.e., 2008 solar minimum conditions) we find a lobe inclination angle of 14°.0 ± 3°.5, whereas for Case 2 (i.e., 2014 solar maximum conditions) we find a lobe inclination of 5°.0 ± 2°.5.We do not expect the observed high-latitude lobe inclinations to match either case exactly, even during periods where ENAs primarily originate from solar minimum/maximum conditions, respectively.Yet, considering the layering effect and the use of a heliotail flooded with solar minimum/maximum conditions, we can place approximate constraints such that, for a given BISM orientation and magnitude (for Cases 1-4, those of Opher et al. 2020, with α BV = 40°and B ISM = 3.2 μG), the lobe inclinations should likely fall within the solar minimum/maximum inclinations, progressing toward the solar minimum case during times reflecting solar minimum conditions and vice versa.Therefore, on average, we would expect the average of our two solar minimum/maximum steady-state cases to reflect the average of the observed lobe inclinations.
While our results serve as a useful tool for indicating the necessity for 0°< α BV < 20°in the HDP, a comparison of the ENA lobe positions with a time-dependent MHD model of the heliosphere is required for a more detailed analysis.Steadystate MHD models are unable to properly capture the latitudinal profile of the high-latitude lobes owing to the lack of fast/slow solar wind layering in the high-latitude tail, so a timedependent model of the heliosphere is required for a more detailed study.Additionally, for all modeled and observed lobe centers included in Table 1, the ecliptic longitudes are higher than the downwind ecliptic longitude of 75°.4.One interesting element that a time-dependent model of the heliosphere may reveal is the source of this systemic positive longitudinal shift of the lobe centers from downwind.
One important method for validation of the B ISM is a direct comparison of the modeled magnetic field magnitude and orientation along the Voyager trajectories with the in situ data.To date, no MHD model has been able to replicate the azimuthal and elevation angles of the observed B ISM and its bend away from the heliopause (Burlaga et al. 2022).Opher et al. (2017) argued that the lack of rotation of the magnetic field along the heliopause following the Voyager crossings is related to the draped B ISM twice reconnecting with the solar magnetic field in the flanks, creating a family of field lines with solar-like orientation ahead of the heliosphere at Voyager 1 and 2, as seen in observations.Turner et al. (2023) further argued that the Voyager probes are likely in a reconnecting boundary layer beyond the heliopause, which would account for the currently observed orientation of the B ISM .This reconnecting boundary region is further suggested by the 40-139 keV radial ion outflow measured to be perpendicular to the magnetic field by Voyager 2 for at least 28 au past the heliopause, which were suggested to be ions escaping from the heliosheath to the ISM (Dialynas et al. 2021).The α BV range presented in this work will need to be validated further with Voyager observations and MHD modeling, which has so far failed to predict the observed trends.
Another limitation of this study was the reliance of a B-V plane as the HDP.Lallement et al. (2005) suggested that the interstellar helium and hydrogen flow vectors constrain the direction of the magnetic field in the HDP.However, as mentioned in Section 1, both Izmodenov et al. (2005) and Pogorelov et al. (2008) demonstrated that an B ISM in the HDP is required to replicate the observed hydrogen deflection angle.For this study, we assumed that the B ISM is in the HDP.For the HDP, Lallement et al. (2010) found the angle between the B ISM and the ecliptic plane to be 52°± 7°.Gurnett et al. (2006) suggested that an angle between the B ISM and the ecliptic plane of 46°± 5°could explain the 2-3 kHz radio emissions observed by Voyager, while Opher et al. (2007) found that an angle of ∼30°-60°best explains the 2-3 kHz emissions.Opher et al. (2009) found that an angle between the B ISM and the ecliptic plane of ∼80°-90°best explains the heliosheath flows measured by Voyager.While most evidence points to the B ISM lying in the HDP, a more rigorous study of B ISM , which includes orientations out of the HDP, has the potential to provide further information on the exact orientation of the B ISM .The HDP itself has significant uncertainty because the plane is defined by two points (the interstellar He and H flow directions) relatively close to each other; therefore, the uncertainty in the HDP increases with distance away from these points on the order of ±10°.
in the HDP.While Zirnstein et al. (2016b) suggested B ISM = 2.94 μG based on model fitting to the IBEX ribbon, Opher et al. (2020) found B ISM = 3.20 μG to better fit with the Voyager 1 and 2 magnetic field data for their multi-ion, multifluid simulations (where PUIs and thermal ions are treated as separate fluids) of the heliosphere, as this presented a better comparison with Voyager magnetic field data.As in Opher et al. (2020), we use B ISM = 3.20 μG.

Figure 1 .
Figure1.Cartoon depicting a tail-centered Mollweide projection showing the inclination angle (ω) of the high-latitude heliotail lobes.The inclination angle is defined here as the angle of the north/south lobe centers relative to the solar meridian.λ and β correspond to the ecliptic longitude and latitudinal coordinates of each lobe center, respectively, with NL and SL denoting the north and south lobes, respectively.

Figure 3 .
Figure 3. Vectors connecting the north/south lobe centers from modeled ENAs and IBEX data in the heliotail extracted from the model/data (left panels) and shifted such that the midpoints of the vectors are centered on an ecliptic longitude and latitude of (0°, 0°) (right panels).Included is the vector for the IBEX data corresponding to the years 2012-2013 (black), 2014-2015 (orange), and 2016-2017 (turquoise) from Dayeh et al. (2022).Top: comparison with changing solar wind profile and solar magnetic field (Cases 1-4).Cases with latitudinally varying solar wind profiles are represented with solid lines, and cases with uniform solar wind profiles are represented with dashed lines.The arrows indicate change in lobe center over time.Arrow colors reflect the line color of the initial lobe center.Middle: comparison with changing B ISM orientation (Cases 5-8).Bottom: comparison of changing B ISM magnitude and n p,ISM (Cases 9-11).

Figure 4 .
Figure 4. Same as in Figure2, but with non-normalized ENA flux in units of (cm 2 s sr keV) −1 .Cases 5-11 all use the same solar wind profile (1995-2017 averaged) and solar magnetic field strength (B r,SW = 37.5 μG).For Cases 5-8, we vary α BV from 60°to 10°.For Cases 9 and 10 we replicate Case 6, but with a lower (B ISM = 2.9 μG) and higher (B ISM = 4.4 μG) interstellar magnetic field intensity, respectively.For Case 11, we also replicate Case 6, but with a higher interstellar plasma density (n p,ISM ) to investigate the effect of higher ISM ram pressure on the lobes.

Table 1
Effect of Solar and Interstellar Conditions on the High-latitude Heliotail Lobe Centers Dayeh et al. (2022).(2022)corresponding to ENA production conditions as during the solar minimum(2012-2013), transition period (2014-2015), and solar maximum (2016-2017), when the ENA propagation time in the heliosphere is accounted for.Listed model solar wind profiles correspond to temporally averaged observations at 1 au during the specified time.Uniform corresponds to identical plasma conditions for all latitudes and longitudes at 1 au.