\section{Introduction}
Pulsar wind nebulae (PWNe) are the bubbles inflated by the relativistic winds of neutron stars as they interact with their surroundings.
Their radio to X-ray emission
is believed to result from synchrotron radiation of the high-energy particles injected by
a rotation-powered neutron star in the presence of a strong magnetic field.
As such, PWNe represent an ideal laboratory to study the physics of neutron stars and particle acceleration.  They are also good indicators for the presence of an active rotation-powered pulsar;
see~\citet{Gaensler2006} for a detailed
review on the structure and evolution of PWNe.
While the Crab Nebula has been viewed to represent the prototype PWN, a growing list of PWNe show properties unlike the Crab: namely, different spectral properties (spectral index, spectral break in the radio) and different morphology.  These are referred to in the literature as ``plerions of the second kind'' or non-Crab-like plerions \citep[see][]{Woltjer1997}, a class that we have been targeting with combined radio and X-ray studies;
see~\citet{SafiHarb2012} for a recent review.
This paper addresses, using X-ray observations, the nature of CTB~87; the least studied object (at least in X-rays) among this growing class of PWNe.
In a subsequent paper (Kothes et al. in preparation), we present a dedicated study of the remnant at radio wavelengths.

CTB~87 (G74.9+1.2)\footnote{Main properties summarized at www.mrao.cam.ac.uk/surveys/snrs/\\snrs.G74.9+1.2.html~\citep{Green2009}}$^,$\footnote{High-energy properties summarized at www.physics.umanitoba.ca/snr/SNRcat/\\SNRrecord.php?id=G074.9p01.2 \citep{Ferrand2012}} has been classified as a plerionic type SNR based on its centrally-filled morphology, lack of evidence of a shell, and linearly polarized radio flux~\citep{Dickel1975, Duin1975}.  In the radio, CTB~87 has a size of 8\arcmin$\times$6\arcmin\ and a flux density of 9~Jy at 1~GHz~\citep{Green2009}.
Early X-ray observations with the \textit{Einstein} satellite also identified CTB~87 as a plerionic SNR, with an X-ray luminosity nearly 100 times weaker than the Crab Nebula in the $0.15-3$~keV band~\citep{Wilson1980}.
Newer sensitive radio studies indicated that the central bright structure has a relatively steep radio spectrum, with a spectral index $\alpha$ = 0.5
($S_{\nu}~\sim~\nu^{-\alpha}$), which is unusual for a plerionic SNR.
In addition, larger diffuse emission, $\sim$16\arcmin\ in size, has been identified and characterized by a much flatter spectrum which is more consistent of a plerionic SNR (Kothes et al. 2003; Kothes et al. in preparation).
\citet{Sitnik2010} suggested that CTB~87 is located in the Cygnus Arm,
and detected a very faint optical shell with 20\arcmin\ diameter in the vicinity of CTB~87.  This optical shell has not been however confirmed and the author concludes that the small amount of data prevents them from drawing a firm conclusion.
Locating CTB~87 in the Cygnus Arm puts the remnant at a much smaller distance than that estimated using HI data from the Canadian Galactic Plane Survey ($6.1 \pm 0.9$~kpc,~\citet{Kothes2003}), a distance that
implies a location within the Perseus spiral arm.
An \textit{EGRET} source has been suggested to be possibly associated with CTB~87, although~\citet{Halpern2001} argued that the SNR is too weak and far away (using an old distance estimate of 12~kpc) to be the most likely candidate.
The revised kinematic distance to CTB~87 of 6.1~kpc
weakened this argument, further increasing interest in CTB~87 as a potential $\gamma$-ray source.
Other $\gamma$-ray studies of the area near CTB~87 with \textit{MILAGRO}  led to the detection of an unresolved source,
MGRO~J2019+37~\citep{Abdo2007},
which was further studied and resolved with \textit{VERITAS}.
The \textit{VERITAS} source, VER~J2016+372, showed an absence of variability throughout the observations and a spectrum similar to other PWNe previously detected in the VHE band~\citep{Aliu2011}.
Finally, \textit{FERMI} detected a high-energy source, 2FGL~2015.6+3709, which was found to be likely variable and therefore a physical association with CTB~87 was deemed improbable~\citep{Nolan2012}.

Here we present a detailed study of the X-ray observations of
CTB~87.
 Preliminary results have been published in~\citet{SafiHarb2011} and~\citet{SafiHarb2012}.
This paper is focused on the observations acquired with the \textit{Chandra X-ray Observatory}, with a total exposure time of $\sim$87~ks.
We also briefly present an unpublished archival observation with the \textit{Advanced Satellite for Cosmology and Astrophysics (ASCA)}, primarily to compare the overall spectrum of the nebula with that obtained with \textit{Chandra}.  Section~\ref{section:obs} discusses the observations and data processing. In Section~\ref{section:images}, we present the observed X-ray morphology and show the first evidence for a compact nebula, with a torus/jet-like structure, as well as a diffuse nebula that is offset from the radio nebula. These structures surround a point source, the putative neutron star powering CTB~87.
 Section~\ref{section:spectra} presents the spectral analysis, including the first spatially resolved spectroscopic analysis. We show that the spectra from the point source and the compact and diffuse nebulae are well described by a power law with a hard photon index, and that the index generally steepens away from the point source.  In Section~\ref{section:disc} we discuss the PWN morphology, properties, and stage of evolution, and infer the properties of the putative neutron star powering the PWN.  We propose that CTB~87 is most likely an evolved PWN, with overall multi-wavelength properties similar to the PWNe in the SNRs Vela and G327.1$-$1.1.

\section{Observations}
\label{section:obs}

