Influence of Solar Disturbances on Galactic Cosmic Rays in the Solar Wind, Heliosheath, and Local Interstellar Medium: Advanced Composition Explorer, New Horizons, and Voyager Observations

We augment the heliospheric network of galactic cosmic ray (GCR) monitors using 2012–2017 penetrating radiation measurements from the New Horizons (NH) Pluto Energetic Particle Spectrometer Science Investigation (PEPSSI), obtaining intensities of ≳75 MeV particles. The new, predominantly GCR observations provide critical links between the Sun and Voyager 2 and Voyager 1 (V2 and V1), in the heliosheath and local interstellar medium (LISM), respectively. We provide NH, Advanced Composition Explorer (ACE), V2, and V1 GCR observations, using them to track solar cycle variations and short-term Forbush decreases from the Sun to the LISM, and to examine the interaction that results in the surprising, previously reported V1 LISM anisotropy episodes. To investigate these episodes and the hitherto unexplained lagging of associated in situ shock features at V1, propagating disturbances seen at ACE, NH, and V2 were compared to V1. We conclude that the region where LISM magnetic field lines drape around the heliopause is likely critical for communicating solar disturbance signals upstream of the heliosheath to V1. We propose that the anisotropy-causing physical process that suppresses intensities at ∼90° pitch angles relies on GCRs escaping from a single compression in the draping region, not on GCRs trapped between two compressions. We also show that NH suprathermal and energetic particle data from PEPSSI are consistent with the interpretation that traveling shocks and corotating interaction region (CIR) remnants can be distinguished by the existence or lack of Forbush decreases, respectively, because turbulent magnetic fields at local shocks inhibit GCR transport while older CIR structures reaching the outer heliosphere do not.


Introduction
Since their discovery by Victor Hess (1912) over a century ago, cosmic rays have been recognized as a fundamental component of Earth's external radiation environment. In contemporary usage, the term "cosmic ray" is usually limited to two components. The first is galactic cosmic rays (GCRs): high-energy (100 MeV nuc −1 ) charged particles originating far outside the heliosphere (e.g., particles accelerated at plasma shocks associated with supernovae or at strong shocks formed by stellar wind collisions; see Hamaguchi et al. (2018)). The second population is anomalous cosmic rays (ACRs): less energetic (∼10-100 MeV nuc −1 ) charged particles accelerated within or in the vicinity of the heliosphere (Klecker et al. 1998). The longest-standing interpretation centers on the heliospheric termination shock (TS) (e.g., ), but ACRs can also potentially be accelerated in the heliosheath (HS), for example via magnetic reconnection (e.g., Drake et al. 2010;Zank et al. 2015;). These particles are influenced by their source populations (interstellar pickup ions) as well as their transport to and through regions local to the observer. Zank et al. (2018) and  have recently reported on their theoretical progress and simulations of pickup ions, demonstrating the importance of pickup ions in mediating the solar wind and finding agreement with NH and V2 plasma observations.
A useful quantity in this context is rigidity, P, particle momentum divided by charge, which amounts to a magnetically normalized gyroradius. Lower-rigidity particles are more easily influenced by local conditions, where the influence is usually observed as a reduction in the count rate signal (typically proportional to differential intensity). Historically, such reductions have been called modulation. Although the term "modulation" should strictly apply to the measured signal, in the literature it is often used to label the process where cosmic ray transport is disrupted, resulting in a reduced signal. One of the most characteristic features of cosmic ray measurements is the Forbush decrease (Forbush 1938;Burlaga et al. 1993), which was first discovered on the ground and later in space. A Forbush decrease occurs when an abruptly reduced GCR intensity is temporally associated with local interplanetary events, such as traveling shock waves, interplanetary coronal mass ejections, corotating high-speed solar wind streams, and globally merged interaction regions (GMIRs). However, while the ultimate source of the GCR variation is the Sun, the challenge of connecting specific events on the Sun to features in the particle measurements is well known. For example, Chowdhury et al. (2013), writes, "This is a very complex phenomenon that occurs throughout the heliosphere and depends on several factors. No single solar parameter can account for the GCR intensity." Likewise, Zhao et al. (2014) have found through simulation that some potential transport parameters do not have an effect on GCR intensity. Our paper focuses on the relation between interplanetary events and GCRs. Because there is ambiguity between the physical processes listed above and their signatures, we will refer to either of them as "disturbances" or "variations" unless stated otherwise, e.g., when we attempt to resolve some of this ambiguity (Section 4.1).
The local activity is associated with magnetic field turbulence on a spatial and temporal scale that results in more scattering of particles in a rigidity-dependent manner; lowerrigidity particles are more easily modulated than higher-rigidity particles. When modeled as a diffusion process, the turbulence is said to produce scattering centers arising from magnetic field irregularities. The study of cosmic ray diffusion shows that such rigidity dependence is seen prominently in parallel diffusion, but only weakly in perpendicular diffusion, suggesting that the relative role of parallel and perpendicular diffusion could be probed by observations, such as we provide in this report, of GCRs at different energies. In particular, the rigidity dependence of the parallel mean free path has been found to be proportional to~P 0.33 from 10 to 1000 MV, but in the distant heliosphere, its dependence increases for higher rigidities and is proportional to ∼P 2 . In contrast, the perpendicular mean free path is only weakly influenced by cosmic ray rigidity (Zank et al. 1998Zhao et al. 2017Zhao et al. , 2018. We are interested in the propagation of disturbances into the outer heliosphere as seen in GCRs, and therefore consider ACE at 1 au to be the inner boundary of our region of interest since we have no GCR measurements closer to the Sun during our period of interest. Only five spacecraft are on escape trajectories from the solar system: Pioneer 10 (P10), Pioneer 11 (P11), Voyager 2 (V2), Voyager 1 (V1), and New Horizons (NH), launched 1972March 3, 1973April 6, 1977August 10, 1977September 5, and 2006 January 19, respectively. In Figure 1, we show the latitudinal and longitudinal positions of these five escaping spacecraft at a selection of times; only P10 is heading down the heliotail, away from the initial interaction with the interstellar medium, whereas the remaining spacecraft are clustered within ∼45°of heliographic inertial longitude of one another. During the period of interest of this study (2012)(2013)(2014)(2015)(2016)(2017), NH ranges from 22 to 40 au, V2 ranges from 96 to 116 au, and V1 ranges from 119 to 141 au. P10 and P11 were at 79 au and 59 au, respectively, when they ceased making observations prior to the time period of interest. P11, V2, V1, and NH are all heading roughly in the direction of the heliospheric nose-that is, the direction of relative motion between the heliosphere and the local interstellar medium (LISM). Note that the LISM includes the region outside the heliopause (HP) but inside of the heliospheric bow wave that is sometimes called the "outer heliosheath" (to differentiate it from the "inner heliosheath" between the TS and HP, which we just call the HS). The P10 and P11 missions have ended, and their calibrated cosmic ray measurements (Teegarden et al. 1973) extend until 1994. There are P10 and P11 Geiger tube data out to 2002, respectively (Van Allen & Randall 2005. Because this study focuses on the time period in which NH has cosmic ray measurements (2012)(2013)(2014)(2015)(2016)(2017), P10 and P11 cannot be used as a point of comparison.
Both V1 ) and V2 have passed through the HP into the nearby LISM (Holzer 1989); therefore, the only operational heliospheric spacecraft remaining beyond the orbit of Jupiter is NH. Effectively, there is no other heliospheric spacecraft outside of Mars because the Juno Jupiter orbiterdespite being outfitted with "Puck" instruments (Mauk et al. 2013) very similar to the Pluto Energetic Particle Spectrometer Science Investigation (PEPSSI; McNutt et al. 2008)-is rarely, if ever, outside Jupiter's massive magnetosphere. While in orbit around Saturn, instruments on the Cassini mission, which ended in 2017, were able to occasionally detect Forbush decreases and solar events (Roussos et al. 2018); however, our study is focused on comparing GCR intensities measured in the solar wind, away from competing sources of energetic particle observations (e.g., radiation belts). Although it is close to the asymptotic heliocentric longitude of V2, NH is the only spacecraft near the plane of the ecliptic in the direction from which the neutral interstellar wind enters the heliosphere ( Figure 1). Thus, NH is in a region of the outer heliosphere being explored for the first time; see Stern et al. (2018) for more information about the NH Kuiper Belt Extended Mission, of which the present study is a part.
In situ cosmic ray measurements began in the LISM on 2012 August 25, when V1 crossed the heliopause Cummings et al. 2016). One of the key surprises was the detection of time-varying anisotropies in >200 MeV GCRs . GCRs are generally isotropic in the heliosphere and there was no expectation that this would be less so in the LISM. What was found was a second-order anisotropy with intensities depressed in the direction perpendicular to the average magnetic field direction. The first quantitative, physical interpretation of the anisotropy was given by Roelof et al. (2013): "[T]ime-dependent depletions in the intensities localized in pitch angle near 90°...would be observed whenever [Voyager 1] is located between (at least) two compressions of the magnetic field on a field line, both (distant) ends of which contain equal isotropic GCR intensities....These compressions would, by their very nature, be time dependent, thus explaining the time dependence of our observed [pitch angle distribution] anisotropies." This twocompression, trapped configuration view has informed subsequent investigations of the V1 GCR anisotropies, including theoretical work by Kóta & Jokipii (2017) and a recent observational study by Rankin et al. (2019), with considerable success. The anisotropic variations have been compared with in situ interstellar shocks and radio measurements of plasma oscillations (Gurnett et al. 2015), but the unusual timing of the cosmic ray anisotropy periods relative to the magnetic field, plasma waves, and locally accelerated particles calls for an explanation. For example, the cosmic ray anisotropy periods may not always correspond with the in situ magnetic field shock signatures (see Section 4.2). To date, there has been no description of the heliospheric interaction with the LISM that comprehensively explains the seemingly inconsistent relationships between the various measurements observed at V1. We present an alternative escaping particle picture in which a disturbance interacting on only one side of an LISM field line results in the observed time-dependent anisotropies. This configuration is required to observationally reconcile the GCR intensity variations observed inside the heliosphere at ACE, NH, and V2 with the V1 interstellar GCR anisotropies.
Although the core of this report is the presentation of new observations from NH, we are also putting forth a framework for interpreting the data that resolves some of the newly established constraints. A large part of the associated analysis involves comparison of GCR intensity in the context of disturbances of solar origin propagating through the heliosphere. This approach has been taken before by many scientists. Observationally, multiple spacecraft have been used to track solar events through the heliosphere (see, e.g., McDonald et al. 1981; Van Allen & Fillius 1992;Webber & Lockwood 1993;Cane et al. 1994;Paulerena et al. 2001;Richardson et al. 2002Richardson et al. , 2005Burlaga et al. 2003aBurlaga et al. , 2003bWitasse et al. 2017, and references therein). Measurements have been used to drive simulations of transient events (see, e.g., Wang & Richardson 2001Luo et al. 2011;Pogorelov et al. 2012;Shen & Qin 2018, and references therein). Observations and modeling have been used to study the effect of traveling disturbances on the LISM (see, e.g., Gurnett et al. 1993;Steinolfson & Gurnett 1995;Wang & Belcher 1999;Zank & Müller 2003;Webber et al. 2007Webber et al. , 2009Richardson et al. 2017;Schwadron & McComas 2017, and references therein). Variability at the termination shock and heliopause due to long (solar-cycle) and short (transient) timescales have been studied with simulations and observations as well (see, e.g., Izmodenov et al. 2005;Washimi et al. 2007Washimi et al. , 2011Washimi et al. , 2015Washimi et al. , 2017McComas et al. 2018;Zirnstein et al. 2018;Burlaga et al. 2019, and references therein). It is our hope that others will compare these newly prepared data sets with sophisticated simulations and theoretical predictions; our interpretation attempts to assimilate the broad interconnected system of observations, provides a potential avenue for subsequent investigation, and highlights some of the most essential constraints asserted by new observations.