\object{CTB~87} (G74.9+1.2) was observed with \textit{Chandra}\ on two occasions.  An 18~ks observation was performed with ACIS-S3 (2001 July 8, obsID 1037)
and a 70~ks observation (2010 January 16, obsID 11092) was performed with ACIS-I3.
The \textit{Chandra}\ data was processed with CIAO v4.3 and CALDB v4.4.1~\citep{Fruscione2006}\footnote{http://cxc.harvard.edu/ciao/}.  The level 1 event files were reprocessed to new level 2 event files, using \textit{acis\_process\_events} (correct for CCD charge transfer inefficiencies\footnote{http://cxc.harvard.edu/ciao/why/cti.html}, remove bad pixels) and \textit{dmcopy} (remove bad grades and status, filter for good time intervals and energy).
The effective exposures after filtering were 17.8~ks (obs 1037) and 69.3~ks (obs 11092), yielding a total effective exposure of 87.1~ks.

For image creation, the two observations were merged after applying the aspect solution to each observation (\textit{reproject\_events}).
\textit{wavdetect} was used to detect point sources and \textit{dmfilth} was used to remove the sources and fill the holes with a Poisson distribution calculated from a background region surrounding each source.  \textit{asphist}, \textit{mkinstmap}, and \textit{mkexpmap} were used to create exposure maps for each CCD, which were then summed.
 \textit{dmimgthresh} was used to remove pixels with very low exposure times and the image was divided by the exposure map
to produce the final exposure corrected image.  This procedure was repeated for various energy bands and the resulting images are presented in Section~\ref{section:images}.

Spectra were extracted from each observation separately, and simultaneously fit using \textit{XSPEC} v12.6.0\footnote{http://xspec.gsfc.nasa.gov}.
\textit{specextract} was used to extract spectra and response files from regions of interest.  The regions used and resulting spectra are discussed in Section~\ref{section:spectra}.  When discussing diffuse emission, each larger region excludes events contained in the interior regions, as well as any counts from point sources detected within the diffuse emission.  Spectra were grouped using \textit{dmgroup}, with the minimum number of counts per bin listed in Table 2.  All uncertainties quoted on model parameters are to 90\% confidence level.  For calculations requiring a distance estimate, we assume $D=6.1D_{6.1}$~kpc~\citep{Kothes2003}.

CTB~87 was also observed on 1995 May 29 with \textit{ASCA} (sequence\# 53041000).  It was imaged with the GIS for 53 ks but was out of the field of view of the SIS.  Data was processed with \textit{XSELECT} (HEASoft v5.3.1).  We present the \textit{ASCA} results only for the purpose of comparing the spectral result for the PWN to the global fit with \textit{Chandra}.
The X-ray observations of CTB~87 presented in this paper are summarized in Table~\ref{table:observations}.

\section{Imaging}
\label{section:images}

\subsection{Nebula Imaging}
\label{section:nebulaimage}

Figure~\ref{figure:pwn} (left panel) shows the combined \textit{Chandra}\ observations 1037 and 11092, after removing point sources, filling the holes, and exposure correcting the image, as described in Section~\ref{section:obs}.  The energy range covered is 0.3$-$7.0~keV
(events below 0.3~keV and above 7.0~keV were excluded as the background dominates at these energies)
and the net exposure time is 87.1~ks.
The image was smoothed with a 2D Gaussian ($\sigma$~=~3\arcsec) and 21~cm contours are overlaid, using data acquired with the \textit{Dominion Radio Astrophysical Observatory} as part of the Canadian Galactic Plane Survey~\citep{Taylor2003}.
 The entire image is 525\arcsec~$\times$~525\arcsec.  The PWN is $\sim$200\arcsec~$\times$~300\arcsec\ overall with brighter emission to the southeast.  The
total observed X-ray flux of the nebula is 1.8~$\times$~10$^{-12}$~erg~cm$^{-2}$~s$^{-1}$ ($\sim$0.3~cts~s$^{-1}$).
The peak of the X-ray emission is seen in the southeast of the nebula, with a point source clearly visible at the peak (Section~\ref{section:ptimage}).  The X-ray peak is located near the southeastern edge of the bright radio nebula, offset from the center of the nebula, with the diffuse emission extending 100\arcsec, 100\arcsec, 250\arcsec, and 60\arcsec\ from the peak to the SE, NE, NW, and SW, respectively.  We observe arcs at the edges of the X-ray nebula (highlighted on Figure~\ref{figure:pwn}), extending from the southeast near the point source, back toward the radio peak.
The nature of these features is discussed in Section~\ref{section:disc}.

Figure~\ref{figure:pwn} (right panel) is a color image created by dividing the data into three bands,  exposure correcting each, and smoothing ($\sigma$~=~8\arcsec) each before combining.  The red, green and blue colors
correspond to the 0.3$-$2.0~keV, 2.0$-$4.0~keV, and 4.0$-$7.0~keV bands, respectively.
The contrast of each image was adjusted to reduce the background and bring out the PWN emission.  We see the hardest (bluest) emission near the bright point source (the putative pulsar, Section~\ref{section:ptimage}).  The arcs along the edge of the PWN (highlighted on Figure~\ref{figure:pwn}, left panel) appear as possible bow shock structures and are real:
the signal-to-noise ratio of the eastern arc is 12.8 while the signal-to-noise ratio of the western arc is 11.5.
In Section~\ref{section:disc} we discuss the nature of these features.

\subsection{Compact Object Imaging}
\label{section:ptimage}

Thanks to \textit{Chandra}'s superb resolution, a point source - the putative pulsar - is found at the southeastern edge of the bright radio nebula  with the following coordinates: $\alpha$(2000)=20$^{\rm{h}}$16$^{\rm{m}}$09\fs2, $\delta$(2000)=+37\arcdeg11\arcmin10\farcs5.
We refer to this new \textit{Chandra}\ point source as CXOU~J201609.2+371110.\footnote{Despite the lack of pulsations from this source, we refer to it hereafter as `pulsar' or `neutron star'.}
A point spread function (PSF) was created using \textit{mkpsf} for an energy of 3~keV (characteristic of the source's energy histogram) and an off-axis angle of 2\farcm38 (the location of the point source in observation 11092) and found to have a FWHM of 1\farcs15.  We note that this is the positional uncertainty of the point source.
Figure~\ref{figure:pointsource} (left) shows a 20\arcsec~$\times$~20\arcsec\ region, centered on the point source and combining both \textit{Chandra} observations (as for Figure~1).  This exposure corrected image was smoothed using a Gaussian with $\sigma$~=~1\arcsec\
and covers the energy range 0.3$-$7.0 keV.  We see evidence of an arc-like (torus) structure to the northeast of the point source with $\sim$10\arcsec\ diameter, and jet-like features extending $\sim$5\arcsec\ to the northeast and southwest.  Such
high-resolution features are expected to form close to the termination shock, due to  the deposition of the neutron star's wind energy into its surroundings~\citep[see e.g.][]{Gaensler2006}.

Figure~\ref{figure:pointsource} (right) shows a radial profile created from the observational data,
along with a radial profile of the PSF.  The data is consistent with a point source to a radius of $<$2\arcsec, beyond which the data clearly has a higher surface brightness than a point source.  This excess emission we term the `compact nebula' and is comprised of the candidate torus and jets, as seen in the left panel of Figure~\ref{figure:pointsource}.  We also tested the source extent using \textit{srcextent}.  The source observed size was 2\farcs05, the PSF observed size was 1\farcs15, and the estimated intrinsic size was 1\farcs69, confirming that there is an extended source surrounding the point source.

The PSF described above was then used as a convolution kernel in $Sherpa$\footnote{http://cxc.harvard.edu/sherpa4.4/index.html} to model the point source emission and determine the morphology of the excess emission.  Figure~\ref{figure:residuals} shows the data, the model, and the residuals resulting from the fit. The model consisted of a constant background level plus a 2D Gaussian with a FWHM of 1\farcs1 and an amplitude of 117 counts/pixel.  We see a large excess
in the residual map
to the northeast of the point source location, seeming again to have the torus-jet morphology as suggested in Figure~\ref{figure:pointsource}.

\section{Spectroscopic Analysis}
\label{section:spectra}

\subsection{Compact Object and Compact Nebula's Spectra}
\label{section:ptspectra}

To extract the spectrum of the point source, we used a circular region centered at CXOU J201609.2+371110 and with a radius~=~2\arcsec.
The background spectrum was extracted using an annulus surrounding the source with radius~=~2\arcsec$-$4\arcsec.
Given the small number of counts (see Table~2), we freeze the column density, $N_{\rm{H}}$, to the best value of 1.4~$\times$~10$^{22}$~cm$^{-2}$ obtained from fitting the total diffuse emission  (see Section~\ref{section:nebulaspectra}).
We attempted fitting with an absorbed power law model and an absorbed blackbody model, as would be expected from a neutron star.
Fitting with the absorbed power law model
yields a hard photon index of $\Gamma$~=~1.1~(0.7$-$1.6)
 and an absorbed flux of $F_{\rm{abs,~0.3-10.0~keV}}$~=~5.2~$\times$~10$^{-14}$~erg~cm$^{-2}$~s$^{-1}$ ($\chi^2_{\nu}$~=~0.91; see Table~2).
 The corresponding X-ray luminosity is $L_{\rm{0.3-10.0~keV}}$~=~3.0~$\times$~10$^{32}$~$D_{6.1}^2$~erg~s$^{-1}$ or $L_{\rm{2-10~keV}}$~=~2.5~$\times$~10$^{32}$~$D_{6.1}^2$~erg~s$^{-1}$.  The left panel of Figure~\ref{figure:pointspectrum} shows the data and fitted model, with the spectra from the two observations grouped together for plotting purposes only.  We note that allowing the column density to vary gives a column density of $N_{\rm{H}}$~=~1.6~(0.2$-$3.0)~$\times$~10$^{22}$~cm$^{-2}$ ($\chi^2_{\nu}$~=~0.93), consistent with that of the diffuse nebula, but due to the low number of counts the parameters are much less constrained.

Fitting with an absorbed blackbody model we find a similarly suitable fit
($\chi^2_{\nu}$~=~0.99,
Figure~\ref{figure:pointspectrum}, right) with a temperature $kT$~=~1.1~(0.8$-$1.5)~keV, an absorbed flux of $F_{\rm{abs,~0.3-10.0~keV}}$~=~3.2~$\times$~10$^{-14}$~erg~cm$^{-2}$~s$^{-1}$, and a luminosity $L_{\rm{0.3-10.0~keV}}$~=~1.8~$\times$~10$^{32}$~$D_{6.1}^2$~erg~s$^{-1}$ or $L_{\rm{2-10~keV}}$~=~1.6~$\times$~10$^{32}$~$D_{6.1}^2$~erg~s$^{-1}$.  However, the temperature is unreasonably high and, when combined with the luminosity, implies a neutron star radius of $\sim$0.03~km which is unreasonably small.

Attempting a two component model (power law plus blackbody) was not possible given the poor statistics.
Therefore we favor the power law model fit, which is also consistent with the X-ray emission seen from other neutron stars powering PWNe, such as PSR~J2022+3842 in G76.9+1.0~\citep[$\Gamma$~=~1.0~$\pm$~0.2,][]{Arzoumanian2011}, and PSR~J1833-1034 in G21.5-0.9~\citep[$\Gamma$~=~1.14~$\pm$~0.07,][]{Matheson2010}.

We then studied the emission from the compact
nebula by extracting a spectrum from an annulus with radius 2\arcsec$-$10\arcsec, using a larger annulus with radius 10\arcsec$-$20\arcsec\
 as the background spectrum.  This compact nebula, like the point source, has a hard spectrum with $\Gamma$~=~1.2~(0.9$-$1.4) and an X-ray luminosity of  $L_{\rm{0.3-10.0~keV}}$~=~8.7~$\times$~10$^{32}$~$D_{6.1}^2$~erg~s$^{-1}$ or $L_{\rm{2-10~keV}}$~=~6.5~$\times$~10$^{32}$~$D_{6.1}^2$~erg~s$^{-1}$.  The best-fit spectral parameters are summarized in Table~\ref{table:spectra} and the compact nebula spectrum is shown in Figure~\ref{figure:compactnebspectrum}.

\subsection{Diffuse Nebula's Spectrum}
\label{section:nebulaspectra}

To study the properties of the extended diffuse emission, we extracted a spectrum of the extended nebula shown in Figure~\ref{figure:pwn} (right image), omitting the inner 10\arcsec\ radius compact nebula.
A power law model provides a good fit, as would be expected from a PWN,
with the following best-fit parameters (see Table~2):   $N_{\rm{H}}$~=~1.38~(1.21$-$1.57)~$\times$~10$^{22}$~cm$^{-2}$,   $\Gamma$~=~1.68~(1.54$-$1.84),
$\chi^2_{\nu}$~=~1.01, $F_{\rm{abs,~0.3-10.0~keV}}$~=~2.1~$\times$~10$^{-12}$~erg~cm$^{-2}$~s$^{-1}$,
$L_{\rm{0.3-10.0~keV}}$~=~1.6~$\times$~10$^{34}$~$D_{6.1}^2$~erg~s$^{-1}$, and $L_{\rm{2-10~keV}}$~=~9.6~$\times$~10$^{33}$~$D_{6.1}^2$~erg~s$^{-1}$.

The above fit is consistent with the results we find from archived \textit{ASCA} GIS data.  The point source in CTB~87 was not resolved in the \textit{ASCA} observation and so we fit the entire nebula, a circular region with a radius of 3\farcm4, with a background spectrum extracted from the same detector, outside the diffuse emission.  Using GIS-2 and GIS-3 data (CTB~87 was not visible in the SIS data), we accumulate 3213 counts after background subtraction (compared to 10136 net counts with \textit{Chandra}).  Fitting GIS-2 and GIS-3 data simultaneously with a power law model in the 0.8$-$8.0 keV energy band, we find that $N_{\rm{H}}$~=~1.14~(0.94$-$1.36)~$\times$~10$^{22}$~cm$^{-2}$, $\Gamma$~=~1.77~(1.63$-$1.92), $\chi^2_{\nu}$~=~0.95, $F_{\rm{abs,~0.8-8.0~keV}}$~=~2.2$\times$~10$^{-12}$~erg~cm$^{-2}$~s$^{-1}$, $L_{\rm{0.3-10.0~keV}}$~=~1.4~$\times$~10$^{34}$~$D_{6.1}^2$~erg~s$^{-1}$, and  $L_{\rm{2-10~keV}}$~=~1.1~$\times$~10$^{34}$~$D_{6.1}^2$~erg~s$^{-1}$ (see Figure~\ref{figure:pwnspectrum}, left).

Next, we perform a spatially resolved spectroscopic study of the extended nebula.
To look for variations in the photon index without an artificially varying column density while fitting smaller regions with fewer counts,
we freeze $N_{\rm{H}}$ to the value from the overall diffuse fit.
Table~\ref{table:spectra}
shows the results of fitting a power law model to the regions shown in Figure~\ref{figure:pwn} (right) and Figure~\ref{figure:pwnspectrum} (right) shows the power law fit to each region (obs. 11092 and 1037 were fit simultaneously and are grouped together for plotting purposes only).  Beginning at the point source and moving outward across the PWN, we observe an overall increase (or steepening) in the photon index.  This spectral steepening is expected as we move away from a neutron star due to
expansion and synchrotron losses~\citep{Gaensler2006}.
We also explored background regions within the diffuse emission for the inner and mid regions.  We find a similar trend in the photon index, however, with a lower net flux and with the parameters being
 less constrained.
The brightening along the edges of the diffuse PWN is well fit by a power law with $\Gamma$~=~1.5~(1.2--1.9) for the eastern arc and $\Gamma$~=~1.7~(1.3--2.1) for the western arc.

We have searched for thermal emission that would arise from hot plasma in the SNR or interaction with the reverse shock.
We do not find evidence for any significant emission.
In particular, the residuals (lower panels of Figure~\ref{figure:pwnspectrum}) do not show any significant excess at low energies that may indicate the presence of thermal emission in the remnant.
However, this may be due to the limited sensitivity of the observation and field coverage.
A deeper, more sensitive, X-ray observation covering the entire radio nebula is needed to search for any faint thermal X-ray emission. This is further discussed below.

\section{Discussion}
\label{section:disc}

\subsection{Age}
\label{section:age}
The SNR's age is normally determined using the thermal emission properties of the X-ray emitting plasma associated with the SNR shell. In particular, knowledge of the shock velocity ($v_s$) and the actual SNR size ($R_s$)
yields a direct measurement of the age as: $t$~=~$\eta$~$\frac{R}{v_s}$, where $\eta$ is a factor that depends on the evolutionary phase of expansion of the SNR.
 For CTB~87, no SNR shell has been detected (as mentioned above).

In the absence of a SNR shell detection, we refer next to the offset between the radio and X-ray emission from the PWN, which can be attributed to the neutron star's motion since its birth following the supernova explosion.
The peak of the 21~cm emission is located at $\alpha$(2000)=20$^{\rm{h}}$16$^{\rm{m}}$04\fs4, $\delta$(2000)=+37\arcdeg12\arcmin31\farcs6,
$\sim$100\arcsec\ from the X-ray peak (or CXOU J201609.2+371110).  At a distance of $D=6.1D_{6.1}$~kpc~\citep{Kothes2003}, this separation translates to $2.9D_{6.1}$~pc perpendicular to the line of sight.
Fig. 4b of \citet{Hobbs2005} shows the majority of pulsars with an age $<$3~Myr have a two-dimensional speed in the range $100-500$~km~s$^{-1}$.  From this range, we estimate an age of $\sim$5$-$28~kyr
for CTB~87, assuming the pulsar was born at the location of the radio peak and is now located at the
 X-ray peak.
This is a reasonable assumption particularly since the X-ray emission from a PWN is associated with synchrotron radiation from the freshly injected particles,
 while the radio emission characterizes the older population of synchrotron emitting electrons over the neutron star's lifetime.  Using the average 2D speed of 307 km~s$^{-1}$ for pulsars with an age $<$3~Myr~\citep{Hobbs2005}, we estimate an age of $\sim$9~kyr for CTB~87.
These age estimates are comparable to the 11~kyr age estimated for Vela~X~\citep{LaMassa2008} and the 18~kyr age estimated for G327.1$-$1.1~\citep{Temim2009}, PWNe which we compare CTB~87 to further below.

The advanced age for CTB~87 (at least an order of magnitude higher than young PWNe like the Crab and G21.5$-$0.9) is also hinted at by the large radio size.
The central 8$^{\prime}$$\times$6$^{\prime}$ radio component implies a physical size of $\sim$(14~$\times$~11)$D_{6.1}$~pc.
The extended 16$^{\prime}$ diffuse component (Kothes et al. 2003; Kothes et al. in preparation) implies an even larger pulsar wind nebula of $\sim$28$D_{6.1}$~pc in extent. This is slightly smaller, but comparable in
size to
the unusual PWN G76.9+1.0 (29~pc $\times$ 35~pc) which has been recently argued to be an evolved PWN with an age $<$ 36~kyr and which harbors a newly discovered energetic pulsar with a 9~kyr characteristic
age (Arzoumanian et al. 2011).

In the remainder of this discussion, we will assume an age of $t = 10 t_{10}$~kyr for CTB~87 where calculations require an age estimate, keeping in mind that the age is uncertain given the estimates above (see also an estimate of the putative pulsar's characteristic age in Section~\ref{section:pulsarproperties}).
Moreover, if we take into account the even larger offset (an additional 30$^{\prime\prime}$) between the peak of the \textit{extended} radio emission and the peak of the X-ray emission,
the age would be 30\% larger for the same assumed velocity. This will be further discussed in the forthcoming radio paper  (Kothes et al. in preparation).
Additionally, if the offset is due to a reverse shock (Section~\ref{section:reverseshock}) and not related to the proper motion, this estimate may not reflect the age of the system.

We should note that there is accumulating evidence for pulsars' spin axes to be roughly aligned with the direction of their motion, albeit with a small number of systems for which 3D information is available~\citep{Ng2007}.
\citet{Noutsos2012} provide further evidence that alignment between projected spin-axes and pulsar proper motion is more likely than orthogonality.
In Figure~\ref{figure:pointsource} (left), we see the possible jets (believed to define the spin-axis) are nearly perpendicular to the apparent proposed proper motion of the putative pulsar.
If our assumed pulsar's direction of motion and jets' direction are correct, then CTB~87 would be a system that doesn't follow the expected alignment from the above-mentioned studies.
A deeper \textit{Chandra} observation and proper motion studies of the X-ray source will confirm this.

\subsection{Putative Pulsar \& Pulsar Wind Properties}
\label{section:pulsarproperties}
To the best of our knowledge, no pulsations have been reported from CTB~87~\citep{Gorham1996, Biggs1996, Lorimer1998}.  We can however infer the properties of the neutron star powering CTB~87 from the X-ray luminosity of the nebula.
The unabsorbed flux of CTB~87 in the $2-10$~keV band is $F_{\rm{unabs,~2-10~keV}}=2.4 \times 10^{-12}$~erg~cm$^{-2}$~s$^{-1}$, implying a luminosity of $L_{\rm{X,~2-10~keV}}=1.1\times10^{34}~D_{6.1}^2$~erg~s$^{-1}$.  We can estimate the spin-down energy loss of the neutron star powering the PWN from the PWN luminosity using the empirical relationship $L_{\rm{X,~2-10~keV}}=10^{-19.6} \dot{E}^{1.45}$~\citep{Li2008}.  This yields an estimate of $\dot{E}=9.6 \times 10^{36}~D_{6.1}^{1.38}$~erg~s$^{-1}$, consistent with the estimate of $\dot{E}=6.9 \times 10^{36}~D_{6.1}^{1.49}$~erg~s$^{-1}$ obtained with the empirical relationship $L_{\rm{X,~2-10~keV}}=1.34 \log \dot{E} -15.34$~\citep{Possenti2002}.  The efficiency in converting the spin-down luminosity to synchrotron emission is then $\eta_X \equiv L_X/\dot{E} \approx 0.001$, which is within the range observed for PWNe powered by rotation-powered pulsars.
The PWN luminosity can also be used to estimate the characteristic age ($\tau$) of the neutron star using $L_{\rm{X,~2-10~keV}}=10^{42.4} \tau^{-2.1}$~\citep{Li2008}.  We estimate an age of $\tau=9800~D_{6.1}^{-0.95}$~years.

Assuming an age of $t~=~10t_{10}$~kyr (keeping in mind the large uncertainty on this estimate, noted in Section~\ref{section:age}), $\dot{E} \sim 10^{37} \dot{E}_{37}$~erg~s$^{-1}$, a braking index of $n$=3 (magneto-dipole braking), a moment of inertia of $I$~=~10$^{45} I_{45}$~g~cm$^2$, and a neutron star radius of $R$~=~10$R_{10}$~km, we can estimate the properties of the putative pulsar~\citep[equations described for e.g. in][]{Gaensler2006}.  The period of the pulsar is estimated as $P=\left( \frac{4 \pi^2 I}{(n-1)\dot{E}\tau} \right)^{1/2}$ $\approx 0.08~I_{45}^{1/2}~\dot{E}_{37}^{-1/2}~t_{10}^{-1/2}$~seconds.  The corresponding period derivative is $\dot{P} = \frac{P}{(n-1)\tau} \approx 1.3 \times 10^{-13}~I_{45}^{1/2}~\dot{E}_{37}^{-1/2}~t_{10}^{-3/2}$ and the corresponding equatorial surface dipole magnetic field is $B = 3.2~\times~10^{19}~I_{45}^{1/2}~R_{10}^{-3}~( P \dot{P} )^{(1/2)} \approx 3.2 \times 10^{12}~I_{45}~\dot{E}_{37}^{-1/2}~R_{10}^{-3}~t_{10}^{-1}$~G,
within the range of the field expected for pulsars with prominent PWNe~\citep[$1\times10^{12} - 5\times10^{13}$~G,][]{Gaensler2006}.

We expect to observe a termination shock around a rotation-powered neutron star, at a radius $r_t=\sqrt{\dot{E}/(4 \pi c P_{\rm{PWN}})}$; where $P_{\rm{PWN}}$ is the nebular pressure.
 Figure~\ref{figure:pointsource} suggests that the arc seen at a radius of $\sim$5\arcsec\  corresponds to the termination shock.
 At an assumed distance of 6.1~kpc, this corresponds to a distance of 0.15~pc from the point source.
 This estimate is in agreement with the typical scale of 0.1~pc observed in other PWNe~\citep{Gaensler2006}.
Using the spin-down energy inferred above, and equating the expression of the termination shock site to 5\arcsec,  we estimate the pressure in the PWN to be $P_{\rm{PWN}}\sim 1.2 \times 10^{-10} D_{6.1}^{-0.62}$~dyne~cm$^{-2}$.  We can then estimate the magnetic field in the nebula as $B_n=(8 \pi P_{\rm{PWN}})^{1/2}\sim 55 D_{6.1}^{-0.31}~\mu$G.
This estimate is in agreement with the preliminary value inferred from the radio study assuming equipartition between particle energy and magnetic field energy (Kothes et al. in preparation).

\subsection{Ambient Density \& the Missing SNR Shell}
The absence of a SNR shell detection implies expansion into a low-density medium.
To infer the upper limit on the
density of the medium that the pulsar wind nebula is expanding into (which would be the shocked ejecta for a Vela-like SNR),
we estimate next the maximum thermal contribution to the diffuse emission allowed by the data by adding a thermal component to the power law model.
 We freeze the best fit power law parameters (Table~\ref{table:spectra}, `Total Diffuse' region), and add an Astrophysical Plasma Emission Code (APEC) thermal component~\citep{Smith2001}.
We increase the emission measure ($EM$) until the model level under 1~keV exceeds the data by 2$\sigma$, finding a conservative upper limit for the normalization on the thermal component of
9~$\times$~10$^{-4}$~cm$^{-5}$ (for $kT$~=~0.5~keV).
This implies that $\int n_e n_{\mathrm H} dV \sim f n_e n_{\mathrm H} V$~=~10$^{14}$ $\left( 4 \pi D^2\right) \left( 9~\times~10^{-4}~\rm{cm}^{-5} \right)$~$<$~4.0~$\times$~10$^{56}$~$D^2_{6.1}$~cm$^{-3}$,
where $f$ is the volume filling factor, $n_e$ is the post-shock
electron density, and $n_{\mathrm H} \sim n_e / 1.2$.  Using a volume
of 5.6~$\times$~10$^{57}$~$D^3_{6.1}$~cm$^3$ for the PWN, $n_e < 0.29 f^{-1/2}~D^{-1/2}_{6.1}$~cm$^{-3}$.
If a SNR shell exists beyond the detected PWN, the shocked ISM must have a density lower than that estimated above for the ejecta, $0.29 f^{-1/2}~D^{-1/2}_{6.1}$~cm$^{-3}$ = $0.92~D^{-1/2}_{6.1}$~cm$^{-3}$ (for a volume filling factor $f \sim 0.1$).  This then implies an upper limit on the upstream ambient density of $n_0 < 0.2~D^{-1/2}_{6.1}$~cm$^{-3}$ ($\rho_0 < 3 \times 10^{-25}~D^{-1/2}_{6.1}$~g~cm$^{-3}$ where $n_0$ includes only hydrogen), assuming a compression ratio of 4.  The inferred
low density supports the lack of detection of the SNR shell or X-ray emitting ejecta,
and is suggestive of expansion into a stellar wind-blown bubble \citep{Weaver1977}.
This upper limit is consistent with estimates obtained by extracting the background spectrum (outside the PWN) and assuming that the thermal emission from a putative shell is buried underneath the background.  From these background spectra, we estimate an upper limit of $9 \times 10^{-14}$~erg~cm$^{-2}$~s$^{-1}$~arcmin$^{-2}$ on the observed X-ray flux of the putative SNR shell.

The above density upper limit is similar to density estimates for G21.5-0.9~\citep{SafiHarb2001, Matheson2010}, whose shell was revealed after accumulating $\sim$0.5~Msec using \textit{Chandra}.  A deep X-ray observation of CTB~87 to search for the missing SNR shell is therefore warranted.
Additionally, future studies of the surrounding environment with CO or HI data may help to further constrain the density.  This scenario will be revisited and further discussed in the follow-up radio paper (Kothes et al. in preparation).

Assuming that the remnant is evolving in a uniform medium, the reverse shock will reach the center of the SNR in a time $t_{\rm{Sedov}} \approx 7 (M_{ej}/10~M_{\sun})^{5/6} (E_{\rm{SN}}/10^{51}~\rm{ergs})^{-1/2} (n_0/1~\rm{cm}^{-3})^{-1/3}$~kyr.  For an ejected mass of
a few solar masses, an explosion energy of $10^{51}$~ergs, and a number density of $n_0 < 0.2$~cm$^{-3}$, this is greater than 4~kyr for CTB~87, consistent with our estimate above for the age (see Section 5.1).

For a remnant in the Sedov stage, the SNR radius can be estimated as $R_{\rm{SNR}} \approx 1.15 \left( \frac{E_{\rm{SNR}}}{\rho_{\rm{ISM}}} \right)^{1/5} t^{2/5}$~\citep{vanderSwaluw2003}.
Again assuming an explosion energy of $E_{\rm{SNR}} = 10^{51} E_{51}$~erg, an ISM density of $n_0 = 0.2~n_{0.2}$~cm$^{-3}$, and an age of $t = 10 t_{10}$~kyr (large uncertainty noted in Section~\ref{section:age}), we estimate the SNR radius to be $R_{\rm{SNR}} \approx 18 \left( \frac{E_{51}}{n_{0.2}} \right)^{1/5} t_{10}^{2/5}$~pc,
which is larger than the revised size of the PWN ($\sim$28~pc diameter at 6.1~kpc; see Section 5.1).
We conclude that all of the above mentioned arguments point to CTB~87 being an evolved PWN expanding into a low density medium, likely due to the stellar wind bubble blown by the progenitor star.

\subsection{Morphology Interpretation}
In Section 5.2 we have interpreted the compact $\sim$5\arcsec\ nebula shown in Figure~2 as a wind-blown bubble associated with the deposition of the neutron star's wind energy into the surroundings.
As shown in Figure~1, in addition to this compact nebula, the X-ray emission from CTB~87 is characterized by a cometary-like structure extending $\sim$250\arcsec\  northwest from CXOU~J201609.2+371110
and bridging the X-ray and radio emission peaks.
As mentioned earlier, this extended emission has some brightening along its edges (see Figure~1).
This morphology may be interpreted in two ways.  First, these arcs may be a bow shock structure due to the motion of the pulsar to the southeast through the SNR interior, leading to the formation of
a PWN with a hybrid morphology (i.e. a wind-blown bubble and a bow shock nebula).
Alternatively, the reverse SNR shock may be running into, and crushing, the PWN.  We discuss next these two scenarios.

\subsubsection{Bow Shock Interpretation}

A pulsar given a kick at birth that has sufficient velocity will move out of the PWN and create a bow shock.  As the pulsar moves outward in the SNR, away from the supernova explosion site, the sound speed in the shocked ejecta decreases and the pulsar's motion becomes supersonic.  The ram-pressure from the pulsar's motion confines the nebula and the PWN evolves to a cometary appearance, with the bow shock forming around the head of the PWN after the pulsar has crossed $\sim$68\% of the SNR shell, at half of the crossing time~\citep{vanderSwaluw2004,Gaensler2006}.

The pulsar wind in a bow shock PWN decelerates at a termination shock where the pressure balance is due to ram pressure from the neutron star's motion.  In the direction of the star's motion, the termination shock radius, $R_{\omega 0}$, is defined by
$\dot{E}/\left(4 \pi \omega R^2_{\omega 0} c\right) = \rho_0 v^2_{\mathrm{PSR}}$ ($\omega = 1$ for an isotropic wind), and is located at approximately half of the observed bow shock radius~\citep{Gaensler2006}.

Ram-pressure confinement in this hybrid morphology scenario requires that the ram pressure due to the pulsar's proper motion in the SNR environment,  $\rho_0 v^2_{\mathrm{PSR}}$, exceeds the nebular pressure, $P_{\mathrm{PWN}}$, estimated in Section~5.2 ($\rho_0$ is the density of the medium in which the PSR is moving).
Assuming that the brightening along the PWN edge in Figure~\ref{figure:pwn} is due to a bow shock, $R_{\omega 0} \sim 1.8 \times 10^{18}~D_{6.1}$~cm.
 Using our inferred value for $\dot{E}$ (Section 5.2), we find that $\dot{E}/\left(4 \pi \omega R^2_{\omega 0} c\right)~\sim~8 \times 10^{-12}$~dyne~cm$^{-2}$.
This is only $\sim$7\% of the pressure $P_{\mathrm{PWN}} = 1.2 \times 10^{-10}$~dyne~cm$^{-2}$ estimated from the termination shock of the more compact nebula (Section~5.2), a ratio that is consistent with $(\frac{r_t}{R_{\omega 0}})^2$ .
 Furthermore, for  $\rho_0 v^2_{\mathrm{PSR}}$  to exceed this nebular pressure of $\sim$10$^{-10}$~dyne~cm$^{-2}$, $v_{\mathrm{PSR}}$ should exceed $\sim$170~km~s$^{-1}$ (for an assumed ambient density of 0.2~cm$^{-3}$, Section 5.3).
As the formation of a prominent bow shock requires the nebular pressure to be exceeded by a large amount, the above suggests a high pulsar velocity is needed for creation of a bow shock and argues against the cometary morphology being caused by ram-pressure confinement.

\subsubsection{Reverse Shock Interaction Interpretation}
\label{section:reverseshock}
In the following, we discuss the observed morphology in the light of an evolved PWN encountering the SNR reverse shock. This is motivated by the age inferred above (Section~5.1) and the offset observed between the radio
and X-ray nebulae, as seen in other evolved PWNe further described below.
A detailed description of the evolution of a PWN in a SNR can be found in \citet{Blondin2001},  \citet{vanderSwaluw2004}, and \citet{Gaensler2006}.
We here briefly summarize.

As the expanding SNR sweeps up sufficient mass from the ISM, it evolves into a Sedov-Taylor phase. A reverse shock develops, running back into the supernova ejecta and eventually crushing the PWN, causing it to rebound and oscillate several times.
The timescale for the reverse shock to interact with the entire PWN is a significant fraction of the PWN life~\citep{vanderSwaluw2004}.  After these oscillations fade, the PWN again expands, now into hot, shocked ejecta.  The pulsar may have traveled a distance comparable or larger than the equivalent spherical PWN around a stationary pulsar.  The pulsar leaves its original wind bubble and generates a new smaller PWN at its current location, off-center within the PWN.
Therefore, after the passage of the reverse shock the PWN consists of (i) a \textit{relic} PWN from the initial energetic stage of the pulsar wind, appearing as a central (possibly distorted) radio PWN with little corresponding X-ray emission, and (ii) a head containing the pulsar, directed toward the SNR shell, and emitting in X-rays.

\Citet{vanderSwaluw2004} express an upper limit on the age, $t_{col}$, at which the reverse shock collides with the entire PWN surface as $t_{col} = 1045 E_{51}^{-1/2} \left( \frac{M_{ej}}{M_{\sun}} \right)^{5/6} n_0^{-1/3}$~yr where $E_{51}$ is the total mechanical energy of the SNR in units of $10^{51}$~erg, $M_{ej}$ is the ejected mass, and $n_0$ is the ambient hydrogen number density.  For CTB~87 we assume (as above) the explosion energy is $10^{51}$~erg, that the ejected mass is
a few solar masses, and that the ambient density is $<$0.2~cm$^{-3}$. We find that $t_{col} > 4$~kyr for CTB~87, at the low end of the estimated age (Section~5.1).
This suggests that the reverse shock has potentially encountered the PWN, and is the cause of the observed morphology.
 The supernova explosion may also have been asymmetric, or the ambient density may be non-uniform, leading to a more rapid evolution in one direction.
 This is to be further discussed in the forthcoming radio paper (Kothes et al., in preparation).

This conclusion is further supported by the TeV emission detected from CTB~87~\citep{Aliu2011}.
Focusing on \textit{Chandra}-detected PWNe near extended TeV sources,  \citet{Kargaltsev2008} find that faint extensions of compact PWNe point toward offset centers of TeV emission, supporting the notion that these TeV sources are crushed PWNe.  Recent \textit{VERITAS} studies of CTB~87 indicate that the TeV emission overlaps with the radio emission, further supporting the picture of a `relic' PWN (E. Aliu and M. Roberts 2012, private communication).
The lack of X-ray emission from the SNR ejecta or SNR shell remains however a question to be addressed with deeper observations.

\subsubsection{Comparison to other evolved PWNe}
G327.1$-$1.1 is a composite remnant with a symmetric radio shell and off-center non-thermal component (PWN), similar in structure to CTB~87.  There is a radio finger extending from the PWN, at the end of which is a compact X-ray source that may be the pulsar~\citep{Temim2009}.  \textit{Chandra}\ images show a compact source embedded in a cometary structure that trails back toward the radio PWN, as well as two prong-like structures extending into a large bubble in the opposite direction from the PWN~\citep{Temim2009}.
We see a similar structure in CTB~87 (Figure 1) with the X-ray nebula offset from the radio nebula and most of the X-ray nebula trailing back toward the radio nebula.  We also see emission ahead of the point source in CTB~87, extending to the southeast away from the radio nebula.
\citet{Temim2009} similarly favor a reverse shock interaction scenario over a bow shock scenario.
G327.1$-$1.1 however has more detailed spectra available, and thermal X-ray emission has been detected with \textit{XMM-Newton}.
The SNR shell detection in G327.1$-$1.1 also allows further constraints to be placed on the SNR properties.
TeV emission from G327.1$-$1.1 has been detected by \textit{H.E.S.S.}, showing the $\gamma$-ray centroid is located between the radio centroid and the pulsar candidate~\citep{Acero2011}.

Vela-X (associated with PSR~B0833-45) is another example of a PWN whose morphology has been affected by its interaction with the reverse shock (Blondin et al. 2001).  The radio nebula extends south of the pulsar, the northern portion possibly having been disturbed by the reverse shock.  Vela-X has a complicated filamentary structure in X-rays, with thermal emission and enhanced O, Ne, and Mg abundances detected, signifying ejecta that has been mixed into the nebula by the propagation of the reverse shock~\citep{LaMassa2008}.
The explanation of Vela-X as a relic PWN was confirmed by \citet{Aharonian2006} with the detection of TeV emission that is offset from the pulsar.

\section{Conclusions}

We have presented the first detailed X-ray study of CTB~87, a large radio nebula classified as a PWN.  Thanks to \textit{Chandra}'s unprecedented resolution, we have detected a point source, CXOU~J201609.2+371110, the likely putative pulsar powering CTB~87,
inside a compact X-ray nebula that is $\sim$100\arcsec\ offset from the peak of the radio nebula.
The compact $\sim$5\arcsec-radius nebula is characterized by a torus/jet-like structure indicating its origin as the putative pulsar's wind-blown bubble.
In addition to the X-ray source and compact nebula, a more extended nebula was found to
trail back toward the peak of the radio emission.
A spatially resolved spectroscopic study of the nebula shows that the PWN is well fit by a power law model, with the overall spectrum steepening away from the point source.  We also observed brightening at the edges of the PWN, and discussed two possible interpretations:
a bow shock caused by the motion of the pulsar or
a PWN crushed by the SNR's reverse shock.
Based on the observed multi-wavelength morphology, timescale for the reverse shock to interact with the PWN, and estimated age of the PWN,
the reverse shock interaction scenario seems more plausible.
A deep X-ray observation is needed to search for thermal emission associated with the shocked ejecta and/or the SNR shell, as well as to confirm the high-resolution structures detected in the compact nebula.
Timing studies in the radio and X-ray will allow the search for pulsations from
CXOU~J201609.2+371110. Finally, in the long term, a proper motion study of the putative pulsar is needed to shed light on its motion within the radio nebula.

\acknowledgments

We thank the referee for useful comments on the manuscript.
The research presented in this paper has used data from the Canadian Galactic Plane Survey, a Canadian project with international partners, supported by the Natural Sciences and Engineering Research Council (NSERC),
and made use of the NASA's Astrophysics Data System and the High-Energy Astrophysics Science Archive Research maintained at NASA's Goddard Space Flight Center.
This research was primarily supported by NSERC. SSH acknowledges support from an NSERC Discovery Grant,  the Canada Research Chairs program, the Canadian Space Agency, the Canadian Institute for
Theoretical Astrophysics, and Canada's Foundation for Innovation.

{\it Facilities:} \facility{ASCA (GIS)}, \facility{CXO (ACIS)}, \facility{DRAO:Synthesis Telescope}.


\section{INTRODUCTION}\label{intro}

The Voyager 1 and 2 spacecraft have been mapping the structure of the
outer heliosphere as they leave the solar system. In 2005, V1 crossed
the termination shock \citep{stone05a,burlaga05a,decker05a}, where the
supersonic solar wind becomes subsonic, and has since been traversing
the HS. The HP, whose location and structure are unknown, separates
the magnetic field and plasma associated with the sun from that of the
LISM \citep{parker63a,baranov79a}.  The magnetic field in the HS has
remained dominantly in the azimuthal (east-west) direction given by
the Parker spiral but could rotate and acquire measurable north-south
and radial components upon crossing the HP. In ideal (non-dissipative)
models of the heliosphere, the local magnetic field is transverse to
the boundary and the HP is a tangential discontinuity
\citep{parker63a,baranov79a}. However, whether the HP is a smooth
interface, or breaks up due to instabilities, has been the subject of
substantial discussion in the literature
\citep{fahr86a,baranov92a,liewer96a,zank96b,swisdak10a}. The structure
of the HP, and in particular whether the boundary is porous to some
classes of particles, is of great importance because of its impact on
the transport of energetic particles into and out of the heliosphere.

Beginning on day 210 of 2012, the V1 spacecraft measured a series of
dropouts in the intensities of energetic particles produced in the
heliosphere: the Anomalous Cosmic Rays (ACRs) and the lower-energy
Termination Shock Particles (TSPs)
\citep{webber13a,stone13a,krimigis13a}. Simultaneous with the dropouts
were abrupt increases in the Galactic Cosmic Ray (GCR) electrons and
protons and an increase in the magnetic field intensity
\citep{burlaga13a}. Finally, on around day 238, the
heliospheric-produced particles dropped to noise levels and the GCRs
underwent a final increase.  Both have since exhibited no significant
variations, which suggests that V1 crossed the HP, with the repeated
dropouts and increases perhaps due to radial fluctuations caused by
changes in the solar wind dynamic pressure. However, during this time
the direction of the magnetic field remained dominantly azimuthal
\citep{burlaga13a}, consistent with the spacecraft remaining in the
HS. While MHD models of the heliosphere suggested that the rotation of
the magnetic field across the HP at the location of V1 would be small
\citep{opher09a,opher09b}, the lack of any significant change in the magnetic
field direction across the final transition on day 238 suggested that
V1 remained within the magnetic domain of the HS.

We present here a model of the magnetic structure of the HP at V1's
location that produces particle and magnetic signatures consistent
with the observations.  By pairing a global MHD simulation with a
local PIC simulation, we show that magnetic reconnection can produce a
complex, nested set of magnetic islands at the HP. Tongues of LISM
plasma penetrate into the HS along reconnected field lines. These
tongues correspond to local depletions of the HS plasma and
enhancements in the magnetic pressure. A key result of the simulations
is that sharp anti-correlated jumps in the HS and LISM number density
can occur across the separatrices emanating from reconnection sites
while the magnetic field undergoes essentially no rotation. Such
behavior undercuts the primary argument suggesting that V1 has not
crossed the HP -- that no field rotation was seen on day 238 where the
final drop in ACRs was measured \citep{burlaga13a}. We therefore
suggest that V1 actually crossed the HP on day 209, the time of the
last reversal in the azimuthal magnetic field $B_T$, and that the
steady values of the normal $B_N\sim 0.12\text{ nT}$ and $B_T\sim
-0.40\text{ nT}$ fields since that time are the draped interstellar field just
outside of the HP.

\section{MHD Simulation}\label{MHD}
To establish local conditions at the HP, we first explore the
heliosphere's large-scale structure with a global MHD simulation that
includes both neutral and ionized components (and both thermal and
pick-up ions in the solar wind) \citep{zieger13a}.  The LISM field,
$B_{\text{ISM}}$, has a magnitude of $0.44\text{ nT}$ and a direction
defined by $\alpha_{Bv}=15.9^{\circ}$ and $\beta_{Bv}=51.5^{\circ}$,
where $\alpha_{Bv}$ and $\beta_{Bv}$ are the angle between
$B_{\text{ISM}}$ and the velocity of the interstellar wind
$v_{\text{ISM}}$ and the angle between the
$B_{\text{ISM}}$-$v_{\text{ISM}}$ plane and the solar equator
\citep[for further discussion of these choices see][]{opher09a}. The
$Z$-axis is along the solar rotation axis and the $X$-axis is chosen
so that $v_{\text{ISM}}$ lies in the $X$-$Z$ plane.  The MHD
simulation did not include the sector zone (where the Parker spiral
field periodically reverses polarity due to the tilt between the solar
magnetic and rotation axes) since this leads to field reversals that
cannot be numerically resolved upstream of the HP and therefore
produces incorrect values of $B=|\mathbf{B}|$
\citep{opher11a,borovikov11a}. The solar field polarity corresponds
to solar cycle 24, with the azimuthal angle $\lambda$ (between the
radial and $T$ directions in heliospheric coordinates) $90^{\circ}$ in
the north and $270^{\circ}$ in the south.

In Fig.~\ref{MHDfig}, $B$ from the simulation reveals the solar wind
compression at the termination shock, the downstream HS, and the
HP. Profiles (solid curves in Fig.~\ref{cutsfig}) along the V1
trajectory of the pick-up ($n_{\text{pui}}$) and thermal
($n_{\text{th}}$) ion densities and the azimuthal ($B_T$) and normal
($B_N$) magnetic fields near the HP are inputs for the PIC
simulations.  $n_{\text{pui}}$ decreases from $\approx 7\times
10^{-4}/\text{cm}^3$ in the HS to zero in the LISM while
$n_{\text{th}}$ rises from $0.003/\text{cm}^3$ to $\approx
0.08/\text{cm}^3$.  $B_N$ (Fig.~\ref{cutsfig}C) is small at V1's
latitude in the LISM. $B_T$ flips direction across the HP, but remains
the dominant component on both sides of the boundary
(Fig.~\ref{cutsfig}D). The polar angle $\delta$ (the
angle between $B_N$ and the equatorial field) in the simulations approaches
$14^{\circ}$ just outside of the HP, which is consistent with the steady values
seen in the V1 data.

The MHD simulation does not match Voyager observations in several
respects --- the sign of the HS azimuthal magnetic field
orientation, the strength of the flows in the HS and the
characteristic scale length of the HP --- none of which is essential
for calculating initial conditions for the PIC simulations. First,
because V1 continued to measure sector boundaries in the HS during
2012, and therefore probably remained in the sector zone, the sign of
$B_T$ in the HS in the MHD model is irrelevant since a ``correct''
model should include the reversals associated with the sector region.
Second, in contrast to the simulation, indirect measurements by the V1
LECP instrument indicate little to no normal flow in the HS
\citep{decker12b}.  No published global model has explained the
observed flows, although simulations that include the sectored field
\citep[e.g.,][]{opher12a} are closer to the observations than those presented
here \citep[see also][for an alternative
explanation]{pogorelov12a}. Finally, since the MHD model does not
include the physics necessary to describe the structure of the HP, the
scale length of this transition is not physical, but is instead a
numerical artifact.  On the other hand, what is essential for input
into the PIC simulations is the strength of the HS field and the
strength and orientation of the field in the LISM.

\section{PIC Simulations}\label{PIC}
The initial profiles of the magnetic field, density, and temperature
for the 2-D PIC simulations (dotted lines in Fig.~\ref{cutsfig};
right-hand scale) were constructed with input from the MHD profiles
although, in keeping with the Voyager 1 observations, there are no
initial flows. The PIC code is written in normalized units based on a
field strength $B_0$ and density $n_0$ (lengths normalized to the ion
inertial length $d_i =c/\omega_{pi}$, with $\omega_{pi}$ the ion
plasma frequency, times to the ion cyclotron time $\Omega_{i0}^{-1}$
and velocities to the Alfv\'en speed $c_{A0}$). In the HS, thermal
($n_{\text{th}}=0.25n_0$, $T_{\text{th}}=0.25 m_ic_{A0}^2$) and
pick-up ions ($n_{\text{pui}}=0.01n_0$,
$T_{\text{pui}}=15.0m_ic_{A0}^2$) were included as independent species
while the LISM only included a thermal component ($n_{\text{th}} =
2.0n_0$, $T_{\text{th}}=0.2m_ic_{A0}^2$).  The simulations were
performed in a domain with dimensions $(L_T, L_R)= (409.6d_i,
204.8d_i)$. The ion-to-electron mass ratio was 25 and the velocity of
light was 15$c_{A0}$.  Not shown in Fig.~\ref{cutsfig} are the three
current sheets, of initial half-width $0.5 d_i$, that produce the
sectored HS field.  This scale reflects satellite measurements at the
Earth's magnetopause that such boundaries collapse to kinetic scales
\citep{sonnerup81a}. Pressure balance across each reversal is achieved
by adjusting the out-of-plane component $B_N$ \citep{smith01a}.

For HS-appropriate values, $n = 10^{-2} \text{
cm}^{-3}$ and $B = 0.3 \text{ nT}$, $d_i \approx 2\times 10^{-5}
\text{ AU}$, $\Omega_{ci}^{-1} \approx 30 \text{ s}$, and $c_{A}
\approx 100 \text{ km/s}$.  Resolving kinetic scales forces the
simulation domain to be much smaller than the actual system.  Despite
this limitation, the important physical processes can still be
understood by appropriately scaling the results \citep{schoeffler12b}.

The simulations are evolved with no initially imposed
perturbations. Because of the lower density, which leads to a locally
higher $c_A$ and effectively thinner current sheets (when normalized
to the local $d_i$), magnetic reconnection first starts in the
sectored HS.  Small magnetic islands grow on individual current layers
in the HS and merge to become larger islands until they are comparable
in size to the sector spacing (Fig.~\ref{PICfig}A). A chain of small
islands grows at the HP. These merge, forming larger islands, and are
compressed by HS islands pushing against the HP
(Fig.~\ref{PICfig}A). By late time, the HS magnetic field has
reconnected with that of the LISM, forming a complex, nested chain of
islands (Fig.~\ref{PICfig}A) at the HP with sizes comparable to the
original sector spacing. Along a cut from the HS to the LISM (dark
line in Fig.~\ref{PICfig}A) the HP is at $\Delta R/d_i\sim 25$ where
$B_T$ reverses sign (Fig.~\ref{PICfig}D). The rotation to the LISM
field direction is complete by $\Delta R/d_i\sim 35$ after which
$\delta$ and $\lambda$ are nearly constant.

We can independently track every particle in the PIC model and
therefore can explore the mixing of the LISM and HS plasmas.  The
overall result is a highly structured distribution in the densities of
the LISM ($n_{\text{LISM}}$) and HS ($n_{\text{HS}}$), with each
experiencing sharp jumps across the separatrices bounding the outflows
ejected from reconnection sites (Fig.~\ref{PICfig}A-C). Particles
initially in the LISM continue to dominate the density on the
un-reconnected LISM field lines, have mixed with HS particles in the
nested islands formed from HS-LISM reconnection, and are largely
excluded from islands formed from reconnection of the HS sectored
field (Fig.~\ref{PICfig}B). Particles initially in the HS dominate
islands resulting from reconnection of the sectored field, are mixed
with LISM particles in the HP islands, and are nearly excluded from
un-reconnected regions of the LISM (Fig.~\ref{PICfig}C).

Radial cuts through the simulation reveal that the increases and
decreases in $n_{\text{LISM}}$ and $n_{\text{HS}}$ are typically
anti-correlated (Fig.~\ref{PICfig}G). Moving from a pure HS magnetic
island into an island or outflow jet where LISM and HS plasma has
mixed reduces $n_{\text{HS}}$ and increases $n_{\text{LISM}}$. Along the
cut the first drop in $n_{\text{HS}}$ occurs downstream of a magnetic
separatrix, where HS particles have an open path to the LISM along
open field lines ($\Delta R/d_i \sim 8$ in
Fig.~\ref{PICfig}G). Similar behavior has been documented in satellite
measurements at the Earth's magnetopause \citep{sonnerup81a} and
echoes V1's observations of the anti-correlated variations in the
fluxes of ACRs/TSPs and galactic electrons/GCRs
\citep{stone13a,krimigis13a,webber13a}. Most important, the cuts
further reveal that, when crossing the last magnetic separatrix on the
LISM side of the HP before finally entering the pristine LISM plasma,
the sharp decrease in $n_{\text{HS}}$ and increase in
$n_{\text{LISM}}$ (in the interval $\Delta R/d_i = 38-50$ in
Fig.~\ref{PICfig}G) occur over an interval where there is no
directional change in the magnetic field (Fig.~\ref{PICfig}D-E). The
absence of a directional change in $\mathbf{B}$ at locations with
strong variations in the particle densities is consistent with one of
the most significant of the V1 observations \citep{burlaga13a}. The
fact that, in our simulation, this occurs on the LISM side of the HP
therefore suggests that it may be incorrect to conclude that V1 has
not crossed the HP.

The simulation cuts also reveal that local decreases in the HS density
typically correspond to increases in the local magnetic field
(Figs.~\ref{PICfig}F-G). The total pressure across the HP is
balanced. While the dominant pressure in the HS is from the plasma,
the dominant pressure in the LISM is magnetic. Thus, when reconnection
opens a path for HS plasma to escape into the LISM and mix with the
lower-pressure LISM plasma, there is nothing to balance the total
pressure and so the region compresses to increase the magnetic field
amplitude. This behavior is primarily seen at separatrix crossings
remote from where reconnection locally reduces the magnetic field
strength (the interval $\Delta R/d_i=38-50$ in Figs.~\ref{PICfig}F-G). In
the V1 data, the magnetic field strength is also observed to increase
where the local flux of HS plasma decreases
\citep{stone13a,krimigis13a,webber13a,burlaga13a}.

Thus, based on our simulations, we suggest that the V1 observations of
simultaneous drops (increases) in HS (LISM) particle fluxes occur at a
series of separatrix crossings outside of the HP that are associated
with a nested set of magnetic islands that form at the HP
(Fig.~\ref{cartoonfig}). At such crossings the magnetic field
direction does not change significantly, while, as seen in the
simulation data, particle fluxes can change sharply. Three active
reconnection sites at the HP, and associated separatrices with two
nested islands, are sufficient to explain the sequence of Voyager
events. On day 166 the spacecraft crossed a current layer, on day 190
the flux of HS electrons began dropping, on day 209 another current
layer was crossed and on days 210, 222 and 238 three successive drops
(increases) in the HS (LISM) particle fluxes occurred. The day 190
drop in the HS electrons suggests that after this time the magnetic
field was no longer laminar so that these electrons, with their small
Larmor radii and large velocities, could leak into the LISM.

Islands and x-lines flowing away from an active x-line (e.g., the
rightmost x-line in Fig.~\ref{cartoonfig}) correspond to reconnection
sites that developed earlier in time. The separatrix field lines
connect to x-lines, which act as bottlenecks to particle transport
across the HP. There are two reasons for this. First, at the x-line
the magnetic field turns into the $N$ direction since the $R$ and $T$
magnetic field components are zero. Thus, the x-line halts the
field-aligned streaming of particles across the HP. Second, to the
extent that $B_N$ is weak compared with $B_T$ energetic particles can
scatter near x-lines, which further limits transport across the HP. In
contrast, downstream of separatrices (to the left in
Fig.~\ref{cartoonfig}) particles can freely stream across the
HP. Thus, the day 210 drop in ACRs (rise in GCRs) occurred at the
separatrix corresponding to the left-most x-line of
Fig.~\ref{cartoonfig} which blocked the transport of ACRs (GCRs)
across the HP. The ACR (GCR) intensity rose (dropped) as the
spacecraft crossed field lines that formed an open corridor across the
HP to the left of the middle x-line in Fig.~\ref{cartoonfig}.  In this
region the flow of GCRs into the HS acts as a sink for the intensity
of these particles. The second ACR drop on day 222 and subsequent
recovery is similar. The final drop of the HS particle fluxes on day
238 occurred at the separatrix of the right-most x-line in
Fig.~\ref{cartoonfig}. LECP measurements of ACR anisotropies
\citep{krimigis13a} show that particles propagating parallel to the
magnetic field dropped more rapidly than those with perpendicular
pitch angles on day 238.  Such behavior is consistent with our
schematic -- in a weakly stochastic field parallel moving particles
will more quickly escape in the LISM.

An inconsistency between our PIC simulations and the observations
concerns the spatial region where sharp jumps in the ACRs and GCRs
take place. In the simulations the anti-correlated jumps occur on both
sides of the HP but if our cartoon is correct, they occur only on the
LISM side in the V1 data. The contradiction is possibly because the
simulations are in a 2-D system. A magnetic field in a real 3-D system
will likely be at least mildly stochastic so that wandering field
lines will smooth the variability of ACR and GCR intensities far from
the HP boundary. More challenging 3-D simulations will be required to
explore this issue.

A second issue is the angle $\delta$ of the magnetic field outside of
the HP, which for the present MHD simulation is
$14^{\circ}$. $B_{ISM}$ for this simulation was chosen to match
heliospheric asymmetries \citep{opher09a,opher09b}. Other MHD
simulations based on fitting the IBEX ribbon yield $\delta\sim
30^{\circ}$ \citep{pogorelov13a}. Such values are considerably smaller
than those implied by cartoons in recent publications
\citep{burlaga13a}. In any case our conclusion that significant
variations in the density of the ACRs and GCRs can occur in regions
with essentially no variation in the field line orientation is not
sensitive to the value of $\delta$ in the LISM. Any rotation in the
magnetic field outside of the HP propagates at the local Alfv\'en
speed, which is well below the velocities of the particles of
interest. Thus, outside of the HP separatrices should retain their
original LISM orientation in locations where there are significant
variations in particle intensity.

Combining a typical aspect ratio of reconnection-produced islands
(0.1), the typical time between dropouts (10 days), and the speed of
Voyager 1 with respect to the HS plasma ($\approx 20\text{ km/s}$),
yields the approximate size of the islands in Fig.~\ref{cartoonfig} as
1 AU. This roughly equals the sector spacing downstream of the TS and
is consistent with previous simulation-derived estimates of magnetic
islands in the HS \citep{schoeffler12b}.

If the schematic with nested magnetic islands (Fig.~\ref{cartoonfig})
is correct, the dropouts in the HS particle fluxes occurred on the
LISM side of the HP on field lines that had a LISM source. Thus,
according to this picture V1 crossed the HP on day 208 and has been
crossing LISM fields since that time.  Our results thus suggest that
$B_T$ in the LISM is negative (i.e., has a polarity of $270^{\circ}$).

\begin{acknowledgments}

The authors acknowledge the support of NSF grant AGS-1202330 to the
University of Maryland, and NSF grant ATM-0747654 and NASA grant
NNX07AH20G to Boston University. The PIC simulations were performed at
the National Energy Research Scientific Computing Center and the MHD
simulations at the Ames NASA Supercomputer Center. We acknowledge
fruitful discussions with L. F. Burlaga, R. B. Decker, M. E. Hill and
E. C. Stone on the Voyager observations. This research benefited
greatly from discussions held at the meetings of the Heliopause
International Team at the International Space Science Institute in
Bern, Switzerland.

\end{acknowledgments}