New Cosmic Ray Monitor: NH/PEPSSI
Here, we present newly analyzed cosmic ray measurements made with the solid-state detector (SSD) system of the PEPSSI instrument on the NH spacecraft (McNutt et al. 2008). The new channels correspond to penetrating ions having energies of 75-120 MeV, 120 MeV-1.4 GeV, and ∼1.4-5 GeV. These channels were used to study Jovian electrons by Haggerty et al. (2009), but in the absence of <1 MeV electrons, they are dominated by cosmic rays. Further specifications of these channels and experimental considerations can be found in Appendix A. An overview of the PEPSSI particle measurements, spanning six orders of magnitude in energy, from ∼2.0 keV nuc −1 to 5 GeV, is provided in Figure 2 for the 2012-2017 period. During this period, the PEPSSI instrument was operating nearly continuously; poweroff periods ranged from several hours to several days and occurred during spacecraft activities such as trajectory correction maneuvers, transitions from spin-stabilization to three-axis stabilization, and reductions in the power load to permit other spacecraft operations.  P10 was launched in 1972P10 was launched in , P11 in 1973P10 was launched in , V1 and V2 in 1977P10 was launched in , and NH in 2006P10 was launched in . P10 and P11 ceased regular data transmissions in 2003 au from the Sun) and 1995, respectively. V1 exited the heliosphere in 2015, and V2, at comparable heliolongitude as NH, exited the heliosphere at the end of 2018, making NH the lone operating spacecraft in the heliosphere outside the orbit of Jupiter.
The pickup and suprathermal ion and energetic particle measurements are described in more detail elsewhere (McNutt et al. 2008;Hill et al. 2009;Kollmann et al. 2019a) and in Appendix A.2. The high-energy particle intensities (Figures 2(a)-(c)) show a minimum around the time of solar maximum and have short-term behavior that is distinct from that of suprathermal ions and energetic particles observed at the same time (Figures 2(d) and (e)). Cosmic rays are known to be anticorrelated with solar cycle, sunspot number (see Figure 3(a)) being a common proxy (e.g., Forbush 1958;Parker 1965;Usoskin et al. 1998;Van Allen 2000), and show Forbush decreases when transient variations in lowenergy intensities occur, as we observe here. Particle activity evidenced by intensity increases in the suprathermal or energetic particles is highlighted in Figure 2, with 30 blue vertical lines for cases in which the cosmic ray intensities do not simultaneously decrease, while the 31 orange lines are for cases in which there is a coincident decrease in the GCR intensity, indicating the existence of a Forbush decrease (see also Section 4.1). However, in principle, these channels respond to other populations as well. One candidate population would be solar energetic particles (SEPs), either ions or electrons. If there were a large event, the spectrum below 1 MeV (observable by PEPSSI) would reveal the heliospheric origin of the event because of the high intensities that would result from the characteristic inverse relationship between intensity and energy. Moreover, PEPSSI measures the particles in the 1-75 MeV energy range, and we observe that there is no significant contribution (less than ∼2%) from such SEP foreground (see Appendix A.2). Once we exclude SEPs and other local populations, we are left with a measurement dominated by cosmic rays. GCRs have been measured throughout the heliosphere, so above ∼100 MeV, it is the only known population that can be appreciably contributing to our measurement. Below this energy, however, there could be a contribution from ACRs, which have been shown to vary more readily with changes in transport conditions (Cummings et al. 1995;Hill et al. 2003) and to have peak intensities in the ∼30-100 MeV nuc −1 energy range, depending on global transport conditions, such as large-scale drifts and magnetic field turbulence (see, e.g., Potgieter 1993). Therefore, our lowest-energy channel could contain both ACRs and GCRs. However, ACRs and GCRs both respond similarly to interplanetary disturbances and we are interested in the relative variations in the cosmic ray to identify telltale features of transient disturbances, so the specific scale of the variations is not central; in fact, we specifically normalize out much of the variation before we compare the features between the propagation-shifted GCR time series. As such, the propagation timing would not be expected to be significantly impacted by a convolution of ACR and GCR counts. High-energy (>75 MeV) electron contribution to the PEPSSI GCR measurements of from ∼4% to 30% is discussed in Appendix A.2 and indicated in Figure 2.
In the Discussion (Section 4.1), we briefly discuss the relationship between solar disturbances in the form of energetic and suprathermal ions and the cosmic ray intensity modulation, in particular short-term Forbush decreases and long-term solar cycle variations. The energy-dependent recovery following a Forbush decrease has been used to distinguish types of solar disturbances (Zhao & Zhang 2016). The cosmic ray intensities observed at NH vary in an energy-dependent manner, with smaller intensity variations seen in the highest-energy, ∼1-5 GeV, channel. Without compositional measurements, we cannot use NH measurements to directly investigate the rigidity dependence of cosmic ray diffusion beyond 1 au mentioned in Section 1, but predictions could be tested under the assumption that protons dominate the intensity variation. . Six years of PEPSSI particle observations from ∼2 keV nuc −1 to ∼5 GeV. Top three grouped panels, (a), (b), and (c), show the intensities of penetrating particles having energies from 75 to 120 MeV, 120 MeV to 1.4 GeV, and 1.4 to 5 GeV. These predominantly GCR protons show the characteristic hallmarks of solar cycle activity, with the local solar cycle #24 maximum (when cosmic rays intensities reach their minimum) occurring in 2015 (see Figure 3) and the recovery period clearly underway by the beginning of 2016. Bottom two grouped panels show lower-energy, local populations, with (d) dominated by helium, mostly pickup ions and suprathermal particles (Kollmann et al. 2019a), down to 2 keV nuc −1 , and (e) showing energetic protons (from 20 to 800 keV). It can be seen that the intensities of these lower-energy populations show no obvious solar cycle dependence but more short-term changes. Vertical lines (in the bottom panels) mark activity observed in energetic and suprathermal particles; blue lines do not correspond with clear Forbush decreases but orange lines (shown in the top three panels too) do. Systematic uncertainties, represented by downward-pointing arrows on the right of the GCR panels, are dominated by potential GCR electrons (see Appendix A.2).

ACE and Voyager Cosmic Ray Measurements
With the addition of the NH cosmic ray monitoring capability to those of the enduring ACE (Stone et al. 1998a) and Voyager (e.g., Stone 2001) missions, we have the opportunity, for the first time, to directly link cosmic ray observations from the inner heliosphere, outer heliosphere, heliosheath, and LISM simultaneously. Whereas Cummings et al. (2016) reported on LISM cosmic ray spectra at V1, including a comparison with 1 au measurements and transport models, we are focusing on linking the modulation of GCR observations at different heliographic radial distances to study how these variations in GCR intensity propagate from the Sun into interstellar space. We have assembled these cosmic ray observations from across the heliosphere for comparison with the PEPSSI measurements, permitting a global examination of cosmic rays.
To anchor the cosmic ray measurement close to the Sun (meaning 1 au, which is near the Sun for the large scales in this study), we use observations from ACE's Cosmic Ray Isotope Spectrometer (CRIS) instrument, which was designed to precisely measure cosmic ray intensities and isotopic composition for heavy ions, primarily Be to Fe, using SSD stacks and a fiber-optic hodoscope (Stone et al. 1998b). An integral >120 MeV proton rate, known as E9, was developed (Mewaldt et al. 2010) using the bottom detectors in CRIS's four solidstate detector (where each nine-element stack is equal to 14.5 mm of silicon), which are housed behind the scintillating-fiber trajectory system. The E9 rate has since been greatly improved with better background corrections and quiet-time selection criteria, and it has very good statistics, averaging between ∼200 and 400 counts per second depending on the phase of the solar cycle, making it an excellent detector of predominantly GCRs near the Sun.
The V1 Figure 6. Statistical error bars are smaller than the symbol sizes, and typical numbers of counts per daily data point are listed on the right below the channel name. Spacecraft, instrument, ion energy range, and approximate helioradius range are provided for reference. Upper (solid, R UpEnv ) and lower (dashed, R LoEnv ) envelopes used in the propagation analysis are shown. (c) Variations between the lower and upper envelope of each time series are normalized to range from 0 to 1 using , and then additive integer offsets (+2, +1.25, +0.6, +0 for ACE, NH, V2, and V1 respectively) are applied for convenience of display. Horizontal axis is the time at which the given disturbance is detected at NH, determined by the propagation analysis. Intervals have been highlighted alternately in pale yellow and pale blue to highlight decreasing and increasing features, respectively, and pale gray and pale pink to highlight peak and valley features, respectively. Black arrows at the bottom of panel (b) and the top of panel (c) represent the dominant variation (increasing, decreasing, peak, or valley) during each interval. Some of the features are subtle and some pronounced, but most are observed at multiple spacecraft at roughly the same propagation reference time (PRT), which aligns with the NH reference frame. Above the panels, the letters from A to N (sans I) label the highlighted regions/features. The 13 labeled highlighted regions are dramatically distorted in panel (b), where no propagation correction is done, but are roughly aligned in panel (c).
spacecraft were designed to measure cosmic rays: the Cosmic Ray Subsystem (CRS; Stone et al. 1977) and the Low Energy Charged Particle (LECP) instrument (Krimigis et al. 1977). Here, we use LECP with the >210 MeV proton rate and continuous 360°scanning, which is achieved with an eightsector stepper motor that pivots 45°every 192 s. The relevant LECP subsystem is an SSD stack composed of five detector layers. The SSD telescope is double-ended, and the cosmic ray channel, known as CH31, derives from event logic on energy signals from two 2450 μm thick SSDs; the event logic selects primarily minimum ionizing protons-that is, the high-energy penetrating ions capable of reaching these detectors by penetrating the instrument body and spacecraft boom structures (note that 200 MeV H can penetrate 10 cm of aluminum and that even a solid meter of aluminum is insufficient to stop all 1 GeV protons, which are numerous). In 2012, we used CH31 to first observe the unexpected cosmic ray anisotropies in the LISM at V1 ; these anisotropies correspond to depressed intensities near 90°pitch angles.

ACE, NH, and V2 Solar Wind Speeds
As discussed in Section 4.1, cosmic rays tend to display Forbush decreases related to local, transient solar-initiated disturbances and longer-term solar cycle intensity changes (often termed modulation). Therefore, to study what is happening to the cosmic rays globally, we need to understand when these transient disturbances reach the spacecraft so we can quantitatively compare the measurements at the same propagation-corrected time. We chose to use the NH time reference for all the comparative analysis herein, using the solar wind speeds measured at ACE, NH, and V2 to make these propagation corrections (see Section 3.2 and Appendix B).
The ACE solar wind speeds V swACE measured by the Solar Wind Electron Proton Alpha Monitor (SWEPAM; McComas et al. 1998), an angularly sensitive electrostatic analyzer, have been used to propagate disturbances, which travel at the supersonic solar wind speed, from ACE to NH. Specifically, we used the daily averaged proton bulk speed obtained from the ACE Science Center. 14 For disturbance propagation beyond NH, moment calculations  from the measurements made by the Solar Wind Around Pluto (SWAP) electrostatic analyzer-based instrument  were used to yield the required solar wind speed V swNH in the outer heliosphere.
For the heliosheath propagation speeds from the TS to V2, V HS2 , we employ the V2 Plasma Science (PLS) instrument (Bridge et al. 1977), which is composed of four Faraday cups enabling directional flow determination. Because of damage during the Jovian encounter, the V1 PLS instrument ceased working in 1980 and is unable to provide plasma measurements in the HS or LISM. Also, in the LISM, the energetic particle densities are too low to enable the Compton-Getting-based speed determination that has been carried out using LECP observations (Decker et al. 2005), so an alternate technique is required for determining propagation speeds to V1 (see Appendix B.2).

Largest Array of Cosmic Ray Monitors
Cosmic rays have been used to study the global heliosphere for many years, but temporal comparisons have never before been conducted for such a large in situ measurement array encompassing such a diversity of regions. At lower energies, the intensities are modulated in response to local and global changes in the magnetic field turbulence spectrum and global drift patterns arising from the 22 yr heliomagnetic polarity cycle (e.g., Potgieter 1993). The intent of the analysis detailed in Section 3.2 is to align the measurements so that, at a given NH propagation reference time (PRT), the same spheroidal shell of plasma is passing all the spacecraft except V1 (because it is outside the heliosphere). This spheroidal shell picture can only be approximately true, given the fact that the solar wind speed has latitudinal, longitudinal, and temporal dependencies. Furthermore, the solar wind is known to slow down as a result of the transfer of energy to interstellar pickup ions (e.g., Richardson et al. 2008;McComas et al. 2017), and the position of the TS also varies with time (e.g., Barnes 1993;Ratkiewicz et al. 1996), further complicating the picture. Moreover, there are stream interactions at the junctions between speed changes that will affect the resulting flow speed and properties of the plasma, and the propagation of solar wind is not a perfect proxy for propagation of transport conditions; for example, global heliospheric drifts are not accounted for. Nonetheless, this technique permits a comparison of measurements across ∼100 au so that similarities and differences can be examined to study the physical processes that control the behavior of cosmic rays throughout the heliosphere.

Propagated Cosmic Ray Observations
In Figure 3(b), the linearly adjusted "raw" cosmic ray counting rates at each spacecraft are plotted at the observed time (see the figure caption for details of the vertical adjustment). Envelopes were defined for ACE, NH, and V2 by conducting a year-long boxcar average of the rate data to capture the long-term variations (e.g., solar-cycle variations) and then shifting the curve above and below until the time profile of the GCR rates were predominantly bounded. Since the V1 LISM GCR rate is known to have no long-term variation (Cummings et al. 2016), the V1 rates were enveloped using constant bounding values, defined by the maximum and minimum daily rates taken measured in the LISM. These envelopes show that the two heliospheric GCR observatories (ACE and NH) observe solar cycle variations (with the minimum of the GCR envelopes occurring at a time-lagged solar cycle maximum; see discussion associated with Figure 3 in Section 4.1). Meanwhile, the V2 GCR envelope is not observed to have a predominant solar cycle variation. This would suggest that the solar cycle variation, which is largely due to changes in solar magnetic topology, does not strongly affect the heliosheath region. The GCR rate observed by V1 in the LISM, as expected, does not observe long-term variations driven by the heliosphere (see Gurnett et al. 2015;Cummings et al. 2016).
As seen in Figure 3(b), transient variations (e.g., Forbush decreases) in the GCR rates (indicated by letters) are seen in each of the observatories as solar wind disturbances, such as shocks and corotating interaction regions (CIRs), propagate away from the Sun. To better understand how these disturbances propagate to the outer heliosphere, we detrend the observations to allow for a comparison of only these transient modulations of the GCR intensity, and shift the observation times to a uniform PRT frame. This PRT is defined such that the NH observations are unshifted. Observations from ACE and V2 are propagated using simple ballistic propagation (i.e., using only the propagation speed and the radial separation of the spacecraft, without additional effects included, such as slowing due to mass loading or stream interactions; see Appendix B.1 for more information on this propagation) as done in other studies (e.g., Elliott et al. 2016). This simple ballistic approach does fairly well in aligning the identified features (Figure 3(c)). For the purposes of this paper, the alignment is only required to be approximate, rather than precise, and so this approach is seen to be sufficient. More sophisticated MHD propagation schemes are beyond the scope of this study. While the ballistic propagation scheme is generally well-constrained, observationally, from ACE to V2, attempts to use this approach for the V1 observations are less well-constrained due to the absence of velocity measurements in the LISM. In applying the same ballistic approach for V1 using a range of velocities (17-1000 km s −1 , or the minimum V1 speed to well beyond the nominal Alfvén velocity of 40 km s −1 in the LISM from Kóta & Jokipii (2017), not shown), no velocity was found to align the features observed in the GCR rates of V1 to the other spacecraft observations. Additionally, no ballistic propagation in this simple manner can explain why features occurring after 2015 appear to be observed by V1 before V2 (see Figure 3(b)). To overcome this, we resorted to a linear propagation between V1 and V2 in which (see Section 3.3 for more information). The result of this propagation for V1 leads to a good alignment of the GCR features between V1 and V2 ( Figure 3(c)).
Undertaking comparisons of temporal variations across multiple in situ spacecraft measurements is notoriously challenging, and has been done here over a longer spatial baseline than has previously been attempted and with spacecraft operating in a wider range of environments (i.e., inner heliosphere, outer heliosphere, heliosheath, and LISM). Elliott et al. (2016Elliott et al. ( , 2018 and Kim et al. (2016) did such comparisons with plasma measurements and modeling, with success, but did not include observations from the heliosheath or LISM. At the outset, our aim was to make this comparison and identify sufficient consistency to have confidence that we understand the propagation throughout the system. Now, we go through the detailed features in the data to establish this consistency.
The features labeled with an A and B seen by ACE at 1 au are a set of smaller Forbush decreases, forming a ∼1 yr depression, terminating with a recovery (marked by a C). These A-B-C variations can be seen in both panels (b) and (c) of Figure 3. The GCR observations at NH began at the end of the B period, and as such there is no NH data for most of the A-B interval, but the C recovery is captured. At V1 and V2, the A-B-C features are very similar to one another, with distinct U-shaped decreases in A-B and smaller depression in B-C (most apparent in Figure 3(c)). This is evidently a merging of the smaller Forbush decreases seen at 1 au. The D-E features describe an inverted V form in the V2 data (Figure 3(c)). This D-E feature tracks from 1 au, where the number of narrow ∼26 days peaks seen at ACE appears to be reduced by the time the disturbance reaches NH (data gaps prevent a conclusive count) and reduced even further when V2 is reached, evidenced by a W shape above the V1 V shape. At ACE and NH, the F trough is seen to have some peaks near the minimum. At V2, there are indications of this feature, but it is not identifiable at V1.
Feature G, taken alone, is a large peak easily seen in either the raw (Figure 3(b)) or normalized (Figure 3(c)) representation at ACE, NH, and V2, but there is no increase at V1 (a sharp kink and small peak toward the end of period G could possibly be associated with these large features, but it is unclear). Following the large G peak is the most distinctive feature, a large Forbush decrease and recovery designated H-J, which is well-defined at all four spacecraft, tying the LISM decrease back to 1 au in the most unambiguous manner.
The decrease and recovery K-L is traceable through all four spacecraft, but is least distinct at V1. The trough, labeled M, is easily observable at ACE, NH, and V2, but is very shallow at V1. The final feature tracked, N, is a falling trend with significant narrow peaks at ACE, with the lowest rate (Figure 3(b)) at the end of the period. See Section 4.1 for a discussion of solar cycle maximum timing.
The results of the preceding analysis show that we are identifying the propagating disturbances properly with the time shifts between ACE, NH, and V2 arising from simple ballistic propagation (see Appendix B.1). It should be emphasized that the expectation is not that these time series would align with one another in most cases; there are radial variations in the solar wind structures (e.g., GMIRS) that contribute to a different response from the GCRs. Moreover, NH, V2, and V1 are all widely spaced in latitude (Figure 1), which would be expected to result in sizable differences since features on the Sun associated with propagating disturbances, such as coronal holes, depend on latitude. Considering this, the good agreement between ACE, NH, and V2 provides additional confidence that the new NH PEPSSI channels are indeed predominantly measurements of GCRs. As such, the addition of NH to the network of heliospheric cosmic ray monitors is invaluable for tying the inner and outer heliospheric observations together because of the evolution of disturbance features over great distances.

Propagation of Solar Disturbances
In Section 3.2, we have taken four roughly comparable cosmic ray counting rate profiles and detrended and timeshifted them to the NH-based PRT. It is important to point out that the ballistic propagation method works well only within the heliosphere and heliosheath, and we found that ballistic propagation cannot be used out to V1 to explain the GCR anisotropies there. A linear propagation scheme is used to compute the time lag in transient features observed in the V1 GCR rate to the other observations in PRT (Section 3.2). In this section, we seek to gain a better physical understanding of the implications of the linear propagation scheme for the V1 GCR features through the development of a more generalized set of equations to describe the propagation from within the heliosphere to V2 in the heliosheath and then to V1 in the LISM (Appendix B.2).
We rely only on piecewise linear equations to analytically model the propagation. The constituent steps in building up this propagation scheme are simple relations between propagation speed, distance, and time. Elliott et al. (2016) conducted such propagation from 1 au to New Horizons using a similar method. The following concepts guided the development of the parameterized propagation.
(1) We follow "signals" (i.e., the transient variations in the GCR rates) in a propagation distance versus time space as they are observed at spacecraft or change at boundaries (e.g., the TS and the HP). Spacecraft trajectories and boundaries are defined by linear equations of time. For spacecraft positions, these equations are linear approximations of the trajectories expressed as radius as a function of time, with a speed and offset parameter. There are two subtleties here to keep in mind (detailed in Appendix B.2). First, boundaries are parameterized by time, but this does not imply physical time dependence. Second, the distance versus time plane can be roughly conceptualized as helioradius versus time, but not strictly so.
(2) We approximate the solar disturbances as one-dimensional (i.e., spherically symmetric) until they reach the TS (the inner boundary of the heliosheath), which itself is not necessarily spherically symmetric (e.g., McComas & Schwadron 2006). This approximation is supported by the success of the simple ballistic propagation used to shift the ACE and V2 GCR observations to the PRT. While the solar wind speed in the heliosheath is known to be time-varying (e.g., Krimigis et al. 2011;Richardson 2011), we use a nominal speed value for these computations. As such, the model can be thought of as a means to understand the nominal conditions due to the general relationship of values. To further generalize the system of equations, we do allow for systematic differences in the heliosheath flow speeds along the path to V2 and along the path of a "signal" propagating to V1. This is due to possible longitudinal and latitudinal variations in the nominal flow speed in the heliosheath, similar to the magnetosheaths of planets. For these calculations, the solar wind speed prior to reaching the TS is assumed to be a nominal 400 km s −1 .
(3) The point along the solar disturbance path at the TS is not required to be along the trajectory of V1 or V2 where those spacecraft crossed the TS shock. We use the term subcontact point to describe the radial location of where the signal reaches the TS. This subcontact point is left as a free parameter in the formation of the generalized propagation scheme, and will be further explored later in this section. This analytical scheme allows for the fact that the TS could be nonspherical, so the distance to the TS along the disturbance path R TS can change as a function of latitude and longitude, which we parameterize by time (see Appendix B.2), although latitude and longitude are not explicitly included in the formulation. (4) To further generalize the model, the path in the LISM that a "signal" would take to reach V1 is not predetermined. This allows for both the possibility that a solaroriginating disturbance reaches the location of V1 to modulate the observed GCR anisotropy locally, as well as for the disturbance to affect V1 observations remotely through changing the magnetic field-aligned energetic particle populations magnetically connected to V1. For example, for every LISM magnetic field line threading the V1 position at a given point along its trajectory, there is a point on the surface of the heliopause that is "closest" to that field line (assuming a nonpathological HP surface, such as one that is locally concave). This point is near the field line draping region around the heliosphere (e.g., McComas et al. 2009;Opher et al. 2017;Schwadron & McComas 2017). The system of equations constructed here allows for both an analysis of the "signal" propagating to V1, as well as to the closest field line connected to V1. Both of these possibilities are further explored later in the section.
We sought the simplest generalized representation of the propagation correction between V2 to V1 and found that a linear relationship was sufficient, namely (see Appendix B.2 for more information). We rewrite the equation for t V1 by subtracting the time of origin t o from both sides as follows: Using this linear scheme, the coefficients were determined empirically to best align the features observed in Figure 3(b). This alignment was determined by optimizing the Pearson linear correlation coefficient between the V2 and V1 deenveloped GCR time series while varying the slope m and temporal offset D ¢ t (the reduced chi-squared measure was checked for consistency, but the Pearson coefficient showed more readily identifiable peaks). The slope was varied from 0.5 to 1.0 in steps of 0.0005 and the temporal offset was varied over±3 yr in half-day steps. (see Section 3.2, Appendix B.2, and Figure 3(c)). We found three local maxima of the correlation coefficient, around which we searched for the absolute maximum. Focused parameter searches were then conducted at the cluster of three peaks. We used the coordinates of the three peaks to estimate an uncertainty, arriving at m=0.8160±0.0035 and t o =2015.136±0.030.
This empirically determined relationship provides-under the analytical, linear description of the V1-to-V2 propagation (Appendix B.2)-the time of arrival at V1 (t V1 ) of a propagating structure observed at V2 at time t V2 . The second equality in Equation (1) is a definition from which it follows that, when the time origin is t o =Δt′/(1 − m), the offset is Δt=0. This implies that t o is the time when a disturbance at V2 is seen simultaneously at V1 (i.e., feature L in Figure 3). These values can now be compared with the physical parameters of the generalized propagation scheme (see Appendix B.2). Thus, the slope and intercept parameters in Equation (1) are:  . Diagram (not to scale) shows propagation of disturbances from 1 au, through NH, into the heliosheath, and into the LISM where the signal travels to V1 by two routes. Further details are described at length in the text (see Section 3.2), but the crux is to note that this is a simple one-dimensional picture and threedimensional effects are handled by two bifurcations when an interface (like the TS) exists only for one signal and not the other. Figure 9 provides a simple version of this diagram for the case of one signal propagating between two trajectories. This simple version is used in Appendix B.2, where the analytic part of the propagation correction from the heliosheath to interstellar medium is derived, resulting in Equations (1)- (3) characterizes the rate at which the interaction point position changes. The propagation speed V HS2 applies from the V2 spacecraft location back to the V2 TS crossing, where !R HS2 =R TS2 -R 2o is the distance from the intercept of the time-dependent V2 spacecraft location, R V2 = R 2o +v 2 (t−t o ), to the V2 TS crossing position (v 2 is the radial speed of V2 and R 2o is the position of V2 at time t o ). Figure 4 shows a schematic representation of the parameterized propagation we employ. Setting = R R int HP for simplicity and defininga = V V HS HS2 , we rewrite Equations (2) and (3) For simplicity, we do not treat the problem three-dimensionally. Instead, we treat all paths the same way, i.e., onedimensionally assuming spherical symmetry, until there is a bifurcation upon reaching a boundary (the TS or HP) at which a signal propagates differently along one or the other path, and then each of the daughter paths are treated separately (illustrated in Figure 4). (Appendix B provides additional details of the propagation technique.) As an example, referencing Figure 4, consider a solar spherical disturbance propagating from the Sun, past ACE (at a time t ACE ), until NH is reached (t NH ), and onward toward V2 in the HS. Then, when the signal propagating out to V2 distances crosses the TS (t 1 ), the first bifurcation between the V2 path and the other paths occurs. Following the V2 path, notice that the disturbance continues until the V2 spacecraft is reached (t V2 ). Note that the "eye" symbols in Figure 4 indicate points of observation of the disturbance signal in the GCR intensity measurements. Also, physical interactions affecting the signal only occur where open circles are drawn in Figure 4; for example, when a propagating disturbance traveling along the V1 path reaches the radial distance of V2, there is no effect on the V1 signal because there is no boundary being crossed there along the V1 path. This linear formulation makes no explicit attempt to represent complications, such as mass loading of the solar wind plasma due to pickup ion populations (Richardson et al. 1995), which could cause continuous radial changes, such as a gradual slowdown of the propagation speed. Returning to t 1 , the signal traveling to V1 can continue in the solar wind until there is another bifurcation (t 2 ) between the solar disturbance path and the V1 path at the TS. Our generalized propagation equations do not require t 2 to equal t S1 (i.e., when the signal would reach the distance V1 crossed the TS), although it does not exclude that as a possibility. The solar disturbance continues, traveling through the heliosheath until the heliopause is reached (t 3 ), where, as done for t 2 , the t 3 parameter allows for HP distances different than the V1 observed distance (t HP1 ). The signal then propagates away from the contact point at the LISM propagation speed until it reaches either the interaction point on the LISM field line that is connected to V1, or to V1 itself, at time t int . It should be noted that this formulation is technically agnostic to the path that the effects of a solar disturbance takes to reach V1 and the conventional approach is taken below, without success. Alternatively, if the speed of propagation in the LISM is assumed (i.e., using the Alfvén velocity) then the distance from the HP to a "connection point" can be established.
We categorize each parameter in the equations as being (1) empirically derived (e.g., m and t o ), (2) known (e.g., locations of the spacecraft), (3) assumed from nominal known values (e.g., V sw ), or (4) as free parameters determined by the equations (e.g., u int and V HS ). Table 1 summarizes the categories and values we have chosen for this nominal configuration that result in the agreement presented in Section 3.2.
Setting the system up in the conventional sense, where the solar-originating disturbance must propagate directly to the location of V1 to affect the observed GCR rate, would require a propagation speed in the LISM of ∼107 au yr −1 (∼508 km s −1 ). This is notable, as the estimated Alfvén speed in the LISM is only ∼40 km s −1 (Kóta & Jokipii 2017). This motivates an alternative picture, namely that the solar-originating disturbance is instead remotely modulating the V1 observations through interactions along the magnetic field line near the HP connected to V1. This is illustrated in Figure 5. In this picture, the lag time of the solar disturbance signal may not correspond to the shortest radial distance from the Sun to the heliopause (taken to be a constant distance, R HP =120 au from the Sun), but rather the point that first transmits the signal to the eventual interaction point through the LISM to V1. As a guiding analogy, note that this bears similarity to the path of a light ray traveling through dissimilar optical media under Fermat's principle of least time, which results in the path with the shortest travel time rather than the path with the shortest distance. showing the nose and near flanks and a single illustrative magnetic field line in the LISM that threads the V1 position beyond the heliopause (HP). V2 is in the heliosheath (HS) beyond the TS, and NH and ACE are in the supersonic solar wind in the outer and inner heliosphere, respectively. A disturbance encloses and originates at the Sun, but we consider only the solar disturbance path (a path of least time; see text) because it is along this path that the signal can first reach the LISM, first at the contact point on the HP and then at an interaction point on the threading field line. Decreasing intensities at near-perpendicular pitch angles result in anisotropies in GCR intensities measured at V1 (see Figure 6). This happens when the particles traveling along the field reach the interaction point where an expanding disturbance traps and adiabatically cools particles that lose energy in a nongyrotropic manner, resulting in the distinctive anisotropies. Because large disturbances, such as GMIRs, will reach much of the HP, the disturbance at the HP near V1 will eventually transmit the signal directly to the V1 position, at which time it might show up, for example, as an in situ shock feature.

Cosmic Rays in the Supersonic Solar Wind and Heliosheath
By comparing the NH measurements with similar cosmic ray data from ACE at 1 au and V2 in the heliosheath (Figure 3), we see temporal variations in the GCR rate propagate into the outer heliosphere and the heliosheath in an approximately ballistic fashion (Figure 3(c)), with the cosmic ray minimum (solar cycle 24 maximum) conditions expanding outward correspondingly. We examine the rate data before de-enveloping (Figure 3(b)). At ACE, the lowest comparable rates of the cycle cluster at three separate times: at the beginnings of intervals L and M, and at the end of N. Complicating issues, there is a local maximum between the M and N intervals, but the overall L-M-N period ( At V2, the motion through the steep GCR intensity gradient makes it difficult to identify the minima. At V1, the behavior of solar maximum conditions in the LISM is unclear, but the GCR features (namely the local maximum between intervals M and N), in the context of the propagation we performed, allow us to identify nominal solar maximum timing at all four spacecraft in this study. We thus identify solar maximum conditions to have reached the various spacecraft, as follows: 2014.8 at 1 au (ACE), 2015.2 at 31 au in the middle heliosphere (NH), 2017.8 at 115 au in the heliosheath (V2), and 2017.5 at 139 au in the interstellar medium (V1). The fact that solar maximum first reaches the more distant V1 before the nearer V2 is a consequence of the multiple routes for propagating disturbances, the nonuniform HS thickness, and the relative positions of the spacecraft relative to the HP, as discussed in Section 3.3 (see also Figures 4 and 5), but because of these complications, the concept of solar maximum conditions might not be useful in the LISM.
Superposed on these large-scale variations are Forbush decreases, which, from Voyager measurements, are often observed in the outer heliosphere and are particularly obvious when multiple interplanetary shocks and/or CIR structures merge to form GMIRs. Smaller features in the GCR modulation can also be associated with shocks. One example (shown in Figure 2 Half (51%) of the lower-energy particle activity events in Figure 2 72% and 20%, respectively, which is consistent with the tendency for recurrent CIR-related events, seen during the latter period, to have less turbulent, eroded magnetic fields that are relatively ineffective at shielding cosmic rays, resulting in fewer Forbush decreases. This distinction between traveling shocks, with their turbulent magnetic field and locally accelerated particles, and CIR-associated recurrent events, with their worn-down magnetic field structure and residual energetic particles, has been pointed out by Decker et al. (2005). The CIR events occur preferentially after solar maximum, consistent with our observations of fewer Forbush decreases after solar maximum and with prior inner heliospheric CIR studies (e.g., Richardson 2018; Allen et al. 2019, and references therein). Recently, Kollmann et al. (2019a) showed that there is an association between the spectral slopes of suprathermal tails observed by PEPSSI and the solar wind speeds from SWAP. This association demonstrates that the events in the later period are residual CIR events that propagate into the outer heliosphere, distinct from the activity before and at solar maximum. We now confirm that GCRs behave as would be expected in this picture. This is important specifically for NH because it enables us to make inferences about the magnetic field, which NH cannot measure; for example, apparent shocks observed with SWAP, which cannot unambiguously be identified as shocks without a measured magnetic field, could now be compared with the cosmic ray variations to see whether there is a Forbush decrease in the PEPSSI GCR observations that would support a shock classification, whereas the lack of a Forbush decrease would suggest that the feature is the residual of a shock formed in the inner heliosphere at CIR interfaces.
The V2 cosmic ray measurements, in comparison with NH and ACE, show that Forbush decreases often evolve once in the heliosheath. Although the GCR features in the heliosheath from V2 do seem to evolve, they also still appear to follow, to zeroth-order, 1D ballistic propagation through the heliosheath. The disturbances propagate at roughly the fast mode speed (Richardson et al. 2017), which changes depending on region but is dominated by an Alfvén speed that is larger than the sound speed in the regions considered here. Our assumptions do admit the possibility that reflections inside the heliosheath (e.g., Washimi et al. 2011) between the TS and heliopause could contribute to the more complicated GCR profiles at V2.

Cosmic Rays in the LISM
We find a simple explanation for the anisotropies in the GCR fluxes measured by V1 in the LISM as shown in Figure 6. A combination of solar disturbances propagating beyond the heliopause and the field lines at V1 magnetically connecting back to the draping region of the heliopause can result in the observed distributions. Such a scenario is shown schematically in the right panel of Figure 7. The prevailing interpretation Kóta & Jokipii 2017) that the spacecraft must be colocated with trapped GCRs and therefore be between two or more magnetic structures for the anisotropy episodes to occur (illustrated in the left panel of Figure 7) does not appear to be consistent with the data presented here. The primary challenges to the trapped GCR interpretation include the following: (1) a comparison of the GCR intensity time profiles at V1 and V2 cannot be reconciled by this explanation using either the simple 1D ballistic propagation or the more generalized equations that allow for variations in the TS and HP locations using physical propagation speeds in the LISM (Appendix B.2); and (2) the observations at V1 of enhancements in the magnetic field relative to the anisotropy episodes (see Figure 6) are inconsistent with a disturbance enveloping the V1 spacecraft. Although our new explanation also involves particle trapping and cooling in a disturbed region, a distinguishing characteristic is that escaping GCRs exhibit the characteristic anisotropies and are also observable far from the disturbed region.
With respect to (1), there is no constant temporal offset that can reconcile the V1 and V2 time series (even approximately). In addition, because in 2016 and 2017 the disturbance features are seen simultaneously or even earlier at the more distant V1 than at the closer V2, it can be inferred that propagation is much faster in the interstellar medium than in the heliosheath, which is at odds with expectations of slow Alfvén-like LISM propagation. In turn, this implies the shock signatures and the anisotropy episodes should be seen at similar times and the shock signature should not increasingly lag behind the anisotropy episodes, as is seen (Figure 6). While several studies (e.g., Washimi et al. 2011) have shown that waves from transient disturbances both transmit and reflect at the heliopause boundary, there is no reason to expect that this would only affect disturbances after 2015. As such, reflection at the heliopause is unable to fully explain this systematic change in the ordering of events between V1 and V2.
To explain these apparent dichotomies, we note that the empirically derived propagation of the signal in the LISM depends weakly (at most) on the distance between the heliopause and V1-that is, once the disturbance reaches the HP, there is very little delay before the signal (the anisotropy episode) is observable at V1. This can occur if V1, despite being up to ∼30 au upstream, is still magnetically connected to field lines that are close to the heliopause. In that case, the GCRs traveling along the field lines move orders of magnitude faster than the disturbances simply propagating radially outward through the LISM, heliosheath, or solar wind. This configuration follows earlier works (e.g., Opher et al. 2017), which indicate that the LISM field lines approach the heliosphere at an oblique angle and drape over the HP asymmetrically such that a magnetic field line in close proximity to the HP along the flank of the heliosphere is also far from the HP at the heliospheric nose. One problem for this scenario is that it only produces a disturbance along the field line on one side of the spacecraft. The field on the far side of the spacecraft, deeper in the LISM, is comparatively undisturbed, so there is no reason to expect strong field compressions  (Burlaga et al. 2018) in the LISM responds to in situ activity, such as transient shocks (Burlaga et al. 2013) and compressions (labeled in red) reaching V1. LECP instrument on V1 measures >211 MeV GCRs using a dual-ended telescope mechanically scanning through  360 . Intensities of the GCRs decrease around  90 pitch angles (sectors 1 and 5, perpendicular to the field), whereas the intensities in other directions remain essentially unchanged. Onsets of three long periods of anisotropy are indicated in green; these onsets are interpreted as times when the GCRs escaping from the disturbance region reach V1. Disturbance is associated with the transient events' entry into the LISM and generation of an expanding structure that changes the energy and pitch angle distribution of GCRs that are processed in this region. Systematically increasing lag between the anisotropy onset and the in situ evidence of the disturbance (shock or compression), from 1 day to 102 days to 620 days, is a natural consequence of the one-sided LISM field interaction we present here.
to be present at just the right time to result in V1 being positioned between two compressions, as would be required.
Thus, this interpretation requires a mechanism in which a disturbance on only one side of the spacecraft can result in a decrease in intensities near 90°pitch angles. (We compare the trapped and escaping explanations in Figure 7.) The trapped configuration picture ) discussed in the introduction relies on magnetic mirroring to reflect large pitch angles and transmit small, field-aligned pitch angles (the loss cone). Kóta & Jokipii (2017) have extended the Roelof mechanism to not rely solely on mirroring (although it does include it); it has the feature that perpendicular GCRs tend to be trapped in the expanding field and become adiabatically cooled. When the shock travels directly over the spacecraft and expands outward, it compresses the LISM magnetic field. This compressed post-shock field then slowly decreases with the shock expansion. Kóta & Jokipii (2017) considered cases where particles are reflected off of or transmitted through the field increases as well as the case where particles enter and become trapped behind the shock in the decreasing field region. Particles gyrating in the expanding region behind the shock and its weakening magnetic field lose energy due to adiabatic cooling that results from conservation of the first adiabatic invariant (Roederer 1970). Particles eventually escape and continue along the field line toward V1. Energy loss for particles traveling more perpendicular than parallel to the magnetic field will result in a portion of the population falling below the lower-energy threshold of the cosmic ray channel. Because the channels that measure cosmic rays are integral channels, responding to all particles above a given energy (∼200 MeV for LECP), there is an observed net loss of particles. This effect will occur regardless of the GCR spectrum, but is more pronounced the steeper the spectrum is falling with increasing energy. Lower-energy particles decrease the count rate when they leave the measured energy range, while higher-energy particles are counted regardless. Adiabatic cooling is able to reduce the intensity of perpendicular GCRs that were traveling toward V1 and were intercepted by a propagating shock. It prevents GCRs traveling in the opposite direction from reflecting back toward V1 and filling in the ∼90°pitch angles, and therefore it explains the observed anisotropy ( Figure 6). Rankin et al. (2019) found this mechanism to be effective. They conducted an analysis of the anisotropies observed at V1 using the CRS instrument. Using magnetic field perturbations on either side of the spacecraft, they then related GCR anisotropies to shocks near V1, and their observations are consistent with the CRS data. Figure 1 of Kóta & Jokipii (2017), Figure KJ1, and Figure  10 of Rankin et al. (2019), Figure R10, illustrate a transient disturbance entering the LISM that is responsible for GCR anisotropy. Figure R10 includes heliospheric structures like the HP and the relevant LISM magnetic field line, and we have done similarly in Figure 7. Figure KJ1 highlights the trapping region and draws a comparison to another figure, Figure KJ2, that highlights trapped and reflecting particle paths. Our panel (a) represents Figures R10a and R10b closely, and Figure R10c is not fundamentally different since, like the other two panels, it shows the V1 spacecraft caught between compression regions. Our Figure 7(b) is different because we show the disturbance coinciding with the draping region and the V1 spacecraft "beyond" the disturbance, i.e., there is a disturbance required on only one side of the spacecraft in our scenario.
While previous studies have investigated individual V1 anisotropy events in depth (e.g., Gurnett et al. 2015), they have not yet connected the V1 observations to other spacecraft measurements or fully deconvolved systematic variations in the anisotropy onset times in relation to other spacecraft observations. The fate of GCRs that subsequently escaped from the post-shock region after being trapped have not been considered. The results of the propagation methods used in this paper also indicate that the disturbances must propagate at unphysical speeds (>500 km s −1 ) in the LISM to reach the V1 spacecraft, suggesting that the disturbances must be acting on the V1 GCR anisotropies remotely, rather than locally. Although other Figure 7. GCR anisotropy in the prevailing (a) trapped particle configuration (Kóta & Jokipii 2017;Rankin et al. 2019) and the (b) new escaping particle configuration. Far away from the heliosphere or before the appearance of the disturbance in cases (a) and (b), the GCRs are isotropic, as depicted in the anisotropy "pie plot" 1 (not shown in panel (b)), which represents, using a green filled circle, the same intensity in all directions and at all pitch angles. Magnetic field line is represented by a brown line. When in the trapped particle case (a), the disturbance enters the LISM and the particles are isotropic outside the disturbed region, as shown in pie plots 1 and 3. Only within a disturbance, such as behind a shock, can the GCRs reflect and lose energy as the magnetic field expands. This expansion preferentially reduces the perpendicular velocities, resulting in suppressed ∼90°pitch angles, which are indicated by white notches in the green circles perpendicular to the field direction. Because V1 is inside the shock region in case (a), the anisotropy episode is observable only as the disturbance passes. For the escaping particle case (b), once the particles enter the trapping region, there is the same sort of energy loss as in case (a); therefore, pie plots 2 and 5 are the same, showing the perpendicular notches. However, outside the disturbed region, the GCRs that escape maintain a reduced intensity near 90°on both sides of the disturbance, as seen in pie plots 4 and 6, so when the GCRs reach V1, an anisotropy episode will commence even though V1 is not surrounded by field compressions. Because of the LISM magnetic field's geometry, the draping region is the most likely location for disturbances to first contact field lines in front of the heliopause. See investigators have compared GCR temporal variations in the heliosphere at greatly separated positions (see Section 1), this study is the first to do so with direct measurements in the inner and outer heliosphere, the heliosheath, and the LISM, simultaneously.
To summarize, straightforward propagation of the observations within the heliosphere to the heliosheath and then into the LISM was not able to bring V1 and V2 observations into any alignment, motivating the construction of the generalized propagation scheme and the possibility of remote connection to V1 on LISM field lines. We then developed the open-field, draping region dominating picture as described in Section 3.3. However, we need to show that perpendicular GCR suppression can arise from a disturbance on only one side of the spacecraft. Like Rankin et al. (2019), we utilize the Kóta & Jokipii (2017) mechanism rather than the original Roelof mirroring mechanism. When the disturbance reaches the heliopause near the draped field and transmits from the contact point to the nearby interaction point, the Kóta & Jokipii (2017) effect rapidly changes the GCR anisotropies because the GCRs travel rapidly along the LISM field line. However, the shock propagates into the LISM itself, and when it reaches V1, the shock is observed in the magnetic field measurements. The in situ shock will lag more and more as V1 pushes farther into the interstellar medium, which is indeed what is observed ( Figure 6), with the lag growing from 1 day, to 102 days, to 620 days for the three major anisotropy episodes beginning in 2012, 2014, and 2015. The LISM propagation speed of ∼8 au yr −1 (∼40 km s −1 ; Kóta & Jokipii 2017) is consistent with all three of these episodes, shock arrivals, time lags, and the different heliosheath traversal distance near the disturbance path versus near V1. This propagation speed is close to the 40 km s −1 Alfvén speed referenced throughout the study.
Our interpretation for the growing lag between the anisotropy episodes and the magnetic field compressions explicitly includes all the features (e.g., small intensity enhancements, plasma oscillations) of the other explanations (see Gurnett et al. 2015;Rankin et al. 2019) because the usual propagation of the disturbance in the LISM is a fundamental part of our picture. The setup described here is a natural consequence of the paradigm in which the draping region is the first point where heliospheric disturbances are communicated into the LISM. It is the usual propagation that communicates the signal from the HP to the contact point in the draping region. Away from the draping region, this propagation happens as well. The only thing that one must assume is that the effect on the anisotropies is at least on the order of, and probably larger than, other mechanisms GCR intensity increases or decreases. Since the usual Forbush decreases do not show any such perpendicular intensity suppression, there is no conflict. Additional to the consistency with the observations and similarity in mechanism, our setup offers a simpler and therefore more likely scenario because it does not require that V1 happens to be within the finite-sized trapping region but can be at any point outside of it

Summary
The development of the NH/PEPSSI cosmic ray monitor and integration of the resulting data into the largest array of in situ cosmic ray observations, extending from 1 au to the interstellar medium (consisting of measurements from ACE, NH, V2, and V1), has led to several conclusions: (1) GCRs are clearly observed at both ACE and NH, revealing that GCR minimum (solar cycle maximum) occurred at 2014.8 and 2015.2 at 1 au and 31 au, respectively. Continued propagation shows that solar maximum conditions arrive at V2 in the heliosheath at 2017.8 at 115 au.
(2) Seventy-two percent of solar disturbances observed by NH before and during solar cycle maximum (2012)(2013)(2014)(2015) are associated with Forbush decreases, whereas only 20% are associated after solar maximum (2016-2017), supporting the interpretation that the earlier events are propagating shocks bringing turbulent magnetic fields (which inhibit GCR transport) and the later events are remnants of CIRs lacking magnetic fields that strongly affect GCRs. This result is consistent with earlier Voyager observations (Decker et al. 2005) and analysis of Cassini and NH suprathermal tail spectral slopes (Kollmann et al. 2019a).
(3) Simple ballistic propagation of solar disturbances brings GCR intensity variations inside the heliosphere (i.e., excluding the V1 LISM measurements) into general agreement between features ( Figure 3). (4) The timing of V1 interstellar observations of GCR intensity variations cannot be brought into agreement with the ACE, NH, and V2 data utilizing traditional disturbance propagation into the LISM using either the 1D ballistic propagation or the more generalized formulation (Appendix B.2). By "traditional disturbance propagation," we mean direct, in situ propagation such that the disturbed parcels of plasma are required to reach the point of observation for there to be a causally related, detectable variation in the GCR intensity.  (1)) leads to the conclusion that the solar disturbances are likely interacting with the LISM magnetic field lines in the draping region. (6) A newly proposed escaping GCR scenario is developed in which shocks first enter the LISM near the field line draping region rather than near the V1 spacecraft trajectory (following formulations of Kóta & Jokipii (2017)). Particles with preferentially perpendicular pitch angles will lose energy while trapped inside the expanding region behind the shock, and then eventually a fraction will escape, continuing along the field line toward V1. (7) This description of the GCR anisotropy episodes along with the traditional propagation of the disturbance to V1 (which does not induce anisotropies) naturally predicts longer time lags between the anisotropy episode onsets and shocks identified in the magnetic field. The shorter observed time lag suggests that solar disturbances first interact with LISM magnetic field lines in the field line draping region that are connected to V1. A propagation speed of ∼40 km s −1 in the LISM is a reasonable Alfvén speed and is consistent with the time lags and all the parameters producing the V1-V2 agreement.
Future avenues for progress could focus on further exploring the escaping particle picture. This could include incorporating more sophisticated modeling efforts as well as investigating how this picture compares to the plasma oscillations observed at V1. Moreover, future work could more comprehensively compare this picture to the idea of precursor events, as well as more in-depth investigations into observed large pressure pulses and how they may play a role in this process. Additionally, further exploring the variations in the propagated GCR profiles may provide additional insight into the evolution of solar disturbances as they propagate into the LISM.
The NH work was supported by NASA NH funding through contract NAS5-97271/Task Order 30 and NASA Voyager Interstellar Mission NNN06AA01C. NH/PEPSSI data are publicly available through NASA's planetary data archive (https://pds.nasa.gov). We thank the ACE SWEPAM instrument team and the ACE Science Center (http://www.srl. caltech.edu/ACE/ASC/level2/new/) for providing the solar wind data. We appreciate the efforts of the Voyager MAG team in providing V1 magnetic field measurements (https:// omniweb.gsfc.nasa.gov). We acknowledge and thank the network of observers who provided the sunspot number measurements, courtesy of the Royal Observatory of Belgium, Brussels (http://www.sidc.be/silso/home). M.E.H. would like to acknowledge helpful conversations with members of the NASA DRIVE science center, Our Heliospheric Shield, NASA grant 18-DRIVE18_2-0029, 80NSSC20K0603, during the revision stage of this paper.

Appendix A PEPSSI Cosmic Ray Monitor: Experimental Consideration
Some details concerning the instrumental hardware and analysis techniques, while unnecessary for the observational and scientific discussion, are nonetheless important to document because this report is the first on measuring GCRs with the PEPSSI detector, which was not designed to make these measurements. To date, three papers represent the instrumental description of the PEPSSI instrument: the original instrument paper of McNutt et al. (2007), Appendix A of Kollmann et al. (2019a), and Appendix A in the present work.

A.1. Hardware
The PEPSSI instrument was not designed to measure cosmic rays, but the ubiquity and high energy of these particles means that they nonetheless produce a strong signal, which is usually considered background to be suppressed with coincidence techniques. However, PEPSSI's nominal electron measurement is based on single SSDs without coincidences behind a 1 μm aluminum flashing (to suppress ions having energies below a few hundred keV nucleon −1 ) and assumes that the electrons are collimated by the instrument aperture (McNutt et al. 2008). This is an effective electron measurement technique in the presence of sufficient electron foreground, such as we observed at Jupiter (McNutt et al. 2007;Haggerty et al. 2009;Kollmann et al. 2014). Cosmic rays are energetic enough to pass through the flashing, the collimator fins, the walls of the instrument, and indeed even the complete body of the spacecraft, providing a large geometric factor. Combined with the high (∼100%) detection efficiency, the SSD-only system becomes a sensitive cosmic ray detector by measuring the deposited energy and relating it to the energy loss for given ranges of incident energies (Figure 8). Although in the heliosphere there are far fewer cosmic rays than low-energy particles, cosmic rays dominate the SSD count rates because these are accepted over 4π steradians incidence whereas low-energy particles need to pass through the narrow apertures. To avoid the confusion between the intended foreground for the PEPSSI experiment and the foreground for this study, we refer to these SSD-only measurements as either collimated when referring to electrons and ions in the intended PEPSSI range below ∼1 MeV or penetrating when referring to the high-energy ions that do not (necessarily) enter through the instrument aperture. The total, uncorrected count rate that is measured is the sum of penetrating and collimated particles.
The study of particles' interaction with matter (e.g., Ziegler 1980) has produced a literature composed of laboratory measurements, theoretical formulations, and computational techniques to quantitatively relate energy loss and deposition with the rate of energy loss through a given range of material (or stopping power). The stopping power (energy loss per distance, normalized to target mass density) for proton projectiles through a silicon target declines from the Bragg curve peak value of 5.7 MeV/(g cm −2 ) at 500 keV to the minimum ionizing energy loss of 1.6 MeV/(g cm −2 ) at 2 GeV. Here, we only consider the energy loss due to ionization that is measured by an SSD. The energy loss by a projectile having a kinetic energy above the minimum ionizing energy increases very slowly-increasing only to 2.1 MeV/(g cm −2 ) at 100 GeV-so it is a good approximation to treat this inverse relationship between particle energy and energy loss as being one to one (in a statistical sense). This approximation is particularly good because the falling cosmic ray spectrum results in a negligible contribution from particles above a few GeV.
For the PEPSSI SSD thickness of 500 μm, we get a minimizing ionizing energy loss of 213 keV for all GCRs 2 GeV. The actual energy loss (corresponding to the Maximum energy bin is an overflow bin that counts all detected particles that deposit energy at or above 1.1 MeV because energy deposits above this energy saturate the electronics and are registered at full scale, data number 1023. High counts within this dump bin therefore do not indicate a separate particle population. measured energy deposit at these energies) for the cosmic ray population is spread out not only because of the incident spectrum but also because of the stochastic nature of energy loss. The actual energy loss is also spread out because the path through the SSD is longer at oblique angles-resulting, for example, in twice the nominal energy deposit at an angle of incidence of 60°. The distribution of energy deposited by over 10 million particles in the PEPSSI SSDs during the period from 2015 DOY 200 to 2017 DOY 199 ( Figure 8) shows a wide peak near 200 keV, fully consistent with expectations. PEPSSI accumulates counts over three energy ranges in each of the three flashed SSDs. These ranges are identified with an "R" (for "rate") and either 00, 01, or 02. The three look directions are identified with an "S" (for "sector") and either 00, 02, or 05 (for three of the six PEPSSI sectors equipped with flashed SSDs); so, for example, R01S05 refers to rate R01 and sector S05 (McNutt et al. 2008). Tables 2 and 3 provide the specifications for these rates; the measured deposited energy is what the SSD detects on the basis of the interaction between the penetrating particle and the silicon crystal; penetrating ion energy is the kinetic energy of the ion as it strikes the SSD (inside the instrument); and ambient particle energy is the energy of the cosmic ray in space before it passes through shielding material (e.g., the walls of the instrument) and is detected by the SSD.

A.2. Instrumental Analysis: Calibration and Background
Energetic particles and pickup and suprathermal ions measured by PEPSSI (McNutt et al. 2008;Hill et al. 2009;Kollmann et al. 2019a) are collimated particles with energies too low to penetrate the sensor from the side. The signal that is dominated by pickup ions and suprathermal ions is associated with "double coincidence" measurements from the time-of-flight-only (TOF) system. When particles trigger the start and stop detection in the TOF system but either miss the SSDs or are not energetic enough to penetrate the flashing (or dead layer), energies down below the pickup ion cutoff are included and dominate. Although both protons and helium ions/nuclei are present, the calibration effort to date, discussed in detail by Kollmann et al. (2019a), indicates that this population is dominated by singly charged helium atoms (He + ) picked up by charge exchange from the neutral helium component of the interstellar medium (Möbius et al. 1999(Möbius et al. , 2004(Möbius et al. , 2009). SWAP has extended those measurements significantly by providing the first in situ observations of interstellar proton pickup ions, produced by charge exchange between solar wind protons and neutral hydrogen, in the outer heliosphere (McComas et al. 2010).
To obtain differential intensities, background corrections are required. We use several techniques to determine the background corrections, a term we are using broadly to include all the conversions and processing required to convert from the instrument-specific data number in the telemetry to fully calibrated physical units. The geometry factor for a single detector in space is G=πA=2.8 cm 2 sr, where A is the detector area of a PEPSSI SSD. This assumes that the instrument walls, structure, and components provide no reduction in detection efficiency, which is a good approximation for cosmic ray energies. We have neglected the "density effect" (Ziegler 1999) because the slight (few percent) decrease in energy loss relative to the Bethe-Bloch formula above ∼1 GeV will tend to lessen the influence of the instrument casing and is a small effect relative to the simplifying assumptions we have made for this first analysis. For the highest-energy channel, we multiply the geometry factor by two to account for the fact that both sides of the detector are accessible to the cosmic rays that can easily penetrate the entire spacecraft bus. Because the cosmic ray distribution function is essentially isotropic and each detector should measure the same thing, we are able to use the counting rate of different sectors at the same energy (e.g., R00S00, R00S02, and R00S05) to determine small effective differences in the geometry factor. For these rates, we found the R00S02 geometry factor to be 97% of the R00S00 and the R00S05 geometry factor to be 117% of R00S00. At these energies, these geometry factors are attributable to different shielding geometries, which have not yet been modeled. Notes. a Includes the effect of energy loss in the instrument housing and the varied angle of incidence of particles on the SSD surface. Uncertainty is based on 70% of the isotropic population having 0°-35°incidence, which corresponds to an energy deposit ranging from 1 to 1.2 times the minimum energy deposit. b This channel also has response from the nominal 5-18 keV range, which is below threshold but contributes some background, which we subtract. c This integral channel has no upper limit, but an estimated passband is allowable because of the known high-energy portion of the cosmic ray spectrum, which is approximately the same throughout the heliosphere and unaffected by local conditions. We also subtracted off estimated background rates because of the collimated particles that are detected in SSD-only mode, radioisotope thermoelectric generator background, and electronic noise. For fine adjustment of the intensities, we are able to use the highest-energy R01 channels and compare them to an independently measured cosmic ray spectrum above 2 GeV, which is essentially unchanging throughout the heliosphere. We have reduced visibility into the energetic particle populations in the ∼1-75 MeV range, but not zero visibility. In Figure 8, the pulse height analysis (PHA) data show the overflow bin, which detects any particles that deposit more than 1.1 MeV. This bin is in both the SSD-only data (used for the cosmic ray) as well as triple coincidence data (associated with the collimated foreground), so we can directly measure whether the integrated spectrum from 1 to 75 MeV is large enough to indicate there is a significant SEP contribution to the GCR signal. Even assuming an extremely hard, E −1 , spectrum above 1 MeV, the maximum SEP contribution to the GCR signal is 2%. The actual 1-75 MeV integrated spectrum (the counts in the overflow bin) is an order of magnitude lower than what it would be in the case of the E −1 spectrum assumption, so we do not have a problem with SEP contaminating the GCR signal.
To determine the species of cosmic rays, we consider the characteristics of energy loss (see Northcliffe & Schilling 1970), noting that the stopping power in the relevant energy range is proportional to the square of the atomic number-the heavier ions lose much more energy, which inhibits heavier ions from penetrating through the side of the instrument. These two facts mean that other prevalent ions would complicate the energy distribution by contributing to a peak that is not centered around 213 keV as seen in Figure 8. Thus, our observations are consistent with the measurements from purpose-built cosmic ray detectors that find the abundance of the heavier ions is down by at least an order of magnitude. Distinguishing penetrating electrons and protons from one another poses an additional challenge because both species have nearly the same minimum stopping power, corresponding to the idealized normal-incidence energy deposit in the SSD of 196 and 213 keV, respectively. What is different is that electrons are minimum ionizing at 1.25 MeV, whereas protons are at 2.5 GeV. However, for electronic stopping power, the relevant mechanism for detection in the SSDs, the higher-energy electrons appear as essentially minimum ionizing particles. Because cosmic ray electrons have significantly lower intensity than protons at the same energy (Blasi 2013), we might assume that the electron contribution is negligible. However, Cummings et al. (2016) have reported GCR electron and proton intensity, which are comparable below 100 MeV, and we used their measurements to estimate the electron contribution ( Figure 2). Fortunately, for the present study of >75 MeV GCRs, our aim is to study the temporal variations in cosmic rays arising from the passage of interplanetary disturbances. If some fraction of the population does not respond according to the usual cosmic ray mechanisms like Forbush decreases, it does not actually negatively impact the conclusions we draw.

Appendix B Propagation-correction Technique
The straightforward propagation inside the heliosphere is handled numerically, as described in Appendix B.1. The physical model associated with propagation from the heliosheath into the interstellar medium is described in detail in the body of this report (Section 3.3) and is based on a new picture of interaction that leads to the unique V1 anisotropy observations. The parameters of the linear function Equation (1) are described in Equations (2) and (3) and derived in Appendix B.2.

B.1. Simple Ballistic Propagation
To determine the propagation-corrected time series, we used ballistic propagation of the solar wind disturbances that affect the cosmic ray intensities. To shift the ACE data to the time at which that solar wind structure would be at the location of NH, the measured solar wind velocity at ACE was taken and smoothed with a 27 days boxcar average (i.e., duration of one solar rotation). This smoothed velocity was then used to ballistically propagate the solar wind until it would arrive at the distance of the NH spacecraft. No propagation adjustments are made to the NH data.
To backpropagate the V2 data to when the solar wind structures would have been at the location of NH, we first took 27 day (solar rotation) smoothed velocities inferred from SWAP measurements at NH (see Elliott et al. 2016Elliott et al. , 2018 to estimate the amount of time the solar wind would take to reach a distance of 79 au (an estimated location of the TS). While this termination shock distance is less than that observed by Voyager 2 (84 au), we use 79 au because it gave the best agreement between the Voyager 2 GCR intensity variations and both ACE and NH. Next, the half-year smoothed flow velocity measured at V2 was used to estimate the time the currently measured plasma would have been at 79 au. A longer smoothing of the plasma velocity was done for Voyager 2 in order to reduce the variation in the estimated time of propagation caused by the slower heliosheath flow speeds (∼100 km s −1 ) leading to a higher percent change in speed between observations due to the large distances between observations. The lag times from NH to the estimated TS were then interpolated to match the times of the backtraced V2 lag times so they could be combined into a total lag time from NH to V2.
To backpropagate the V1 data to match the NH time stamps, we used a simple linear model to compute when the V1 data would align with the V2 data (see Section 3.2). This approach was motivated by the observed change in the order of feature observations in the data (Figure 3). Once time shifted to V2, we used the same time lags to shift the V1 data to NH as we did with the V2 data. The model used to shift the V1 data to match V2 assumes the TS position was at 85 au.

B.2. Heliosheath and Interstellar Analytical Propagation
The propagation of disturbances from the heliosheath into the interstellar medium is an exploratory effort; therefore, we modeled the signal propagation analytically, allowing flexibility to model the mechanisms with a combination of constrained and free parameters. Unlike the method in Appendix B.1, in which a linear equation without free physical parameters is used, here we first determined the physically parameterized linear equation that brought the V1 and V2 time series into agreement and then used the well-constrained parameters to calculate the values of the free parameter that satisfy the equation.
We describe the signal transmission using a sequence of propagation steps between V2 and V1. Each step involves a number of analytic curves, which we will term traces, that describe distance relative to the Sun as a function of time, e.g., a spacecraft trajectory or the position of a heliospheric boundary.
For trajectories, the interpretation of a trace is physical-the function of the form = + x t Vt L ( ) is the linear approximation of the spacecraft position x as a function of time t, where the spacecraft speed is V and the position at time zero is L. x is identical to the heliographic radius if the motion is radial. If there are nonradial excursions (which are significant for  90 pitch angle particles on the LISM field line), x is the distance along the propagation path. Also, to keep the mathematics uniform, we not only describe trajectories but also boundaries (interfaces where signal characteristics can change) that are constant in time as functions parameterized by time. This parameterization is done to account for the nonspherical nature of the system. Thus, instead of requiring more than one spatial dimension in the model, we just allow boundaries to have a positions that depend implicitly on a set of angles J j , that determine latitude and longitude, J j R , ( ). Because, for a given boundary, the latitude and longitude are specified by the intersecting trace, which is only dependent on time, we can rewrite the angularly dependent boundary position in terms of the time parameter, J j J j = = R R t t R t , , ( ) ( ( ) ( )) ( ). We require that this be able to be expressed as a linear function of t with constants u and r, where u has dimensions of speed and r is the boundary position at t=0, thence = + R t ut r ( ) . The path of the disturbance itself is also represented by a trace.
Let us start with two general traces, R i (t) and R j (t), and the path, R ij (t), of the signal propagating between the two trajectories: Figure 9 shows these traces graphically. Given a propagating disturbance signal observable along trajectory R j (t) at time t j , we seek to determine the time t i when the signal propagating at velocity V ij crosses trajectory R i (t). The relationship between t i and t j is found by analyzing the intersection of the trajectories and signal path at both times. Thus, we require = R t R t i i ij i ( ) ( ) and = R t R t j j ij j ( ) ( ), and eliminating r ij , we arrive at the following relation: Now we generalize this. Instead of stepping from i to j and j to k (replace i with j and j with k in (B2) and so on, we step from trajectory 1 to trajectory N by setting = + j i 1 and proceeding recursively to determine the time t N 1 ( ) when the signal reaching trajectory R t N ( ) at time t N is observable at trajectory R t 1 ( ), . Relationship between times t i and t j is used to compare the passage of a propagating disturbance signal between points of observation or transition on the trajectories, and is given in Equation (B2). Figure 4 represents the physical model built up of use of Equation (B2), as described in Appendix B.2. Figure 3 shows the application of the propagation correction. Boundary: disturbance crosses the TS along the path between the Sun and the interaction point

( )
Note that (B1)-(B5) do not require the r i values to be steadily increasing or decreasing. Our propagation model of the disturbance described in Section 3.2 includes N =5 steps that we sketch in Figure 4 and quantify in Table 4. We calculate the offset between the Voyager observations by using a path that mathematically connects the interaction point (which the disturbance reaches just before it is detected by Voyager 1) with the contact points and termination shock located inward, and to Voyager 2 located outward of the termination shock. We recover Equations (2) and (3), inserting the model-specific parameters as given in Table 4 into Equations (B4) and (B5), expanding the result, and defining m≡M 5 and D º t T 5 . The above means that we start out by considering the trajectories of the interaction and contact points on the heliopause (see Figure 5) because even with a steady-state LISM field configuration and heliopause, the interaction points can change position because, as V1 moves across different field lines, it becomes magnetically connected to different interaction points in the draping region.