DISCOVERING THE PROGENITORS TO TYPICAL FIELD MILLISECOND PULSAR BINARIES By Samuel John Swihart A DISSERTATION Michigan State University in partial fulfillment of the requirements Submitted to for the degree of Astrophysics and Astronomy – Doctor of Philosophy 2020 ABSTRACT DISCOVERING THE PROGENITORS TO TYPICAL FIELD MILLISECOND PULSAR BINARIES By Samuel John Swihart Neutron stars in low-mass binaries can accrete matter and angular momentum from a non-degenerate companion star and be spun-up to rapid spin periods, making them detectable as millisecond radio pulsars (MSPs). Although the swollen stellar companion transfers mass onto the neutron star for hundreds of millions of years or more, most MSPs in the Galactic field have a compact white dwarf companion in a wide orbit, representing the end product of the spin-up process. Following the launch of the Fermi Space Telescope, it became clear MSPs were nearly ubiquitous emitters of high-energy -rays. Using multiwavelength follow-up observations of unidentified Fermi sources, the rarely observed binary MSPs with non-degenerate companions started being discovered in greater numbers. In this dissertation, I expand this population of rare MSPs with new discoveries and present a novel method for identifying and characterizing new Galactic compact neutron star binaries: cross-matching unassociated Fermi sources with short-period optical variables from high-cadence all-sky surveys. This technique has allowed for the discovery of a number of new systems with non-degenerate companions in orbits around previously unknown MSPs. Among these are a few systems with red giant companions in relatively wide orbits ( (cid:38) 1 d) that will naturally evolve into the MSP–white dwarfs that represent the bulk of known MSP binaries. The long orbital periods, giant companions, and inferred evolutionary tracks of these systems indicate they are entirely different from the close-orbit MSP binaries being found recently with Fermi, and suggests a new subclass of compact neutron star systems that are the progenitors of typical field MSP binaries. Copyright by SAMUEL JOHN SWIHART 2020 To my mom and dad – you made all this possible iv TABLE OF CONTENTS LIST OF TABLES . LIST OF FIGURES . CHAPTER 1 . . . . . . . . . . . . . . . . . . . . ix . . . . . . . . . INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Pulsars . 1.2 Pulsar Evolution . . 1.3 Millisecond Pulsars (MSPs) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Classes of Millisecond Pulsars 1 1 4 7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.4.1 The Canonical MSP Binaries . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.4.2 Spider Millisecond Pulsars . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.4.3 Transitional Millisecond Pulsars . . . . . . . . . . . . . . . . . . . . . . . 18 1.4.4 A New Spider Subclass: MSP binaries with Giant Companions . . . . . . . 21 . 22 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 1.5.1 Gamma-rays . 1.5.2 X-rays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 . 1.5.3 Optical Variability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 1.5.4 Optical Spectroscopy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 1.5.4.1 Radial Velocity Tracking . . . . . . . . . . . . . . . . . . . . . . 29 1.5.4.2 Determining the Mass Ratio . . . . . . . . . . . . . . . . . . . . 30 Spectroscopy if Accretion Disk is Present . . . . . . . . . . . . . 31 1.5.4.3 . . . . 33 1.5 Observations to Find and Characterize New Millisecond Pulsars . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 This Thesis . . . . . . . . . . . . . . . . . . . . CHAPTER 2 2.1 Introduction . 2.2 Observations . 2.4 Results . 2.3 Optical Spectropscopy . . . . . . . . . . . . . . . 2FGL J0846.0+2820: A NEW NEUTRON STAR BINARY WITH A GIANT SECONDARY AND VARIABLE -RAY EMISSION . . . . . . . . 37 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.2.1 Fermi-LAT Source . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.2.2 Optical and Near-IR Counterpart . . . . . . . . . . . . . . . . . . . . . . . 40 2.2.2.1 Catalina Sky Survey (CSS) . . . . . . . . . . . . . . . . . . . . . 40 PROMPT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.2.2.2 . 42 2.3.1 SOAR and MDM Spectroscopy . . . . . . . . . . . . . . . . . . . . . . . 42 2.3.2 Keck HIRES Spectroscopy . . . . . . . . . . . . . . . . . . . . . . . . . . 44 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 2.4.1 Keplerian Orbit Fitting and Mass Ratio . . . . . . . . . . . . . . . . . . . 46 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 2.4.2 . . . . . . . . . . . . . . . . . . . . . . . . . 53 2.4.3 Optical Light Curve Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 2.4.4 Component Masses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 2.4.5 Long-Term Optical Brightness . . . . . . . . . . . . . . . . . . . . . . . . 59 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fermi-LAT Analysis 2.4.3.1 Distance . . . . . . . . . v 2.5 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 . . . . . . . . . 3.2.2 Introduction . 3.3.1 VLA . 3.3.2 ATCA . 3.1 . 3.2 X-ray Observations . . . 3.2.1 Chandra ACIS . . . . . 3.4 Optical Observations CHAPTER 3 A MULTI-WAVELENGTH VIEW OF THE NEUTRON STAR BINARY 1FGL J1417.7–4402: A PROGENITOR TO CANONICAL MILLISEC- OND PULSARS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 3.2.1.1 X-ray Variability . . . . . . . . . . . . . . . . . . . . . . . . . . 70 Spectral Fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.2.1.2 Swift XRT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 3.3 Radio Continuum Observations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 3.4.1 SOAR Spectroscopy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 3.4.2 Optical/Near-IR Photometry . . . . . . . . . . . . . . . . . . . . . . . . . 82 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 3.6.1 The Spectra Themselves . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 3.6.2 Understanding the H . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 Stellar Chromosphere . . . . . . . . . . . . . . . . . . . . . . . 93 Stellar Wind . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 The Source of H . . . . . . . . . . . . . . . . . . . . . . . . . 95 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 3.6.2.1 A Disk . 3.6.2.2 3.6.2.3 3.6.2.4 3.7 Conclusions . 3.5 Light Curve Fitting . . . . 3.6 H Spectroscopy . 3.5.1 Distance . . . . . . . . . . . . . . . . . . . . . . CHAPTER 4 PSR J1306–40: AN X-RAY LUMINOUS REDBACK WITH AN EVOLVED . . . . . . 4.2.1 Introduction . 4.1 . 4.2 Optical Observations COMPANION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 . 102 SOAR Spectroscopy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 4.2.1.1 Determining the Mass Ratio . . . . . . . . . . . . . . . . . . . . 103 4.2.1.2 The Minimum Neutron Star Mass . . . . . . . . . . . . . . . . . 105 4.2.2 Optical Photometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 4.3 Light Curve Fitting . . . . . 4.4 H Spectroscopy . . 4.5 Discussion . . 4.3.1 Distance . . . . . . . . . . . . . . . CHAPTER 5 A NEW LIKELY REDBACK MILLISECOND PULSAR BINARY WITH 5.1 Introduction . 5.2 Observations . A MASSIVE NEUTRON STAR: 4FGL J2333.1–5527 . . . . . . . . . . . . 119 . 119 . 121 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi . . . 5.2.3.1 5.2.3.2 5.3 Results & Analysis 5.3.1.1 Minimum Neutron Star Mass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.1 The -ray Source & Optical Discovery . . . . . . . . . . . . . . . . . . . . 121 5.2.2 X-ray Observations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 5.2.3 Optical Observations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 SOAR Photometry . . . . . . . . . . . . . . . . . . . . . . . . . 124 SOAR Spectroscopy . . . . . . . . . . . . . . . . . . . . . . . . 125 . 126 5.3.1 Modeling the Radial Velocities . . . . . . . . . . . . . . . . . . . . . . . . 126 . . . . . . . . . . . . . . . . . . . 128 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 5.3.2 X-ray Spectrum . . 5.3.3 X-ray Light Curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 5.3.4 Optical Light Curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 5.3.4.1 Modeling the Light Curves . . . . . . . . . . . . . . . . . . . . . 134 5.3.4.2 Distance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 CHAPTER 6 CONCLUSION AND FUTURE WORK . . . . . . . . . . . . . . . . . . . 141 . 142 . 145 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1 Future Prospects . BIBLIOGRAPHY . 5.4 Conclusions . . . . . . . . . . . . . . . . . . . vii LIST OF TABLES Table 2.1. PROMPT Photometry of 2FGL J0846.0+2820 . . . . . . . . . . . . . . . . . . 41 Table 2.2. Barycentric Radial Velocities of J0846 . . . . . . . . . . . . . . . . . . . . . . . 43 Table 2.3. Best fit ELC parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 Table 3.1. Summary of Chandra Spectral fits for J1417 . . . . . . . . . . . . . . . . . . . . 74 Table 3.2. Swift XRT Observations of J1417 . . . . . . . . . . . . . . . . . . . . . . . . . 74 Table 3.3. VLA radio continuum data for J1417 . . . . . . . . . . . . . . . . . . . . . . . 79 Table 3.4. ATCA data for J1417 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 Table 3.5. SMARTS Photometry of J1417 . . . . . . . . . . . . . . . . . . . . . . . . . . 84 Table 3.6. Summary of ELC fits for J1417 . . . . . . . . . . . . . . . . . . . . . . . . . . 88 Table 4.1. Barycentric Radial Velocities of PSR J1306–40 . . . . . . . . . . . . . . . . . . 104 Table 4.2. SMARTS Photometry of PSR J1306–40 . . . . . . . . . . . . . . . . . . . . . . 107 Table 4.3. Summary of ELC fits for PSR J1306–40 . . . . . . . . . . . . . . . . . . . . . . 109 Table 5.1. SOAR Velocities of J2333 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 viii LIST OF FIGURES Figure 1.1: Toy model of a typical radio pulsar Figure 1.2: The  − (cid:164) diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 6 8 Figure 1.3: Schematic diagram demonstrating the concept of a Roche lobe . . . . . . . . . Figure 1.4: The key phases of the MSP recycling process . . . . . . . . . . . . . . . . . . . 10 Figure 1.5: Orbital period versus companion mass for recycled field MSPs with known companion star types . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 Figure 1.6: Artist rendition of an eclipsing spider MSP . . . . . . . . . . . . . . . . . . . . 15 Figure 1.7: Artist’s renditions of the two main states of a transitional MSP . . . . . . . . . 19 Figure 1.8: Diagram of an intrabinary shock . . . . . . . . . . . . . . . . . . . . . . . . . 25 Figure 1.9: Example X-ray light curves due to strong intrabinary shocks . . . . . . . . . . . 26 Figure 1.10: Diagram of ellipsoidal variations . . . . . . . . . . . . . . . . . . . . . . . . . 28 Figure 1.11: Example double-peaked H emission line due to an accretion disk . . . . . . . 32 Figure 2.1: Optical spectra of 2FGL J0846.0+2820 . . . . . . . . . . . . . . . . . . . . . . 44 Figure 2.2: Radial velocity curve for the optical counterpart to 2FGL J0846.0+2820 . . . . 47 Figure 2.3: -ray light curve of 2FGL J0846.0+2820 . . . . . . . . . . . . . . . . . . . . . 50 Figure 2.4: Map of -ray emission near the location of 2FGL J0846.0+2820 . . . . . . . . 51 Figure 2.5: -ray light curve of 2FGL J0846.0+2820 near the time the flux dropped . . . . . 53 Figure 2.6: PROMPT optical light curve of 2FGL J0846.0+2820 . . . . . . . . . . . . . . . 54 Figure 2.7: Long-term optical light curve of 2FGL J0846.0+2820 . . . . . . . . . . . . . . 60 Figure 2.8: Long-term optical light curve of 2FGL J0846.0+2820 split into separate epochs 61 Figure 3.1: Chandra X-ray light curve of 1FGL J1417.7–4407 . . . . . . . . . . . . . . . . 71 ix Figure 3.2: Chandra X-ray spectra of 1FGL J1417.7–4407 . . . . . . . . . . . . . . . . . . 73 Figure 3.3: Long-term Swift X-ray light curve of 1FGL J1417.7–4407 . . . . . . . . . . . . 76 Figure 3.4: Swift X-ray light curve of 1FGL J1417.7–4407 folded on the binary period . . . 77 Figure 3.5: Optical/near-IR light curves of 1FGL J1417.7–4407 . . . . . . . . . . . . . . . 85 Figure 3.6: Phase-resolved optical spectra of 1FGL J1417.7–4407 near the H region . . . 90 Figure 3.7: Orbital position of the secondary in 1FGL J1417.7–4407 corresponding to the varying H phenomenology . . . . . . . . . . . . . . . . . . . . . . . . . . 91 Figure 4.1: Radial velocity curve of PSR J1306–40 . . . . . . . . . . . . . . . . . . . . . . 103 Figure 4.2: Optical light curve of PSR J1306–40 . . . . . . . . . . . . . . . . . . . . . . . 108 Figure 4.3: Phase-resolved optical spectra of PSR J1306–40 around the H region . . . . . 114 Figure 5.1: X-ray and optical finder images for 4FGL J2333.1–5527 . . . . . . . . . . . . . 122 Figure 5.2: Optical spectra of 4FGL J2333.1–5527 . . . . . . . . . . . . . . . . . . . . . . 126 Figure 5.3: Radial velocity curve of 4FGL J2333.1–5527 . . . . . . . . . . . . . . . . . . . 128 Figure 5.4: X-ray spectrum of 4FGL J2333.1–5527 . . . . . . . . . . . . . . . . . . . . . . 131 Figure 5.5: X-ray light curve of 4FGL J2333.1–5527 . . . . . . . . . . . . . . . . . . . . . 132 Figure 5.6: Optical light curve of 4FGL J2333.1–5527 . . . . . . . . . . . . . . . . . . . . 135 Figure 6.1: Component masses for the redback and huntsmen systems that have well- constrained mass estimates for the neutron star . . . . . . . . . . . . . . . . . . 143 x CHAPTER 1 INTRODUCTION 1.1 Pulsars Neutron stars are the most extreme stars in the Universe. These extraordinary compact objects represent one of the possible end stages of massive stellar evolution and can reach densities that exceed those of atomic nuclei. Many observed neutron stars are highly magnetic, endowed with magnetic fields billions of times stronger than Earth’s, and spin at high speeds, similar to that of a household blender. Stars release energy by fusing their hydrogen and helium into increasingly heavier isotopes. For stars with a birth mass (cid:38) 8–25 (cid:12), their central densities and temperatures eventually become sufficiently high that matter in the center is converted to iron. Since the binding energy per nu- cleon is greatest for isotopes about iron, once the star reaches this point no more nuclear energy is available. The core grows in mass until it reaches about 1.4 (cid:12), at which point it becomes unstable and implodes. As the core collapses, the extreme densities induce electrons and protons to combine, form neutrons, and release neutrinos. The neutrinos heat the envelope of the star and blow it away, which results in a violent event know as a core-collapse supernovae. The resulting object, composed now of mostly neutrons, is known as a neutron star1, and is typically on the order of the Chandrasekhar mass (∼ 1.4 (cid:12)) with a radius of ∼10–15 km (e.g., Steiner et al., 2013; Özel & Freire, 2016). Prior to the collapse, the progenitor star had some rotational energy. Since angular momentum is (mostly) conserved after the supernova explosion and removal of the outer envelope, the leftover 1The above mass range should be taken as a very rough approximation, since the precise range of initial masses that determine whether a massive star will become a neutron star or a black hole is still an active field of research, and depends sensitively on the stellar mass loss rate, which can be highly unpredictable, especially if the star is in a binary system. 1 neutron star core will attain a rapid rotation rate, on the order of tens to hundreds of times per second at the time of its birth (e.g., Migliazzo et al., 2002; Kramer et al., 2003; Lorimer & Kramer, 2004). The precise mass and radii of these rapidly rotating objects is difficult to estimate, and relies heavily on our theoretical understanding of the equation of state of matter2 at these extreme densities. Although a description of the different equations of state for neutron stars is beyond the scope of this thesis, an upper limit of ∼2.5 (cid:12) (known as the “"Tolman-Oppenheimer-Volkov limit”; Lattimer & Prakash, 2001) has been proposed for the theoretical maximum mass for a neutron star. The current observationally confirmed record holders are MSP J0740+6620 with a mass of 2.14+0.10−0.09 (cid:12) (Cromartie et al., 2020), as well as an indirect constraint on the mass of the final merger remnant from the binary neutron star merger GW170817 (NS (cid:46) 2.17 (cid:12), Margalit & Metzger, 2017). In addition to rapid rotation, many neutron stars also have very strong magnetic fields (on the order of 107 − 1015 G), making them ideal laboratories for accelerating charged particles. As the magnetized neutron star rotates, a strong electric field is generated that accelerates charged particles (electrons and positrons) along the magnetic field lines. This dense plasma of charged particles corotates with the neutron star and its surrounding magnetic/electric fields and is known as the “pulsar magnetosphere” (Goldreich & Julian, 1969). The outer limit at which the plasma can still corotate with the neutron star and not exceed the speed of light sets the boundary of the “light cylinder radius.” For a recent review on the structure and properties of the pulsar magnetosphere, see Spitkovsky (2011), and references therein. A simple, toy model of a typical pulsar is shown in Figure 1.1. The light cylinder divides the pulsar magnetosphere into two main regions, one where the magnetic field lines can fully close within the light cylinder, and another where the field lines remain open. As charged particles are accelerated along the open field lines, a cascade of low-energy photons are emitted in a direction tangential to the field lines. As the neutron star rotates, this results in conical beam shaped emission 2The equation of state provides a relation between the density and pressure of the star, and so is tightly linked to the macroscopic properties of the neutron star, such as its mass and radius 2 Figure 1.1: Toy model of a typical radio pulsar, taken from Lorimer & Kramer (2012). In this simple model, the neutron star is approximated as a rotating magnetic dipole along with its surrounding magnetosphere. Charged particles are extracted from near the neutron star surface and are accelerated along the open field lines, resulting conical beam-shaped radio emission. The misalignment of the rotational axis of the neutron star with respect to the magnetic axis results in the “lighthouse effect,” where an observer at Earth observes pulses of radio emission each time the emission cone crosses their line of sight. emanating from a region centered on the magnetic poles. As the emission cone rotates into our line of sight, we see a brief “pulse” of radio emission, analogous to a rotating lighthouse beam. When 3 pulsations are observed, these rotating neutron stars are known as “pulsars,” where the radio pulses we detect will be periodic, modulated on the rotational period of the neutron star. Pulsars are typically observed with sensitive radio telescopes capable of detecting very faint signals. However, typically the individual pulses are so weak that they are swamped by the back- ground noise, and are not easily detectable. Therefore, in order to boost the signal, hundreds or thousands of individual pulses are “folded” or “stacked” using advanced Fourier techniques in order to create an “integrated pulse profile.” The shape and strength of this profile will be different for different pulsars, and depends mainly on the orientation of the pulsar’s magnetic field with respect to our line of sight and on the geometry and size of the radio emission beam. The observed flux density from pulsars are usually described by a basic power-law with a steep spectral index 0 <  < −4, for  ∝  (Maron et al., 2000), and so observations are usually performed at some specific observing frequency that will optimize the detection. This simple observing scenario can become much more complicated when intervening material from a binary companion causes frequency-dependent eclipses of the radio emission (see Section 1.4.2). 1.2 Pulsar Evolution As a pulsar rotates and emits radiation as described above, its rotation rate will gradually slow down as it slowly loses rotational kinetic energy over time. This energy loss is remarkably stable, with some pulsars showing timing stability that rivals those of the best terrestrial atomic clocks. The normal evolution of a pulsar can be well described through a diagram showing the spin period  versus its spin period derivative (cid:164) = / (the rate at which the pulsar rotation is slowing down; Figure 1.2). This diagram also provides a convenient way to classify pulsars based on their intrinsic properties. Assuming the neutron star radiates as a pure magnetic dipole and the initial spin period at birth was much faster than the current one, and by making some basic assumptions about the relative alignment of the magnetic field and rotational axis of the pulsar, the approximate surface 4 magnetic field strength () and the characteristic age of the pulsar () can be approximated by: (cid:19)−1 (cid:19)(cid:18) (cid:18)  (cid:19)1/2 (cid:18)  (cid:164)  (cid:164) 10−15  =  2 (cid:164) ≈ 15.8 Myr  = 3.2 × 1019 G.  (1.1) (1.2) Both of these equations rely on assumptions that do not always hold true for each individual pulsar (for instance, the presice value for  is sensitive to the assumed mass and radius of the neutron star), and so these only provide rough estimates for the true age and surface magnetic field strength of a pulsar. Nevertheless, these approximate values can be used to roughly differentiate between pulsar sub-populations and their typical evolutionary tracks. In general, pulsars are born with spin periods on the order of  ∼ 10−1 − 10−2.5 seconds and with spin-down rates (cid:164) (cid:46) 10−15 (e.g., Migliazzo et al., 2002; Lorimer & Kramer, 2004, 2012), placing them in the upper left portion of Figure 1.2. It is clear from this figure that, ordinarily, older pulsars will have slower rotational periods than their younger, faster-spinning progenitors. As young pulsars lose rotational energy, their spin period increases while the spin-down rate decreases, perhaps due to a slowly weakening magnetic field (e.g., Aguilera et al., 2008; Viganò et al., 2013), although the timescale for this magnetic field decay is still an open question (e.g., Pons et al., 2007; Viganò et al., 2012). These effects cause a typical pulsar to quickly move down and right on the  − (cid:164) diagram over time, towards the bulk of the pulsar population, which have spin periods of ∼0.5 s and spin-down rates of about 10−15. Eventually, the rotational energy of the neutron star will be so low that the pulsar emission mechanism is not powerful enough to produce observable radio pulsations, and the pulsar crosses the so-called “pulsar death line” into the “pulsar graveyard.” Most young pulsars emit energy very efficiently and so they rapidly reach this region of the dia- gram ∼105−7 years after their birth. At this point the evolution slows down substantially due to the weakening magnetic field and slower spin-down rate and so the bulk of the known pulsars live near this “death line.” The vast majority of these neutron stars are in isolated systems. However, there 3https://www.atnf.csiro.au/research/pulsar/psrcat/ 5 Figure 1.2: The  − (cid:164) diagram, showing the measured spin period and period derivative for all 2811 known pulsars from the ATNF pulsar catalogue3. Lines of constant surface magnetic field strength and characteristic pulsar age (see text for definitions) are shown. The “pulsar graveyard,” (see text) is represented by the shaded region in the lower right. This diagram is generally used as a convenient way to classify the different sub-populations of pulsars and track their basic evolution. The dark grey circles highlight the bulk of the pulsar population, with spin periods of ∼0.5 s, ages from ∼ 105 − 108 years, and moderate (1011 − 1013 G) magnetic fields. The top-right of the diagram hosts a population of young neutron stars with extremely strong magnetic fields, known as “magnetars,” or sometimes referred to as soft gamma-ray repeaters (SGRs) or anomalous X-ray pulsars (AXPs). The lower-left part of the diagram shows the objects that are the subject of this thesis. These are the rapidly-rotating millisecond pulsars, which attained their fast spin periods by being spun-up via accretion from a binary companion. Noteworthy features (e.g., presence of a binary companion) and associations (e.g., located in a supernova remnant (SNR)) for some of the pulsars are noted by the legend. Figure produced using the Python package psrqpy (Pitkin, 2018). 6 105yr106yr107yr108yr109yr1010G1011G1012G1013G1014G10−310−210−1100101Period(s)10−2210−2110−2010−1910−1810−1710−1610−1510−1410−1310−1210−1110−1010−9PeriodDerivativeSGR/AXPBinaryRadio-IREmission”Radio-Quiet”RRATPulsedThermalX-raySNR exists a subclass of pulsars on the lower-left part of Figure 1.2 that rotate very rapidly (>30 times per second) but also appear very old. These systems, which are often found in binaries, are known as millisecond pulsars and are described in detail in the following sections. A final subpopulation is visible in the upper-right portion of Figure 1.2, consisting of pulsars with long spin periods but very strong magnetic fields ( (cid:38) 1013.5 G), resulting in extremely high spin-down rates ((cid:38) 10−12). These objects, which are beyond the scope of this thesis, are young pulsars commonly referred to as “magnetars” (Kaspi & Beloborodov, 2017) and experience non-standard evolution, likely due to the presence of a high-mass companion. 1.3 Millisecond Pulsars (MSPs) Normal, isolated pulsars quickly evolve to have relatively long spin periods (∼0.5 sec) and eventually cross the pulsar death line on the order of ∼107−9 years after their birth, and become undetectable. Therefore, there must be some external source to facilitate the “rejuvenation” of a pulsar to rapid spin periods in order to create the population of very old, rapidly spinning pulsars seen in the bottom left of Figure 1.2. This can only be done via accretion from a binary companion. Most stars in the Universe form in binaries. If one star is much more massive than the other, it will evolve faster and reach its evolutionary endpoint sooner than its lower mass companion. If the more massive star is (cid:38) 8 (cid:12), this endpoint will typically manifest in a core-collapse supernova, leaving behind a rotating pulsar (see Section 1.1). Assuming the less massive star is not destroyed by the supernova explosion, and the binary is not disrupted, the secondary companion will evolve normally as a low-mass star orbiting the pulsar primary. This companion star will undergo normal stellar evolution; as its envelope gradually expands while on the main sequence, and especially once it evolves off the main-sequence into the red-giant phase, the star will begin to fill its Roche lobe (e.g., Bhattacharya & van den Heuvel, 1991). At this point the gravitational pull of the neutron star will dominate and pull material down onto its surface via the stellar wind or through the inner Lagrange point (Figure 1.3). This material (which is mostly Hydrogen-rich gas) comes in not 7 directly, but at an angle, imparting angular momentum (in addition to accreting mass) onto the neutron star. Figure 1.3: Schematic diagram demonstrating the concept of a Roche lobe in a binary system. Top: A neutron star (dark grey circle, size not to scale) and a companion star (red circle) are in a binary system. The light grey tear-drop-shaped regions denote the Roche lobes of the two stars, within which material is gravitational bound to each respective object. If the envelope of the companion expands to fill its Roche lobe, material can stream towards the neutron star. Bottom: Two stars orbit around a common center of mass (CM), which lies closer to the more massive object (left black dot). Thin grey lines denote surfaces of equal gravitational potential. Thick black lines are the Roche lobes of the two respective stars. If the companion (right black dot) fills its Roche lobe, material streams towards the massive object through inner Lagrange point 1. Figure credit: Ed Brown; License: Creative Commons Attribution–NonCommercial–ShareAlike 4.0 International (CC BY-NC-SA 4.0) This process (known as “spin-up”) rapidly speeds up the rotation rate of the neutron star (on a timescale of at least ∼ 106−7 years, though it may be much longer depending of the mass and 8 CML4L5L2L3L1 orbital period of the companion (see Section 1.4); Urpin & Konenkov, 1997), analogous to hitting the side of a basketball to make it spin on your finger, resulting in a neutron star that is now spinning with periods of a few milliseconds (i.e., a millisecond pulsar; Tauris & van den Heuvel, 2006). The limiting spin period that a spun-up neutron star could theoretically reach is still an open question4, although it mainly depends on the balance between the magnetic pressure of the accreting neutron star and the inward pressure of the accreting material (e.g., Bhattacharya & van den Heuvel, 1991). This theoretical value can be complicated by gravitational wave emission losses or a complex and changing magnetosphere in the close binary environment. This accretion spin-up process is known as pulsar “recycling” since the old, once slowly rotating, neutron star is spun-up (“rejuvenated”) to rapid spin periods (Alpar et al., 1982a; Radhakrishnan & Srinivasan, 1982). A schematic diagram showing some of the key phases of the MSP recycling process is displayed in Figure 1.4. Millisecond pulsars are extremely stable rotators; some have spin periods that are constant to ∼1 part in 1014 (Ransom, 2013), making them the most precise and stable clocks in the Universe. This is because the same accretion processes that spun-up the pulsar to millisecond spin periods are also likely responsible for decreasing the strength of the dipole magnetic field of the neutron star (from ∼1012 Gauss down to ∼108 Gauss; e.g., Bisnovatyi-Kogan & Komberg, 1974; Bhattacharya, 1995; Cumming, 2003). This relatively weak magnetic field results in extremely low spin-down rates, and hence extraordinarily stable rotation periods over very long timescales, placing them in the bottom-left portion of Figure 1.2. Millisecond pulsars are also ideal laboratories for stringent tests of fundamental physics, in- cluding tests of general relativity, and observational evidence of gravitational waves. In addition, many MSPs are extremely heavy, pushing the theoretical boundaries of neutron star masses, and setting important constraints on the equation of state for ultra dense matter. Mass transfer from a companion must occur in order to complete the spin-up, however, many binary evolution models 4The fastest spinning pulsar known to date is PSR J1748–2446ad, with a spin period of ∼1.39 ms (Hessels et al., 2006). 9 Figure 1.4: The key phases of the MSP recycling process. In the first panel, a massive star is being orbited by a Sun-like companion. Eventually the massive star will explode in a supernova, leaving behind a rotating neutron star visible as a radio pulsar (panel 2). The spin frequency of the pulsar will gradually slow down over time as it loses rotational kinetic energy. In the third panel, the lower-mass star has undergone normal stellar evolution and expands into its red giant phase. At this point material will fall down towards the pulsar forming an accretion disk, making the system visible as a low-mass X-ray binary, and transferring matter and angular momentum to the neutron star, spinning it up to rapid spin periods and likely increasing its mass. Finally, accretion eventually ends and the spun-up neutron star will become visible again, now as a millisecond radio pulsar (panel 4). The strong, relativistic pulsar wind impinges on the companion, ablating material from its surface, causing radio eclipses for substantial fractions of the orbit (see Section 1.4.2). Recent observations have shown that some systems can switch back and forth between the phases shown in panels 3 and 4 (see Section 1.4.3). Figure credit: Bill Saxton, NRAO predict that only (cid:46)0.1–0.2 (cid:12) of material needs to be accreted onto a neutron star in order to spin it up to millisecond periods (e.g., Urpin & Konenkov, 1997; Tauris, 2012; Chen et al., 2013). This implies many of the neutron stars that became millisecond pulsars were actually born more massive than the bulk of the neutron star distribution, which suggests that finding new, massive MSPs can be important for our understanding of fundamental physics, as well as accretion onto compact objects. 1.4 Classes of Millisecond Pulsars 1.4.1 The Canonical MSP Binaries Although accretion from a binary companion is necessary for MSPs to attain their rapid rotation, the vast majority of MSPs are not actively accreting, but have instead already reached the end stage of the canonical recycling process. In the canonical scenario, the secondary star is of relatively 10 low-mass (∼ 1 (cid:12)). After the neutron star has been fully spun-up to millisecond spin periods, the envelope of the companion star is completely stripped (due to a combination of losing material through accretion and via the relativistic pulsar wind impinging on its surface (see Section 1.4.2)), leaving behind a low-mass (∼ 0.1 − 0.3 (cid:12)) Helium white dwarf companion in a wide ((cid:38)3 days), circular orbit. There is strong theoretical support for this canonical scenario, for which a low-mass companion starts mass transfer once it evolves off the main sequence, into a (sub)giant phase (see Tauris & Savonije, 1999, and references therein). On the red giant branch, the luminosity of the star comes almost exclusively from hydrogen shell-burning, and this burning rate depends only on the mass of the degenerate helium core. The central core mass grows as shell-burning deposits “ashes” onto the core. In this phase the star moves up the red giant branch, where luminosity increases with only a small change in effective temperature at the stellar surface. Therefore, the radius of the giant star (star) is tightly correlated with the mass of the degenerate helium core. During mass transfer, the companion fills its Roche lobe (RL), such that star ≈ RL. The Roche lobe radius only depends on the relative masses of the neutron star and its companion, and their or- bital separation (and hence, their orbital period). Thus, from the time the star begins transferring mass through Roche lobe overflow, the mass of the degenerate helium core is directly associ- ated with the orbital period of the system, suggesting there must be a tight correlation between the mass of the resulting white dwarf and the final orbital period of the binary (e.g., Savonije, 1987). Not surprisingly, there is also strong observational evidence for this scenario. In the past ∼20 years, numerous MSP binaries have been discovered, and the vast majority of them are in end-stage systems like this, with white dwarf companions in wide-orbits. The masses and orbital periods of these degenerate companions have been shown to agree very well with the theoretical expectations, and so these systems are often referred to as the canonical MSP binaries (see Figure 1.5, blue circles). The progenitors to these canonical systems are likely MSPs with red giant companions in long orbits ( (cid:38) 1 d, see Section 1.4.4, Chapter 2, Chapter 3). 11 It should be noted that there are other MSP formation scenarios, and they depend primarily on the mass of the companion. If the companion has an intermediate mass (in the approximate range 2–6 (cid:12)), a similar MSP binary as described above will be produced after the recycling process is complete, except with a slightly more massive C-O white dwarf as opposed to a He white dwarf. In the final scenario, the secondary is a high-mass star ((cid:38) 8 (cid:12)). After this companion spins up the neutron star through accretion, it may also undergo a supernova explosion, forming a second neutron star. A few of these double neutron star systems are currently known, and measuring the pulsar motion and evolution in highly relativistic systems like these provide stringent tests of general relativity (e.g., Hulse & Taylor, 1975; Stairs, 2003; Lorimer & Kramer, 2004; Damour, 2009). A further description and analysis of these systems with intermediate and high-mass companions is outside of the goals of this thesis, but some good references are Pfahl et al. (2002) and van den Heuvel (2019), and references therein. MSPs are also present in large numbers in globular clusters. While MSPs in the Galactic field are likely the results of natural evolution in binary systems, MSPs in clusters are more likely formed as a result of a dynamical encounter in the more dense, kinematically hotter environment. Although these systems are important for a number of science cases related to MSPs, they were likely not formed in the canonical evolution scenario, and so what they can tell us about the MSP recycling process is limited, and are not the focus of this dissertation. 1.4.2 Spider Millisecond Pulsars As described above, in the canonical evolution scenario for MSPs, the neutron star accretes matter and angular momentum from a low-mass, main-sequence or giant companion star, and is spun up to millisecond spin periods. Although this mass transfer is expected to be lengthy (e.g., lasting ∼1 Gyr for a 1 (cid:12) star with an initial period of 3 days, and ∼150 Myr for an initial period of 6 days (Tauris & Savonije, 1999)), the vast majority of MSPs in the field of the Galaxy are fully recycled, such that accretion has permanently ended, and have degenerate, white dwarf companions in relatively wide orbits (Tauris & van den Heuvel, 2006). Therefore, most MSP binaries offer only 12 Figure 1.5: Orbital period versus companion mass for recycled field MSPs with known com- panion star types. Most MSPs have He white dwarfs (open blue circles) from the canonical binary evolution scenario, well-represented by models from Tauris & Savonije (1999). These models (black line) assume solar metallicity and an initial secondary mass of 1.0 (cid:12), and denote the endpoints of an ensemble of systems with varying initial period, not the evolution of a single binary. The few CO white dwarfs (filled cyan circles) likely had more massive secondaries and underwent close common-envelope evolution. The field redbacks (red squares) and black widows (open black squares) are visible at short orbital period, and it is clear that these systems will not have normal field MSPs as their progeny. Instead the likely progenitors to the canonical field MSPs are systems with red giant companions in relatively wide orbits, the aptly-named “huntsmen” subclass (Section 1.4.4). The systems that are the subjects of this thesis (Chapters 2–5) are labeled with stars. Open red squares and stars refer to candidate systems for which radio pulsations have not yet been detected. Figure adapted from Strader et al. (2019). 13 this thesis:J1417J0846J1306J2333 indirect constraints on the spin-up process. The observation that non-degenerate companions are rare among MSPs suggests that the pulsed radio emission is quenched during even low-level accre- tion onto the neutron star (e.g., Illarionov & Sunyaev, 1975), such that the systems are undetectable as radio pulsars but visible as low-mass X-ray binaries (LMXBs). Fortunately, in the past twelve years, a breakthrough has been enabled with follow-up obser- vations of -ray sources (0.1–100 GeV energies) detected with the Fermi Gamma-ray Space Tele- scope5, revealing pulsars difficult to detect with radio observations, and leading to an enormous increase in the number of known MSP binaries with non-degenerate, hydrogen-rich companions. These systems show substantial radio eclipses for large fractions of the binary orbit. These eclipses are likely due to the radio emission being absorbed or scattered by intervening ionized material associated with the companion, but this material apparently does not affect the detectability of the GeV emission. This excess material is produced by highly relativistic particles (the “pulsar wind”) impinging on the non-degenerate companion’s surface, in a process known as ablation (similar to evaporation). The relativistic pulsar wind works to blow material off the companion star, which may stay bound to the system, and affects the detectablitiy of the MSP, making it difficult (and in some cases impossible) to detect radio pulsations. A cartoon of what an eclipsing system like this might look like is shown in Figure 1.6. These eclipsing MSP sources are commonly referred to as “spiders” because of the way the neutron star partially (and sometimes completely) “cannibalizes” the secondary. This is also one of the leading theories behind the existence of isolated MSPs, where the pulsar wind completely ablates its companion until there is nothing left (e.g., Fruchter et al., 1988; Kluzniak et al., 1988). These spider systems are typically classified as “black widows” or “redbacks” (after the cannibal- istic spiders), based on the mass and stellar type of the companion. Black widows are systems that have semi-degenerate, very low-mass companions ( (cid:46) 0.08 (cid:12)) that are slowly being ablated by the pulsar wind. Fruchter et al. (1988) discovered the 5https://www.nasa.gov/content/fermi-gamma-ray-space-telescope 14 Figure 1.6: Artist rendition of an eclipsing spider MSP. The pulsar wind blows material off the non-degenerate companion star and prevents material from reaching the neutron star surface. This material may stay bound to the system, making it difficult to detect radio pulsations; Image credit: European Space Agency & Francesco Ferraro (Bologna Astronomical Observatory) first black widow system, PSR B1957+20. This MSP had a binary companion in a 9.2 hour orbit, and the radio pulsar was eclipsed for ∼50 minutes during each orbit. These eclipses were due to the ∼0.02(cid:12) companion being highly irradiated on its tidally-locked inner face by the pulsar wind, which evaporates material from the companion surface and blocks the pulsed radio emission at regular intervals. This was the first system discovered that suggested that pulsar recycling may result in isolated MSPs; however it is unlikely that this particular MSP would fully ablate its com- panion within the age of the Universe. Additional systems have since been found where the level of irradiation (i.e., heating) is stronger than in PSR B1957+20 and the companion has even lower mass, and so this process may indeed be responsible for forming isolated MSPs (e.g., Ginzburg & Quataert, 2020). At the time of this writing, there are currently ∼40 confirmed black widow MSPs, with a handful of candidates that are awaiting confirmation. Redbacks, on the other hand, have significantly more massive ( (cid:38) 0.1 (cid:12)), nearly Roche 15 lobe filling, non-degenerate secondaries (Roberts, 2011), and are usually in slightly wider orbits ( (cid:38) 0.25d). Contrary to the much lighter, semi-degenerate black widows, these companions are much more akin to normal stars like the Sun, and the more massive secondaries and moderately wider orbits among redbacks ensures the companions will survive the wind erosion from the pulsar. The first of these redback systems was discovered by D’Amico et al. (2001a) as the result of tim- ing observations of the MSP PSR J1740–5340 in the globular cluster NGC 6397. The radio pulsar in this binary is eclipsed for ∼40% of its 32 hour orbit, and subsequent observations revealed the companion is consistent with a main-sequence star (or perhaps slightly more evolved) with a mass of ∼0.3 (cid:12) (Orosz & van Kerkwijk, 2003), and complex H absorption features, likely associated with a wind from the companion (Sabbi et al., 2003). These observations are broadly consistent with the observed properties of numerous other redback-like MSPs found in recent years, including 2FGL J0846.0+2820, 1FGL J1417.7–4402, PSR J1306–40, and 4FGL J2333.1–5527 described in detail in the following chapters. The redback population has seen a large increase in recent years as a result of follow-up obser- vations of unidentified Fermi -ray sources (see Section 1.5), with over 10 discoveries published since 2016 (Strader et al., 2019). Currently there are 15 confirmed redbacks outside of globular clusters, with approximately 12 additional candidates where the multiwavelength data suggests a redback classification but for which radio pulsations have not yet been found. The companion masses and orbital periods for all the field MSPs for which the companion star type is known is displayed in Figure 1.5. The companions in both black widows and redbacks are tidally-locked (similar to the Earth- moon system) and are usually in very circular orbits as a result of strong tidal interactions during the mass transfer process. While the companion is being irradiated and stellar material is being evaporated by the pulsar, the size of the orbit changes due to the competition between gravitational wave emission and magnetic braking working to shrink the orbit, and evaporated companion mate- rial being expelled from the system, which effectively adds angular momentum and widens the orbit. 16 It is still an open question whether all neutron star LMXBs with initially short orbital periods will become redbacks and black widows. In the binary models of Podsiadlowski et al. (2002), they found that if the orbital period at the onset of mass transfer is (cid:46)18 hr, then the orbit will shrink due to the dominant effects of gravitational wave emission and magnetic braking extracting angular momentum from the system, resulting in systems with orbital periods similar to the redbacks and black widows. If mass transfer begins when the orbital period is above this bifurcation period, then the orbit will widen over time, leading to a canonical MSP system with a white dwarf companion. Although these models are fairly robust, the evolutionary fate of an individual system depends sensitively on the efficiency of magnetic braking and on the mass of the companion at the start of mass transfer. For nearly all redbacks, the systems are below this bifurcation period such that the orbit is slowly shrinking, and so it is still an open question whether all redbacks will transition to be black widows, although the current evidence suggests that the two subclasses evolve independently (Chen et al., 2013). The level of efficiency with which the companion is being evaporated will have an important effect on this final fate. The pulsar wind energy flux is not isotropic, but will instead follow the complex geometry of the pulsar magnetosphere. Therefore, in determining the evaporation efficiency of the companion star, there is likely a beaming effect that depends on the orientation of the pulsar magnetic field axis with respect to the direction of the companion star. Chen et al. (2013) demonstrated that in general, the companions in redbacks are evaporated more efficiently than black widows, causing the system to shrink less quickly. The correct combination of evaporation and magnetic braking efficiency has been shown to reproduce the strongly bimodal distribution of the spider MSPs fairly well; however, these quantities are still not well known, are sensitive to the underlying assumptions about the system parameters, and likely vary from system to system, so establishing the relationships between all normal neutron star LMXBs, canonical MSPs, and the redbacks and black widows is still an open question. The constant irradiation of the tidally-locked inner face of these spider companions by the pulsar wind creates extreme temperature asymmetries on the stellar surface. These large differences in temperature on different sides of the star can lead to large optical flux variations depending on 17 the orbital phase when the system is observed. The strong temperature gradient also drives rapid circulation in the photosphere, which can create a strong dynamo that may boost the companion magnetic fields. These enhanced magnetic fields can act as ducts, efficiently channeling relativistic pulsar wind particles to the companion surface, possibly leading to starspots (see Chapter 4). 1.4.3 Transitional Millisecond Pulsars When millisecond pulsars are still spinning up, they are typically visible as low-mass X-ray binaries (LMXBs), where an accretion disk forms around the neutron star, and emit high-energy thermal radiation visible in X-rays due to their high temperature. If the material from the disk is funneled to the neutron star surface by strong magnetic fields, hot-spots can form at the magnetic poles that are bright in X-rays. As the neutron star rotates, these spots can sweep by an observer’s line-of-sight, resulting in millisecond X-ray pulsations, similar to what is commonly seen in the radio when no disk is present. These accreting MSPs and systems that share some of their multiwavelength characteristics are vital for understanding the connections between neutron star LMXBs and normal radio MSPs. Redbacks in particular have drawn extensive interest in recent years, as at least three of these sys- tems have been found to switch between this accretion-powered X-ray binary state and a rotationally- powered radio pulsar state, on timescales of days to months (Archibald et al., 2009; Papitto et al., 2013; Bassa et al., 2014; Roy et al., 2015; Bogdanov et al., 2015; Johnson et al., 2015). Since accretion is needed in order to spin-up a pulsar to millisecond rotational periods, it is reasonable to predict that there must be an evolutionary stage where the accretion disk becomes devoid of material and disappears, leaving a system which looks like a redback-type MSP. In fact, recent observations have shown that this process can actually proceed in both directions. In these systems, known as “transitional millisecond pulsars” (tMSPs), the binary system actually switches back-and-forth between two states: 1) a LMXB state, where an accretion disk is formed and material is actively accreting onto the neutron star from the non-degenerate companion, and 2) a redback-like radio MSP state, where the disk disappears and the radio pulsar re-emerges, being once again powered 18 by the spin-down energy of the neutron star, while being orbited by a (partially) Roche lobe filling companion. These two states are characterized by distinct multiwavelength phenomenology as described in detail in Figure 1.7. Figure 1.7: Artist’s renditions of the two main states of a transitional MSP, showing a model tMSP in the radio pulsar state (top) versus the accretion-powered disk state (bottom). In the pulsar state, the radio MSP (green) is observed, although it is typically eclipsed by ionized material blown off the companion. The optical emission is dominated by the tidally-distorted companion resulting in periodic variability (Section 1.5.3), and stellar absorption lines can be measured to track the secondary’s motion around the invisible MSP (Section 1.5.4.1). X-ray emission in this state is often modulated on the binary orbital period due to an intrabinary shock (Section 1.5.2). In the disk state, the pulsar wind can no longer hold off the mass inflow from the companion and an accretion disk forms around the MSP. The radio pulsar becomes undetectable, and the system becomes brighter in the optical/X-rays. The optical emission is dominated by contributions from the disk; the secondary’s absorption lines are often “filled-in” by disk emission and double-peaked Balmer emission lines, a hallmark of Doppler motion in the outer edges of an accretion disk, are observed in the optical spectra (Section 1.5.4.3). Occasionally, a flat-spectrum radio jet (purple) is observed in this state as material is ejected by the rotating pulsar magnetosphere. Figure credit: NASA’s Goddard Space Flight Center So far, three tMSPs have been observed to undergo this transition, proving the long-suspected evolutionary link between radio MSPs and their LMXB progenitors. This special subclass of MSPs is remarkable because they give one view of the end of the MSP recycling process, but also because 19 they provide excellent laboratories to study the physics of low-level accretion onto magnetized compact objects, the interactions between energetic pulsars and their stellar and gaseous environ- ments, and the typical evolutionary paths leading to the general MSP population. It is not entirely clear how and when these systems switch between states, although it is likely related to interactions between the infalling material and the rotating pulsar magnetosphere. These interactions can be better understood with a brief overview of what determines the accretion state around a rotating pulsar. In general, there are three main regimes that depend on where the inner edge of the accretion disk is truncated with respect to the rotating pulsar magnetosphere (e.g., Lipunov, 1987). If the inner disk radius (in) is larger than the light cylinder radius (i.e., the radius at which Keplerian corotation with the spinning neutron star equals the speed of light, see Figure 1.1) lc =  2 , (1.3) then the pulsar’s magneto-dipole radiation can continue to emit unimpeded and the radio pulsar will shine normally, being powered by the spin-down energy of the neutron star. The only way for disk material to be accreted onto the neutorn star surface is if the inner disk is moving faster than the pulsar’s rotating magnetic field at that radius. This implies that the mass inflow rate from the disk is high enough that in is smaller than the corotation radius , (1.4) (cid:16)   2 (cid:17)1/3 co = 42 the radius at which matter in Keplerian orbit corotates with the spinning neutron star of mass  and spin period . Only when the disk can get within this boundary can material flow freely onto the neutron surface along the magnetic field lines. The last regime occurs when the mass inflow rate from the disk is high enough that the inner disk punches down within lc, but not far enough that it is inside the corotation radius (i.e., co < in < lc). In this situation, the inflowing disk matter is actually expelled from the system since the pulsar magnetic field lines are spinning faster than the gas at the inner edge of the disk. This is commonly referred to as the “propeller” phase, where the disk may become entirely disrupted by the rotating magnetosphere (Illarionov & Sunyaev, 1975). 20 1.4.4 A New Spider Subclass: MSP binaries with Giant Companions Although tMSPs provide compelling science cases for studying accretion onto neutron stars, this subclass does not provide valuable insights when trying to understand the formation of all millisec- ond pulsars. In fact, all the known tMSPs have short orbital periods (<0.5 day), and evolutionary models suggest their lives will not end as the typical MSP binaries, which have white dwarf compan- ions with long orbital periods (e.g., Podsiadlowski et al., 2002, see also Section 1.4.1). Likewise, nearly all the confirmed redbacks and black widows have periods (cid:46)1 day, and so they too are unlikely to evolve into wide-orbit MSP-white dwarf binaries (see Section 1.4.2), which dominate the observed population of MSPs in the field (Figure 1.5). Instead, the progenitors to the canonical MSP binaries are likely neutron star binaries that have red giant companions in relatively wide orbits ( (cid:38)1 day, Tauris & Savonije, 1999). These long-period binaries have companions that have evolved off the main-sequence, where accretion onto the neutron star began once the envelope of the secondary expanded during its red giant phase. Unlike redbacks and black widows, the orbital periods in these systems will increase over time, and so they likely represent the late phases of typical MSP binary formation in the Galactic field. Two of these systems with giant companions are presented in detail in Chapters 2 and 3, present- ing a subclass which had previously been unobserved. Given their wide-orbits, giant companions, and inferred evolutionary tracks, these systems do not fit into the existing subclasses of black widow or redback MSP binaries. Instead they suggest a new subclass of neutron star binaries, aptly named “huntsmen,” after the giant Australian spider, and are the likely progenitors to the canonical MSP binaries. These huntsmen systems have been elusive, and proven very difficult to detect as radio MSPs, even compared to the redbacks and black widows (e.g., Camilo et al., 2016; Keane et al., 2018). This difficulty is likely due to strong winds from the evolved companions, which may eclipse the radio pulsar even more readily than for redbacks with main sequence-like companions. Although detect- ing radio pulsations in these huntsmen systems is challenging, the giant nature of the secondaries 21 mean they will be intrinsically more luminous in the optical. Observing these brighter systems not only enables much easier spectroscopic follow-up of the binaries, but also helps circumvent a major hurdle of studying neutron star binaries with non-degenerate companions: it is difficult to model the inclination (a key parameter for fully characterizing the binary (see Sections 1.5.3 and 1.5.4)) of systems for which the light curve is dominated by irradiation of the secondary. For both giants and more luminous main sequence stars, the lower ratio of high-energy to optical flux at the stellar surface means the effects of irradiation are simply lesser than for lighter, less evolved companions, allowing for more reliable constraints on system parameters. Further characterizations of similar systems would provide new insights on the relationship between MSPs and their low-mass X-ray binary progenitors. It is therefore important to continue discovering and studying new neutron star binary systems showing a wide range of phenomenology in order to address some of the key questions related to the population and evolution of MSPs. Fortunately, MSP binaries are excellent factories for accelerating particles to high-energies, and so they are usually strong emitters of GeV -rays. Concerted multiwavelength follow-up observations of these high-energy sources have made it possible to reveal numerous pulsars difficult to detect with radio observations and, following the techniques described in the following section, have lead to an enormous increase in the number of known MSP binaries with hydrogen-rich companions. These systems have been shown to host pulsars with unusually large neutron star masses with respect to the bulk of the millisecond pulsar distribution (Chapter 6), and so present compelling case studies for constraining the physical limits of the neutron star spin-up process. 1.5 Observations to Find and Characterize New Millisecond Pulsars 1.5.1 Gamma-rays Observations with the Fermi Gamma-ray Space Telescope have revealed that pulsars are nearly ubiquitous emitters of GeV -rays (e.g., Abdo et al., 2013a). In fact, recent multiwavelength follow-up observations of previously unidentified Fermi sources have been a boon for discovering 22 new systems, particularly new binary MSPs with non-degenerate companions (see Section 1.4.2) that were not visible in large radio pulsar surveys due to eclipses from companion material (e.g., Kong et al., 2012; Strader et al., 2015; Li et al., 2016; Halpern et al., 2017; Li et al., 2018). However, observing in -rays is not an easy task, because in contrast to visible light, high-energy -rays cannot be refracted by a lens or focused by a mirror. Instead, the workhorse instrument on Fermi, the Large Area Telescope (LAT), uses technologies similar to those in terrestrial particle accelerators, in which a high-energy photon interacts with detector material, producing an elec- tron/positron cascade. In order to perform astronomy with the LAT, the propagation direction of incoming -ray photons are only inferred from models of the resulting electromagnetic cascade. Therefore, the arrival direction takes the form a probability distribution on the sky referred to as the Point Spread Function (PSF). A larger PSF means the possible locations of a specific -ray point source take up a larger area on the sky than if the PSF were smaller. This area on the sky, known as an “error ellipse,” can be anywhere from ∼0.001–1 square degrees in size depending on the source brightness and location on the sky, and may contain hundreds of sources at other wavelengths. In the most-recently released 4FGL catalog6 (the fourth Fermi-LAT -ray source catalog based on 10 years of survey data since August 2008), ∼1600 sources (∼29% of the full catalog) are still not yet associated with known astrophysical phenomena at other wavelengths. Although the majority of these unidentified sources are likely to be Active Galactic Nuclei (AGN) associated with the supermassive black holes at the centers of very distant galaxies, many of these mysterious sources are likely Galactic neutron star binaries in our own Milky Way, just waiting to be discovered. The most straightforward way to discover a new MSP among an unassociated Fermi source is to search for pulsations in the radio. Large single-dish instruments like the Green Bank and Arecibo telescopes have helped discover many new young pulsars and MSPs in this way among the hundreds of unidentified Fermi sources (e.g., Stovall et al., 2014; Bates et al., 2015; Keane et al., 2018). By slowly and methodically scanning over a Fermi error region with a sensitive radio telescope, new 6https://fermi.gsfc.nasa.gov/ssc/data/access/lat/10yr_catalog/ 23 pulsar systems can be discovered “blindly” (Section 1.1). Although pointed radio observations are the most straightforward way to discover a MSP, as described in Section 1.4.2 MSP binaries with non-degenerate companions are much more difficult to find in blind radio pulsation searches due to the extensive eclipses from ionized material ablated from the companion. Therefore, concerted multiwavelength follow-up observations (mainly in radio, X-ray, optical, and infrared regimes) of previously unidentified -ray source regions have been promising, since MSP binaries show distinct phenomenology across the electromagnetic spectrum. In addition, if the orbit of a presumed companion can be measured (e.g., Section 1.5.4.1), deep radio pulsation searches can be performed at precisely localized positions on the sky (as opposed to the blind searches mentioned above), and at favorable binary phases when the MSP is less likely to be eclipsed by intervening material. 1.5.2 X-rays One way to link the high-energy -ray emission to a more localized source is to observe in X-rays, since the spatial resolution of most X-ray telescopes is on the order of arcseconds (as opposed to several arcminutes with Fermi). As mentioned in Section 1.4, if accretion is ongoing X-rays can be emitted via accretion onto the neutron star surface or from thermal emission from the hot inner disk. This is the case for normal neutron star LMXBs and is a relatively “low-cost” way to identify likely candidate MSP binaries or their progenitor systems among unidentified Fermi sources. If a system is not actively accreting, characteristic X-ray emission is still present and can be used to quickly localize sources. But, instead of emission coming from a disk, these X-rays likely come from non-thermal emission from an intrabinary shock. At the interface between the rela- tivistic pulsar wind and the irradiation-driven wind from the companion, a strong shock can form between the neutron star and its companion (see Figure 1.8). Emission from this shock is typically modulated on the orbital period of the binary, resulting in a range of phenomenologies based on the relative strengths of the winds and the properties of the orbit. Deep X-ray observations, combined 24 with sophisticated models for these shocks, can constrain important properties of the system, such as the shock geometry and inclination of the binary. Figure 1.8: Diagram of an intrabinary shock. A cartoon showing the formation of an intrabinary shock, formed at the collision interface between the pulsar and companion winds. The relative strengths of the winds determine whether the shock wraps around the pulsar or the companion. Orbitally modulated X-ray emission is often observed due to the presence of the shock (see text and Figure 1.9). Figure from Venter et al. (2015). For black widow and redback systems, the X-ray light curves show distinct phenomenology depending on the orientation of the shock and the orbital properties of the system. In the simplified intrabinary shock models of Romani & Sanchez (2016), the controlling physical parameter that shapes the X-ray light curve is the ratio between the pulsar and companion wind momentum fluxes. If the pulsar wind dominates, the shock wraps around the companion, resulting in two peaks in the X-ray light curve near companion inferior conjunction (i.e., when the companion is between Earth and the pulsar), where the peaks are due to Doppler boosting of the shock when it points along the line of sight to the observer. If the companion wind dominates, the shock instead wraps 25 around the pulsar and the phases of the peaks shift by half an orbital phase. The relative amplitudes and symmetry of the two peaks depend on the velocity of the companion wind with respect to the orbital velocity, while the system inclination affects the overall amplitude of the peaks (i.e., more face-on systems will show smaller variations than systems that are edge-on). It should be noted that because the secondaries in these systems are tidally-locked, they rotate rapidly when compared to similar, isolated stars. This rapid rotation, together with a strong temperature gradient induced by heating the tidally-locked inner face, results in a strong dynamo effect that can greatly enhance the magnetic field on the companion surface, leading to strong magnetically driven stellar winds and facilitating the formation of a powerful shock. Example X-ray light curves from two spider MSPs that have strong intrabinary shocks are presented and described in Figure 1.9. Figure 1.9: Example X-ray light curves due to strong intrabinary shocks. Example X-ray light curves showing the effects of the intrabinary shock for a black widow MSP binary with a weak companion wind (left), and a redback with a larger companion and a stronger wind (right). The data is repeated over two orbital cycles for visual clarity. The two peaks in the light curves are due to Doppler boosting of the shock when it is pointed towards the observer. In the black widow system, the shock wraps around the companion such that the peaks occur just before and just after the companion is in front of the MSP (pulsar superior conjunction). For the redback system, the companion wind is dominant so the shock wraps around the pulsar and the peaks are shifted by half an orbit (see text). Figure from Romani & Sanchez (2016). 26 1.5.3 Optical Variability Another way to infer the presence of a candidate spider MSP binary is to look for optical variability consistent with tidal deformation of the secondary star. If a non-degenerate companion closely orbits a MSP, the strong gravity from the massive neutron star can pull on and stretch the stellar en- velope, giving it a teardrop shape (see description for Roche lobe in Section 1.3). As the companion orbits, the projected area of the star as seen by an observer at Earth will change depending on the phase of the orbit. If the orbit is edge-on with respect to our line of sight, then we would expect to see two roughly equal maxima in the optical light curves per orbital period, corresponding to when the maximum projected area of the star is observed at quadrature phases, and two minima corre- sponding to the conjunctions. These minima are typically unequal due to the effects of gravity and limb darkening when viewing toward the inner Lagrange point of the tidally locked companion (see Figure 1.10). These effects are known as “ellipsoidal variations” and can be used in practice to find new MSP binaries, since most -ray emitting compact binaries with Roche lobe filling secondaries should have this detectable periodic optical variability. Since a limited number of variable stars are present within the error regions of unidentified Fermi sources, following-up these optical variables presents another “low-cost” method of identifying new candidate spider MSPs. It should be noted that contributions from an accretion disk, or heating of the tidally-locked secondary by the pulsar, can significantly complicate the variability we expect due to ellipsoidal variations (see below). Perhaps most importantly, modeling these ellipsoidal variations using sophisticated light curve modeling tools can constrain the inclination of the binary system, where the maximum effect is ob- served if the system is edge-on ( = 90◦), and no modulations are expected for face-on inclinations ( = 0◦, Figure 1.10). This simple picture of a teardrop-shaped star manifesting in predictable double-peaked light curves is complicated by a process alluded to in Section 1.4.2, known as irradiation (i.e., heating from the pulsar). The companion star is tidally-synchronized on short timescales ((cid:46) 104 years for the orbital periods and mass ratios typical of spider MSPs, Zahn, 1977), so its inner side always 27 Figure 1.10: Diagram of ellipsoidal variations. Schematic demonstration of determining the binary inclination with ellipsoidal variations of a Roche lobe filling secondary (teardrop) and a compact object (small black circle). There are two equal maxima per orbit when the tidally distorted secondary has the largest projected surface area, and two unequal minima per orbit. The “dayside” of the secondary is nominally predicted to have the faintest optical flux due to gravity darkening (a lower effective temperature), but this effect can be reduced or even reversed if this side is also irradiated by the compact object (see text). The amplitude of the variations are largest when the system is edge-on ( = 90◦) but are undetectable if face-on ( = 0◦), with intermediate inclinations producing intermediate amplitude variations. Figure adapted from Greene et al. (2001). faces the same direction towards the neutron star. This means that direct heating by the relativistic pulsar wind will always be impinging on the same cross-section of the star. This will heat the inner face (often referred to as the companion’s “dayside”) to a temperature much hotter than the star’s cooler “nightside.” The presence of an intrabinary shock likely facilitates this heating by channeling the pulsar wind indirectly to the stellar surface along the enhanced magnetic field lines of the star and of the shock itself (e.g., Sanchez & Romani, 2017). In some MSP binaries with non-degenerate companions, particularly ones with short orbital periods (and therefore small orbital separations), this irradiative heating can dominate the effects of ellipsoidal variations. Instead of a double-peaked optical light curve with peaks at quadrature phases, the light curve would show a broad single-peaked maximum, centered on superior con- junction of the companion (i.e., the neutron star is between Earth and the secondary), when we observe the heated dayside of the companion. Half an orbit later, as the secondary crosses in front of the neutron star, we would see the cooler nightside (backside) of the companion, corresponding to a minimum in the light curve. This effect is seen clearly in the two redback-like binary systems 28 i=90°° presented in detail in Chapter 4 and Chapter 5. Combining the effects of irradiative heating and ellipsoidal modulations in order to model the optical/near-IR light curves is one of the few ways to constrain the inclination of the binary. This measurement is vital for interpreting, and ultimately, fully characterizing new MSP binaries, since the system inclination, together with spectroscopic orbital information (see Section 1.5.4), allows a direct estimate for the mass of the neutron star and its orbiting companion. 1.5.4 Optical Spectroscopy 1.5.4.1 Radial Velocity Tracking Once a candidate MSP binary is found via the radio, X-ray and/or optical observations introduced above, the best way to confirm the presence of the massive invisible primary is to track the secondary in its orbit. This is done with optical spectroscopy. By measuring the positions (i.e., wavelengths) of the stellar absorption lines intrinsic to the photosphere of the companion, and comparing these to the laboratory values measured on Earth, a radial velocity can be calculated. Once multiple measurements are performed, a radial velocity curve can be built, corresponding to the secondary’s radial motion around the center of mass of the system. Typically, a Keplerian model is fit to these radial velocity data in order to estimate key orbital parameters of the binary, such as the orbital period orb, eccentricity , and the radial velocity semi-amplitude 2 (half the peak-to-trough amplitude of the radial velocity curve). The semi-amplitude and orbital period of the binary immediately yield a key quantity for single-line spectroscopic binaries, known as the binary mass function  (): (cid:16)(1 − 2)3/2(cid:17)  () = 3 2 orb 2 . (1.5) Following Kepler’s Third Law, if one assumes a circular orbit ( = 0), this equation can be rewritten as , (1.6)  () = 1sin3(i) (1 + )2 29 where  is the gravitational constant, 1 is the mass of the invisible primary,  is the system inclination, and  = 2/1 is the binary mass ratio. The fact that this value can be measured solely from the observable quantities 2, orb, and  is powerful, and sets a lower limit on the mass of the invisible primary object by assuming  = 90◦ and  → 0. 1.5.4.2 Determining the Mass Ratio Armed with the binary mass function, the only remaining quantities that are needed in order to completely solve for the component masses of the system are the mass ratio  and the orbital inclination . While the inclination can be estimated by modeling the optical/near-IR light curves (Section 1.5.3), the mass ratio can be a more difficult quantity to measure. If the suspected radio pulsar is already known and can be precisely timed over many epochs, then the mass ratio can be measured directly by inspecting the pulsar timing residuals, allowing an observer to infer the relative masses of the two components based on the pulsar’s motion around the common center of mass. However, as is the case for most MSP binaries with non-degenerate companions, the radio pulsar is elusive, and if no radio timing solution is possible, other routes must be taken to estimate the mass ratio. Fortunately, optical spectroscopic observations provide another avenue for measuring this quan- tity. As a consequence of Roche geometry, measuring the projected rotational velocity of the secondary, rot sin(), combined with the secondary semi-amplitude 2 immediately returns an es- timate for the mass ratio. The reason this is possible is due to the assumption that the secondary is tidally-locked to the suspected neutron star. As described above, this tidal synchronization happens on a short timescales (see Section 1.5.3), and typically will be much less than the time it takes to spin-up the MSP through accretion. By definition, a tidally-locked star will have a rotational period equal to the binary orbital period, and for typical MSP binaries with non-degenerate companions, this implies the star is rotating rapidly with respect to similar, isolated stars (as an example, the Sun has a spin period of 30 ∼25 days, whereas similar type stars orbiting MSPs have periods of <1 day). This rapid rotation will cause the stellar absorption lines to appear much wider than they normally would, a phenomenon known as Doppler broadening (Shajn & Struve, 1929). Using high-resolution spectroscopy, the width of the secondary’s absorption features can be compared to those of a standard star with similar spectral type to infer the secondary’s projected rotational velocity(cid:0)rot sin()(cid:1). Combining this result with the radial velocity semi-amplitude 2, one can directly determine the mass ratio  = 2/1 by solving the following equation for close binaries using the Roche lobe approximation of Eggleton (1983): (1.7) rot sin(i) = K2 0.49 q2/3 (1 + q) 0.6 q2/3 + ln(1 + q1/3) . This value for , together with a spectroscopic orbital solution and an estimate for the binary inclination (see Sections 1.5.3 and 1.5.4.1), finally allow the component masses of the system to be determined. 1.5.4.3 Spectroscopy if Accretion Disk is Present In addition to tracking the secondary’s orbital motion, optical spectroscopic observations can allow one to infer other characteristics of a presumed neutron star binary, including the presence of an accretion disk. If an accretion disk is present, typically this component dominates the optical light from the system, making it difficult to detect regular, periodic variations in the photometry (e.g., ellipsoidal variations, Section 1.5.3). The presence of a disk also complicates measurements of absorption lines from the secondary, as the disk emission often “fills in” these weak photospheric absorption features. One hallmark of an accretion disk is a double-peaked H profile (as well as additional Balmer or Helium lines in some systems). This emission structure is consistent with Doppler motions in the outer edges of the disk, where the receding half of the disk is responsible for the redder peak of the profile while the approaching side results in the bluer peak (Figure 1.11). 31 Figure 1.11: Example double-peaked H emission line due to an accretion disk. A high- resolution spectrum showing a double-peaked H emission line coming from an accretion disk. Due to the Doppler effect, the redward peak comes from the receding half of the disk, while the blue peak comes from the approaching side. This double-peaked Balmer emission is generally interpreted as arising from an accretion disk. However, an intrabinary shock, and material blown from a companion star and being heated by a pulsar wind, could provide an additional natural physical explanation for complex H profiles, as shown in Chapter 3 and Chapter 4. Cartoon credit: ESA; Spectrum taken with the PST1 telescope (Baranowski et al., 2009). Since a number of accreting MSPs and tMSPs in their disk-state have shown these double- peaked features, it was believed that nearly all neutron star binaries that showed similar structures in their H profiles must also contain accretion disks. However, some redback-like systems show similar phenomenology but are unlikely due to disk emission, particularly in the MSP binary pre- sented in Chapter 3. Instead, an intrabinary shock could provide a natural physical explanation for complex H profiles. As the shock orbits along with the secondary, parts of the shock may trail in the orbit as the pulsar wind pushes it in the opposite direction (see Figures 1.6 and 3.7). This can 32 Ω lead to H profiles that appear double-peaked, but are not due to a disk. 1.6 This Thesis In this thesis, I present an in-depth look at four neutron star binaries that represent new dis- coveries among the unassociated Fermi -ray source population. Among these are three systems that are the likely progenitors to the canonical MSP binaries. These field MSP binaries with giant companions in wide orbits represent a new subclass of spider MSPs and provide new insights into the late stages of the MSP recycling process. In the following chapters, I describe the discovery of these systems and also provide an overview of the range of multiwavelength phenomenology we expect from MSP binaries with non-degenerate companions. Uncovering these new redback and huntsmen MSP binaries among bright unassociated Fermi sources suggests that identifying and characterizing new and interesting systems like these will continue to be fruitful in the months and years to come. In Chapter 2 we present the discovery of the likely optical stellar counterpart to the Fermi-LAT -ray source 2FGL J0846.0+2820. This source is an optical variable coincident with a variable - ray source. Using optical photometry and spectroscopy, we found this star is a partially-stripped red giant companion with a mass of ∼ 0.8(cid:12) in a wide (8.1 d) orbit around a heavy, invisible primary object of ∼ 2(cid:12). The inferred mass and association with the -ray source strongly suggests the primary is a MSP in the late-stages of pulsar recycling, however no radio pulsar has been found to date despite continued monitoring. The giant nature of the companion may explain the difficulty in detecting radio pulsations, since the tidally-locked red giant likely has a strong magnetically driven stellar wind that is efficiently absorbing or scattering the pulsar emission. Discovery of this source shows that such positional coincidence searches between optical variables and unassociated Fermi sources can be productive and provide a relatively simple way to identify new compact binaries in the Milky Way. This source was detected with high significance by Fermi in the first year of the mission, but the -ray flux dropped substantially in mid-2009 and since then it has been unde- 33 tectable by Fermi. Since two of the three tMSP systems have shown large -ray flux changes when transitioning between states, it is conceivable that 2FGL J0846.0+2820 transitioned from a disk to a pulsar state during this time. This flux change was accompanied by an increased variation in the optical brightness, perhaps due to a changes in the accretion flow. H emission is also observed from this source, perhaps consistent with a faint accretion disk or an intrabinary shock. Although the evidence for an accretion disk in this system is mixed, the evolutionary fate of the system is certain; eventually the low-mass companion will become a white dwarf in a wide orbit around the presumed MSP primary, and so this system represents a previously unobserved subclass of neutron star binaries that are the progenitors to the canonical field MSPs: the huntsmen. In Chapter 3 we describe the second system in the emerging huntsmen subclass of MSP bi- naries: 1FGL J1417.7–4407. Using new radio, infrared, optical, and X-ray data, we present the unique multiwavelength phenomenology we observe that allow us to better understand this unusual system. The counterpart to this persistent -ray source was initially discovered using optical and X-ray data to infer the presence of a compact X-ray binary with a neutron star primary and a red giant companion in a ∼5.4 day orbit (Strader et al., 2015). The MSP nature of the source was confirmed when a 2.66 ms radio pulsar was detected at the same location as the optical/X-ray source (Camilo et al., 2016). Initial spectroscopic measurements of complex, double-peaked H emission suggested the presence of an accretion disk. However, the simultaneous discovery of the radio pulsar and our subsequent finding that the radio spectral index (based on VLA/ATCA obser- vations) is broadly consistent with a MSP disfavors a scenario where an accretion disk is present. Instead, we use optical spectroscopy to show that the complex morphology of the persistent H emission line can be better explained by the presence of a strong, magnetically driven stellar wind from the giant secondary and its interaction with the pulsar wind. This scenario is consistent with our Chandra and Swift X-ray observations, which show evidence for a double-peaked orbital light curve and a hard spectrum, similar to that observed in some redback MSP binaries and likely due to an intrabinary shock. The X-ray luminosity we infer is also slightly more luminous than all known redback systems in the rotational-powered pulsar state, perhaps due to the wind from the giant 34 companion. In this context, the source of the H emission likely comes from some combination of material flowing to the inner Lagrange point, hot gas swept out beyond the companion’s orbit by the pulsar’s radiation pressure, and the intrabinary shock. A substantial magnetic field generated by the tidally-locked giant may facilitate the formation of such a luminous shock. Similar to 2FGL J0846.0+2820, the giant companion and inferred evolutionary track makes this system different from the redback MSPs, presenting another system in the newly discovered subclass of MSP bina- ries with giant companions in wide orbits, the progenitors to the canonical field MSP binaries. In Chapter 4, we present new optical photometry and spectroscopy of the MSP binary PSR J1306–40, initially discovered in a blind radio survey for MSPs among unassociated Fermi sources (Keane et al., 2018). Our optical light curve models in conjunction with our spectroscopic results imply a minimum neutron star mass of 1.75 ± 0.09 (cid:12) (1-sigma) and a high, nearly edge-on incli- nation. These data also suggest that the companion is filling a substantial fraction of its Roche lobe and that it has a somewhat evolved, subgiant-like radius, luminosity, and gravity. H line profiles seen in the optical spectra switch quickly from persistent emission to absorption during pulsar su- perior conjunction, suggesting that the H emitting region is eclipsed when the companion crosses between Earth and the pulsar, consistent with emission from an intrabinary shock. Additional evidence for a strong shock comes from the X-ray light curve, which shows one maximum and one minimum per orbital period modulated on the binary orbital period (Linares, 2018). This can be attributed to emission from an intrabinary shock that becomes at least partially eclipsed at the same phases as the H emission. Given that PSR J1306–40 has been spun up to obtain its rapid spin period (2.2 ms), this would put PSR J1306–40 on a standard evolutionary track of low-mass X-ray binaries that started the mass transfer process after leaving the main sequence. This means this object, similar to 2FGL J0846.0+2820 and 1FGL J1417.7–4407, is in the late stages of the MSP recycling process that will terminate with a wide orbit MSP–white dwarf binary. In Chapter 5, we use archival X-ray data and new optical photometric and spectroscopic observations to present the discovery of a likely new redback millisecond pulsar binary associated 35 with the Fermi -ray source 4FGL J2333.1–5527. This system hosts a low-mass main sequence- like companion in a 6.9-hr, highly inclined orbit around a suspected massive neutron star primary. Although the radio pulsar has not yet been found, these data tentatively suggest the mass of the neutron star is likely well in excess of ∼ 1.4 (cid:12) assuming the secondary mass is typical of the known redbacks. The single peak in the phased optical light curves imply the companion star is being heated substantially on its tidally-locked dayside and further supports the idea that optical variables inside unassociated Fermi regions present compelling sources to follow-up with additional multiwavelength observations in an effort to discover new compact neutron star binaries. 36 CHAPTER 2 2FGL J0846.0+2820: A NEW NEUTRON STAR BINARY WITH A GIANT SECONDARY AND VARIABLE -RAY EMISSION 1We present optical photometric and spectroscopic observations of the likely stellar counterpart to the unassociated Fermi-Large Area Telescope (LAT) -ray source 2FGL J0846.0+2820, selected for study based on positional coincidences of optical variables with unassociated LAT sources. Using optical spectroscopy from the SOAR telescope, we have identified a late-G giant in an eccentric ( = 0.06) 8.133 day orbit with an invisible primary. Modeling the spectroscopy and photometry together lead us to infer a heavy neutron star primary of ∼ 2(cid:12) and a partially stripped giant secondary of ∼ 0.8(cid:12). H emission is observed in some of the spectra, perhaps consistent with the presence of a faint accretion disk. We find the -ray flux of 2FGL J0846.0+2820 dropped substantially in mid-2009, accompanied by an increased variation in the optical brightness, and since then it has not been detected by Fermi. The long period and giant secondary are reminiscent of the -ray bright binary 1FGL J1417.7–4407, which hosts a millisecond pulsar apparently in the final stages of the pulsar recycling process. The discovery of 2FGL J0846.0+2820 suggests the identification of a new subclass of millisecond pulsar binaries that are the likely progenitors of typical field millisecond pulsars. 2.1 Introduction Neutron stars in low-mass binaries can accrete matter and angular momentum from a non- degenerate companion star and be recycled to very fast spin periods, making them detectable as millisecond pulsars (MSPs; Alpar et al., 1982b). During active accretion, the system is observable as a low-mass X-ray binary. As the orbital period grows and accretion eventually ends, the neutron star turns on as a radio MSP powered by the spindown energy of the neutron star. 1This chapter is taken mostly word-for-word from Swihart et al. (2017), as published in the Astrophysical Journal 37 Most MSPs in the Galactic field are binaries with a degenerate white dwarf companion of 0.2- 0.3 (cid:12) and orbital periods in the range of days to weeks (Tauris & van den Heuvel, 2006). These are the end products of the recycling process. Recent advances—especially follow-up of sources discovered with the Fermi Large Area Telescope (LAT) in GeV -rays—have allowed the discovery of a subclass of MSP binaries with non-degenerate companions in which recycling is apparently not yet complete. These sources, which typically show radio eclipses, are categorized based on their companion’s mass. Black widow systems have very light companions ( (cid:46) 0.08(cid:12)) that are being actively ablated by the wind of the pulsar, while redbacks have non-degenerate, nearly Roche-lobe filling, main sequence companions of mass  (cid:38) 0.2(cid:12) (Roberts, 2011). Recently, at least three of these redbacks have been found to switch between accretion- powered X-ray binary states and rotationally-powered pulsar states on timescales of days to months (Archibald et al., 2009; Papitto et al., 2013; Bassa et al., 2014; Roy et al., 2015; Bogdanov et al., 2015; Johnson et al., 2015). It is not clear how and when these systems switch on or off as a radio MSPs, nor the cause of the rich phenomenology observed when an accretion disk is present, though both are likely related to the interaction between the inner accretion flow and the pulsar magnetosphere. These transitional MSPs are key systems for understanding the physics of the pulsar recycling process, the interactions between energetic pulsars and their binary companions and surroundings, and the typical evolutionary paths leading to the general MSP population. Cur- rent evolutionary models predict such transitions on timescales comparable to the billion-year-long evolution of the binary, not in weeks to months as is observed in the known transitioning systems (Benvenuto et al., 2015). Furthermore, these transitional MSPs all have short orbital periods ((cid:46) 0.5 days); while the evolutionary endpoint of these systems is uncertain, they will likely not be typ- ical of MSP binaries in the field, which have degenerate companions with periods of days to months. In this work, we report observations of a bright -ray source, 2FGL J0846.0+2820, studied as part of an ongoing survey of unassociated Fermi-LAT sources (Strader et al., 2014, 2015). Using follow-up optical observations of the presumed companion, including photometry and spec- 38 troscopy, we find that 2FGL J0846.0+2820 is likely a Galactic compact binary with a massive neutron star primary and a giant secondary in a relatively wide orbit. The long orbital period, giant companion, and component masses are remarkably similar to those of the recently discovered -ray bright binary 1FGL J1417.7–4407 (Strader et al., 2015), which was independently found to host a radio MSP (Camilo et al., 2016). Although no pulsar has yet been discovered to be associated with 2FGL J0846.0+2820, our observations provide evidence of a second system in a new subclass of long-period -ray bright binaries with heavy neutron star primaries and giant secondaries. These systems are plausible progenitors for typical MSP–white dwarf binaries observed in the field. 2.2 Observations 2.2.1 Fermi-LAT Source The -ray source was first detected as 2FGL J0846.0+2820 (Nolan et al., 2012a), listed in the second full catalog of Fermi-LAT sources, based on the first two years of LAT data obtained from 2008 August to 2010 August using the earlier P7V6 instrument response functions (IRFs). 2FGL J0846.0+2820 was one of 774 new -ray sources that did not appear in the 1FGL catalog of LAT sources (Abdo et al., 2010) based on overlapping 95% source location confidence contours. The overall significance of 2FGL J0846.0+2820 was 4.1, just above the 2FGL catalog detection threshold of test statistic, TS > 25 (Mattox et al. 1996, for four total degrees of freedom – two positional and two spectral; see §3.2 in Nolan et al. 2012a for a description). With the limited statistics, the average LAT spectrum was best fit as a single power law with photon index Ɖ = 2.51 ± 0.20 and a 0.1–100 GeV flux,  = (1.1 ± 0.3) × 10−8 photons cm−2 s−1, corresponding to a luminosity of  ≈ (4.0 ± 1.0) × 1034 ( / 8.1 kpc)2 erg s−1 (see § 2.4.3.1 for distance estimate), which is comparable to the luminosity of 1FGL J1417.7–4407 (∼3 × 1034erg s−1, Strader et al., 2015) The 2FGL catalog 1-month binned light curve showed only 95% confidence upper limits indicating each point has TS < 10 or flux error uncertainty > 50% of the flux value. We present an updated analysis of the recent Pass 8 Fermi-LAT data in §2.4.2. 39 2.2.2 Optical and Near-IR Counterpart The -ray source 2FGL J0846.0+2820 was selected for follow-up study based on a search of positional coincidences of periodic optical variables found in the Catalina Sky Surveys Data Release-1 (CSDR1; Drake et al., 2014) catalog with Fermi LAT (Atwood et al., 2009) sources. We identified the  ∼ 15.7 Catalina source CSS J084621.9+280839, with a USNO B1.0 catalog (Monet et al., 2003) J2000 sexagesimal position of (R.A., DEC.) = (08:46:21.89, +28:08:41.0), as a periodic variable 0.◦215 offset from the 2FGL centroid. This object is also associated with a 2MASS point source with  = 14.21± 0.02,  = 13.59± 0.03,  = 13.50± 0.03 mag (Cutri et al., 2003). The CSS source is also listed in the Sloan Digital Sky Survey as J084621.87+280840.8 with magnitudes:  = 18.24,  = 16.51,  = 15.74,  = 15.42,  = 15.22 (Ahn et al., 2012). 2.2.2.1 Catalina Sky Survey (CSS) To analyze the CSS photometry of this variable, we retrieved 471 CSS photometric measurements of the variable optical counterpart taken between 2005 April 10 and 2013 September 25. We searched for periodicity using a Lomb-Scargle periodigram (Scargle, 1982), finding peak power at a period of ∼16.2 d, agreeing with the 16.1866 d period found by Drake et al. (2014). When phased on this period, small periodic flux modulation is evident. However, we have found this period to be an alias of the real orbital period determined via spectroscopy (§2.4.1). The spectroscopic orbital period is 8.1328 d. Augmented by additional photometry, we analyze the phased light curve of the system in §2.4.3. The long-term CSS light curve offers the opportunity to study whether the optical flux from the system has changed significantly since 2005. After removing five measurements that were >3 outliers, there are a total of 466 measurements. The median magnitude is equiv = 15.72 and the median uncertainty 0.06 mag. The long-term CSS light curve reveals a monotonic increase in the mean system brightness over the past ∼ decade. A least-square linear regression fit to the data shows the system has brightened at a rate of 0.011± 0.001 mag yr−1. We discuss this trend as well 40 Table 2.1. PROMPT Photometry of 2FGL J0846.0+2820 JD-2450000 Band Mag (d) 7057.65342 7057.65421 7057.65500 7057.65579 7057.65658 7057.65579 7057.65658 · · · B B B B B B B · · · 17.001 16.366 16.141 16.683 16.755 16.683 16.755 · · · Err 0.289 0.157 0.153 0.214 0.225 0.214 0.225 · · · Note. — This table is published in its en- tirety in machine-readable format. A por- tion is shown here for guidance regarding its form and content. These magnitudes are not corrected for extinction. as the strange phenomenology present in the phased CSS light curves in § 2.4.5. 2.2.2.2 PROMPT We obtained time series photometry of the optical source in , , and  bands with the 16-inch PROMPT-5 telescope (Reichart et al., 2005) at Cerro Tololo International Observatory between 2015 February 4 and 2015 June 7. Each observing night consisted of multiple 60-second exposures of the target field, which included the candidate optical counterpart to 2FGL J0846.0+2820 as well as five nearby comparison stars. We performed differential aperture photometry to obtain instrumental magnitudes of the target source using the five (non-variable) comparison stars as a reference. We calibrated our instrumental magnitudes using observations of the Landolt (1992) standard star field RU149. Our final sample includes 1643 photometric measurements in , 1038 in , and 500 in , with mean magnitudes  = 16.73,  = 15.95, and  = 15.40. All the PROMPT photometry can be found in Table 2.1. 41 2.3 Optical Spectropscopy 2.3.1 SOAR and MDM Spectroscopy We began spectroscopic monitoring of the source with the Goodman Spectrograph (Clemens et al., 2004) on the SOAR 4.1-m telescope on 2014 Dec 8, continuing through 2016 Dec 31. For the initial epochs we used a 2400 l mm−1 grating in the region of the Mg triplet and a 1.03” slit, yielding a resolution of ∼ 0.7 Å. The observations from 2015 December onward used a 2100 l mm−1 grating, covering a similar region of the spectrum, but with a slightly lower resolution of ∼ 0.9 Å. Two or three 600 s exposures were taken per epoch. All spectra were reduced in the standard manner, with wavelength calibration performed using FeAr arcs obtained after each set of spectra. A small number of spectra were obtained with a low-resolution 400 l mm−1 grating (resolution ∼ 5.8 Å) to check for evidence of emission. We measured barycentric radial velocities for the medium-resolution spectra through cross- correlation with bright template stars taken with the same setup. Given the long period of the system, at each epoch we determined the radial velocity as a weighted average of the two or three measurements from the individual 600 s exposures. This gave 19 independent SOAR epochs of velocities (with a 20th coming from the Keck/HIRES observation discussed below). These mea- surements are listed in Table 2.2. In some of our low-resolution spectra we see evidence for H in emission, though the line is generally weak. To better study the emission, we also obtained some low-resolution spectra with OSMOS on the Hiltner 2.4-m telescope at the MDM Observatory at Kitt Peak.These spectra were obtained in 2–3 exposures of 20 min each on eight epochs from 2016 October 20 to 2017 January 8. Reduced in the usual manner, they cover a usable wavelength range of ∼ 3960–6840 Å at a resolution of about 3.9 Å. These spectra demonstrate a wider range of H morphology than observed in the small sample of SOAR spectra, from an epoch with very deep H absorption (1.7 Å equivalent width on 2017 42 Table 2.2. Barycentric Radial Velocities of J0846 BJD (d) 2456999.8170356 2457003.8053888 2457003.8295105 2457012.8609441 2457037.7599617 2457071.6414565 2457119.5191800 2457120.5289206 2457158.4547783 2457166.4751144 2457170.4641598 2457378.8366717 2457417.7571464 2457463.5367120 2457473.6121242 2457495.5359481 2457508.5254554 2457657.1257481a 2457744.7835625 2457752.8115122 RV (km s−1) 52.5 32.5 33.2 63.7 83.0 97.0 92.4 96.3 36.9 37.1 58.6 66.0 3.4 44.1 -12.4 64.9 47.7 95.9 66.5 58.4 Err. (km s−1) 1.6 1.6 1.5 2.4 1.7 1.7 2.0 1.7 1.9 2.0 1.9 2.0 3.4 1.9 2.1 1.9 1.9 0.9 2.1 1.8 aThis epoch comes from the Keck/HIRES spectrum. Jan 5) to one with obvious H emission (about 1.0 Å equivalent width, and substantially larger if corrected for absorption, on 2016 Oct 21). These changes are accompanied by substantial overall changes in the spectral morphology, with a lower effective temperature implied in the Oct 21 spec- trum and a higher effective temperature in the Jan 5 spectrum. In Figure 2.1 we show a comparison of two OSMOS spectra taken several months apart that illustrate the extremes of varying effective temperature and H emission. In the 2016 Oct 21 spec- trum the effective temperature is lower, with more deeper atomic lines and a hint of the emergence of molecular features, along with clear H emission and fill-in of H. In the 2017 Jan 5 spectrum the effective temperature is clearly higher and there is no longer any H emission visible; the line appears to entirely be in absorption. Other low-resolution spectra are intermediate between these 43 extremes, though unfortunately we do not have enough spectra to track the temperature or emission effectively as a function of orbital phase. Figure 2.1: Optical spectra of 2FGL J0846.0+2820. Two of our OSMOS spectra showing the most extreme examples of varying effective temperature and H emission (rest = 6562.8 Å). Our other low-resolution spectra lie somewhere between these extremes. Analysis of the low-resolution spectra using the spectral classification program MKCLASS (Gray & Corbally, 2014) suggests a late-G spectral classification with a subsolar metallicity of –1.1 to –0.6. 2.3.2 Keck HIRES Spectroscopy If the secondary is tidally synchronized with the presumed central compact object, we would expect to see evidence of this through the broadening of spectral lines due to rapid rotation. Using the projected rotational velocity  sin  and orbital semi-amplitude, the mass ratio  = 2/1 can be 44 estimated. We found evidence for line broadening in our SOAR spectra, but the resolution was too low to precisely measure  sin . To address this, we obtained a high-resolution spectrum with HIRES (Vogt et al., 1994) on the Keck I telescope on 2016 Sep 25. The spectrum was taken using the C5 decker, yielding a nominal resolution of ∼ 36000, and covered a wavelength range of ∼ 3900–8100 Å. The single 1200-sec exposure was reduced using HIRedux (Bernstein et al., 2015). We additionally obtained a number of spectra of bright late-G to mid-K giant stars with the same resolution and binning to use as templates. We created a set of rotational convolution kernels (assuming a standard limb darkening law) for projected rotational velocities ( sin ) 10-100 km s-1. Then, on an order-by-order basis, we convolved the spectra of our template stars with the set of kernels and cross-correlated these broadened templates with each of the original, unbroadened spectra. We then fit relations between  sin  and the full-width at half-maximum (FWHM) for each pair of cross-correlations, excluding regions with very strong, broad lines. Finally, we cross- correlated the 2FGL J0846.0+2820 spectrum with the unconvolved template stars and used the FWHM values derived from these cross-correlations to estimate the projected rotational velocity. For  templates, this procedure produces 2 estimates of  sin  per order. In practice, we find that the dispersion in  sin  among templates is much smaller than the dispersion among orders, echoing a similar result found for star cluster velocity dispersions (Strader et al., 2009). The final value derived in this manner is  sin  = 23.2 ± 1.0 km s−1, where the uncertainty is the standard deviation of the measurements among all of the templates and orders. Further examination of the HIRES spectrum showed some evidence of chromospheric activity, and also provided us with an estimate of the foreground reddening from the strength of the NaD lines (Munari & Zwitter, 1997). The ( − ) estimate is 0.10, which is slightly higher than the 0.04 predicted from the Schlafly & Finkbeiner (2011a) all-sky map. 45 2.4 Results 2.4.1 Keplerian Orbit Fitting and Mass Ratio After correcting the observation epochs to Barycentric Julian Date (BJD) on the Barycentric Dy- namical Time system (Eastman et al., 2010), we performed a Keplerian fit to our radial velocity data using the custom Markov Chain Monte Carlo sampler TheJoker (Price-Whelan et al., 2017). We fit for the period , BJD time of periastron passage , eccentricity , argument of the pe- riastron , systemic velocity , and the semi-amplitude 2. The posterior distributions were generally close to normal, with equivalent best-fit orbital elements:  = 8.13284 ± 0.00043 d,  = 2457007.9589+0.2378 −0.2907 d,  = 0.061 ± 0.017,  = 81.◦8+10.7−12.9,  = 42.6 ± 0.7 km s−1, and 2 = 54.4 ± 1.0 km s−1. This fit is remarkably good, with an rms for the median values of 1.8 km s−1. The phased radial velocity curve is shown in Figure 2.2. We use the posterior samples from this fit to derive the mass function  ()  () = 3 2 2 = 1 (sin )3 (1 + )2 (2.1) for mass ratio  and inclination . We find  () = 0.136+0.008 −0.007 (cid:12). The mass ratio  can be directly determined using our estimate of  sin  and the semi-amplitude of the secondary using the standard formula  sin  = 0.462 2 1/3 (1 + )2/3 (Casares, 2001). Using the values presented above, this gives  = 0.402+0.034 −0.037. If the secondary fills its Roche lobe, the orbital period and mass ratio are used to calculate the mean density (Eggleton, 1983), yielding ¯ = 0.003 g cm-3. Together with an analysis of our low-resolution spectra, this result leads us to interpret the secondary as a late-G type giant. We defer a detailed discussion of the inclination to § 2.4.3 below, but note that for typical neutron star masses in the range 1.4–2.0 (cid:12), these measurements suggest an inclination in the approximate range 30–35◦. 46 Figure 2.2: Radial velocity curve for the optical counterpart to 2FGL J0846.0+2820 obtained in 20 epochs between 2014 December 8 and 2016 December 31. The high quality fit yields excellent constraints on the orbital ephemeris. 2.4.2 Fermi-LAT Analysis To examine the -ray source properties with an improved Fermi-LAT event reconstruction and instrument characterization, we analyzed ∼8.4 years of Pass 8 data from 2008 Aug 04 (15:43:36.0; all UTC times) to 2017 Jan 01 (00:00:00.0), selecting 0.1–100 GeV events within a circular region of interest (ROI) with radius 15◦ centered on the 2FGL position. The photons belonged to the SOURCE class (front and back converted events) as defined under the P8R2_SOURCE_V6 IRFs and were within a zenith angle of 90◦ of the LAT instrument to minimize contamination from Earth limb photons. We filtered the events to make sure the data were flagged as good and filtered out times corresponding to occurrences of bright, LAT-detected -ray bursts and Solar flares. Analysis of the -ray data was performed using Fermi Science Tools (STs) v10r00p05 and 47 binned likelihood analysis. The background model included all 3FGL catalog sources (Acero et al., 2015) and the diffuse -ray backgrounds (Acero et al., 2016) using the respective Galactic and isotropic templates files2, gll_iem_v06.fit and iso_P8R2SOURCE_V6_v06.txt. The quoted uncertainties are statistical only and are larger than expected LAT systematic uncertainties of ∼8% for fluxes and ∼0.1 on the spectral slopes (Ackermann et al., 2012). The normalization parame- ters of 3FGL sources with average significance ≥ 5 in 4 years and within 6◦ of the ROI center were left free to vary. Additionally, the normalization parameters of sources flagged as signifi- cantly variable and within 8◦ of the ROI center, and of the diffuse components were also free to vary. We first fit the entire ∼8.4-year data set centered at the 2FGL J0846.0+2820 position, modeled as a single power law (normalization and index free). The best-fit parameters of this initial fit were  = (5.2 ± 1.7) × 10−9 photons cm−2 s−1, Ɖ = 3.2 ± 0.4, and point-source TS = 12 (∼3 fitting only two spectral parameters). The flux is significantly less than the 2FGL value (see § 2.2.1), implying a variable source. To quantify the variability, we constructed a ∼ 0.5-year binned light curve with 17 bins, the first 16 of which had a 180 day duration and the last one 191 days. For the fit in each time bin, we allowed the same parameters to vary as in the fit of the entire data set but kept the photon index of 2FGL J0846.0+2820 fixed to the initially preliminary fit value. The -ray emission was found to be concentrated within the first two 180-day bins with TS = 6.2 and 29.5 respectively, coinciding with the first half of the time interval analyzed in the 2FGL catalog. As a next step, we checked for any nearby point sources new to the 8.4-year dataset and not in the 3FGL 4-year catalog. This was done by generating a large TS map of the field, spanning 5◦ × 5◦, with 0.◦25 per pixel. We found two candidate sources, and used gtfindsrc to localize them, then refit their spectral parameters with gtlike adopting single power-law models. Only one of the candidates appeared significant and was included in the model, with  = (4.1 ± 1.3) × 10−9 photons cm−2 s−1, Ɖ = 2.3 ± 0.4, and TS = 43.0. The  = 1.283 blazar, B2 0849+28 (Hewett & Wild, 2010), is just 0.(cid:48)7 offset from its best-fit LAT position, R.A. = 133.◦028, Decl. = +28.◦556, 2http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html 48 and 95% confidence error radius, 95 = 4.(cid:48)6, and is likely the source of the -ray emission. With the new point source in the background model, we fit the LAT position of the target source with only the first two 180-day bins of data, resulting in R.A = 131.◦83, Decl. = 28.◦17 (0.◦212 from the CSS source), 95 = 0.◦14 and average flux,  = (2.2 ± 0.7) × 10−8 photons cm−2 s−1 (∼2x the flux measured when fitting the entire dataset (§ 2.2.1)), Ɖ = 2.7 ± 0.2, and TS = 35. For the same 360-day dataset, we compared the best-fit likelihood for a exponentially cutoff power law (normalization, index, and cutoff energy free) and found no significant evidence for curvature in the spectrum (TScut = 0.0; see Abdo et al., 2013b). For completeness, we also repeated the analysis with the full 8.4-year dataset, and confirmed the initial finding of a low-significance source (TS = 6.4) averaged over the larger dataset. We utilized the fitted index = 2.7 that was based on the first 360-days of data, and the new best-fit LAT position to generate the final ∼ 0.5-year lightcurves as presented in Figure 2.3 where we only report points for bins in which 2FGL J0846.0+2820 was found with both TS ≥ 4 (∼2) and ≥ 4 predicted counts, else a 95% confidence upper limit is shown. We confirmed the initial result that the -ray emission is concentrated within the first two 180-day bins with TS = 5.3 and 35.4, respectively. We quantify the significance of variability as TSvar = 161.6 following the Fermi-LAT catalog analysis (see §3.6 in Nolan et al., 2012a, for the definition), which is well above the threshold of 32 needed to flag this source as variable at 99% confidence level. We produced a TS map of the region surrounding 2FGL J0846.0+2820 in a grid of 0.◦1 per pixel side using the Fermi ST gttsmap and data from the first 360 days of the Fermi mission, and a model which did not include a source corresponding to 2FGL J0846.0+2820. The resulting TS map is shown in Figure 2.4 with the localization contours, 2FGL ellipse, and optical position. The optical position is just outside both 95% contours, but within the new 99% localization. The new Pass 8 localization and 2FGL 95% contours overlap indicating these are the same source, however, the shift in the position is notable despite this being located at a fairly high Galactic latitude,  = 36.◦3, and our work in identifying all possible significant point sources in the field using the 8.4-year dataset. The shift could be due to several factors, including using different event reconstructions 49 Figure 2.3: -ray light curve of 2FGL J0846.0+2820. Fermi-LAT 0.1–100 GeV -ray light curve (Top) and corresponding TS values (Bottom) in ∼ 0.5 year bins. In the light curve, down arrows indicate 95% confidence upper limits on the fluxes, and the upper limit derived from the data after the first two bins is indicated with the horizontal dashed line. (P7V6 IRF) and time range compared to the 2FGL catalog analysis, diffuse models for the diffuse background models, and could also indicate a region with poorly modeled diffuse emission seen as extended regions of signal in the residual TS map. To investigate the regions of excess signal visible in Figure 2.4 that are not associated with 2FGL J0846.0+2820, we constructed another TS map with our target source in the model at the new position. This resulted in no excesses above TS = 13.5 and confidence contours which indicated possible extension, suggesting a positive fluctuation of unmodeled diffuse emission (there are no known features in the adopted diffuse emission model in the region of this source; see Fig. 5 in Acero et al. 2016). This is further supported by the fact 50 )(cid:173)1 s2 cm(cid:173)9Flux (0.1 to 100 GeV) (10102030MJD550005550056000565005700057500TS0102030 that no significant TS excesses were observed at these positions in the TS map covering the full time range. Figure 2.4: Map of -ray emission near the location of 2FGL J0846.0+2820. The LAT TS map centered on 2FGL J0846.0+2820 in J2000.0 coordinates (bottom bar indicates TS values). The green contours are, from inside out, 68, 95, 99% confidence on the localization, newly determined from the Pass 8 data analysis. The 2FGL 95% confidence ellipse (white ellipse) is shown as well as the optical CSS position (white cross). In the interest of better constraining when the flux of 2FGL J0846.0+2820 dropped, we analyzed the second 180 day time bin (2009 Jan 31 to 2009 July 31) in more detail. We constructed flux light curves with 30-day bins spanning ±60 days of our 180 day time bin. We then constructed a second 30-day flux light curve with the start times shifted by 15 days (Figure 2.5). Due to the rather faint nature of this source, we can not break the flux light curve into smaller intervals, but we estimate that the flux dropped sometime between 2009 May 31.7 (MJD 54982.7) and 2009 Jul 30.7 (MJD 55042.7). 51 03691215182124273027.5 28.0 28.5 29.0 29.5133.0 132.0 131.0 130.0Decl. (deg)RA (deg)CSS2FGL J0846.0+2820 (95%) An analysis fitting the entire LAT data set except the first 360 days, with the photon index fixed (Ɖ = 2.7), resulted in a non-detection with TS = 0 and a 95% flux upper limit of <2.2 × 10−9 pho- tons cm−2 s−1. This amounts to a drop in flux of ∼10x with respect to the second 180-day detection. Using the spectroscopic optical period and the second (and most significant) 180 day time bin, we searched the LAT data for evidence of modulation with the orbit. In order to do so, we calculated spectral weights, calculated with the best-fit model for this time range and the Fermi ST gtsrcprob, which have been shown to enhance sensitivity to periodic signals (e.g., Kerr, 2011). We then calculated the exposure at the position of 2FGL J0846.0+2820 in 30 second intervals over this time period in order to correct for exposure differences with orbital phase (as described in, for instance, Johnson et al., 2015). We used the Fermi ST gtpphase to assign orbital phase values to the LAT events and, correcting for exposure and using the spectral weights, tested for modulation using the H-test (de Jager et al., 1989; de Jager & Büsching, 2010) and the Z2 m test with two harmonics. Both tests showed no significant modulation, returning 0.1 and 0.3, respectively. 2FGL J0846.0+2820 is one of 234 unassociated sources not present in the 3FGL catalog, despite a detection at high significance in the earlier (2FGL) catalog. Many of the sources lost between 2FGL and 3FGL were at low Galactic latitudes where the Galactic diffuse emission is strongest, such that improvements in modelling this emission was expected to have the most influence in detecting sources. However, 2FGL J0846.0+2820 is out of the Galactic Plane. A number of 2FGL sources out of the Plane were also spurious, some apparently due to contamination from the Moon (Corbet et al., 2012). However, at the declination of 2FGL J0846.0+2820 (ecliptic coordinates: (, ) = (126.◦4, 9.◦8)), the effects of the Moon should be minimal, and given the refined event reconstruction and characterization presented here, there is no compelling evidence that it was a spurious source. Instead, our analysis shows that the most likely explanation for the absence of the source from the 3FGL catalog is the source variability over the lifetime of Fermi. 52 Figure 2.5: -ray light curve of 2FGL J0846.0+2820 near the time the flux dropped. The LAT light curve for the most significant 180-day interval (±60 days) in shorter 30-day time bins (black points). The data were shifted by 15 days (blue points) to find the last significant time bin. The upper limits are indicated when TS < 4 with fewer than four predicted photons. The dashed line marks the best fit average flux. 2.4.3 Optical Light Curve Models We show the phased light curves of the binary in Figure 2.6. There are two maxima and minima per orbit, consistent with ellipsoidal variations due to a tidally deformed secondary orbiting the presumed neutron star primary. Motivated by the H emission and change in the -ray flux of the binary, we start from the assumption that the secondary is filling its Roche lobe, and then explore models where the Roche lobe filling factor is free to vary. The baseline expectation for ellipsoidal variations are two equal maxima when the projected area of the tidally distorted star is largest, and two unequal minima due to varying effects of gravity darkening when the system is viewed along the axis connecting the primary and secondary. Mod- elling these ellipsoidal modulations can constrain the inclination of the system between the extremes where the maximum effect is observed if the system is edge-on ( = 90◦) and no modulations are expected for face-on inclinations ( = 0◦). 53 MJD5490054950550005505055100 )(cid:173)1 s(cid:173)2 cm(cid:173)8Flux (0.1 to 100 GeV) (1012345 Figure 2.6: PROMPT optical light curve of 2FGL J0846.0+2820. PROMPT , , and  light curves folded on the spectroscopic period. The best fit ELC models that included (right) and excluded (left) a steady accretion disk component are shown with black lines. Uncertainties for each measurement are not shown in the upper panels, but are displayed in the fit residuals (lower panels). Two orbital phases are shown for clarity. The broad filter used by CSS and its relatively large photometric uncertainties do not make these data ideal for modeling ellipsoidal variations in 2FGL J0846.0+2820. Instead, to model the light curve we use the more precise and well-sampled PROMPT photometry in , , and  bands. 54 0.00.51.01.5Phase−0.30−0.150.000.150.3017.116.916.716.5B0.00.51.01.5Phase−0.2−0.10.00.10.216.216.116.015.915.815.7V−0.15−0.080.000.080.1515.615.515.415.315.2RNoDiskDisk We model these light curves using the Eclipsing Light Curve (ELC; Orosz & Hauschildt, 2000a) code. We try fits that include contributions from a Roche-lobe filling secondary and a steady ac- cretion disk component, as well as fits that exclude the disk contribution and allow the filling factor to vary. Throughout, we assume the primary object is invisible. In addition, we set the scale of the system by fixing the value of the orbital period, semi-amplitude, binary mass function, mass ratio, and eccentricity to those determined via spectroscopy. These values imply an average orbital separation of  ∼ 26 (cid:12). We also require the observed value of  sin  be matched. From an analysis of our low-resolution spectra and the photometric colors, we fix the effective temperature of the secondary (2) to 5250 K. We note this value is not well-constrained as the ELC code fits normalized light curves in all bands and is thus insensitive to the effective temperature of the com- panion. The assumed effective temperature does affect the inferred distance, an issue we revisit in the following section. Given the long period of the system and its luminous secondary, irradiation of the companion by a central X-ray source is likely to be less important than in typical black widows or redbacks. While Chandra X-ray observations of 2FGL J0846.0+2820 are forthcoming, we have only upper limits on its X-ray flux from all-sky monitors (no clear detections from the Monitor of All Sky X-ray Image (MAXI, ∼1-20 keV; Matsuoka et al., 2009) or the Swift Burst Alert Telescope (BAT, ∼15-150 keV; Barthelmy et al., 2005)), which is a comparably weak constraint. For an X-ray luminosity of 1032-1033 erg s−1, the ratio of optical to X-ray flux at the surface of the star is >100, suggesting irradiation is not very important. On the other hand, recent detailed modeling of intra- binary shocks in redbacks have suggested the pulsar wind flow is channeled onto the stellar surface due to the magnetic field of the secondary (Romani & Sanchez, 2016; Wadiasingh et al., 2017), which may explain the unmodeled features seen in the light curve, such as a phase shift, in the con- text of irradiation. This issue should be explored in more depth when future X-ray data are available. We fit the folded light curves for , , and  bands simultaneously. For models that included an accretion disk component, we fixed the filling factor of the secondary to 1 while fitting the following parameters: the inclination , the inner disk temperature disk, the opening angle of the disk rim , the inner (in) and outer (out) radii of the disk, the power-law index of the disk temperature profile 55 Table 2.3. Best fit ELC parameters Free Parameter R.L. filling factora Inclination (◦) Inner disk temp. (K) Disk opening angle (◦) Inner disk radiusb Outer disk radiusb Disk powerlaw index () Phase shift (Ɗ ) 2(dof) No Disk 0.86 ± 0.03 27.1+1.1−1.0 · · · · · · · · · · · · · · · –0.022 ± 0.006 3987.0 (3178) Disk 1.0 30.7+3.0−1.7 3600+4600 −2500 15.0 0.05 ± 0.04 0.70 ± 0.13 –0.62 ± 0.2 –0.024 ± 0.005 3989.6 (3174) aThe Roche-lobe filling factor of the secondary was fixed to 1.0 for all models that included a disk. bExpressed as a fraction of the effective Roche-lobe radius of the primary. , and a small phase shift Ɗ . For models that did not include a disk, we fit the inclination, the filling factor of the secondary 2, and a phase shift. The best fit parameters with 1 confidence levels are listed in Table 2.3. Parameters without listed errors were not well constrained by the models (i.e. flat 2 distributions) so we refrain from quoting their exact uncertainties, and instead only report the values associated with the best fit. The first clear result from the fits is that a small phase shift (Ɗ ) is required: models with no shift can be clearly ruled out at the 4–5 level. Approximately the same Ɗ  was found for both the disk and no-disk case. Due to the long period, the small value of Ɗ  actually corresponds to a substantial offset of ∼ 4.3–4.7 hr, depending on the exact model. Similar phase offsets have been observed in some black widow and redback systems, likely due to asymmetries in the heating of the secondary or in the disk (e.g., Romani & Sanchez, 2016; Li et al., 2016). One possibility is that these asymmetries may be caused by hot spots in a disk or star spots on the companion. Although the ELC code allows us to fit for a variety of disk or star spots, 56 we find that adding such features do not significantly change the fit quality, probably due to the relatively face-on nature of the system and our typical measurement uncertainties. Nonetheless, starspots may be important for understanding the long-term variability in the system (see § 2.4.5). For models that include a disk component and a Roche-lobe filling secondary, we find our best fits to the data occur when the accretion disk contributes a small fraction (<3%) of the total light from the system. This relatively low veiling is not necessarily unexpected, given the weak H emission and the dominance of the photometric light curve by ellipsoidal variations. For fits with a disk we find an inclination of 30.◦7+3.0−1.7. For models without a disk,the best fit 2 is slightly lower than our best “with disk” model, but given the large number of data points the difference is negligible. Our best fit occurs when the companion is slightly underfilling its Roche lobe, with a filling factor 2 = 0.86. This is very similar to the results of McConnell et al. (2015) for the transitional MSP J1023+0023, who found 2 = 0.83+0.03−0.02 while the system was in the radio MSP state. For 2FGL J0846.0+2820, the inclination implied by the no-disk fit is 27.◦1+1.1−1.0. This value is slightly lower than inferred for the fit including a disk, consistent with the general result that disk veiling leads to an underestimate of the inclination from ellipsoidal variations (e.g., Kreidberg et al., 2012; McConnell et al., 2015; Wu et al., 2016). 2.4.3.1 Distance We estimate the distance to the system by comparing the luminosity inferred from the best-fit secondary radius and assumed temperature to the observed magnitude of the system. We determined bolometric corrections using 10 Gyr isochrones assuming [Fe/H] = –1 (Marigo et al., 2008). The best fit model assuming the star fills its Roche lobe has an effective radius of 7.1 (cid:12). Using an assumed effective temperature of 5250K, the secondary has a predicted  = 1.1 mag. Compared to the observed 0 = 15.63 mag (corrected for extinction using ( − ) = 0.10), the distance is ∼ 8.1 kpc. For an effective temperature of 5000 K, the distance would instead be ∼ 7.0 kpc. 57 We performed similar calculations for the case of an underfilled Roche lobe, finding very similar results: ∼ 8.3 and 7.2 kpc for temperatures of 5250 K and 5000 K, respectively. We emphasize that these distance calculations are uncertain and dominated by systematic effects. Even given the large uncertainty in the distance, the location of the system in the Galactic anti-center ( = 197◦) and its substantial distance above the plane (∼ 4 kpc for  = 36◦) makes it unlikely that the binary was formed in the thin disk. 2.4.4 Component Masses Using the posterior samples for the orbital period, semi-amplitude, and mass ratio obtained from our spectroscopic analysis (§ 2.4.1), along with our estimates of the inclination from fitting the photometry, we can infer the primary and secondary masses of the system. For the fits that excluded a disk, we find 1 = 2.81 ± 0.36 (cid:12) and 2 = 1.12 ± 0.21 (cid:12). The primary mass implied by these models is larger than any known neutron star and approaches the maximum theoretical mass (Chamel et al., 2013). For fits that included a disk, our models imply masses 1 = 1.96 ± 0.41 (cid:12) and 2 = 0.77 ± 0.20 (cid:12). In contrast to the fits with no disk, the primary mass is more typical of those inferred for neutron stars in binaries associated with Fermi sources (e.g., Ransom et al., 2011; Romani et al., 2015b; Strader et al., 2015, 2016), and we think these masses are more likely to be accurate. The uncertainties in the component masses are large due to the relatively face-on inclination. Furthermore, the potential for indirect, asymmetric heating on the face of the secondary may be adding an additional source of systematic uncertainty to our mass estimates. A few previous publications that describe the light curve modeling of similar black widow and redback systems have inferred unusually large neutron star masses (>2(cid:12)), which were later found to be unreliable because of inadequacies in the assumed model (e.g. not accounting for indirect heating via pulsar spin-down power reprocessed in an intra-binary shock, see Romani et al., 2015a,b, and references 58 therein). However, given the clear differences between this system and typical MSP spiders, includ- ing the luminous secondary, the wide orbit, and the face-on orientation of the binary, unmodeled irradiation is expected to make a less important contribution to the systematic uncertainties in our mass estimates than for most MSP binaries with hydrogen-rich companions (see discussion in § 2.4.3). Nonetheless, other interpretations of the system are unlikely: neither stellar-mass black holes nor white dwarfs in quiescence have been observed to emit GeV gamma rays. The inclination would need to be considerably lower ( = 21◦) to be consistent with even a low-mass 5 (cid:12) black hole, and would need to be somewhat higher (35–40◦) to accommodate a white dwarf. Furthermore, an accreting white dwarf with a giant secondary in an 8.1 d orbit would be extraordinarily unusual; most such “symbiotic" systems have periods of hundreds of days and accrete from a wind (Belczyński et al., 2000). Future observations can help improve constraints on the presence of a disk and hence on the inclination. 2.4.5 Long-Term Optical Brightness As discussed above, there is compelling evidence that the system has brightened at a rate of 0.011 ± 0.001 mag yr−1 over the last decade, with only a few small interruptions in this monotonic trend (Figure 2.7). There are a number of possible explanations for this trend: the secondary could be slowly increasing in radius; the mean effective temperature could be increasing (either globally or due to a change in the properties of starspots); or an accretion disk could be getting brighter. The last of these was proposed as an explanation for similar secular brightness trends in the qui- escent black hole binary systems Nova Muscae 1991 (Wu et al., 2016) and A0620-00 (Bailyn, 2017). The trend for 2FGL J0846.0+2820 is not confined to the CSS data: we find that our recent PROMPT data fits the trend, after correcting for a small zeropoint offset of -0.30 mag in the listed Catalina magnitudes using non-variable comparison stars in common between the PROMPT and CSS data. 59 Figure 2.7: Long-term optical light curve of 2FGL J0846.0+2820. Long-term optical light curve showing a monotonic increase in the system brightness over the past decade. The Catalina data are shown in red and our recent PROMPT -band data are in green (far-right on figure). The dashed line and rate of brightening are the result of a linear fit to the Catalina data only. A small zeropoint offset has been applied to the PROMPT data as described in § 2.4.5 In addition to the overall trend, the data taken between 2009 November and 2010 June shows a larger scatter than observations before or after this. Above (§2.4.2), we showed the -ray emission disappeared around the start of 2009 July, and these high-scatter observations are the first CSS data taken after this change in the -ray flux. Given the long time baseline for comparison, the photometric behavior around this epoch is unusual, and we consider it unlikely that these two events are coincidental: something happened in mid-2009 that affected both the -ray flux and optical brightness of the binary. To show the CSS photometry in more detail, we break the data into eight chunks corresponding roughly to the separate epochs visible in Figure 2.7. When we phase the data on the spectrocopic period, strange behavior is apparent (Figure 2.8). The earliest and latest data show variations more consistent with the ellipsoidal variations found in the PROMPT data (black lines, adjusted for the linear brightening trend), though the phase coverage is not ideal. However, the epochs that bracket the decrease in gamma rays do not show the same variations, with large changes in either the phase of the variations, or no clear variations at all. For instance, the worst fit to the best ELC model occurs for the data taken between 2009 November and 2010 June (2/dof = 82.3/36), while the data from 2005 October to 2006 May clearly matches the model much better (2/dof = 23.2/71). 60 20052006200720082009201020112012201320142015201605001000150020002500300035004000MJD-5300015.415.515.615.715.815.916.016.1CSSV-equivalent(mag)0.011±0.001mag/yr Figure 2.8: Long-term optical light curve of 2FGL J0846.0+2820 split into separate epochs. The long-term CSS data split into eight separate epochs (first 4 on left, last 4 on right) and phased on the spectroscopic period. The best fit ELC model (with a disk) from § 2.4.3 is shown in black, adjusted for the linear trend seen in Figure 2.7. The shape of the light curve does not appear consistent with typical ellipsoidal variations across all epochs. 61 15.915.815.715.6Oct’05-May’0615.915.815.715.6Nov’06-Jun’0715.915.815.715.6Oct’07-May’08CSSV-equivalent(mag)0.00.51.0Phase15.915.815.715.6Oct’08-May’09Nov’09-Jun’10Nov’10-May’11Nov’11-May’120.00.51.0PhaseOct’12-May’13BestfitmodeltoPROMPTdata This strange phenomenology as a function of epoch is similar to that observed by van Staden & Antoniadis (2016) for the optical light curves of the redback PSR J1723–2837, which they explain with distinct groups of starspots that vary in number and lifespan over time. 2.5 Discussion and Conclusions Continuing our program to follow up on unassociated Fermi sources using photometry and spectroscopy, we have shown that the previously unidentified -ray source 2FGL J0846.0+2820 is likely associated with a heavy neutron star of roughly 2 (cid:12) with a giant secondary companion (2 ∼ 0.8 (cid:12)) in an 8.133 d orbit. Unrelated variable stars will exist in the error ellipses of some unassociated Fermi sources. For instance, the space density of CRTS variable stars is about 2.7 per deg2 at the latitude of 2FGL J0846.0+2820. It is therefore not surprising that, in addition to the source that is the subject of this paper, there is another periodic optical variable (CSS J084632.3+282524) present inside the newly determined Pass 8 99% confidence ellipse. This other source is classified as a contact binary with an optical period of 0.36 days (Marsh et al., 2017). Based on the massive, invisible primary inferred for the original CSS source (CSS J084621.9+280839), and the rarity of compact binaries in the field, the evidence for an association between this source and the Fermi -ray emission is strong. An analogous example is the recent discovery of the likely redback MSP counterpart to the Fermi source 3FGL J0838.8–2829. Halpern et al. (2017) observed the 3FGL field in the optical and X-rays, motivated by the presence of a cataclysmic variable (CV) that showed variability in these bands. However, it was quickly determined that this CV was unlikely to be the -ray source. Instead, optical photometric and spectroscopic observations suggested a previously unknown X- ray/optical variable source, which showed properties consistent with redback MSPs in their pulsar state, was the correct counterpart to the -ray source (see also, Halpern et al., 2017). Forthcoming X-ray observations will help confirm the connection between the Fermi -ray source 2FGL J0846.0+2820 and the optical binary CSS J084621.9+280839. 62 Nonetheless, this work shows that such positional coincidence searches between optical vari- ables and unassociated Fermi sources can be productive and may provide a relatively “low-cost” method of identifying new compact binaries in the Galaxy. Overall, 2FGL J0846.0+2820 is very similar to the recently discovered -ray bright binary 1FGL J1417.7–4407 (Strader et al., 2015), which was independently found to host a radio MSP (Camilo et al., 2016). This system has a heavy neutron star primary and giant secondary in a long 5.4 d orbit. The high primary mass and long period of 2FGL J0846.0+2820 are also similar to the famous MSP J1614–2230, though this system has a white dwarf secondary and has ceased mass transfer (Demorest et al., 2010). Similarly to the conclusions reached by Strader et al. (2015) and Camilo et al. (2016) for 1FGL J1417.7–4407, 2FGL J0846.0+2820 has properties consistent with a system on the standard evo- lutionary track of low-mass X-ray binaries that started Case B mass transfer after leaving the main sequence and whose orbital periods will increase over time (e.g., Tauris & Savonije, 1999). These systems do not fit into the existing classes of black widow or redback neutron star binaries, and we have suggested that an apt name for this new subclass is “huntsman", after a large spider that does not engage in sexual cannibalism. One relevant difference between 1FGL J1417.7–4407 and 2FGL J0846.0+2820 is the non-zero eccentricity of the latter. Nearly all binary MSPs in the field have circular orbits due to tidal interactions that occur when the system undergoes mass transfer, spinning up the pulsar. It is possible that 2FGL J0846.0+2820 has only partially completed the recycling process and not yet circularized, or that a non-zero eccentricity has emerged due to the interaction of the binary with a circumbinary disk, as has been theorized for some MSP systems (Antoniadis, 2014). It is now clear that both rotation-powered millisecond pulsars and accreting neutron stars can be associated with -ray emission. Such emission is generally weak for millisecond pulsars but is apparently ubiquitous in these systems, while it has been observed in only a small number of low-mass X-ray binaries. These latter systems all belong to the class of transitional millisecond 63 pulsars or have other phenomenological similarities with this class. For two of the three systems that have shown actual transitions, the -ray flux is observed to be higher in the accreting state. Thus, it would be natural to associate the decrease in -ray flux for 2FGL J0846.0+2820 with a transition from a disk to a pulsar state in mid-2009. As we discussed in § 2.4.5, there are some unusual features in the optical light curve around that possible transition time. The issue with a direct analogy to the other transitional millisecond pulsars is that the higher variations in the optical flux were transient, with no long-term evidence of a state change. For example, the 2013 state change in PSR J1023+0038 was associated with an optical brightening of nearly 1 mag (Bogdanov et al., 2015). Hence, these optical data are more consistent with a “glitch" in the binary than a full-fledged state change. However, the relatively nearby main-sequence companion of PSR J1023+0038 is much less luminous than the more distant giant secondary in 2FGL J0846.0+2820. If the companion to PSR J1023+0038 were replaced with the giant in 2FGL J0846.0+2820, the secondary would swamp the disk, resulting in an inferred disk fraction that falls well below 10%. Therefore, similar state changes in 2FGL J0846.0+2820 might not be associated with such obvious changes in the optical brightness. Given that the mechanism for these transitions is not understood, and the notable physical differences between these hunstmen and the redback transitional millisecond pulsars, a difference in phenomenology of the optical and -ray emission may not be unexpected. Some evidence for this can be found in 1FGL J1417.7–4407. In this system there is strong, double-peaked H observed at nearly all epochs from early 2013 to the present (with the latter statement from continuing spectroscopic monitoring with SOAR). Strader et al. (2015) took this as evidence for an accretion disk. This system also has a hard X-ray spectrum and a high ratio of -ray to X-ray flux, consistent with transitional MSPs in their disk states. However, this simple interpretation was challenged by the discovery of an MSP associated with this source by Camilo et al. (2016). The overlap between the observation epochs of these two studies show that the source was active as an MSP while it also appeared to have an accretion disk. It is possible that a disk is 64 present but that the accreted material is not reaching the surface of the neutron star, perhaps due to ejection in a propeller (Strader et al., 2015). Alternatively, the double-peaked H could be due to complex line profiles associated with a wind (Camilo et al., 2016) or an intra-binary shock that forms due to the interaction of the pulsar wind outflows with the disk and/or material from the companion (e.g. Bogdanov et al., 2011; Rivera Sandoval et al., 2017). At present the evidence for an accretion disk in 2FGL J0846.0+2820 is mixed. The gradual optical brightening and H emission suggests a disk could be present. However, the current ab- sence of gamma rays, the dominance of the photometric light curve by ellipsoidal variations, and our best-fit ELC models favor a low-mass disk at the most. Given that a millisecond pulsar was detected in 1FGL J1417.7–4407 (Camilo et al., 2016) despite the presence of strong double-peaked H suggesting a more substantial accretion disk (Strader et al., 2015), a search for a radio pulsar in 2FGL J0846.0+2820 is well-motivated. The apparent variability in the Fermi data along with the peculiar long-term optical light curve suggest we are far from fully understanding 2FGL J0846.0+2820. However, this work provides convincing evidence for the presence of a heavy neutron star primary with a giant secondary in a relatively wide orbit. These long-period -ray bright binaries with massive neutron stars and giant companions likely represent the late phases of typical MSP binary formation in the Galactic field, phases which until recently had been unobserved. Further characterizations of similar systems would also provide new insight on the relationship between MSPs and low-mass X-ray binaries. 65 CHAPTER 3 A MULTI-WAVELENGTH VIEW OF THE NEUTRON STAR BINARY 1FGL J1417.7–4402: A PROGENITOR TO CANONICAL MILLISECOND PULSARS 1The Fermi -ray source 1FGL J1417.7–4407 (J1417) is a compact X-ray binary with a neutron star primary and a red giant companion in a ∼5.4 day orbit. This initial conclusion, based on optical and X-ray data, was confirmed when a 2.66 ms radio pulsar was found at the same location (and with the same orbital properties) as the optical/X-ray source. However, these initial studies found conflicting evidence about the accretion state and other properties of the binary. We present new optical, radio, and X-ray observations of J1417 that allow us to better understand this unusual system. We show that one of the main pieces of evidence previously put forward for an accretion disk—the complex morphology of the persistent H emission line—can be better explained by the presence of a strong, magnetically driven stellar wind from the secondary and its interaction with the pulsar wind. The radio spectral index derived from VLA/ATCA observations is broadly consistent with that expected from a millisecond pulsar, further disfavoring an accretion disk scenario. X-ray observations show evidence for a double-peaked orbital light curve, similar to that observed in some redback millisecond pulsar binaries and likely due to an intrabinary shock. Refined optical light curve fitting gives a distance of 3.1±0.6 kpc, confirmed by a Gaia DR2 parallax measurement. At this distance the X-ray luminosity of J1417 is (1.0+0.4−0.3) ×1033 erg s−1, which is more luminous than all known redback systems in the rotational-powered pulsar state, perhaps due to the wind from the giant companion. The unusual phenomenology of this system and its differing evolutionary path from redback millisecond pulsar binaries points to a new eclipsing pulsar “spider" subclass that is a possible progenitor of normal field millisecond pulsar binaries. 1This chapter is taken mostly word-for-word from Swihart et al. (2018), as published in the Astrophysical Journal 66 3.1 Introduction Millisecond pulsars (MSPs) form when matter and angular momentum are accreted onto a neutron star from a non-degenerate companion, recycling them to very rapid spin periods. These systems are generally observable as low-mass X-ray binaries (LMXBs). The typical MSP recycling process involves accretion from a giant donor star overfilling its Roche-lobe, and ends as the orbital period grows and accretion onto the neutron star eventually stops. The resulting system has a rotationally-powered pulsar primary and a low-mass (0.2–0.3 (cid:12)) white dwarf companion, with orbital periods ranging from days to weeks (Tauris & van den Heuvel, 2006). These types of systems constitute the bulk of the known MSPs in the Galactic field. The most important advance in expanding the known population of MSP binaries has been multi-wavelength (X-ray, optical, and radio) follow-up observations of Fermi-LAT GeV -ray sources, which have led to the discovery of many MSPs in binaries with non-white dwarf com- panions. With improved statistics, these systems have been categorized based on their companion mass: black widows have very light ( (cid:46) 0.10(cid:12)), highly ablated companions, while redbacks have non-degenerate, main sequence-like companions ( (cid:38) 0.2(cid:12)) that typically fill a substan- tial fraction of their Roche lobe (Roberts, 2011). These systems frequently show lengthy, irregular radio eclipses, making them challenging to discover in normal radio pulsar surveys. These new binaries show interesting phenomenology, but have also drawn extensive interest as the long sought-after links between MSPs and their LMXB progenitors. In particular, three redbacks have been observed to undergo transitions between an accretion-powered “disk state” and a rotationally-powered “pulsar state” on timescales of days to months (Archibald et al., 2009; Papitto et al., 2013; Bassa et al., 2014; Roy et al., 2015; Bogdanov et al., 2015; Johnson et al., 2015). These systems, known as transitional millisecond pulsars (tMSPs), give one view of the end of the MSP recycling process, but also provide essential insights into the physics of low-level accretion onto magnetized compact objects and the interactions between MSPs and their stellar and gaseous environment. 67 One challenge in extending these insights to understanding the formation of all MSPs is that all the known tMSPs have short orbital periods ((cid:46) 0.5 days), and hence will likely not end their lives as the typical field MSP binaries, which have degenerate companions and longer orbital periods. The progenitors to such systems are thought to be red giant–neutron star binaries with > 1 day periods (e.g., Tauris & Savonije, 1999). One of the few known systems fitting the latter description is a recent discovery. 1FGL J1417.7– 4407 / 3FGL J1417.5–4402 (hereafter J1417) is a Galactic compact binary whose stellar counterpart was first discovered by Strader et al. (2015) in an optical survey of unassociated Fermi-LAT -ray sources. A variable X-ray source near the center of the Fermi error ellipse matches a bright optical source. Optical photometry and spectroscopy revealed that the optical source is a red giant, and modeling its ellipsoidal variations, radial velocities, and rotational velocity showed that the system is an evolved ∼0.35 (cid:12) late-G/early-K giant in a 5.37 day orbit with a ∼2 (cid:12), suspected neutron star primary. This compact object interpretation was confirmed when Camilo et al. (2016) found a 2.66 ms radio pulsar at the same location and with the same orbital period and phase as the optical binary discovered by Strader et al. (2015). The long orbital period, giant secondary, and short spin period suggest this system is likely in the late stages of the standard MSP recycling process that will end with a white dwarf companion in a (cid:38) 8 day orbit (Podsiadlowski et al., 2002), making it the first typical MSP binary progenitor identified in the Galaxy. The associated -ray emission, non-degenerate secondary, and MSP primary make J1417 sim- ilar to other redbacks, while its wider orbit, giant companion, and inferred evolutionary track are unique, leading us to suggest it as the first member of the “huntsman" subclass of MSP binaries. Furthermore, despite the fact that J1417 hosts an MSP, this system shows phenomenology dissim- ilar to redbacks/tMSPs in their pulsar states: in particular, the optical spectra show double-peaked H emission lines at nearly all epochs (Strader et al., 2015), which is a classic signature of an accretion disk. This emission was observed even during times when the system was observed by 68 Camilo et al. (2016) as a radio pulsar, precluding the possibility that a pulsar to disk state change had occurred. Another peculiarity is that, at the distance inferred by Strader et al. (2015, 4.4 kpc), the implied X-ray luminosity was ∼ 1.4 × 1033 erg s−1, typical for a redback/tMSP in the disk state but higher than usually observed for a system in the pulsar state. Camilo et al. (2016) suggest a nearer distance, which would produce a more typical X-ray luminosity of ∼ 1032 erg s−1, but which still leaves the spectroscopic evidence for an accretion disk unexplained. To better constrain the properties of J1417 and its connection to redbacks/tMSPs, here we report on new multi-wavelength observations of the system. We describe our newly acquired X-ray, radio, and optical/near-IR observations of the source in §3.2, §3.3, and §3.4, respectively. In §3.5 we perform detailed modeling of the ellipsoidal light curve, while we attempt to explain the source of the complex H structures in §3.6. Concluding remarks are presented in §3.7. 3.2 X-ray Observations 3.2.1 Chandra ACIS Strader et al. (2015) analyzed a 2 ks Chandra/ACIS-I exposure of the 1FGL field taken in 2011 Feb (ObsID 12843; PI Ricci). The neutron star binary discussed in this work (J1417) was the only significant X-ray source in the field. They found that the X-ray spectrum was well-fit by an absorbed single power-law with photon index Ɖ = 1.32 ± 0.40, but that this could not be distinguished from an absorbed blackbody spectrum (kT = 0.78+0.16−0.12 keV) due to the small number of counts. No significant limits could be placed on the source variability. These data sampled orbital phases  = 0.981 − 0.986. To achieve improved constraints on the short-term X-ray variability, as well as an improved spectrum, we obtained a single, uninterrupted ∼50 ks Chandra exposure of J1417 (ObsID 17786; PI Strader) on 2016 Jun 6 using the ACIS-S3 CCD configured in FAINT mode. To avoid pileup, 69 the CCD was used in a custom timed exposure subarray mode with 256 pixel rows, starting from CCD row 385. The data were reprocessed, reduced, and analyzed with CIAO 4.9 using the CalDB 4.7.4 calibration files (Fruscione et al., 2006). We found no evidence for background flares (these are thought to come from low-energy solar or geomagnetic protons; Grant et al., 2002). The data gave a net exposure of 49771.0 s, covering ∼11% of the orbit from  = 0.86 − 0.97. We note that here and throughout this paper, we adopt the ephemeris from Camilo et al. (2016), where  = 0 is the ascending node of the pulsar. We extracted spectra of the source with specextract using a circular region of radius 2.5”, with the background determined from a nearby circular, source-free region with radius 25”, yielding a net count rate of 60.3 ± 1.1 cts ks-1 (2886 photons). To analyze the spectra using Xspec (Arnaud, 1996), we binned the data to have at least 30 counts per bin, and assumed Wilms et al. (2000) abundances for the photoelectric absorption models throughout. 3.2.1.1 X-ray Variability To analyze the light curve, we applied a barycenter correction to the observation times using the CIAO tool axbary and grouped the data into 0.5 and 2 ks time bins. The corresponding 0.5–7 keV light curves are shown in Figure 3.1. Visual inspection of the light curves reveal the observed 0.5–7 keV count rates vary by a factor of two over this relatively small range of orbital phases. To quantify the significance of variability, we fit the 2 ks bins to a constant count rate model set to the weighted average of the dataset (Figure 3.1, dashed line). From a 2 fit, we find a probability of ∼0.5% that these data arise from a constant flux distribution. Data over a longer time period would be necessary to verify this variability and to determine whether the X-ray flux is modulated on the orbital period, as observed for many redbacks in the pulsar state, likely due to an intrabinary shock (Bogdanov et al., 2015). 70 Figure 3.1: Chandra X-ray light curve of 1FGL J1417.7–4407. The Chandra ACIS light curve between 0.5–7 keV. The ∼50 ks exposure covers ∼11% of the binary orbit. Photon arrival times have been barycenter corrected and are grouped into 500s (black) and 2ks (blue) bins. The black dashed line marks the weighted average of the 2ks dataset. 3.2.1.2 Spectral Fitting We first attempted to fit our newly acquired Chandra data with a purely thermal model, however, both an absorbed single blackbody (specified in Xspec as tbabs*bbodyrad) as well as a neutron star hydrogen atmosphere model (tbabs*nsatmos) yielded unsatisfactory fits (2/dof = 271.9/76 and 814.1/76, respectively). Next, we fit the spectrum with an absorbed single power-law model (tbabs*pegpwrlw), finding an excellent fit (2/dof = 63.7/76, for a null hypothesis probability of 0.84 in Xspec) with photon index Ɖ=1.41 ± 0.09 (90% confidence). This value is consistent with the 2011 data, but with higher precision. In fitting for the absorption component, we found a neutral hydrogen column density () of (1.5 ± 0.6) × 1021 cm-2, which is about twice the Galactic column density in this direction 71 0.880.900.920.940.96BinaryPhase0102030405060Timesince2016Jun0622:38:51.816UTC(ks)0.020.040.060.080.100.12ChandraACISCountRate(cts/s)0.5ks2ks inferred from H1 data (6.54×1020 cm-2, Kalberla et al., 2005), and is consistent with the findings of Camilo et al. (2016), who inferred  = 2.2+1.7−1.3 × 1021 cm-2 from an analysis of 2015 Swift data (see §3.2.2). Holding the absorption component fixed to the Galactic value yields a slightly worse fit (2/dof = 69.5/77). The extra absorbing material could be Galactic or could be intrinsic to the system, a possibility we discuss below. Given that the absolute value of  is low and not all that well-constrained, for our spectral analysis we only quote fits in which  was free to vary. The unabsorbed 1–10 keV flux of this model is 8.75× 10−13 erg cm-2 s-1, consistent with the 2011 Chandra observations, which had large ((cid:38)30%) uncertainties (Strader et al., 2015). While a simple power-law provides an acceptable fit to the data, we test for the presence of a thermal component in the spectra by adding a blackbody component (tbabs*(pegpwrlw + bbodyrad)). This model provides a slight improvement in the fit over the pure power-law model, with 2/dof = 57.1/74, for a null hypothesis probability of 0.93. The power-law component of this model has a photon index Ɖ = 1.08+0.25−0.15 while the blackbody has  = 0.60+0.15−0.09. An  test gives a probability of 1.7% that the improvement in the fit that includes a thermal component occurs by chance if the pure power-law model was correct. Although this is not strong evidence for a thermal component in the spectrum, it is consistent with previous X-ray studies of tMSPs in their pulsar states (e.g., Archibald et al., 2010; Bogdanov et al., 2011), which also show marginal evidence for a thermal component coming from a small “hot spot” on the surface of the neutron star, likely from the heated magnetic polar caps. We also note the spectral index in this model is harder than those typically observed in tMSPs during their accretion states (Ɖ ∼ 1.7; e.g., de Martino et al., 2010; Bogdanov et al., 2015), and is broadly consistent with synchrotron emission from an intrabinary shock. Both of these models are shown in Figure 3.2 and summarized in Table 3.1. All errors are quoted at 90% confidence. Assuming our new optical distance of 3.1 kpc (see §3.5.1), the unabsorbed flux of the pure power-law model corresponds to a 1–10 keV X-ray luminosity of (1.0+0.06−0.04)× 1033 (d/3.1 kpc)2 72 erg s−1. At the geometric distance derived from the Gaia parallax (3.5 kpc, §3.5.1), the X-ray luminosity is x∼1.3× 1033 erg s−1. At either distance, J1417 is brighter in X-rays than all known redbacks in the pulsar state, likely due to a strong intrabinary shock. We return to this in §3.6 and §3.7. Figure 3.2: Chandra X-ray spectra of 1FGL J1417.7–4407 with best fit models in red. In the upper figure, the data is fit with a pure absorbed power-law model, while the lower figure shows an absorbed power-law (dashed) plus the addition of a blackbody component (dotted). Fit residuals are shown in the bottom panels. 73 10−410−30.01normalized counts s−1 keV−111025−101(data−model)/errorEnergy (keV)Absorbed powerlaw10−410−30.01normalized counts s−1 keV−111025−1012(data−model)/errorEnergy (keV)Absorbed powerlaw + blackbody Table 3.1. Summary of Chandra Spectral fits for J1417 Parameters  (1021 cm-2) Photon Index (Ɖ)  (10−13 erg cm-2 s-1)a  (keV) Blackbody Fraction Blackbody Emission Radius (km)b 2 (dof) power-law power-law + Thermal Component 1.5 ± 0.6 1.41 ± 0.09 8.75 ± 0.44 · · · · · · · · · 0.84 (76) <1.0 1.08+0.25−0.15 8.66 ± 0.50 0.60+0.15−0.09 0.13+0.05−0.09 0.09 ± 0.06 0.77 (74) aUnabsorbed X-ray flux between 1–10 keV bAssuming a distance of 3.1 kpc (§3.5.1) Table 3.2. Swift XRT Observations of J1417 Date (UT) 2016 May 07 2016 Jun 05 2016 Jul 05 2016 Aug 02 2016 Aug 31 2016 Sep 29 Exposure Time (ks) 1.14 0.94 1.05 0.88 0.95 0.98 XRT Count Rate (0.3–10 keV) (10−2 s−1) 2.0 ± 0.5 1.7+0.6−0.5 1.9 ± 0.5 2.7 ± 0.6 1.5+0.5−0.4 1.6+0.5−0.4 Phase 0.152 0.564 0.153 0.297 0.751 0.218 Note. — Uncertainties are listed at 90% confidence. Orbital phase convention is that  = 0 is the ascending node of the pulsar. 74 3.2.2 Swift XRT Camilo et al. (2016) analyzed nine Swift X-ray Telescope (XRT, Burrows et al., 2005) observations of J1417 totalling 12.5 ks between 2015 Mar 13 and Jun 18. In addition to finding slight evidence for orbital variability in the light curve, the X-ray flux and spectral index they find are consistent with our new, deep Chandra data. We also observed this source with Swift/XRT in photon counting mode on six epochs between 2016 May 7 and Sept 29 for a total of 5.9 ks. Similar to the analysis presented by Camilo et al. (2016), we extracted source counts in the 0.3–10 keV energy range from a circular region with radius 47”, with the background determined from a nearby circular, source-free region with radius 200”. We summarize these observations in Table 3.2. Overall, we see that the mean 0.3–10 keV count rate of J1417 has remained relatively constant since the 2015 data were taken (Figure 3.3). These data also allow us to investigate the possibility of orbital variability. This is challenging for J1417, given its long period and the visibility constraints imposed by Swift’s low Earth orbit. We chose to analyze the folded X-ray light curve per individual “snapshot” so that each data point corresponds to a relatively narrow range of binary phases. The accuracy on the binary period derived from radio pulsar timing is 0.00003 d (Camilo et al., 2016), such that the absolute phase uncertainty of the XRT observations is only ∼0.0005 (about 4 minutes) after building a phase-coherent solution forward to the most recent XRT epochs. This is negligible compared to the exposure times of each XRT snapshot. We show the folded XRT light curve in Figure 3.4. Before we discuss the phased light curve, it is worth a brief review of theoretical expectations for the orbital variations in the X-ray flux associated with an intrabinary shock. In the relatively simple model of Romani & Sanchez (2016), the shock is assumed to occur between the (irradiation-driven) wind from the companion and the pulsar wind (see also the more sophisticated model of Wadiasingh et al. (2017)). The controlling physical parameter is the ratio between the companion and pulsar wind momentum fluxes. If the pulsar wind dominates, the X-ray light curve has a double-peaked 75 Figure 3.3: Long-term Swift X-ray light curve of 1FGL J1417.7–4407. The 0.3–10 keV Swift XRT count rates, extracted from 15 observations between 2015 Mar and 2016 Sept. maximum around  = 0.25 (when the companion is located between the MSP and the observer) and a minimum at  = 0.75. If the companion dominates, the minimum and maximum are shifted by a phase of 0.5. The velocity of the companion wind with respect to the orbital velocity affects the symmetry of the X-ray light curves, and more face-on orbits mute the orbital variations. Overall, the morphology of the J1417 X-ray light curve shows marginal evidence for features expected to arise due to an intrabinary shock, such as the “peaks” around  ∼ 0.15 and 0.35. If real, the double-peaked emission in this model is due to Doppler boosting along the line of sight, with a central dip due to the companion eclipsing part of the shock. However, we emphasize that the low count rate and incomplete phase coverage precludes a detailed comparison of the X-ray light curve with possible models. For example, the simple intrabinary shock model does not consider the possibility that the magnetic field of the secondary may be energetically relevant for shaping the interaction between the winds (e.g., Roberts et al., 2014; Sanchez & Romani, 2017). 76 570505710057150572005725057300Time(MJD)0.000.010.020.030.040.05SwiftXRTCountRate(cts/s)Camiloetal.(2016)ThisWork57500575505760057650 Figure 3.4: Swift X-ray light curve of 1FGL J1417.7–4407 folded on the binary period. Snapshots from the 15 Swift XRT observations in Figure 3.3 folded on the orbital period of the binary. The times spanned by each exposure are smaller than the symbol sizes. The dashed line shows the weighted average of all snapshots. Two orbital cycles are shown for clarity. 3.3 Radio Continuum Observations In all the known tMSPs that have been observed in their disk states, bright, flat-spectrum radio emission is present. This has been explained by partially self-absorbed synchrotron radiation com- ing from a compact jet-like outflow (e.g., Hill et al., 2011; Deller et al., 2015; Papitto & Torres, 2015), though the short timescale anti-correlation between X-ray and radio emission in the tMSP J1023+0038 in its disk state clouds a simple jet interpretation for the radio emission (Bogdanov et al., 2018). In a number of black hole LMXBs, the connection between outflow processes and the inflow 77 0.00.51.01.52.0BinaryPhase0.000.010.020.030.040.05SwiftXRTCountRate(cts/s)Camiloetal.(2016)ThisWork from an accretion disk have resulted in an empirical radio-X-ray luminosity (  −  ) correlation, which ties the X-ray luminosity (a probe of the inner regions of the accretion flow) to the radio luminosity, which is set by the power of the jet (Corbel et al., 2000; Gallo et al., 2003; Corbel et al., 2013; Gallo et al., 2014; Plotkin et al., 2017a,b). A similar disk-jet coupling has also been observed (to a lesser extent) in a number of “normal” accreting neutron star systems, albeit with radio emission that is at least a factor of 10 dimmer than observed for black holes at the same X-ray luminosity (Migliari et al., 2003; Migliari & Fender, 2006; Tudose et al., 2009; Tetarenko et al., 2016). However, this similarity is obscured by the fact that different neutron star subsamples have shown different correlations in the   −   plane, suggesting there may not be a “universal”  −  correlation for neutron stars as there is for black holes (Tudor et al., 2017; Gallo et al., 2018). Due to the less luminous jets in “ordinary” neutron star LMXBs compared to black holes (by a factor ∼22, Gallo et al., 2018), it is difficult to detect these systems in the radio at the relatively low X-ray luminosities we see in J1417. The vast majority of neutron star jets that have been explored to date have X-ray luminosities (cid:38) 1036 erg s-1. In fact, only upper limits on the 5-GHz radio luminosity of a few neutron star LMXBs have been reported in the literature at X-ray luminosities (cid:46) 1034 erg s-1 (e.g., Tudor et al., 2017; Gallo et al., 2018). By contrast, while in their disk states, two of the three tMSPs have been firmly detected at 5-GHz with X-ray luminosities between 1032– 1034 erg s-1 (Hill et al., 2011; Deller et al., 2015) (the third tMSP is at much higher  , Papitto et al., 2013). Despite the limited statistics, tMSPs appear to be more radio bright than typical neutron star LMXBs at the same X-ray luminosity, possibly due to the ejection of disk matter by a propeller mechanism (Deller et al., 2015), though the existing radio upper limits for such systems at   ∼ 1032–1034 erg s-1 are not particularly constraining. Furthermore, the anti-correlation between the radio/X-ray luminosity in the high and low X-ray modes of the accreting state of the tMSP J1023+0038 induces movement in the   −   plane that is perpendicular to the slope of any mooted correlation (Bogdanov et al., 2018), suggesting a crucial difference between tMSPs and “normal” accreting neutron star binaries. 78 Table 3.3. VLA radio continuum data for J1417 Date 5.1 GHz flux density (Jy) 7.1 GHz flux density (Jy) spec. index phase 2015 Jun 15 2015 Jun 16 2015 Jun 17 2015 Jul 15 Combined < 31.0 56.5 ± 9.8 < 27.7 < 29.9 22.4 ± 4.8 < 30.3 40.9 ± 9.5 < 26.7 < 31.4 13.5 ± 5.1 −1.0 ± 0.9 · · · · · · · · · −1.5 ± 1.3 0.269 0.454 0.640 0.843 · · · Note. — Radio upper limits are 3 values. Orbital phase convention is that  = 0 is the ascending node of the pulsar. Table 3.4. ATCA data for J1417 Date configuration central freq. (GHz) 2015 Oct 24 6A 2016 Jun 06 1.5B 1.43 1.89 2.59 1.41 1.86 2.37 2.83 flux density (Jy) 294 ± 42 88 ± 27 145 ± 21 258 ± 90 135 ± 71 86 ± 47 131 ± 51 spec. index −1.1 ± 0.3 −1.2 ± 0.8 3.3.1 VLA Motivated by the evidence that J1417 was in a disk state based on the persistent double-peaked H emission seen by Strader et al. (2015), and prior to the publication of a pulsar detection, we obtained VLA data to search for the presence of radio continuum jets. We observed J1417 with the Karl G. Jansky Very Large Array (VLA) under program VLA/15A-491 (PI Strader). Observations were obtained over four epochs during 2015 June/July, while the VLA was moving from BnA configuration into A configuration (which has a maximum baseline of max = 36.4 km). Each epoch was one hour in duration, yielding 38 minutes on source. The observations were obtained 79 with 4096 MHz of bandwidth covering the entire C band (4–8 GHz), divided into 32 spectral windows each of 128 MHz bandwidth and sampled with 64 frequency channels. The VLA data were calibrated using J1352-4412 for complex gains, and 3C147 for abso- lute gain and bandpass. Data were edited, reduced, and imaged using standard routines in AIPS (Greisen, 2003) and CASA (version 4.5.0; McMullin et al., 2007). To minimize artifacts from wide fractional bandwidths and to constrain the spectral index, we divided the data into two frequency chunks (4–6 GHz, and 6–8 GHz), and imaged each chunk separately. The central frequencies of these subbands are 5.1 and 7.1 GHz, with beams of 1.91” × 0.32” and 1.40” × 0.27”, respec- tively. We also note that these observations are all taken at very low elevation ((cid:46) 12◦), as J1417 is near the southern declination limit of the VLA; these conditions somewhat hamper our sensitivity. J1417 is detected in the combined 4 hr image at 5.1 GHz (22.4 ±4.8 Jy) and shows marginal evidence of radio emission at 7.1 GHz (13.5 ± 5.1 Jy). However, after imaging the individual 1-hr blocks separately, it is clear that this combined flux density does not accurately represent the system: the source is not detected at either 5.1 or 7.1 GHz in three of the four blocks, while it is detected at high significance in both subbands during the second block (on UT 2015 Jun 16), with flux densities of 56.5 ± 9.8 (5.1 GHz) and 40.9 ± 9.5 Jy, respectively (Table 3.3). For flux density  ∝  with frequency , we find a spectral index  = –1.0 ± 0.9 at this epoch. This spectral index does not constrain the nature of the source: it is consistent with either the steep-spectrum slope of an MSP ( ∼ -1.6, e.g., Lorimer et al., 1995; Maron et al., 2000) or the flat spectrum emission observed for tMSPs in the disk state (Papitto et al., 2013; Bassa et al., 2014; Archibald et al., 2015; Deller et al., 2015). We summarize our VLA continuum observations in Table 3.3. We also separately imaged the individual ∼8 min scans of the block with the strongest 5.1 GHz detection (2015 Jun 16) to search for variability on these shorter timescales. The source was detected in some scans and not in others, consistent with its modest flux density and the higher rms noise of the individual scans. We found no evidence for statistically significant variations in flux density on timescales of minutes. 80 Although radio timing measurements do not always give an accurate flux density scale, Camilo et al. (2016) report a 4.8 GHz flux density of ∼ 40 Jy for J1417 at one epoch of GBT observations, generally consistent with the flux density we observe during the epoch in which the source is detected. We discuss the observed variations in the flux density below in §3.7. 3.3.2 ATCA We observed J1417 with the Australia Telescope Compact Array (ATCA) on 2015 Oct 24, a few months after our VLA observations, under program ID CX336 (PI Miller-Jones). We observed in the relatively extended 6A configuration (max ∼ 5.9 km) with the 16cm receiver, providing coverage from 1–3 GHz. The bandwidth is sampled by 2048 frequency channels, each 1 MHz wide. We obtained 9.4 hours of data on-source, and calibrated the complex gains using PKS B1424–418 and the absolute flux and bandpass using PKS 1934–638. These data covered an orbital phase range of  = 0.782 to 0.824. We also observed J1417 with ATCA on 2016 June 6 (program ID CX360, PI Miller-Jones), simultaneous with our Chandra X-ray observations. The array was in the 1.5B configuration (max ∼ 4.3 km), and we again observed with the 16cm receiver and the same correlator set-up and calibrators as the 2015 Oct observations. We obtained 5.0 hours of data on-source, over an orbital phase range of  = 0.737 to 0.760. For both observations, the ATCA data were edited and reduced using standard routines in Miriad (Sault et al., 1995). To minimize imaging artifacts caused by wide bandwidths, we split the data into frequency chunks (three for the 2015 Oct data and four for the 2016 June data). Each frequency chunk was self-calibrated and imaged in AIPS. We summarize our ATCA observations in Table 3.4. For the 2015 Oct 24 data, the source was well-detected in each chunk, with flux densities of 294 ± 42 Jy (1.43 GHz), 88 ± 27 Jy (1.89 GHz), and 145 ± 21 Jy (2.59 GHz). For a simple power-law fit these measurements imply a spectral index of –1.1 ± 0.3, in excellent agreement with the VLA measurement, though the flux densities together are not clearly well-described by 81 this model. This may suggest the presence of frequency-dependent refractive scintillation (e.g., Stinebring & Condon, 1990), frequency-dependent eclipses (e.g., Broderick et al., 2016), or both. For the 2016 June 6 data, the lower on-source time and more compact configuration led to a higher rms noise, and the source was only marginally detected at 2–3 in each chunk: 258 ± 90 Jy (1.41 GHz), 135 ± 71 Jy (1.86 GHz), 86 ± 47 Jy (2.37 GHz), 131 ± 51Jy (2.83 GHz). Nevertheless, the spectral index implied by a power-law fit to these values is –1.2 ± 0.8, again consistent with the previous ATCA and VLA measurements. In fact, if we fit a single power-law with either the 2015 October or 2016 June ATCA data together with the 2015 Jun 16 VLA data, the best-fit spectral index is –1.0 ± 0.2 or –1.0 ± 0.3, broadly consistent with emission from a MSP and not from a LMXB jet. 3.4 Optical Observations 3.4.1 SOAR Spectroscopy We have continued occasional spectroscopic monitoring of J1417 with the Goodman spectrograph (Clemens et al., 2004) on the SOAR telescope, obtaining 10 epochs of spectroscopy from 2015 Jun 23 to 2017 Sep 01 (UT). The setup and reduction was identical to that from Strader et al. (2015), giving calibrated spectra with wavelength coverage of ∼ 5375–6640 Å at a resolution of about 1.7 Å. We present a subset of these data and analyze them, along with a number of spectra obtained by Strader et al. (2015) between 2014 Mar and 2015 Feb, in §3.6. 3.4.2 Optical/Near-IR Photometry We obtained optical   and near-IR  photometry of J1417 using ANDICAM on the SMARTS 1.3-m telescope at CTIO roughly every night between 25 Jul and 19 Sep 2017 (UT). On each night, we took 200, 90, and 50 sec exposures in , , and  bands, respectively. We simultaneously obtained dithered near-IR observations in . After excluding poor frames, the on-source exposure 82 time in  was typically 120-200 sec per visit. Data were reduced following the procedures in Walter et al. (2012). We performed differential aperture photometry to obtain instrumental magnitudes of the target source using fifteen ( ) or seven () comparison stars as a reference. Absolute calibration was done with respect to the Landolt (1992) standard field TPheD ( ) or to the 2MASS catalog (). Our final dataset includes 36 measurements each in   and 28 in , with mean observed magni- tudes  = 17.18,  = 16.16,  = 15.00, and  = 13.38. Median errors in each band are ∼ 0.03 in   and ∼ 0.07 in . Our full sample of SMARTS photometry is listed in Table 3.5. Observation epochs have been corrected to Barycentric Julian Date (BJD) on the Barycentric Dynamical Time (TDB) system (Eastman et al., 2010). Strader et al. (2015) obtained and analyzed time series photometry of J1417 in   taken in mid-2014 using the PROMPT-5 telescope (Reichart et al., 2005). This dataset included 29, 49, and 40 photometric measurements in , , and  bands respectively. After applying a small zeropoint offset to PROMPT  and  to match the ANDICAM filters, we reanalyze this sample in conjunction with the SMARTS data in the following section. 3.5 Light Curve Fitting Strader et al. (2015) analyzed the PROMPT   photometry under the assumption that the light curves were dominated by ellipsoidal variations due to the tidal deformation of a Roche lobe filling secondary, finding that the inclination was  = 58 ± 2◦. Here we relax these assumptions and provide more comprehensive light curve fits to the full photometric dataset discussed above. Although the overall scale of the binary system was already well constrained by the spectro- scopic results presented in Strader et al. (2015), in our light curve fits we adopt the pulsar ephemeris from Camilo et al. (2016), which provide improved precision on the orbital period and mass ratio. Throughout our analysis, we assume a circular orbit (consistent with the optical radial velocity curve), and fix the orbital period , radial velocity semi-amplitude of the secondary 2, mass 83 Table 3.5. SMARTS Photometry of J1417 BJD–2450000 (d) 7960.480563 7960.483566 7960.485358 7960.481234 7964.485152 7964.488149 7964.489947 7964.485814 7968.496390 ... Band mag err B V I H B V I H B ... 17.143 16.105 14.979 13.382 17.167 16.196 15.045 13.502 17.190 ... 0.022 0.039 0.056 0.073 0.149 0.113 0.080 0.111 0.039 ... Note. — The full SMARTS dataset is avail- able in machine-readable format. We show a portion of the table here as a preview of its form and content. Magnitudes have not been corrected for extinction. ratio , and projected semimajor axis of the pulsar orbit 1, to the values presented by Strader et al. (2015) or Camilo et al. (2016):  = 5.37372 d, 2 = 115.7 km s-1,  = 2/1 = 0.171, and 1 = 1 sin  = 4.876(9) lt-s. We model the light curves with the Eclipsing Light Curve (ELC; Orosz & Hauschildt, 2000a) code, which fits normalized, filter-specific light curves in each band independently. We fit two main classes of models. The first assumes no accretion disk is present. In these “No Disk” fits, the free parameters are the inclination , the Roche lobe filling factor of the secondary 2, and its intensity-weighted mean temperature 2. Given the possible evidence for a disk from the H spectroscopy, we also explore fits in which a non-negligible accretion disk component was present. Under these circumstances, the secondary was assumed to completely fill its Roche lobe ( 2 = 1), while the disk is characterized by five parameters: the inner disk temperature disk, the inner (in) and outer (out) radii of the disk, the power-law index of the disk temperature profile , and the opening angle of the disk rim . 84 Figure 3.5: Optical/near-IR light curves of 1FGL J1417.7–4407. PROMPT (diamonds) and SMARTS (circles) optical/near-IR photometry with best fit ELC models. Two orbital phases are shown for clarity. The black lines show 8 random samples from the posterior marginalized over all the free parameters in the SMARTS+PROMPT “No Disk” MCMC run (Table 3.6, top panel, column 1). We present fits using both the SMARTS data only and with both datasets together; given the consistency of the two datasets, using both together provides the best time and phase coverage. These ELC fits were done using its included Markov Chain Monte Carlo (MCMC) sampler, and we summarize the best-fit values as the median of the posterior distribution, with 1 equivalent uncertainties given as the 15.9 and 84.1% quantiles of the posterior. These values are given in Table 3.6, and sample models are plotted in Figure 3.5. Following the recommendations of Hogg & Foreman-Mackey (2018) for reporting results of an MCMC run, we refrain from showing only a “best-fit” model, taken as the collective medians 85 0.00.51.01.52.0Phase1314151617magBVRIH of each individual posterior distribution. Instead, we have selected a few randomly drawn example samples from the posterior of our MCMC run marginalized over all parameters, and plot these with the data in Figure 3.5. This ensures we retain probabilistic information about the sampling results while also showing models that are satisfactory fits to the data. The first main result from the fits is that there is no evidence for a disk in the optical/near-IR light curves: its inclusion does not improve the quality of the model fits, and in the best-fitting case it makes up a tiny fraction of the light (∼ 1% in ). Either a disk is not present or it is so faint compared to the giant secondary that its presence is irrelevant in the broadband photometry. Considering the fits with no disk, the main result is the well-known degeneracy between the inclination and the filling factor. In the absence of significant irradiation, the relative amplitude of ellipsoidal modulations depends mainly on the observable projected area of the secondary. There- fore, light curve models with lower inclinations (more face-on) and larger filling factors will look essentially identical to models of more inclined systems whose secondary fills a smaller fraction of its Roche lobe. We find a filling factor 2 = 0.83+0.05−0.07 and an associated inclination 64.◦2+5.0−7.8. The slightly higher (and more uncertain) inclination than found in Strader et al. (2015) leads to a slightly lower primary mass: 1 = 1.62+0.43−0.17(cid:12), but with a larger uncertainty. If instead we fix the secondary to fill its Roche lobe (Table 3.6, top panel, column 2), then  = 55.◦2+0.7−1.0 and 1 = 2.13+0.08−0.06 (cid:12), consistent with the neutron star mass found by Strader et al. (2015). In addition to the above models, we also explored ones that included irradiation. This would appear in the optical/near-IR light curves as extra luminosity on the “day” side (the side of the tidally locked companion facing the pulsar) as the star is heated directly by high-energy photons from the pulsar (e.g., Romani & Sanchez, 2016), or indirectly through reprocessing of the pulsar spin-down power in an intra-binary shock (e.g., Romani et al., 2015a,b; Rivera Sandoval et al., 2017; Wadiasingh et al., 2017). 86 It is evident from our light curves that their behavior is well-described by standard ellipsoidal variations, with no obvious evidence for irradiation. This was reflected in our ELC fits that included irradiation as a free parameter, which found a negligible contribution from an irradiating source. Obviously, high-energy emission is present in this system, but we find no sign that it affects the optical/near-IR light curves. Overall, we have improved on the light curve fitting presented by Strader et al. (2015) by exploring a much larger parameter space in conjunction with additional data and more precise dynamical constraints. The model with the fewest components that best represents the data is the SMARTS+PROMPT “No Disk” model with the filling factor left to vary (Table 3.6, top panel, column 1), which we focus on for the remainder of the paper. The component masses associated with this fit are 1 = 1.62+0.43−0.17 (cid:12) and 2 = 0.28+0.07−0.03 (cid:12), well within the range of typical values inferred for other neutron star binaries associated with Fermi sources (e.g., Ransom et al., 2011; Romani et al., 2015b; Strader et al., 2016; Halpern et al., 2017; Shahbaz et al., 2017). 3.5.1 Distance Following similar steps as detailed in Strader et al. (2015) and Swihart et al. (2017), we estimate the distance to J1417 by comparing its observed magnitude to its intrinsic luminosity derived from the best-fit radius and effective temperature of the secondary from our optical/near-IR light curve models. After applying bolometric corrections to each band as a function of temperature and metallicity from fitting 10 Gyr, [Fe/H] = –0.65 isochrones (Marigo et al., 2008), we compared the extinction-corrected magnitudes (corrected using the Schlafly & Finkbeiner (2011b) reddening maps) to the predicted absolute magnitude in each band. Given the inferred effective temperature and Roche lobe filling factor and their uncertainties (Table 3.6, top panel, column 1), we find a distance of 3.1 ± 0.6 kpc. For a Roche lobe filling secondary ( 2=1.0) at the same temperature, the distance would instead be 3.6 ± 0.5 kpc. From consideration of the mean colors and two low-resolution spectra of J1417, Strader et al. 87 Table 3.6. Summary of ELC fits for J1417 SMARTS+PROMPT SMARTS+PROMPT ( 2 = 1.0 fixed) SMARTS+PROMPT ( 2 = 0.83 fixed) SMARTS Only Fitted Parameters incl (◦) k R.L filling factor ( 2) s i D o N 2 (K) 2 (dof) Derived Parameters 1 ((cid:12)) 2 ((cid:12)) 2 ((cid:12)) log(g) (cgs) Fitted Parameters incl (◦) 2 (K) disk (K) inb outb  (◦)  2 (dof) Derived Parameters 1 ((cid:12)) 2 ((cid:12)) 2 ((cid:12)) log(g) (cgs) k s i D 64.2+5.0−7.8 0.83+0.05−0.07 4560+460−336 386.1(250) 1.62+0.43−0.17 0.28+0.07−0.03 3.66+0.30−0.33 2.76 ± 0.01 56.0+1.4−2.1 4633+259−536 3191+1480 −2159 0.22+0.20−0.18 0.58+0.20−0.17 15.1+5.7−14.3 -0.75+0.18−0.23 388.0(246) 2.08 ± 0.15 0.36 ± 0.03 4.16 ± 0.10 2.75 ± 0.01 55.2+0.7−1.0 1.0a 4523+297−317 390.8(251) 2.13+0.08−0.06 0.36 ± 0.01 4.20+0.05−0.04 2.75 ± 0.01 64.0+1.1−1.4 0.83a 4559+270−344 386.1(251) 1.63+0.06−0.05 0.28 ± 0.01 3.66+0.04−0.03 2.75 ± 0.01 · · · · · · · · · · · · · · · · · · · · · ... · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ... · · · · · · · · · · · · 62.3+6.7−9.6 0.83+0.06−0.10 4509+388−415 127.1(132) 1.70+0.65−0.25 0.29+0.11−0.04 3.72+0.52−0.64 2.76 ± 0.01 54.7+2.2−3.4 4679+1118 −684 3261+1919 −3095 0.25+0.22−0.18 0.66+0.18−0.22 14.6+4.5−14.1 -0.74+0.39−0.27 127.9(128) 2.18+0.28−0.29 0.37 ± 0.05 4.23+0.18−0.20 2.76 ± 0.02 aFilling factor held constant bInner and outer disk radii are expressed as a fraction of the effective Roche-lobe of the primary 88 (2015) assumed the effective temperature of the secondary is around 5000K, about 450 K warmer than the effective temperature we find in the new light curve analysis. This warmer temperature, combined with the assumption of a Roche Lobe-filling secondary, led to a larger, brighter star and a greater inferred distance of 4.4 kpc. As mentioned above, these fits were performed at the subsolar metallicity inferred from the optical spectrum (Strader et al., 2015), but we note that assuming solar metallicity instead produces consistent distances to within 0.1 kpc. J1417 is present in Gaia DR22 with parallax 0.22 ± 0.07 mas. Given the low significance of this measurement, we use a Bayesian inference procedure recommended by Bailer-Jones et al. (2018) to estimate the distance to J1417. Using the weak distance prior of Gandhi et al. (2018) based on the expected distribution of LMXBs in the Galaxy, the geometric distance estimate is 4.2+2.2−1.1 kpc (1). Jennings et al. (2018) adopt a distance prior based on a MSP-specific model for the space density, estimating 3.50+2.09−0.38 kpc. Both these values are fully consistent with our distance estimates from the light curve model, but not with the DM-derived distance (1.6 kpc) from Camilo et al. (2016). For the remainder of the paper, we adopt the distance inferred from our light curve fits (3.1± 0.6 kpc). 3.6 H Spectroscopy 3.6.1 The Spectra Themselves In tMSPs, the presence of a broad, double-peaked H emission profile has generally been accepted as evidence for an accretion disk (e.g., Wang et al., 2009; Bogdanov et al., 2015). Yet the observation of J1417 as a radio pulsar, and the other evidence that this system is in the pulsar state (i.e. the X-ray light curve and radio spectral index), is in opposition to this evidence for a disk. Here we carefully consider the morphology of the H emission in J1417 and discuss models that can explain its origin. First, we echo the basic statement in Strader et al. (2015) that the majority of the H profiles 2https://www.cosmos.esa.int/web/gaia/dr2 89 Figure 3.6: Phase-resolved optical spectra of 1FGL J1417.7–4407 near the H region. Example optical spectra of J1417 around the H region. All spectra have been shifted to the rest frame of the secondary (dashed line). Each panel shows spectra that are representative of nearly all the spectra obtained at that orbital phase. 90 φ=0.438φ=0.5162015Feb172014Mar15I-1000-50005001000RadialVelocity(km/s)φ=0.043φ=0.1412017Jul102017Jul31II-1000-50005001000RadialVelocity(km/s)Flux(Arbitrary)φ=0.7172014Aug25IIIφ=0.158φ=0.1852017Jul212014May1IV65406550656065706580Wavelength(˚A)φ=0.366φ=0.4412017Feb112014Jun14V65406550656065706580Wavelength(˚A)φ=0.5372015Feb17VI Figure 3.7: Orbital position of the secondary in 1FGL J1417.7–4407 corresponding to the varying H phenomenology. A schematic diagram of J1417 (as seen from Earth) showing the orbital position of the secondary corresponding to the range of H phenomenology we show in the panels of Figure 3.6. show two peaks. However, we also find that the components vary in both shape and amplitude over time, and that these variations appear to be associated with the orbital phase of the binary in a consistent manner over the 3.5-year time baseline of our optical spectra. In Figure 3.6 we show example spectra of J1417 around the H region. These have all been shifted to the rest frame of the secondary. We show six panels, each with one or two spectra, that are representative of nearly all the spectra obtained at that phase. In panels I & II, the emission is double-peaked but asymmetric, with stronger blue or red components, respectively. Here the phases are near 0.5 and 0 (at quadrature). In III & IV (taken near conjunction) one of the components has almost completely vanished. In general, when the secondary is receding from Earth along our line-of-sight (=0.25–0.75), the H profile has a stronger blue component when compared to the 91 =0.438=0.5162015Feb172014Mar15I-1000-50005001000RadialVelocity(km/s)=0.043=0.1412017Jul102017Jul31II-1000-50005001000RadialVelocity(km/s)Flux(Arbitrary)=0.7172014Aug25III=0.158=0.1852017Jul212014May1IV65406550656065706580Wavelength(˚A)=0.5372015Feb17V65406550656065706580Wavelength(˚A)=0.366=0.4412017Feb112014Jun14VI=0.438=0.5162015Feb172014Mar15I-1000-50005001000RadialVelocity(km/s)=0.043=0.1412017Jul102017Jul31II-1000-50005001000RadialVelocity(km/s)Flux(Arbitrary)=0.7172014Aug25III=0.158=0.1852017Jul212014May1IV65406550656065706580Wavelength(˚A)=0.5372015Feb17V65406550656065706580Wavelength(˚A)=0.366=0.4412017Feb112014Jun14VI=0.438=0.5162015Feb172014Mar15I-1000-50005001000RadialVelocity(km/s)=0.043=0.1412017Jul102017Jul31II-1000-50005001000RadialVelocity(km/s)Flux(Arbitrary)=0.7172014Aug25III=0.158=0.1852017Jul212014May1IV65406550656065706580Wavelength(˚A)=0.5372015Feb17V65406550656065706580Wavelength(˚A)=0.366=0.4412017Feb112014Jun14VI=0.438=0.5162015Feb172014Mar15I-1000-50005001000RadialVelocity(km/s)=0.043=0.1412017Jul102017Jul31II-1000-50005001000RadialVelocity(km/s)Flux(Arbitrary)=0.7172014Aug25III=0.158=0.1852017Jul212014May1IV65406550656065706580Wavelength(˚A)=0.5372015Feb17V65406550656065706580Wavelength(˚A)=0.366=0.4412017Feb112014Jun14VI=0.438=0.5162015Feb172014Mar15I-1000-50005001000RadialVelocity(km/s)=0.043=0.1412017Jul102017Jul31II-1000-50005001000RadialVelocity(km/s)Flux(Arbitrary)=0.7172014Aug25III=0.158=0.1852017Jul212014May1IV65406550656065706580Wavelength(˚A)=0.366=0.4412017Feb112014Jun14V65406550656065706580Wavelength(˚A)=0.5372015Feb17VI=0.438=0.5162015Feb172014Mar15I-1000-50005001000RadialVelocity(km/s)=0.043=0.1412017Jul102017Jul31II-1000-50005001000RadialVelocity(km/s)Flux(Arbitrary)=0.7172014Aug25III=0.158=0.1852017Jul212014May1IV65406550656065706580Wavelength(˚A)=0.366=0.4412017Feb112014Jun14V65406550656065706580Wavelength(˚A)=0.5372015Feb17VIOrbital Motion(i ~ 64°)orbital plane red (Figure 3.6, panels I, III), and vice versa when the companion approaches Earth (=0.0–0.25 & 0.75–1.0; panels II, IV). Panel V shows two spectra taken between these extremes, which shows nearly symmetric H emission. Panel VI shows the most outlying spectrum observed, with unique phenomenology in our sample: it has a bright, narrow main peak superimposed on a broader, asymmetric component. In Figure 3.7, we show a schematic diagram of the binary at the orbital phases corresponding to the six panels in Figure 3.6. For a typical H equivalent width of about 5Å and the range of distances from the last section, we find that the H luminosity is ∼ 1031 erg s−1. 3.6.2 Understanding the H Now we consider possible origins of the H emission and how it would vary with orbital phase. H emission might conceivably originate (1) from the outer regions of an accretion disk, (2) from the chromosphere of the secondary, (3) in material lost from the secondary, such as an outflow or wind, possibly enhanced by an intrabinary shock between the pulsar and the secondary. We consider these in turn. 3.6.2.1 A Disk Double-peaked Balmer emission is generally interpreted as arising from an accretion disk, and has been observed for the tMSPs PSR J1023+0023 and XSS J12270–4859 in their disk states (e.g., Wang et al., 2009; de Martino et al., 2014a). Such emission is not observed during the pulsar state of these systems, presumably due to the disappearance of the disk. The spectra obtained of PSR J1023+0023 in its disk state in 2013 Oct (Bogdanov et al., 2015) show slightly asymmetric peaks, with a stronger blue component, and the overall shape and binary phase of these profiles are roughly consistent with the spectra we show in panels I & V of Figure 3.6. 92 Coti Zelati et al. (2014) also show the H profile of PSR J1023+0038 in its disk state taken a few months after Halpern et al. (2013). The emission is still double-peaked, but with a slightly larger red peak than the blue, similar to Figure 3.6, panel II. However, the spectrum they present is a co-addition of four exposures taken at a number of binary phases (∼0.25–0.6, Coti Zelati, private communication) so we cannot make a direct phase-resolved comparison to J1417. However, several pieces of evidence suggest that the J1417 H emission is unlikely to originate from a standard LMXB accretion disk. First, the minimum of the double-peaked profile is station- ary in the rest frame of the secondary. This might conceivably be explained if the secondary has strong H absorption that sweeps through an essentially fixed emission profile (the primary motion is only ∼ 20 km s−1), but this does not explain why the emission components also appear to lie at essentially the same velocity as a function of orbital phase in the secondary rest frame. Second, as pointed out by Strader et al. (2015), the peak separation when the “classic" double- peaked profile is evident (Figure 3.6, panel V) is comparable to the orbital velocity, which would unphysically suggest that the disk fills the entire binary. (At other phases the peak separation is broader and more easily accommodated.) Finally, there is no independent evidence of a disk: one is not necessary to fit the optical/near- IR light curves, and Camilo et al. (2016) find no evidence for UV emission from a disk using Swift/UVOT observations. We do note that these constraints are weaker than for typical redbacks owing to the higher luminosity of the giant. 3.6.2.2 Stellar Chromosphere If the emission were coming directly from the atmosphere of the giant, as is the case in chromospher- ically active stars such as FK Com and RS CVn systems (Drake, 2006), the H emission would be relatively narrow and track the secondary’s motion. It would not be expected to vary substantially with orbital phase. Instead, the H emission we observe is broad and variable in velocity, and is never observed at the systemic velocity of the secondary, which instead shows (likely intrinsic) H 93 absorption. Thus the chromosphere of the secondary is not the main contributor to the emission. 3.6.2.3 Stellar Wind An outflow or wind from the secondary would form a natural medium for the production of an H profile that varies, but is centered on the velocity of the secondary. The observed rotational velocity of the star is consistent with the star being tidally locked to the neutron star (Strader et al., 2015), as expected for a close compact binary in which extensive mass transfer has occurred. For low-mass stars like J1417, this rapid rotation can result in a strong dynamo that can boost the surface magnetic fields by a factor of ∼ 102–103 compared to similar isolated stars (Morin, 2012). This enhanced surface field can lead to strong magnetically driven winds, since mass loss from low-mass red giants is caused primarily by Alfvén wave pressure (Cranmer & Saar, 2011). For J1417 specifically, using the values in Table 3.6, the escape velocity from the secondary is no larger than 170–180 km s−1. For all of the emission profiles a substantial (and sometimes dominant) amount of the emission is present at velocities larger than the escape velocity of the secondary, showing that the gas is not simply an outflow but a wind that has escaped the red giant. The wind from J1417 could emit H photons directly given that a source of ionization is present, or indirectly as it interacts with the pulsar wind, forming an intrabinary shock near the companion surface (Romani & Sanchez, 2016; Wadiasingh et al., 2017). Such shocks are likely responsible for at least some of the observed X-ray variability in redbacks and the known tMSPs (e.g., Roberts et al., 2014) and provide a natural location for a region of gas with high enough temperature to emit H photons. This latter scenario can at least partially explain the complicated morphology of the emission: the red giant wind would not flow away from the star unimpeded, but instead would be halted and swept away by the pulsar wind, perhaps wrapping around the binary and extending out beyond the orbit of the secondary. 94 This “radio ejection” mechanism (Burderi et al., 2002) has been used successfully to explain the H emission in the globular cluster binary PSR J1740-5340A, which hosts a stripped subgiant in a 1.3 day orbit (D’Amico et al., 2001a; Bogdanov et al., 2010; Mucciarelli et al., 2013). Sabbi et al. (2003) analyzed the phase-resolved H emission profiles of the companion to PSR J1740-5340A, finding that a combination of material flowing to the inner Lagrange point and gas swept out beyond the companion’s orbit could explain the broad components not associated with the stellar atmosphere. The overall shape and binary phase of some of our spectra are remarkably similar to the profiles seen in the optical spectra of PSR J1740-5340A. In particular, the spectrum we present in panel VI of Figure 3.6 shows a main, narrow peak as well as a broad, fainter component that is analogous to the spectrum analyzed in detail by Sabbi et al. (2003). An important difference is that, contrary to PSR J1740-5340A, the offset of the main peak by >100 km s-1 from the radial velocity of the companion for J1417 implies this emission feature does not likely originate in the stellar atmosphere, echoing our previous discussion. Furthermore, in the spectra we show nearest to opposition and conjunction of the secondary (panels III & IV, respectively), the profiles show markedly less pronounced double-peak morphology; while one of the H peaks is very broad and enhanced, the other virtually disappears. This odd Doppler-shifted “flaring” has not been observed in tMSPs in either state, but similar profiles in the black hole binary V404 Cyg have been explained as being due to an outflow from the central regions of the binary that absorbs one of the components (Hynes et al., 2002). 3.6.2.4 The Source of H Given all the evidence, we suggest that the strange phenomenology in H is likely not due to emission from an accretion disk, but instead comes from some combination of material flowing to the inner Lagrange point, hot gas swept out beyond the companion’s orbit by the pulsar’s radiation pressure, and an intrabinary shock. H photons coming from an intrabinary shock or giant winds would be slightly offset from the orbital motion of the secondary, and display a very broad range of velocities, as we observe in nearly all our spectra. In both these circumstances, the emitting 95 gas would slightly trail the secondary in its orbit. Matter streaming to the inner Lagrange point, which is then pushed away by the pulsar wind, would also generally track the orbital motion of the secondary, but always be slightly shifted, trailing the companion and moving roughly perpendicular to the line connecting the primary and secondary (see Figure 2 in Burderi et al., 2002). In general, at binary phases  ∼ 0.25–0.75, these H emitting photons would be blueshifted with respect to the secondary, and vice versa at phases =0.0–0.25 & 0.75–1.0. As described at the beginning of this section, this is precisely the global trend we see in the profiles across all epochs (Figure 3.7). Whether this qualitative scenario can reproduce the detailed profiles at all epochs, as well as the X-ray light curve of the binary, remains to be tested with improved data and future detailed modeling. 3.7 Conclusions We have presented multiwavelength follow-up observations of the MSP binary J1417, a system likely in the late stages of the standard MSP recycling process that began after the secondary evolved off the main sequence and filled its Roche Lobe. Recently, Swihart et al. (2017) reported on the dis- covery of a compact binary assocated with the unidentified Fermi source 2FGL J0846.0+2820. This system, similar to J1417, also has an inferred heavy neutron star primary with a giant secondary in a wide 8.1 d orbit. Although no radio pulsations have yet been discovered in 2FGL J0846.0+2820, given the clear differences between these systems and the “redback" and “black widow" subclasses of MSP binaries, we have termed these systems with giant companions huntsman binaries after the large huntsman spider. From the perspective of the pulsar, the solid angle subtended by the secondary in J1417 is 1.3–1.4%, comparable to the value for typical redback systems with unevolved donors (Roberts et al., 2015). However, if some of the X-ray emission comes from a shock close to the secondary, then the effective solid angle is larger. While the light curve itself shows no evidence for irradiation (§ 3.5), consistent with the high luminosity of the evolved secondary, the evidence for a wind from the optical spectroscopy shows that a substantial amount of material is leaving the system. This 96 material is likely responsible for the difficulty in detecting pulsations from the MSP, even compared to typical redbacks (Camilo et al., 2016). It is possible that similar huntsman systems remain undetected for this same reason. The radio flux densities for J1417 are consistent with those expected for a MSP with modest scintillation, excepting the few VLA epochs where the pulsar is not detected. As it happens, these epochs were at superior or inferior conjunction, which is when the pulsar was also not detected in timing observations (Camilo et al., 2016). This would be consistent with a scenario in which the radio emission was absorbed rather than scattered by the material in the system. On the other hand, both ATCA detections were obtained at or just after inferior conjunction. Hence there may be both stochastic and orbital components to the variations in eclipsing material, consistent with the conclusions we draw for the H variations in §3.6. Strader et al. (2015) pointed out that the 3FGL Fermi-LAT spectrum of J1417 was a power-law with no evidence for curvature. This is unusual among MSPs, and is more similar to the -ray bright low-mass X-ray binaries 3FGL J1544.6–1125 and 3FGL J0427.9–6704 (Acero et al., 2015; Bogdanov, 2016; Strader et al., 2016). If this power-law spectrum is confirmed in the forthcoming 4FGL Fermi-LAT catalog, then it could suggest another source of high-energy emission in the system beyond simply the pulsar magnetosphere. The new X-ray observations have allowed a better comparison of the X-ray properties of J1417 to known redbacks and tMSPs. At our new optically inferred distance (3.1 ± 0.6 kpc),   = 0.7− 1.4× 1033 erg s−1, which is more typical of tMSPs in their accretion-powered state and brighter than all known redbacks in their pulsar state. If there is no accretion disk, as we suspect, it is likely that the bulk of this X-ray emission is coming from a strong intrabinary shock near the surface of the secondary. The tidally locked giant could generate a substantial magnetic field that may facilitate the formation of a luminous shock. The most important new data to interpret J1417 would be an X-ray light curve complete in orbital phase, elucidating the nature of an intrabinary shock and allowing modeling of the H 97 emission in this context. Forthcoming Gaia data releases should address systematics in the parallax measurements and definitively confirm the unusually high X-ray luminosity of the system. The unusual properties of this huntsman binary motivate ongoing searches for similar pulsar binaries. 98 PSR J1306–40: AN X-RAY LUMINOUS REDBACK WITH AN EVOLVED COMPANION CHAPTER 4 1PSR J1306–40 is a millisecond pulsar binary with a non-degenerate companion in an unusually long ∼1.097 day orbit. We present new optical photometry and spectroscopy of this system, and model these data to constrain fundamental properties of the binary such as the component masses and distance. The optical data imply a minimum neutron star mass of 1.75 ± 0.09 (cid:12) (1-sigma) and a high, nearly edge-on inclination. The light curves suggest a large hot spot on the companion, suggestive of a portion of the pulsar wind being channeled to the stellar surface by the magnetic field of the secondary, mediated via an intrabinary shock. The H line profiles switch rapidly from emission to absorption near companion inferior conjunction, consistent with an eclipse of the compact emission region at these phases. At our optically-inferred distance of 4.7 ± 0.5 kpc, the X-ray luminosity is ∼1033 erg s-1, brighter than nearly all known redbacks in the pulsar state. The long period, subgiant-like secondary, and luminous X-ray emission suggest this system may be part of the expanding class of millisecond pulsar binaries that are progenitors to typical field pulsar–white dwarf binaries. 4.1 Introduction Millisecond pulsars (MSPs) form in binary systems, where the central neutron star accretes matter and angular momentum from a non-degenerate companion, spinning it up to very rapid spin periods. Although this mass transfer process is expected to last hundreds of Myrs or more (Tauris & Savonije, 1999), most MSPs in the Galactic field have degenerate white dwarf companions in wide orbits (orb (cid:38) 5 d), which represent the end stage of the recycling process (Tauris & van den Heuvel, 2006). Recent multiwavelength (optical, radio, and X-ray) follow-up observations of unidentified Fermi 1This chapter is taken mostly word-for-word from Swihart et al. (2019), as published in the Astrophysical Journal 99 Large Area Telescope (LAT) -ray sources have revealed a substantial number of close neutron star binary systems that host hydrogen-rich secondaries rather than He white dwarfs. Unlike most field MSPs (as described above), these appear to be systems in which the standard recycling process is not yet complete (e.g., Benvenuto et al., 2014), providing valuable insights into the spin-up process. These binaries are classified by the mass of the secondary: “redbacks" have non-degenerate, main sequence-like companions ( (cid:38) 0.1 (cid:12)), compared to the less massive ( (cid:46) 0.05 (cid:12)), highly ablated, semi-degenerate “black widows” (Roberts, 2011). Both systems show radio eclipses due to ionized material from the secondary; these eclipses are typically more extensive for redbacks (e.g., D’Amico et al., 2001b; Camilo et al., 2015; Cromartie et al., 2016). A few redbacks have been observed to transition back and forth between an accretion-powered disk state and a rotationally-powered pulsar state. These “transitional millisecond pulsars” proved the suspected evolutionary link between some recycled MSPs and neutron star low-mass X-ray binaries (e.g., Archibald et al., 2009; Papitto et al., 2013; Bassa et al., 2014). However, like the bulk of the redback population, the short orbital periods in these systems ((cid:46) 0.5 d) mean they are unlikely to end their lives as the wide-orbit MSP–He white dwarf binaries that dominate the observed population of MSPs in the field (Chen et al., 2013). Instead, the progenitors to these canonical MSP–white dwarf systems are neutron stars with evolved red giant companions whose orbital periods are > 1 day (Tauris & Savonije, 1999). In this context, the MSP binary 1FGL J1417.7–4407 (PSR J1417–4402) was the first MSP binary discovered that is a likely progenitor to the typical MSP-white dwarf binaries (Strader et al., 2015; Camilo et al., 2016). Due to its long period, red giant companion, and inferred evolutionary track, a new “huntsman” subclass was coined to distinguish unique systems like these from the more common redbacks (Strader et al., 2015; Swihart et al., 2018). Optical spectra of 1FGL J1417.7–4407 show a strong, persistent double-peaked H emission profile that is unusual among redbacks in their pulsar states and reminiscent of a classic accretion disk. However, due to the small separation between the peak components, and the stationary profiles 100 as a function of orbital phase in the rest frame of the secondary, this H phenomenology is likely not due to a disk. Instead, the emission likely comes from a combination of material swept off the companion and a strong intrabinary shock (Camilo et al., 2016; Swihart et al., 2018). Another piece of evidence for an unusually luminous shock is the X-ray luminosity in this system (∼ 1033 erg s-1), which is higher than other redbacks in the pulsar state. Swihart et al. (2018) discuss the possibility that the shock luminosity in 1FGL J1417.7–4407 is enhanced over typical redbacks by a strong, magnetically driven wind from the tidally-locked red giant secondary. The candidate MSP binary 2FGL J0846.0+2820 shows many similarities to 1FGL J1417.7–4407, making it a possible second member of the huntsman subclass, though no radio pulsar has yet been confirmed in this system (Swihart et al., 2017). Pending future evolution studies, it might be reasonable to include other long period redback-like systems with evolved companions in this class (e.g., PSR NGC 6397A). The subject of this paper, PSR J1306–40, was discovered in the SUPERB survey (Keane et al., 2018), where it was identified as a candidate redback binary. No pulsar timing solution was pre- sented owing to the frequent eclipses. Linares (2018) found the X-ray spectrum of this source is a hard powerlaw typical of redbacks, and that the X-ray and optical light curves are modulated on the same orbital period, confirming the source as a likely redback. The orbital period is one of the longest known among compact MSP binaries (orb ∼ 1.097 d). Here we present the first optical spectroscopy and multi-band optical photometry of PSR J1306– 40, and show that it has properties more similar to huntsman-type MSP binaries than classic redback systems. We detail our new optical spectroscopic and photometric observations of the system in §4.2. We model the optical light curves in §4.3, and examine the phase-resolved H profiles in §4.4. Finally, we make concluding remarks and discuss this system in the context of known redbacks and its connection to MSP–white dwarf progenitors in §4.5. 101 4.2 Optical Observations 4.2.1 SOAR Spectroscopy We obtained spectra of the companion to PSR J1306–40 using SOAR/Goodman (Clemens et al., 2004) from 2017 Jul 11 to 2018 Jun 9. All data used a 1200 l mm−1 grating and a 0.95” slit, giving a resolution of about 1.7 Å. Exposure times ranged from 15 to 30 min per spectrum, depending on weather conditions. The spectra covered a wavelength range of ∼ 5485–6740 Å. Barycentric radial velocities were determined through cross-correlation with bright standards taken with the same setup, primarily in the wavelength range 6050–6250 Å. The resulting 43 velocities are listed in Table 2.2. Observation epochs have been corrected to Barycentric Julian Date (BJD) on the Barycentric Dynamical Time system (Eastman et al., 2010). Even though the SUPERB survey has detected a pulsar in this system (Keane et al., 2018), it has not yet been timed, so the ephemerides must be determined from our data. Using the cus- tom Markov Chain Monte Carlo sampler TheJoker (Price-Whelan et al., 2017), we fit a circular Keplerian model to the data in Table 4.1 in order to determine the orbital period , BJD time of the ascending node 0, systemic velocity , and the semi-amplitude 2. Here, and throughout the paper, uncertainties are given at 1-sigma. We find  = 1.097195(161) d, 0 = 2457780.8373(19) d,  = 32.0 ± 1.8 km s−1, 2 = 210 ± 2 km s−1. (We note that the spectroscopic period we find here is fully consistent with the photometric period found in the discovery paper Linares (2018)). A fit using these median values is shown in Figure 4.1. For the remainder of this work, we assume the period derived from our spectroscopy. This fit has an rms scatter of 11 km s−1 (compared to a median uncertainty among the velocities of about 8.4 km s−1) and a 2/d.o.f of 74/39, perhaps suggesting the fit could be improved. How- ever, we see no clear trends in the residuals. A fit with the eccentricity left free does not find a value significantly different from 0, so there is no evidence that the orbit is non-circular. It may just be that the velocity uncertainties are slightly underestimated. It would be useful to obtain additional 102 radial velocities in the range  = 0.7− 1.0; it is challenging to get complete phase coverage for this system since its period is close to 1 d. Future pulsar timing should improve the ephemerides by orders of magnitude. Our SOAR spectra also show resolved H in emission in most, but not all epochs. We present an analysis of the H morphology in §4.4. Figure 4.1: Radial velocity curve of PSR J1306–40. Circular Keplerian fit to the SOAR/Goodman barycentric radial velocities of PSR J1306–40 listed in Table 4.1. 4.2.1.1 Determining the Mass Ratio By assuming the secondary is tidally synchronized with the pulsar, we can estimate its projected rotational velocity ( sin ) by comparing the rotationally broadened spectra with non-rotating template stars of similar spectral type. The assumption of synchronization is reasonable since this timescale is (cid:46) 1 Myr for a system with the period of PSR J1306–40 and a typical redback mass 103 Table 4.1. Barycentric Radial Velocities of PSR J1306–40 BJD (d) 2457945.5419639 2457945.6301728 2457956.5212405 2457956.5352468 2457956.5546942 2457956.5687017 2457956.5948590 2457956.6088705 2457966.4981876 2457966.5122271 2457966.5399208 2457966.5608956 2457966.5914866 2457967.5998999 2457996.4786139 2457996.4931540 2457996.5105239 2458001.4795192 2458001.4935388 2458139.7564130 2458139.7704090 2458140.7400266 2458140.7540222 2458140.7770050 2458140.7910007 2458161.6721983 2458161.6872977 2458202.5822163 2458202.5962046 2458202.8350917 2458202.8490793 2458223.5928781 2458223.8317433 2458223.8456941 2458243.6779211 2458243.6919415 2458243.7115353 2458243.7263475 2458243.7459698 2458247.5062024 2458247.5435496 2458278.5472619 2458278.5613912 RV (km s−1) –119.5 –16.9 –123.1 –107.9 –95.8 –65.6 –57.2 –38.0 –27.2 –18.4 38.7 49.6 115.1 –1.5 231.8 238.7 228.8 –144.0 –124.6 –136.8 –117.1 –194.5 –166.4 –150.8 –158.5 –121.2 –107.7 177.8 197.7 190.9 176.1 231.9 41.2 32.1 –87.6 –112.3 –125.5 –141.9 –130.3 131.9 160.5 197.1 183.7 err. (km s−1) 6.5 7.2 6.8 6.4 8.1 8.2 7.0 6.0 7.6 20.1 7.2 8.7 10.2 12.8 8.4 8.0 9.4 11.4 10.2 7.7 6.8 9.8 8.3 7.4 8.2 7.9 7.3 8.5 9.6 7.4 7.6 7.7 8.4 14.9 8.6 10.5 9.5 10.2 10.9 11.5 14.2 15.9 15.7 104 ratio (Zahn, 1977). Combined with our measurement of the orbital semi-amplitude, we use the  sin  value to constrain the mass ratio of the binary. Following similar procedures as described in Strader et al. (2014) and Swihart et al. (2017), we obtained spectra of bright late-G to mid-K giant stars to use as templates and convolved these with a set of rotational convolution kernels reflecting a range of  sin  values (including limb dark- ening). After cross-correlating the broadened templates with the original, unbroadened spectra, we fit a relation between the full-width at half-maximum (FWHM) values and the input value of  sin . We then cross-correlated the spectra of the companion to PSR J1306–40 with that of the unconvolved standard stars and used the FWHM values to estimate the projected rotational velocity. Our final estimate of  sin  derived in this manner is 75.5 ± 3.8 km s−1, where the uncertainty represents the standard deviation of the measurements among all templates and does not account for systematic uncertainties. Along with our measured value of the semi-amplitude, we use this equation, which uses the Roche lobe approximation of Eggleton (1983):  sin  = 2 0.49 2/3 (1 + ) 0.6 2/3 + ln(1 + 1/3) (4.1) where  = 2/  is the mass ratio. This equation is valid assuming that the secondary fills its Roche lobe and is synchronized. We find that  = 0.290±0.031. This measurement can be directly tested once a timing solution for the pulsar is available. 4.2.1.2 The Minimum Neutron Star Mass 2/2 =   sin3()/(1 + )2. We have determined The binary mass function  () = 3 all of these quantities from the spectroscopy alone, except for the inclination. Propagating the sin3() = uncertainties appropriately, we can determine the minimum neutron star mass   1.75 ± 0.09 (cid:12). This quantity is independent of the light curve modeling we carry out in §4.3 and is a robust lower limit. Hence we see that, consistent with many other redbacks (Strader et al., 2019), the mass of PSR J1306–40 is well in excess of the canonical 1.4 (cid:12). 105 4.2.2 Optical Photometry We obtained optical   photometry of the companion to PSR J1306–40 using ANDICAM on the SMARTS 1.3-m telescope at CTIO over 101 nights between 2017 Jul 17 and 2018 Apr 29 (UT). On each night, we took single 340 and 250 sec exposures in  and  bands, respectively, and two 250 sec exposures in  that were merged during the reduction process. Data were reduced following the procedures in Walter et al. (2012). We performed differential aperture photometry to obtain instrumental magnitudes of the target source using fifteen nearby comparison stars as a reference. Absolute calibration was done with respect to the Landolt (1992) standard fields TPheD (2017 Jul 17 – Sep 03) and RU149 (2018 Jan 24 – Apr 29). After excluding measurements with large errors, our final dataset includes 92, 93, and 95 measurements in , , and  bands, respectively. The mean observed magnitudes (not corrected for extinction) are  = 19.21,  = 18.35, and  = 17.20, with median errors of ∼ 0.03 mag in  and  and ∼ 0.06 in . The SMARTS photometry is listed in Table 4.2. To determine the orbital phase of these observations, we folded the data on the orbital period and ephemeris derived from our spectroscopic observations, where  = 0 corresponds to the ascending node of the pulsar (in this convention, the secondary is between the pulsar and Earth at =0.25). We show the folded light curves in Figure 4.2. The most obvious features of the light curves are the broad maxima centered near  ∼ 0.75 and narrower minima at  ∼ 0.25. A maximum at this phase suggests the tidally locked “day” side is being heated substantially, while the minimum corresponds to the “night” side of the secondary as we view it at inferior conjunction. As pointed out by Linares (2018), this heating plays a dominant role in shaping the optical light curve over the tidal deformation of the secondary (i.e., ellipsoidal modulations). If heating were not dominating the light curve, ellipsoidal modulations would result in two roughly equal maxima at  = 0 and  = 0.5, and two unequal minima at  = 0.25 and  = 0.75. The minimum at  = 0.75 would be dimmer due to the effects of gravity and limb darkening when viewing along the axis connecting the primary and secondary. Overall, the shape and amplitude 106 Table 4.2. SMARTS Photometry of PSR J1306–40 BJD (d) 2457962.46060 2457962.46769 2457962.47235 2457967.46273 2457967.46981 2457967.47448 2457968.46610 2457968.47319 2457968.47785 2457969.47164 · · · Band mag err B V I B V I B V I B · · · 19.104 18.212 17.007 19.357 18.631 17.461 19.334 18.426 17.289 19.140 · · · 0.104 0.029 0.062 0.099 0.044 0.095 0.073 0.037 0.033 0.062 · · · Note. — The full SMARTS dataset is avail- able in machine-readable format. We show a portion of the table here as a preview of its form and content. Magnitudes have not been corrected for extinction. of the light curves appear broadly consistent with the unfiltered Catalina Sky Survey (Drake et al., 2009) data presented by Linares (2018), suggesting no significant state change has occurred in this system since late-2005. Another notable feature of the light curves is that the maxima are not symmetric about the expected  = 0.75. Instead, the light curves slope downward between 0.5 <  < 1.0, suggesting that the source of heating is asymmetric with respect to the rotational axis of the secondary. Similar asymmetric heating has been observed in a number of other redbacks and black widows, and may be common in these systems due to the effects of an asymmetric intrabinary shock, heating mediated by the magnetic field of the secondary, or other magnetic activity such as star spots (e.g., Stappers et al., 2001; Breton et al., 2013; Schroeder & Halpern, 2014; Romani et al., 2015a,b; Deneva et al., 2016; Romani & Sanchez, 2016; van Staden & Antoniadis, 2016; Sanchez & Romani, 2017; Cho et al., 2018). We show below that asymmetric heating is a promising (though not unique) mechanism to explain the observed light curves of PSR J1306–40. 107 Figure 4.2: Optical light curve of PSR J1306–40. SMARTS optical photometry of PSR J1306– 40. Black lines show the best fit “No Spot” (dashed) and “Hot Spot” (solid) ELC models from Table 4.3. Two full orbital phase cycles are shown for clarity. 4.3 Light Curve Fitting We modeled the   light curves of PSR J1306–40 using the Eclipsing Light Curve (ELC; Orosz & Hauschildt, 2000b) code. Given the recent pulsar detection, we assume there is no ac- cretion disk; the light curves are dominated by a tidally distorted secondary that is being heated on its tidally locked day side. We also assume a circular orbit and fix the orbital period , semi- amplitude 2, and mass ratio  to the values derived from our spectroscopy (§4.2.1). These values immediately set the scale for the system, giving an orbital separation of ∼ 6 R(cid:12). In our most basic model, we fit for the orbital inclination , the Roche lobe filling factor 2, the intensity-weighted mean surface temperature of the secondary eff, the central isotropic irradiat- 108 0.00.51.01.52.0Phase17.017.518.018.519.019.5magBVI Table 4.3. Summary of ELC fits for PSR J1306–40 Parameters incl (◦) filling factor ( 2) eff day Ɗ  spot(◦)a spot(◦)b spot (K) spot(◦) 2(d.o.f.) 1 ((cid:12)) 2 ((cid:12)) log  (cgs) No Spot 70.4+16.7−7.0 0.88+0.06−0.02 4728+77−67 5542+275−165 -0.0354+0.0016 −0.0019 · · · · · · · · · · · · 370.3(275) 2.08+0.35−0.34 0.59 ± 0.10 3.745+0.025 −0.010 Hot Spot 83.5+3.3−4.8 0.94+0.04−0.05 4913+198−191 6187+397−451 2+5−2 27+14−19 11063+306−414 52+3−4 -0.0228 ± 0.0039 220.7(271) 1.77+0.07−0.03 0.51+0.02−0.01 3.725 ± 0.005 aSpot latitude. spot = 0◦ and spot = 180◦ represent the North and South pole, respectively bSpot longitude. spot = 0◦ and spot = 180◦ denote the inner Lagrange point and night side of the star, respectively ing luminosity from the pulsar, and a phase shift Ɗ . We characterize the irradiating luminosity indirectly through the most directly observed quantity: the maximum day side temperature of the heated secondary day. In Table 4.3, we summarize this model (“No Spot”) with the medians of the posterior distributions and corresponding 1 uncertainties for each parameter. We also plot this model along with the data in Figure 4.2 (dashed line). The best-fit inclination from these models is not all that well constrained ( = 70.4+16.7−7.0 ), but implies a massive neutron star (1 = 2.08+0.35−0.34 (cid:12)) that would be among the heaviest known. For the best-fit parameters, the formal 2/d.o.f. = 370/275, suggesting the fit can be improved. In particular, this model overestimates the -band flux at  ∼ 0.25 and underestimates the -band flux at this same phase, indicating the temperature profile of the night side is not well-matched. This model also does a poor job of reproducing the asymmetric, downward sloping curves between 0.5 <  < 1.0, as the only source of heating is symmetric around conjunction. 109 A promising solution for matching the shape of the downward sloping light curves around com- panion superior conjunction is by invoking an indirect, anisotropic heating model. The interaction between the pulsar and companion winds can form an intrabinary shock, generating asymmetric heating on the secondary as the spin-down power of the pulsar is reprocessed and indirectly illumi- nates the day side of the companion (Romani & Sanchez, 2016; Wadiasingh et al., 2017). In fact, a skewed intrabinary shock that slightly trails the companion during its orbit has been used (and is perhaps required) to explain the optical light curves in a number of black widow and redback MSP systems. (e.g. Romani et al., 2015a,b; Deneva et al., 2016; Cho et al., 2018). Such a shock could also be a natural location for producing H emission (see §4.4). In addition, strong magnetic fields on the companion’s surface may play an important role in the location of the anisotropic heating of the secondary. The rapid rotation and large temperature asymmetry of the tidally-locked companion can produce a strong dynamo, enhancing the surface magnetic fields (Morin, 2012). These magnetic field lines can duct the particles accelerated in the shock, channeling them directly to the magnetic poles and creating intense localized heating (i.e., one or more hot spots) (Tang et al., 2014; Sanchez & Romani, 2017). ELC does not allow us to directly model an offset intrabinary shock or a channeled pulsar wind. But we do have the ability to model its effect on the companion: we can add a hot spot to our underlying model of a heated, tidally distorted secondary. The hot spot is modeled as a circle with a temperature structure that falls off linearly towards the edges. We fit for the central temperature (spot), size (spot), and location (spot, spot) of the hot spot (see Table 4.3). Our main result is that adding a single hot spot to the heated, distorted secondary provides an excellent fit to the light curve (2/d.o.f. = 221/271). The reduced 2 in this model is < 1 likely due to the overestimation of a subset of the photometric uncertainties. In the best fitting hot spot model, the system is highly inclined ( ∼ 83◦), with a large hot spot near the north pole of the secondary. We show this best-fit model in Figure 4.2 (solid line) and summarize the posteriors in Table 4.3. Using our best hot spot model, the inferred mass of the primary is 1 = 1.77+0.07−0.03, fairly 110 massive for a typical neutron star, but fully consistent with the neutron star masses found in many redback MSPs (e.g., Kaplan et al., 2013; Romani et al., 2015b; Bellm et al., 2016; Strader et al., 2016; Sanchez & Romani, 2017; Shahbaz et al., 2017; Linares et al., 2018; Swihart et al., 2018). The secondary has a mass of 2 = 0.51+0.02−0.01, placing it in the high-mass tail of the redback companion mass distribution (Strader et al., 2019). The inferred gravity of the companion is log g = 3.725 ± 0.005 (cgs), suggesting it is slightly evolved off the main sequence, which may be relevant for interpreting the high energy emission from this system in the context of magnetically driven winds/outflows. For the hot spot model, we require a small phase shift to obtain adequate fits to the light curves. This is primarily due to the light curve minimum occurring ∼30 minutes earlier than what we expect from our spectroscopic ephemeris. This shift is larger than what we can account for from our formal uncertainties on the orbital period and 0, and the maxima appear roughly centered on the expected  = 0.75, so this feature of the light curves is likely to be real. One interpretation of this feature is that in addition to the hot spot, the intrabinary shock could be slightly trailing the companion, leading to off-center heating. For the hot spot itself, echoing our previous discussion, its best-fit location is centered near the companion’s north pole, consistent with magnetic ducting of intrabinary shock particles. The spot covers roughly ∼ 19% of the surface area of the star, but contributes about a third of the observed -band flux owing to its high temperature. This value is broadly consistent with the contribution from a similar hot spot modeled for the optical light curves of the redback binary 3FGL J0212.1+5320 (Shahbaz et al., 2017). 4.3.1 Distance Using the results from our light curve models, we estimate the distance to the binary by comparing its apparent magnitude to its intrinsic luminosity following similar procedures outlined in Strader et al. (2015) and Swihart et al. (2018). We estimate the bolometric luminosity by assuming eff 111 and the inferred radius of the secondary ( ∼ 1.6 (cid:12)) from our best-fit model (Table 4.3, column 3). We fit 10 Gyr solar metallicity isochrones (Marigo et al., 2008) to estimate bolometric correc- tions and used these to obtain the predicted absolute magnitude in each band ( ). Finally, we compared these values to the mean apparent magnitudes after applying extinction corrections using the Schlafly & Finkbeiner (2011a) reddening maps to get a distance to the binary: 4.7 ± 0.5 kpc. For the remainder of this paper, we adopt this value for the distance. We note that this method has proven successful in the past at estimating reliable distances to similar systems (e.g., Swihart et al., 2018). However, due to systematic effects such as the unknown metallicity and precise evolutionary state of the star, our estimate is likely uncertain by at least 20%. Future Gaia releases will help us address these systematics. PSR J1306–40 is listed in Gaia DR22 with parallax 0.106 ± 0.278 mas. Although the current parallax measurement is not significant, we estimate the geometric distance to the binary by using a weak distance prior based on an exponentially decreasing space density Galaxy model with a scale length of 0.94 kpc (Bailer-Jones et al., 2018). We also calculated the distance using a larger scale length (1.35 kpc), suggested by Astraatmadja & Bailer-Jones (2016). The resulting distances are 2.83+1.65−0.99 kpc and 3.47+2.37−1.36 kpc, respectively, somewhat smaller than our optically-inferred distance but within the uncertainties. Although the above estimates are uncertain, they are all wholly inconsistent with the dispersion measure based distance estimates made by using the measured dispersion measure and the Cordes & Lazio (2002, CL02) or Yao et al. (2017, Y17) electron density models (CL02∼1.2 kpc and Y17∼1.4 kpc, respectively). The discrepancies between the light curve/parallax distances and the dispersion measure-based estimates are consistent with the results of Jennings et al. (2018), who use Gaia DR2 parallaxes to measure the distances to a number of black widows and redbacks, and find that some dispersion measure distances, particularly those at high Galactic latitude, are underestimated. This also agrees 2https://www.cosmos.esa.int/web/gaia/dr2 112 with previous results for the distances of normal MSPs outside the Galactic Plane (Gaensler et al., 2008; Roberts, 2011). 4.4 H Spectroscopy A number of MSP binaries have shown strong H emission lines. In some of these systems, the existence of double-peaked Balmer emission indicates the presence of an accretion disk (e.g., de Martino et al., 2014b; Bogdanov et al., 2015), although an intrabinary shock and/or material streaming off the companion have also been used to explain the complex morphology (Sabbi et al., 2003; Swihart et al., 2018). Here we present the phase-resolved H profiles and suggest an intra- binary shock near the companion’s surface can explain the origin of the emission. Throughout most of our spectroscopic monitoring of PSR J1306–40, a moderate H emission line was persistent, occasionally showing complex, double-peaked morphology. At other times, the H appears only in absorption. We show the phase-resolved H profiles in Figure 4.3. Each spectrum has been corrected to the rest frame of the secondary (dotted line) and arbitrarily shifted in flux to display changes in the profile shape. When in emission, the H profile always peaks directly at, or very near, the velocity of the secondary. This suggests the emission is either coming directly from the star or from a region relatively close to the companion’s surface. Since the emission tracks the orbital motion, the line could arise from the stellar chromosphere as is the case in RS CVn systems (e.g., Drake, 2006). However, the emission we observe is broad, with a significant amount of the emission present at velocities (cid:38)200 km s-1 from the radial velocity of the companion. Furthermore, near companion inferior conjunction (=0.25), the emission line disappears and instead only appears in absorption. Since H absorption is typical in low-mass stars like the companion to PSR J1306–40 (Cram & Mullan, 1985; Pickles, 1998), it is likely that the emission region is being eclipsed at these phases and we are only seeing H absorption intrinsic to the star. Such a scenario can be explained by invoking an H emitting intrabinary shock near the com- 113 Figure 4.3: Phase-resolved optical spectra of PSR J1306–40 around the H region. Spectra have been corrected to the rest frame of the secondary (dashed line) and arbitrarily shifted in flux to display changes in the profile shape. The switch from emission to absorption near companion inferior conjunction ( = 0.25) is apparent, likely due to an eclipse of the emission region. 114 65306540655065606570658065906600Wavelength(˚A)Flux(Arbitrary)=0.097=0.099=0.114=0.121=0.227=0.252=0.271=0.299=0.534=0.615=0.633=0.902-1500-1000-500050010001500RadialVelocity(km/s) panion’s surface. Given our evidence for a rapidly-rotating Roche-lobe filling companion, it is likely the magnetic field on the surface is enhanced, which can lead to strong winds from the likely subgiant secondary. These winds can then be heated directly by the pulsar’s radiation pressure or through an interaction with the pulsar wind as it forms an intrabinary shock (see §4.3). Given the orbitally-modulated X-ray emission found by Linares (2018), such a shock almost certainly exists in this system and can provide a natural explanation for the observed X-ray and H phenomenology. 4.5 Discussion We have presented new optical photometry and spectroscopy of the source PSR J1306–40, a redback-like MSP binary with one of the longest known orbital periods among MSPs with non- degenerate companions. Our modeling suggests the companion is somewhat evolved and massive compared to typical redbacks, which may be relevant for interpreting the high-energy emission in the context of an intrabinary shock. As pointed out by Linares (2018), PSR J1306–40 lies in the error region of the Fermi-LAT -ray source 3FGL J1306.8–4031 (Acero et al., 2015), about 3.4’ from its center. In the preliminary LAT 8-year point source catalog (FL8Y3), which includes four additional years of survey data, the updated error region corresponding to the source FL8Y J1306.8–4035 is now ∼60% smaller in area than in 3FGL and has been shifted ∼3.2’ South, such that PSR J1306–40 lies only ∼0.86’ away from its center. Hence there is compelling evidence that PSR J1306–40 is indeed associated with a Fermi-LAT -ray source. Although the -ray spectra of many MSPs are typically highly curved, the relatively flat 3FGL spectrum shows no significant evidence for curvature with no cutoff up to (cid:38)10 GeV, reminiscent of huntsman MSP 1FGL J1417.7–4407 (Strader et al., 2015; Camilo et al., 2016), which has a red giant secondary in a 5.4 day orbit and also shows a flat -ray spectrum. Similar to PSR J1306–40, 1FGL J1417.7–4407 likely has a strong intrabinary shock and a significant wind from 3https://fermi.gsfc.nasa.gov/ssc/data/access/lat/fl8y/ 115 the companion that may be influencing -ray production (Swihart et al., 2018). Linares (2018) attributes the unusual -ray spectrum of 3FGL J1306.8–4031 to contamination from two nearby, active galaxies. We do not rule out this possibility, but given the improved source localization in FL8Y, contamination from other sources may not be important. The upcoming official 4FGL catalog will allow a reassessment of the spectrum of the -ray source. The X-ray light curve shows one maximum and one minimum per orbital period (Linares, 2018), which can be attributed to emission from an intrabinary shock that becomes at least partially eclipsed during inferior conjunction of the companion (e.g., Bogdanov et al., 2011; Rivera Sandoval et al., 2018). Although the minimum X-ray count rate presented by Linares (2018) approaches zero, there is no clear evidence for a total eclipse of the X-ray emitting region. At our best fit inclination ( = 83.5◦), an X-ray emitting point source would be completely eclipsed between  ∼ 0.21− 0.29, suggesting the shock is somewhat extended. Additional data will be needed to probe the precise geometry of the shock. Linares (2018) estimated the epoch of inferior conjunction of the companion from fitting the time of minimum of the X-ray light curve (0,X). Comparing this value with our spectroscopic ephemeris, we find the X-ray light curve lags behind the companion orbit by about one hour; the X-ray flux reaches a minimum at  ∼ 0.29, with a similar offset needed to match the X-ray maxima. However, 0, and our spectroscopic 0 are separated in time by 1471 d, approximately 1340 orbital periods. Assuming our uncertainty on the binary period (0.000161 d), building a phase-coherent solution for the radial velocity data backward to the epoch of the X-ray observation results in an absolute uncertainty of ∼12.9 hours, nearly half an orbital period. Therefore, the lag we find is not reliable with the current data. In the future, a precise pulsar timing solution would enable us to check the relative phase alignments between the X-ray and optical light curves. PSR J1306–40 has a long orbital period and a relatively massive secondary compared to the population of known redbacks in the Galactic field (Strader et al., 2019). In the binary evolu- tion models of Podsiadlowski et al. (2002), PSR J1306–40 lies above the bifurcation period that 116 distinguishes systems whose orbits will shrink to become black widows/redbacks and those that will grow to become MSP–white dwarf binaries. The light curves of PSR J1306–40 also show strong signs of irradiation, which can further contribute to an increase in the orbital period if the companion undergoes even low levels of evaporation (Chen et al., 2013). Our light curve models also suggest that the companion is filling a substantial fraction of its Roche lobe and that it has a somewhat evolved, subgiant-like radius, luminosity, and gravity (log g = 3.725 ± 0.005). Given that PSR J1306–40 has been spun up to obtain its rapid spin period, this would put PSR J1306–40 on a standard Case B evolutionary track of low-mass X-ray binaries that started the mass transfer process after leaving the main sequence, placing it in the late stages of the MSP recycling process that will terminate with a wide orbit MSP–white dwarf binary (Tauris & Savonije, 1999). The relatively long period and subgiant-like secondary of PSR 1306–40 resemble the huntsman systems 1FGL J1417.7–4407 and 2FGL J0846.0+2820 (Strader et al., 2015; Camilo et al., 2016; Swihart et al., 2017, 2018), which are likely progenitors to the typical MSP–He white dwarf bina- ries observed in the Galactic field. Another similarity between PSR 1306–40 and at least 1FGL J1417.7–4407 is the X-ray luminosity: Linares (2018) inferred a 0.5–10 keV X-ray luminosity of 8.8 × 1031 erg s-1 for PSR 1306–40 when using the dispersion measure-based distance of 1.2 kpc. But at our optical light curve-inferred distance estimate of 4.7 kpc, the X-ray luminosity of PSR 1306–40 is ∼1033 erg s-1, brighter than nearly all known redbacks in the pulsar state—except 1FGL J1417.7–4407. Detecting radio pulsations in these huntsman systems has proven difficult even compared to typical redbacks (Camilo et al., 2016; Keane et al., 2018). This difficulty may be related to the strong winds from the evolved companions, which when ionized could eclipse the radio emission more readily than for redbacks with main sequence-like companions. This high inferred eclipse fraction, and the rarity of these systems due to the shorter length of the red giant phase of evolution compared to the main sequence phase, suggests discovering other huntsman systems may be even more reliant on multiwavelength follow-up of unassociated 117 Fermi-LAT error regions than for typical redbacks. Fortunately, the companions in these systems are intrinsically brighter than redbacks, and so more readily observable in ground-based optical variability studies. A precise pulsar timing solution would be a significant step for fully understanding the huntsman PSR J1306–40 and its connection to known redbacks, placing tighter constraints on the orbital dynamics and permitting a search for spin and/or orbitally modulated -ray pulsations in the Fermi data. Additional deep X-ray observations would also be useful to compare with the most recent optical light curves and with detailed intrabinary shock models. The subgiant nature of the nearly Roche lobe filling, highly irradiated companion, and the clear evidence for intrabinary material makes this system a good candidate for future multiwavelength monitoring. 118 CHAPTER 5 A NEW LIKELY REDBACK MILLISECOND PULSAR BINARY WITH A MASSIVE NEUTRON STAR: 4FGL J2333.1–5527 1We present the discovery of a likely new redback millisecond pulsar binary associated with the Fermi -ray source 4FGL J2333.1–5527. Using optical photometric and spectroscopic observa- tions from the SOAR telescope, we identify a low-mass, main sequence-like companion in a 6.9-hr, highly inclined orbit around a suspected massive neutron star primary. Archival XMM-Newton X- ray observations show this system has a hard power-law spectrum Ɖ = 1.6± 0.3 and   ∼ 5× 1031 erg s−1, consistent with redback millisecond pulsar binaries. Our data suggest that for secondary masses typical of redbacks, the mass of the neutron star is likely well in excess of ∼ 1.4 (cid:12), but future timing of the radio pulsar is necessary to bolster this tentative conclusion. This work shows that a bevy of nearby compact binaries still await discovery, and that unusually massive neutron stars continue to be common in redbacks. 5.1 Introduction Prior to the launch of the Fermi Gamma-ray Space Telescope (Atwood et al., 2009) in 2008, most known millisecond pulsar (MSP) binaries had white dwarf companions in relatively wide or- bits (orb (cid:38) 5 d). Systems like these represent the end point of the MSP recycling process (Tauris & van den Heuvel, 2006), in which a stellar companion accretes matter and angular momentum onto the neutron star surface, spinning it up to rapid spin periods. However, recent multiwavelength follow-up observations of newly-discovered -ray sources from the Fermi Large Area Telescope (LAT) have revealed a substantial population of MSPs with non-degenerate companions. These systems generally display radio eclipses due to material from the secondary, making pulsations difficult to detect in blind radio timing searches. They are usu- 1This chapter is taken mostly word-for-word from Swihart et al. (2020), as published in the Astrophysical Journal 119 ally classified by the mass of their secondaries: black widows have very lightweight ( (cid:46) 0.05 (cid:12)), semi-degenerate companions that are highly ablated by the pulsar, while redbacks have more massive ( (cid:38) 0.1 (cid:12)) main sequence-like companions (Roberts, 2011; Strader et al., 2019) and typically show more extensive radio eclipses (e.g., Camilo et al., 2015; Cromartie et al., 2016; Keane et al., 2018). The redback class has drawn considerable attention in recent years due to three systems that have been observed to transition between an accretion-powered disk state and a rotationally-powered pulsar state on short ((cid:46) month) timescales, proving the suspected evolutionary link between low-mass X-ray binaries and MSPs (Archibald et al., 2009; Papitto et al., 2013; Bassa et al., 2014) and showing that the MSP recycling process in redbacks may be ongoing. Concerted multiwavelength observations of binary MSPs have enabled mass measurements of a number of neutron stars, placing tight constraints on the equation of state of dense matter and improving our understanding of fundamental physics at supranuclear densities (Steiner et al., 2013; Özel & Freire, 2016). Using 32 precise neutron star mass measurements, Antoniadis et al. (2016) found that the MSP mass distribution is strongly bimodal, with a narrow peak near the canonical ∼1.4 (cid:12) and a second, broad peak around ∼1.8 (cid:12). Intriguingly, the neutron star masses in the approximately two dozen known redbacks are consistent with being drawn almost exclusively from the more massive second peak of this proposed neutron star mass distribution (Strader et al., 2019). For MSP binaries with non-degenerate companions, the neutron star mass estimates are typi- cally dependent on the binary inclinations inferred from fitting the orbital variations in the optical light curves. When the companion is significantly heated by the pulsar (true for all black widows and many redbacks), this modeling is complex, and hence the inclinations can be plagued by sub- stantial systematic uncertainties. This leads, in turn, to conflicting mass estimates (e.g., Schroeder & Halpern, 2014; Romani et al., 2015b; Sanchez & Romani, 2017; Linares et al., 2018). However, in systems that are nearly edge-on, valuable constraints on the neutron star mass are possible that are practically independent of the less certain light curve modeling. In this paper, we present the discovery of the compact binary associated with the Fermi -ray 120 source 4FGL J2333.1–5527, showing that it is likely a redback MSP binary with a low-mass main- sequence companion in a highly inclined orbit around a massive neutron star. 5.2 Observations 5.2.1 The -ray Source & Optical Discovery 4FGL J2333.1–5527 is a bright unassociated Fermi-LAT -ray source with an overall detection significance of 19.4 in the 0.1–100 GeV energy range based on eight years of survey data (The Fermi-LAT collaboration, 2019). The source is persistent in -rays and has appeared in the previous 1FGL (Abdo et al., 2010), 2FGL (Nolan et al., 2012b), and 3FGL (Acero et al., 2015) catalogs. The source has a significantly curved -ray spectrum with no evidence for variability, consistent with most other -ray MSPs (Abdo et al., 2013a). The 4FGL 95% positional error ellipse lies entirely within the 3FGL region and is ∼70% smaller in area (Figure 5.1). The 0.1–100 GeV flux corresponds to a luminosity of 4.4× 1033 (/3.1 kpc)2 erg s−1. This distance is derived from the photometric light curve fitting in § 5.3.4.1 and is a bit more distant than the median redback distance (1.8 kpc; Strader et al. 2019). In any case, none of our central results depend on the exact distance, and given the faintness of the secondary, it is very unlikely this distance is off by more than a factor of ∼ 2. The subject of this paper is the strong candidate optical counterpart to 4FGL J2333.1–5527. We first identified this source as a possible variable star within the 4FGL error ellipse due to its unusually large Gaia  mag uncertainty compared to nearby sources of similar brightness. This source is listed in Gaia DR22 with a brightness  ∼ 20.3 mag and is located at a J2000 sexagesimal position of (R.A., Dec.) = (23:33:15.967, –55:26:21.06). For the remainder of this work, we refer to the optical and X-ray source (see below) as J2333. 2https://www.cosmos.esa.int/web/gaia/dr2 121 Figure 5.1: X-ray and optical finder images for 4FGL J2333.1–5527. Top: XMM-Newton/EPIC X-ray image of the field of 4FGL J2333.1–5527 with the 95% error regions from the 4FGL (red) and 3FGL (orange) catalogs overlaid. The subject of this work is marked with the blue circle. Bottom: Optical DSS image of the field with the location of the X-ray source marked in blue. The inset displays a zoomed-in view of one of our SOAR (cid:48) images showing the variable optical counterpart associated with the X-ray source. 122 3FGL4FGLNE3 arcminX-ray: XMM4FGL3FGL3 arcminOptical: DSSNEX-ray 5.2.2 X-ray Observations Prior to the launch of Fermi, the field containing 4FGL J2333.1–5527 was observed briefly and serendipitously with XMM-Newton on 2007 Nov 29 (ObsID 0505381001; PI Boehringer). We downloaded the publicly-available European Photon Imaging Camera (EPIC) observation from NASA’s High Energy Astrophysics Science Archive Research Centre (HEASARC) archive3. Data were reprocessed using standard tasks in the Science Analysis System (SAS) version 18.0.0 software package. We used standard flagging criteria FLAG=0, in addition to #XMMEA_EP and #XMMEA_EM for pn and MOS, respectively. Patterns 0-4 were used for pn and 0-12 for MOS. We found no evidence of strong, extended background flares, giving a total exposure time of ∼8.7 ksec. The optical variable discussed above matches the position of the brightest X-ray source within the error ellipse of 4FGL J2333.1–5527 (Figure 5.1), and this is the source we analyze. We extracted background-subtracted spectra and light curves using circular source extraction regions of radius 15” and local background regions three times larger. For our timing analysis, barycentric corrections were applied to all event files using the SAS task barycen. Individual MOS1, MOS2 and pn light curves were extracted using the SAS tasks evselect and epclccorr before being combined into a single EPIC light curve using the FTOOLS (Blackburn, 1995) task lcmath. Similarly, for our spectral analysis, individual MOS1, MOS2 and pn spectra were extracted using standard xmmselect tasks and then combined using epicspeccombine in order to create a weighted-average EPIC spectrum. Spectral fitting was performed using XSPEC (Arnaud, 1996) version 12.10.1 (§ 5.3.2). Due to the limited number of source counts, XSPEC’s implementation of Cash statistics (Cash, 1979) modified for a background-subtracted spectrum, W-statistics, was used to test the goodness of our spectral fits. 3https://heasarc.gsfc.nasa.gov/docs/archive.html 123 5.2.3 Optical Observations 5.2.3.1 SOAR Photometry We observed J2333 on 2018 Aug 27, Sep 19, and Oct 23 and 2019 Aug 19 and Sep 15 using SOAR/Goodman in imaging mode. On each of our 2018 nights we took a series of 180 sec ex- posures with the SDSS (cid:48) filter (Fukugita et al., 1996); during the 2019 nights we also included alternating exposures in (cid:48). Each image covered a circular 7.2’ diameter field and was binned 2×2 giving a final plate scale of 0.3”/pixel. Typical seeing was 0.8”, 1.4”, and 1.0” on the 2018 August, September, and October nights, respectively, and 1.1” and 1.0” on the 2019 August and September nights. The raw images were corrected for bias and then flat-fielded with a combination of dome and sky flats using standard packages in IRAF (Tody, 1986). We performed aperture photometry to obtain the instrumental magnitudes of the target, which were then calibrated using the (cid:48)- and (cid:48)-band magnitudes of a number of nearby comparison stars taken from The Blanco Cosmology Survey (Desai et al., 2012). To ensure that any variability we observed from our target was real, we examined the light curves of our comparison stars to choose only the most stable stars, finding 14 reference stars in (cid:48) and 24 in (cid:48). Images that were affected by clouds or that were taken under very poor seeing conditions ((cid:38)1.8”) were removed. The final photometric sample consists of 121 and 272 measurements in (cid:48) and (cid:48), respectively. The mean observed magnitudes of the target source (not corrected for extinction) are (cid:48) = 21.123 mag and (cid:48) = 19.868 mag, with median per epoch errors of ∼0.039 mag and ∼0.013 mag in (cid:48) and (cid:48), respectively. For the remainder of the paper, we refer to the bright phases of the optical light curve as the “dayside”, corresponding to when we observe the inner heated face of the tidally locked companion, while the “nightside” corresponds the companion’s non-illuminated backside (see § 5.3.4). 124 5.2.3.2 SOAR Spectroscopy We obtained spectra of J2333 using the Goodman spectrograph (Clemens et al., 2004) on the SOAR telescope on Cerro Pachon, Chile from 2018 Aug 20 to 2019 Aug 6. All data used a 400 l mm−1 grating with an approximate wavelength range of either ∼ 4800 to 8800 Å or a setup about 1000 Å bluer. Depending on the seeing, we used slit widths of 0.95” or 1.2” yielding resolutions of 5.7 or 7.3 Å FWHM, respectively. Owing to the faintness of the star, the exposure times were relatively long, 25 or 30 min per spectrum. The spectra were reduced and optimally extracted in the normal manner. Radial velocities were measured through cross-correlation with bright templates of similar spectral type around the Mg region of the spectra. 42 spectra had high enough signal-to-noise for velocities to be measurable. The mid-exposure observation times are reported as Modified Barycentric Julian Dates (BJD - 2400000.5 d) on the TDB system (Eastman et al., 2010). The nightside spectra are consistent with having spectral types as low as mid-K, while most of the dayside spectra are best matched as early-G stars (Figure 5.2). One or two spectra appear even warmer, into the range of an early-F star with weak metal lines, perhaps reflecting a transient episode of increased heating. We show below that the mean colors of the star on its night and day sides are fully consistent with the typical spectra observed. Nearly all of the spectra cover both H and H, and while there is Balmer absorption present in most of the warmer spectra, there is no clear evidence for emission in any of the spectra. This is unlike some other redbacks, where the Balmer emission is attributed to a wind or shock (Swihart et al., 2018). In some redbacks that show strong evidence for irradiation, the radial velocities depends on the absorption lines used to measure them. For example, Balmer lines will tend to trace the brighter (warmer) portions of the companion, which are closer to the binary center of mass and therefore move slower, whereas metallic lines track the dark (cooler) side, which is moving faster (e.g., 125 Figure 5.2: Optical spectra of 4FGL J2333.1–5527. Two of our SOAR spectra showing examples of the varying effective temperature of the companion. The dayside spectrum (top) is consistent with an early-G star, while the nightside spectrum (bottom) is markedly cooler, matching that of a mid-K dwarf. Our other low-resolution spectra typically lie somewhere between these extremes. Linares et al., 2018). For the few spectra that showed clear H absorption, we compared the radial velocities of these lines to the Mg velocities, but did not find evidence for a clear offset between these velocities at any phase. 5.3 Results & Analysis 5.3.1 Modeling the Radial Velocities We fit Keplerian models to the radial velocities using the Monte Carlo sampler TheJoker (Price- Initially we fit a four-parameter circular model, finding period  = Whelan et al., 2017). 126 0.2876450(14) d, epoch of ascending node 0 = 58463.46881(65) d, semi-amplitude 2 = 360±5 km s−1, and  = 44 ± 4 km s−1. These parameters represent a good fit, with 2 = 36.4 for 38 d.o.f. and an rms of 20.5 km s−1. This model is plotted with the original velocity data in Figure 5.3 (here and throughout this paper, we use the orbital phase convention where  = 0.25 is the inferior conjunction of the companion). However, this model of the data may not be complete: the substantial heating of the secondary implied by the light curves (see § 5.3.4 and Figure 5.6) suggests that heating could also be affecting the measured radial velocities of the secondary. Such heating can displace the center of light outward from the center of mass, inflating the measured velocity semi-amplitude or causing a false eccentricity (Davey & Smith, 1992). There is no statistical support for an eccentric model: again fitting with TheJoker, the posterior for the eccentricity is peaked at 0, with a median posterior value of  = 0.009. Using median values of this fit reduces the 2 value by only 0.2 despite two additional free parameters. We also performed circular fits for relevant subsets of the data (most/least heated phases, or closest to quadrature), but found this did not substantially affect the measured 2 within its uncertainty. Hence there is little direct evidence from the velocities themselves that they are affected by irradiation. Nonetheless, to make a first-order estimate of the effect heating might have on our velocities, we use the ELC models calculated to match the light curve in § 5.3.4.1. At each phase of the model, we take the intensity-weighted velocity offset of each grid element from the center of mass as viewed by an observer, and combine it with an effective temperature–equivalent width relation for Mg taken from Johansson et al. (2010). This method is similar to that used by Shahbaz et al. (2017) except that we do not set the equivalent widths of the elements heated beyond a certain point to zero, since this is not demanded by the data. Future spectroscopy with a larger telescope would allow improvements in the quality of the spectra that could show the need for additional modeling. Table 5.1 shows the original velocities as well as those with this correction applied. The strength of the correction is highest near quadrature. Refitting a circular model to the corrected velocities 127 Figure 5.3: Radial velocity curve of 4FGL J2333.1–5527. Circular Keplerian fit to the observed radial velocities. As discussed in the text, the fit to the “corrected" velocities has a lower semi- amplitude but is otherwise extremely similar. gives parameters identical to the original fit within the uncertainties, excepting the semi-amplitude, which is lower at 2 = 338(4) km s−1. The overall fit is of slightly lower quality (2 = 38.2 for 38 d.o.f.; rms of 20.8 km s−1), and indeed there is no particular evidence from the data that this model provides a better fit than the original one. In subsequent sections we note results from both the original and “corrected" velocities to indicate how the conclusions depend on these effects. 5.3.1.1 Minimum Neutron Star Mass Using the posterior samples from our Keplerian fits to the original radial velocities, we derive the 2 (1 − 2)3/2/(2) =   sin3/(1 + )2 = 1.39 ± 0.05 (cid:12), for mass function  () =  3 mass ratio  = 2/  and inclination . For the heating-corrected velocities, the 2 is lower by 128 Table 5.1. SOAR Velocities of J2333 RVcorr (km s−1) 181.8 386.4 397.9 340.9 238.4 –295.9 –175.1 –106.2 70.0 316.6 253.8 109.0 –85.7 –208.3 –217.0 –260.7 –263.1 144.7 22.5 –90.8 –309.7 222.3 71.0 –248.0 –259.6 325.8 295.2 369.9 364.5 337.6 281.2 158.1 –241.1 –269.9 –257.5 102.2 363.8 251.1 –46.3 –186.1 –188.3 –288.2 err (km s−1) 18.2 19.3 20.1 19.1 18.1 23.5 21.5 18.3 21.9 21.8 20.8 25.1 31.2 24.5 28.4 20.1 18.2 21.3 20.9 25.3 21.4 16.5 22.0 20.2 29.1 25.8 27.5 24.0 17.8 16.2 22.8 18.9 23.3 16.5 19.0 22.2 19.8 21.0 20.2 19.0 20.4 23.1 BMJD (d) 58345.3394130 58367.2301482 58367.2476498 58367.2763716 58367.2940899 58368.2446582 58368.2942262 58368.3121900 58368.3340083 58372.1703291 58372.1879465 58372.2087764 58372.2263937 58372.3114531 58380.3103683 58380.3278272 58462.0531349 58463.0939823 58463.1114704 58463.1327180 58463.1701621 58484.0827981 58484.1002837 58491.0546756 58491.1070246 58503.0410349 58637.3104413 58637.3279199 58637.3499781 58637.3674571 58637.3893140 58637.4067921 58664.2250896 58664.2426674 58664.2703763 58664.3298441 58664.3956801 58665.2948716 58665.3356255 58665.3605401 58701.3098865 58701.3323067 RV (km s−1) 184.9 401.0 421.3 361.2 250.2 –320.0 –182.1 –108.9 70.6 335.1 263.8 110.7 –90.3 –220.2 –235.3 –285.2 –277.1 150.3 21.7 –99.8 –334.3 231.8 73.6 –266.7 –272.6 347.6 300.8 384.0 388.9 360.7 294.6 163.7 –263.3 –294.7 –269.8 104.2 388.8 262.6 –50.6 –201.5 –200.6 –310.8 129 about 6%, leading to a mass function of  () = 1.15 ± 0.04 (cid:12). Both values are large compared to typical redbacks, suggesting a relatively edge-on orientation, and in the former case a massive neutron star primary. Assuming the original velocities and the median mass ratio of known red- backs ( = 0.16), the minimum mass of the presumed neutron star for an edge-on ( = 90◦) orbit is   = 1.87 ± 0.07 (cid:12). The mass ratio is not well-constrained from these data and hence the detection and timing of the neutron star as a pulsar will be necessary to improve our estimate of the primary mass. If we instead assume the minimum possible mass ratio ( = 0), the minimum neutron star mass is   = 1.39±0.05 (cid:12) for an edge-on orbit, with smaller inclinations resulting in a heavier neutron star. We revisit the mass estimate in the context of the light curve fitting in § 5.3.4.1. 5.3.2 X-ray Spectrum We initially fit the spectrum of J2333 with simple absorbed power-law (TBabs*powerlaw) and blackbody (TBabs*bbody) models. It is immediately apparent that the spectrum is too broad to be fit with a single blackbody (cstat/dof = 171/217), and the absorbed power-law provides an appreciably better fit (cstat/dof = 132/217; Figure 5.4). Using the power-law model, we find no evidence of additional intrinsic absorption and thus fix the TBabs parameter to the line-of-sight absorption column density H = 1.1 × 1020 cm−2 (Dickey & Lockman, 1990; Kalberla et al., 2005; HI4PI Collaboration et al., 2016). The best-fitting photon index is Ɖ = 1.6 ± 0.3, with an unabsorbed 0.3 − 10 keV flux of X = (cid:17) × 10−14 erg s−1 cm−2. This equates to a luminosity of (cid:16) (cid:17) × 1031 erg s−1, using the reference distance of 3.1 kpc. 4.3+1.4−1.1 (cid:16) 5.0+1.5−1.3 We also tried a more complex TBabs*(powerlaw + bbody) model in an attempt to account for the fit residuals found at energies (cid:38) 5 keV (Figure 5.4). A number of binary MSPs have been shown to have evidence for a thermal component in their spectrum (e.g., Heinke et al., 2009; Archibald et al., 2010; Bogdanov et al., 2011; Roberts et al., 2014, 2015) likely coming from the heated magnetic polar caps on the surface of the neutron star. The introduction of a soft 130 ( = 0.09±0.05 keV) blackbody component causes the power-law to become flatter (Ɖ = 1.1±0.5). In this case, we find a slightly higher total unabsorbed 0.3 − 10 keV flux of s−1 cm−2, of which the power-law contributes ≈ 90%. However, this model (cstat/dof = 127/215) is comparable statistically to the pure absorbed power-law model and thus we prefer the simpler model. 5.7+2.3−1.7 (cid:16) (cid:17) × 10−14 erg Figure 5.4: X-ray spectrum of 4FGL J2333.1–5527. XMM-Newton/EPIC MOS+pn spectrum of J2333, with model and model/data ratios. Due to the low number of source counts, the unbinned data were fitted with W-statistics and later rebinned to a minimum significance of 2.2 for visualization purposes only. The model (red) shows our best-fit absorbed power-law model (TBabs*powerlaw) with photon index Ɖ = 1.6 ± 0.3. 5.3.3 X-ray Light Curve The background-subtracted X-ray light curve of J2333 is shown in Figure 5.5. After applying barycentric corrections, the data were grouped into 150s and 500s time bins. The 8.7 ksec exposure covers approximately 40% of the orbit from  ∼ 0.65 − 1.05, and shows considerable variability on short timescales. The black dashed line in Figure 5.5 represents the weighted average count rate of the 0.5 ks binned dataset and is shown for visualization purposes only. 131 10−510−410−30.01normalized counts s−1 keV−110.5250246ratioEnergy (keV) Although the 0.2–10 keV count rates vary by at least a factor of 3 over this short interval, a 2 fit suggests a probability of ∼3% that this data could be produced from a constant flux distri- bution. Together with the hard X-ray spectrum presented above, this emission can be attributed to an intrabinary shock that occurs at the interface between the wind driven from the companion and the relativistic pulsar wind (e.g., Gentile et al., 2014; Romani & Sanchez, 2016; Al Noori et al., 2018). A longer observation would allow us to determine if this variability were real and possibly modulated on the orbital period of the binary. Consistent with the low flux, this source was not detected in a short Swift/XRT observation performed as a part of a program monitoring unidentified Fermi sources (Stroh & Falcone, 2013). This highlights the need for deeper X-ray observations of other unassociated Fermi -ray sources. Figure 5.5: X-ray light curve of 4FGL J2333.1–5527. The XMM-Newton/EPIC light curve between 0.2–10 keV. The observation covers ∼40% of the binary orbit. Barycenter corrections have been applied to the photon arrival times and the data are grouped into 150s (black) and 500s (blue) bins. The weighted average of the 500s binned dataset is shown with a dashed line. 132 0.650.750.850.951.05BinaryPhase024681012Timesince2007Nov2917:48:54.816UTC(ks)0.000.020.040.060.080.100.12XMMEPICCountRate(cts/s)0.15ks0.5ks 5.3.4 Optical Light Curve The optical light curves show clear evidence of variability, with (cid:48) and (cid:48) magnitudes ranging from ∼20.4–21.8 and∼19.55–20.35 mag, respectively. Our 2019 Sep data covers slightly more than a full orbital cycle, suggesting the period is approximately 6.9 hours. In order to confirm this estimate, we performed a period search on the full photometric dataset (five epochs) using a Lomb-Scargle periodogram. The highest power peak was found at a period of orb = 0.28764467 d ≈ 6.90347 hrs. Since this value is fully consistent with the spectroscopic period, for the remainder of the paper we assume the orbital period and ephemerides (and corresponding uncertainties) derived from the spectroscopy (§ 5.3.1). We calculated the midpoint of each observation and then folded these data on our best fit orbital period and 0. Consistent with our spectroscopic convention,  = 0 represents the ascending node of the suspected pulsar, such that the secondary lies along the line connecting the primary and Earth at  = 0.25. The folded light curves are displayed in Figure 5.6. Similar to a number of other redback MSPs, the light curves show a broad maximum and narrower minimum per orbital period, consistent with the secondary being heated on its tidally locked “day” side by the suspected pulsar primary. This heating is dominating the shape of the light curve over ellipsoidal modulations due to tidal deformation of the secondary, which would manifest in two maxima/minima per orbital cycle. In redbacks, heating patterns like these typically occur when the pulsar spin-down power is reprocessed in an intrabinary shock formed from the interaction between the pulsar and companion winds (Romani & Sanchez, 2016; Wadiasingh et al., 2017; Kandel et al., 2019). Such a shock would also produce a characteristic hard power-law X-ray spectrum (e.g., Roberts et al., 2015; Werner et al., 2016), consistent with our observations (§ 5.3.2). Overall, the light curves are fairly well-behaved with a few notable exceptions. First, in (cid:48) the scatter in the data appears more pronounced after  ∼ 0.65, suggesting increased intrinsic stochastic 133 variability when viewing the inner heated face of the companion. This increased variability also appears in (cid:48), so that the observed maximum appears slightly later than the expected  = 0.75, while the shape and slope of the (cid:48) data around  = 0.85 − 0.95 suggests the light curve is not well-behaved at these phases. Lastly, the scatter in the (cid:48) data when viewing the companion’s nightside ( ∼ 0.1 − 0.4) is appreciably larger than at any other phases, with the brightness varying by nearly 0.2 mag over timescales of a few minutes. This phenomenology suggests that there is variable heating ongoing on the companion surface similar to what has been observed in a number of other redback MSPs (e.g., Cho et al., 2018). This variability could be explained by a number of physical processes such as magnetic reconnection events between the intrabinary shock and stellar surface, a rapidly varying shock geometry, inhomogeneities, asymmetries, and/or magnetic channeling of the pulsar wind, or migrating star spots on the companion (e.g., Deneva et al., 2016; van Staden & Antoniadis, 2016; Halpern et al., 2017; Sanchez & Romani, 2017; Zharikov et al., 2019). Together with our observation in §5.2.3.2 that a few optical spectra have noticeably higher effective temperatures than the others, this variability reinforces the idea that this system undergoes irregular transient heating episodes that affect the tidally locked companion. 5.3.4.1 Modeling the Light Curves Armed with the SOAR (cid:48) and (cid:48) light curves, we model the system using the Eclipsing Light Curve (ELC; Orosz & Hauschildt, 2000b) code. Throughout our modeling we assume a compact, invisible primary with no accretion disk, and a tidally distorted secondary in a circular orbit. We hold the orbital period  and semi-amplitude 2 to the values found from our best radial velocity curve fits (§ 5.3.1). We also assume the secondary is tidally locked, a valid assumption since the synchronization timescale is (cid:46) 104 yr at the orbital period of J2333 (Zahn, 1977). Since we do not yet have an estimate for the mass ratio of the binary, we assume the median value for the known redbacks:  = 2/NS = 0.16 (Strader et al., 2019). Leaving the mass ratio 134 Figure 5.6: Optical light curve of 4FGL J2333.1–5527. SOAR (cid:48) (bottom) and (cid:48) (top) photometry of the companion to J2333 along with the best fit ELC model. The data described in § 5.2.3.1 have been folded on the best-fit period and ephemeris from § 5.3.1, and the uncertainties have been rescaled as described in the text (§ 5.3.4.1). Two orbital phases are shown for clarity. . as a free parameter results in a wide range of well-fitting models mainly due to the degeneracy between this parameter and the filling factor of the secondary (i.e., a smaller star can accommodate smaller mass ratios and vice versa). Future observations, such as high-resolution spectroscopy to measure the companion’s rotational broadening, or a measurement of the projected semi-major axis of the pulsar orbit would allow for a direct determination of the mass ratio and break this degeneracy. The overall scale of the system is already set by , 2, and , giving an orbital separation for the secondary of ∼2.3 (cid:12). Finally, the controlling physical parameters of our light curve models 135 0.00.51.01.52.0Phase(P=0.287644d)22.021.521.020.520.019.5magi´g´ are the effective temperature of the companion’s nightside 2, the isotropic irradiating luminosity of the suspected pulsar (which we characterize here as the maximum temperature on the day side of the heated secondary) day, the Roche lobe filling factor of the companion 2, the binary inclination , and a phase shift Ɗ . This last parameter can account for any small offsets that might be present between the light and radial velocity curves, which are commonly found in black widow or redback binaries dominated by irradiation (e.g., Schroeder & Halpern, 2014; Draghis et al., 2019; Swihart et al., 2019). As mentioned in § 5.2.3.2, we analyzed a number of our low/medium resolution spectra of J2333 that were taken when we were viewing the nightside of the secondary ( = 0.1 − 0.4) to constrain 2. These spectra are consistent with a ∼K5 main sequence star with an effective temperature of ∼4400 K. To confirm this value, we also examined the nightside colors of the companion using our SOAR photometry. After applying extinction corrections using the Schlafly & Finkbeiner (2011a) red- dening maps, the de-reddened color at  ∼ 0.25 is ((cid:48) − (cid:48))0 = 1.428. We compare this value to the theoretical colors of main sequence stars using a 10 Gyr solar metallicity isochrone (Bressan et al., 2012), finding 2 ∼ 4470K. If we instead assume a more metal-poor star with [Fe/H] = –1, we find a lower 2 ∼ 4350K. These values are broadly consistent with our estimate from the spectral type so for the remainder of the analysis, we fix the nightside effective temperature of the secondary to 2 = 4400K. We initially found models that appeared superficially to provide good fits to the data but had a large formal 2/dof = 1702/389, mostly dominated by the small errors on the (cid:48) photometry. Furthermore, there are more than twice as many (cid:48)-band measurements as there are in (cid:48), so when modeling the light curves together ELC preferentially finds models that fit the (cid:48) data well, while the (cid:48) data receives very little weight. In addition, a small fraction of data points have substantially smaller uncertainties than the typical data at that phase, which can dominate the model fitting, essentially forcing the model to go thru (“overfit”) these small error datapoints. This is especially 136 problematic when trying to fit the nightside in (cid:48) and the (cid:48) data near maximum, where the increased scatter is most pronounced and only a couple datapoints with small errors can dominate the model at these phases. Therefore, we rescaled the uncertainties so that all intraband errors were equal and adjusted the values so that the total reduced 2 of the final model was ∼ 1.0. This process ensured the photometry in both bands contributed approximately equally to the total 2 while also insulating against “overfitting” datapoints that had substantially smaller uncertainties than the typical data at that phase. This method has been used in the literature for modeling the light curves of redbacks and achieving fits with reliable uncertainties (e.g., McConnell et al., 2015; Linares et al., 2018; Strader et al., 2019). The best fit ELC model along with the photometry (with rescaled uncertainties) is shown in Fig- ure 5.6. The best fit model suggests the star is only partially filling its Roche lobe ( 2 = 0.75±0.01), has a dayside temperature of ∼5700 K (consistent with our optical spectra; § 5.2.3.2), and is in a nearly edge-on orbit ( ∼ 86◦). Using these values together with our spectroscopic results, we infer the mass of the primary associated with a range of reasonable secondary masses. The distribution of redback companion masses can be modeled with a normal distribution: 2 = 0.36±0.16 (Strader et al., 2019). Utilizing this distribution along with our measurement of the binary mass function  () = 1.39 ± 0.05 (cid:12) (§ 5.3.1.1), the primary mass is   = 1.96+0.25−0.27 (cid:12). In these models the secondary has a radius of 2 ∼ 0.51 (cid:12). If instead we use the smaller mass function corresponding to our “heating- corrected” radial velocity model, the primary mass is   = 1.70+0.23−0.25 (cid:12), though as we discuss in § 5.3.1 there is no evidence that this “heating-corrected” 2 fits the data better than the original 2 value. Regardless, both these estimates for the primary mass suggest the presumed neutron star is well in excess of the canonical ∼ 1.4 (cid:12). The inclination is not precisely constrained owing to the slight degeneracies between this parameter, the filling factor, and the level of heating. We explored models where 2 (and therefore 137 also the relative heating on the dayside) was allowed to vary within the constraints set by our optical spectra, but found that models that both fit the data well and were consistent with our spectra fully agree with our main results within a few percent. The best fit value for the inclination and 1 uncertainty is  = 85.◦8+4.2−19.1. Since the currently assumed value is close to edge-on, allowing a wider range of inclinations would allow a more massive neutron star than discussed above. That said, using the original 2 measurement and the median secondary mass for known redbacks, models with  (cid:46) 74◦ imply neutron star masses larger than the current Shapiro delay record holder (2.14+0.10−0.09(cid:12); Cromartie et al. 2020) and hence are less likely. Lower values of 2 can accommodate a wider range of inclinations. 5.3.4.2 Distance Similar to the steps described in Strader et al. (2015) and Swihart et al. (2019), we use our best light curve models to infer the distance to the system. We first derive the intrinsic luminosity of the system by assuming 2 and 2 from our best fit model, then apply bolometric corrections in each band as a function of temperature using the same solar metallicity 10 Gyr isochrone listed above. Finally, we compare these predicted absolute magnitudes to the mean de-reddened apparent magnitudes in each band to infer a distance of 3.1 ± 0.3 kpc. The uncertainty in this estimate includes systematic effects due to the unknown metallicity and Roche lobe filling factor of the star as well as the dispersion between filters. A future Gaia data release should provide a parallax distance for this source to test the pho- tometric estimate; for the handful of sources with accurate parallax distance measurements, these photometric light curve distances appear to give reasonable results (Jennings et al., 2018; Strader et al., 2019). 138 5.4 Conclusions We have presented the discovery of the suspected optical and X-ray counterparts to the unas- sociated Fermi -ray source 4FGL J2333.1–5527, showing that it is likely a redback MSP binary harboring a massive neutron star with a heated companion in a nearly edge-on orbit. Our Keplerian fit to the radial velocity data suggests a primary mass well in excess of the canonical ∼ 1.4 (cid:12), although tighter constraints on the binary mass ratio, from high-resolution optical spectroscopy to determine the projected rotational velocity of the secondary or future timing of the radio pulsar, are necessary to bolster this tentative conclusion. If the radio pulsar is found and timed, the likely edge-on geometry of this system would be ideal for measuring the relativistic Shapiro delay to precisely measure a massive neutron star (e.g., Demorest et al., 2010; Antoniadis et al., 2013; Cromartie et al., 2020), though unfortunately this is likely to be challenging due to radio eclipses. The hard power-law X-ray spectrum in J2333 is similar to that of other recently discovered redback MSPs and is consistent with emission from an intrabinary shock. The inferred X-ray luminosity of ∼ 5 × 1031 erg s−1 is also in line with known redbacks in the pulsar state, even considering the uncertain distance (Linares, 2014; Strader et al., 2019). Future X-ray observations can better constrain the properties of the shock, especially given the reasonably well-understood inclination of the binary. Although our optical and X-ray data provide compelling evidence that this object is a redback radio pulsar, the most significant advance towards fully understanding and characterizing the binary in 4FGL J2333.1–5527 would be the detection of radio pulsations at the position of the optical/X- ray source, followed by a precise pulsar timing solution. This would confirm our interpretation of the source and would place important constraints on the size of the orbit and mass ratio of the binary that are only weakly constrained at present by our modeling of the light and radial velocity curves. Such a detection would also enable a search for high energy pulsations modulated at the spin period of the binary as has been found in a number of other redback MSPs (e.g., Johnson et al., 2015; Smith et al., 2017). In the absence of radio pulsations, given the brightness of the source in 139 -rays, a brute force -ray pulsation search could also be performed in order to definitively link the optical/X-ray object to the -ray source. The unveiling of a new redback system among the brighter unassociated 4FGL sources suggests an exciting continuing discovery space for compact binaries with Fermi in the months and years to come. 140 CHAPTER 6 CONCLUSION AND FUTURE WORK In this dissertation, we presented the results obtained by analyzing multiwavelength observations of confirmed and candidate field MSP binaries among unassociated Fermi sources. Radio, X-ray, optical, and infrared observations were used to identify and classify four MSP binaries with non- degenerate companions in the late stages of MSP recycling. In Chapter 2 and Chapter 3, we presented the first multiwavelength observations of a new subclass of spider neutron star binaries: MSPs with giant companions in wide orbits, aptly named “huntsmen” after the large Australian spider. These systems are the likely progenitors to the canon- ical MSP binaries, MSPs with white dwarf companions ( (cid:38)3 d), which represent the bulk of the MSP population in the Galactic field. In Chapter 3 we also showed that double-peaked H emission lines, assumed in many previous works as the classic hallmark of an accretion disk, can actually occur as the result of a strong, magnetically driven stellar wind from the giant secondary and its interaction with the pulsar wind. As hot gas from the companion flows towards the inner Lagrange point, it is swept out beyond the companion’s orbit by the pulsar’s radiation pressure and the intrabinary shock. This causes the emission region to slightly trail the orbit of the companion, resulting in the characteristic emission profiles seen in Figures 3.6 and 3.7. This extended material from the secondary can also explain the difficulty in detecting radio pulsations in these huntsmen systems, even when compared to typical redbacks. In Chapter 4, we used newly obtained optical photometry and spectroscopy of the MSP binary PSR J1306–40 to estimate the mass of the MSP and its companion, finding the neutron star is mas- sive (1.75 ± 0.09 (cid:12)) and its companion evolved off the main sequence. The nearly edge-on orbit makes this system ideal for measuring the relative masses of the two components if the pulsar can be precisely timed. However, this will be challenging given the (sub)giant nature of the companion and its strong wind, which likely eclipses the radio pulsar at irregular intervals. This system also 141 highlights the need for precise photometric modeling when determining the binary inclination in irradiated compact binaries, since without the inclusion of starspots in the model, the light curve fits are inadequate (Figure 4.2). Similar to 2FGL J0846.0+2820 and 1FGL J1417.7–4407, this system is likely in the late stages of MSP recycling that will terminate with a wide orbit MSP–white dwarf binary, so may also be included in the huntsmen subclass. Lastly, in Chapter 5 we brought attention to a newly discovered compact binary, 4FGL J2333.1– 5527, that resembles a “normal” redback more so than the three systems described in the previous chapters. Although the radio pulsar has not yet been found, optical photometric and spectroscopic observations confirm the presence of a low-mass main-sequence companions in a a 6.9-hr, highly inclined orbit around a suspected massive neutron star primary that is likely well in excess of ∼ 1.4 (cid:12). The companion shows a distinct optical light curve dominated by irradiative heating of the tidally-locked inner face, and archival XMM-Newton data show a highly variable X-ray light curve, consistent with an intrabinary shock between the pulsar and the companion. Even though this system is a persistent -ray source (present in the first Fermi-LAT catalog release and in each subsequent catalog), and was among the brighter unassociated Fermi sources in the 4FGL catalog, it went undiscovered until this work, and therefore highlights the need for continued multiwavelength follow-up of unidentified Fermi sources in order to find new neutron star binaries. 6.1 Future Prospects There are only about two dozen known MSPs with main-sequence or giant companions (Fig- ure 1.5), over half of which have been discovered since 2017. One intriguing aspect of this popu- lation of MSP binaries is that most of these systems harbor massive neutron stars when compared to the bulk of the neutron star mass distribution. In fact, there is increasing evidence that neutron star masses are bimodal (Özel & Freire, 2016), with a narrow peak at the canonical ∼ 1.4 (cid:12) and a broad secondary peak at ∼ 1.8 (cid:12) with a tail extending to > 2 (cid:12). The neutron star masses in redbacks/huntsmen are consistent with being drawn almost entirely from this proposed second 142 peak of massive neutron stars (see Figure 6.1, Strader et al., 2019). Figure 6.1: Component masses for the redback and huntsmen systems that have well- constrained mass estimates for the neutron star, sorted by neutron star mass. These systems have direct measurements (or meaningful lower limits) of the neutron star mass via a combination of optical spectroscopy, optical photometry, and in some cases the projected motion of the pulsar. Most non-recycled neutron stars with dynamically measured masses are near ∼1.4 (cid:12) (red line), however there is increasing evidence for a broad secondary peak in the neutron star mass distribu- tion that peaks at ∼1.8 (cid:12) (blue line, Özel & Freire, 2016). Nearly all redbacks/huntsmen have neutron star masses consistent with this more massive peak, and therefore, discovering new systems will provide promising opportunities to find the most massive neutron stars. Figure adapted from Strader et al. (2019). The basic expectation is that neutron stars should have initial masses around the Chandrasekhar mass (∼ 1.4 (cid:12), see Section 1.1). However, Figure 6.1 shows that MSPs with non-degenerate companions appear to deviate from this canonical value. These deviations are due both to varia- tions in the initial neutron star masses and in the subsequent accretion (e.g., Özel & Freire, 2016). Determining the mass spectrum of neutron stars constrains the details of supernova explosions 143 and their compact object remnants, including whether there is a mass “gap” (Bailyn et al., 1998) between neutron stars and black holes (currently there are no confirmed black holes of (cid:46) 5 (cid:12), with the possible exception of the recently observed gravitational wave enent GW190814 (Abbott et al., 2020b), and no neutron stars with masses (cid:38) 2.2 (cid:12)). Understanding this gap is especially important in the context of distinguishing between massive neutron stars and lightweight black holes in the new era of gravitational wave events from binary mergers (e.g., Abbott et al., 2020a; The LIGO Scientific Collaboration et al., 2020). Furthermore, the maximum mass of a neutron star has implications far beyond astrophysics, placing tight constraints on the equation of state of dense matter (Section 1.1) and improving our understanding of fundamental physics at supranuclear densities (e.g., Steiner et al., 2013; Lattimer & Steiner, 2014). Although numerous Galactic MSP binaries with non-degenerate companions have been con- firmed in recent years, it is evident from their measured distances that there must be many more out there waiting to be discovered. The median distance of confirmed redbacks is 1.8 kpc, and only a few redbacks or candidates have distances >3 kpc, including the systems presented in this thesis. In addition, new candidates with distances of 1–2 kpc have been discovered in just the past two years. Hence, it is safe to conclude that the redback (and likely MSPs in general) census is highly incomplete beyond 2 kpc, and optical, X-ray, and radio follow-up of new and existing Fermi-LAT sources will continue to yield a substantial return. Even if the limited sample of redbacks that we know of represents 100% of the population within 3 kpc (an extraordinarily conservative estimate), assuming the current fraction of MSP binaries among Fermi sources, and the fact that there are still >1600 -ray sources that are still unidentified, this suggests that hundreds of MSP binaries with non-degenerate companions remain to be discovered in our Galaxy. 144 BIBLIOGRAPHY 145 BIBLIOGRAPHY Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2020a, ApJ, 892, L3 Abbott, R., Abbott, T. D., Abraham, S., et al. 2020b, ApJ, 896, L44 Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJS, 188, 405 Abdo, A. A., Ajello, M., Allafort, A., et al. 2013a, ApJS, 208, 17 Abdo, A. A., Ajello, M., Allafort, A., et al. 2013b, ApJS, 208, 17 Acero, F., Ackermann, M., Ajello, M., et al. 2015, ApJS, 218, 23 Acero, F., Ackermann, M., Ajello, M., et al. 2016, ApJS, 223, 26 Ackermann, M., Ajello, M., Albert, A., et al. 2012, ApJS, 203, 4 Aguilera, D. N., Pons, J. A., & Miralles, J. A. 2008, A&A, 486, 255 Ahn, C. P., Alexandroff, R., Allende Prieto, C., et al. 2012, ApJS, 203, 21 Al Noori, H., Roberts, M. S. E., Torres, R. A., et al. 2018, ApJ, 861, 89 Alpar, M. A., Cheng, A. F., Ruderman, M. A., & Shaham, J. 1982a, Nature, 300, 728 Alpar, M. A., Cheng, A. F., Ruderman, M. A., & Shaham, J. 1982b, Nature, 300, 728 Antoniadis, J. 2014, ApJ, 797, L24 Antoniadis, J., Tauris, T. M., Ozel, F., et al. 2016, arXiv e-prints, arXiv:1605.01665 Antoniadis, J., Freire, P. C. C., Wex, N., et al. 2013, Science, 340, 448 Archibald, A. M., Kaspi, V. M., Bogdanov, S., et al. 2010, ApJ, 722, 88 Archibald, A. M., Stairs, I. H., Ransom, S. M., et al. 2009, Science, 324, 1411 Archibald, A. M., Bogdanov, S., Patruno, A., et al. 2015, ApJ, 807, 62 Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astro- nomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 17 Astraatmadja, T. L. & Bailer-Jones, C. A. L. 2016, ApJ, 832, 137 Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071 Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G., & Andrae, R. 2018, AJ, 156, 58 Bailyn, C. D. 2017, in IAU Symposium, Vol. 324, IAU Symposium, 35 146 Bailyn, C. D., Jain, R. K., Coppi, P., & Orosz, J. A. 1998, ApJ, 499, 367 Baranowski, R., Smolec, R., Dimitrov, W., et al. 2009, MNRAS, 396, 2194 Barthelmy, S. D., Barbier, L. M., Cummings, J. R., et al. 2005, Space Sci. Rev., 120, 143 Bassa, C. G., Patruno, A., Hessels, J. W. T., et al. 2014, MNRAS, 441, 1825 Bates, S. D., Thornton, D., Bailes, M., et al. 2015, MNRAS, 446, 4019 Belczyński, K., Mikołajewska, J., Munari, U., Ivison, R. J., & Friedjung, M. 2000, A&AS, 146, 407 Bellm, E. C., Kaplan, D. L., Breton, R. P., et al. 2016, ApJ, 816, 74 Benvenuto, O. G., De Vito, M. A., & Horvath, J. E. 2014, ApJ, 786, L7 Benvenuto, O. G., De Vito, M. A., & Horvath, J. E. 2015, ApJ, 798, 44 Bernstein, R. M., Burles, S. M., & Prochaska, J. X. 2015, PASP, 127, 911 Bhattacharya, D. 1995, Journal of Astrophysics and Astronomy, 16, 217 Bhattacharya, D. & van den Heuvel, E. P. J. 1991, Phys. Rep., 203, 1 Bisnovatyi-Kogan, G. S. & Komberg, B. V. 1974, Soviet Ast., 18, 217 Blackburn, J. K. 1995, in Astronomical Society of the Pacific Conference Series, Vol. 77, Astro- nomical Data Analysis Software and Systems IV, ed. R. A. Shaw, H. E. Payne, & J. J. E. Hayes, 367 Bogdanov, S. 2016, ApJ, 826, 28 Bogdanov, S., Archibald, A. M., Hessels, J. W. T., et al. 2011, ApJ, 742, 97 Bogdanov, S., van den Berg, M., Heinke, C. O., et al. 2010, ApJ, 709, 241 Bogdanov, S., Archibald, A. M., Bassa, C., et al. 2015, ApJ, 806, 148 Bogdanov, S., Deller, A. T., Miller-Jones, J. C. A., et al. 2018, ApJ, 856, 54 Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127 Breton, R. P., van Kerkwijk, M. H., Roberts, M. S. E., et al. 2013, ApJ, 769, 108 Broderick, J. W., Fender, R. P., Breton, R. P., et al. 2016, MNRAS, 459, 2681 Burderi, L., D’Antona, F., Di Salvo, T., & Burgay, M. 2002, arXiv e-prints, astro Burrows, D. N., Hill, J. E., Nousek, J. A., et al. 2005, Space Sci. Rev., 120, 165 Camilo, F., Kerr, M., Ray, P. S., et al. 2015, ApJ, 810, 85 147 Camilo, F., Reynolds, J. E., Ransom, S. M., et al. 2016, ApJ, 820, 6 Casares, J. 2001, in Lecture Notes in Physics, Berlin Springer Verlag, Vol. 563, Binary Stars: Selected Topics on Observations and Physical Processes, ed. F. C. Lázaro & M. J. Arévalo, 277 Cash, W. 1979, ApJ, 228, 939 Chamel, N., Haensel, P., Zdunik, J. L., & Fantina, A. F. 2013, International Journal of Modern Physics E, 22, 1330018 Chen, H.-L., Chen, X., Tauris, T. M., & Han, Z. 2013, ApJ, 775, 27 Cho, P. B., Halpern, J. P., & Bogdanov, S. 2018, ApJ, 866, 71 Clemens, J. C., Crain, J. A., & Anderson, R. 2004, in Proc. SPIE, Vol. 5492, Ground-based Instrumentation for Astronomy, ed. A. F. M. Moorwood & M. Iye, 331 Corbel, S., Coriat, M., Brocksopp, C., et al. 2013, MNRAS, 428, 2500 Corbel, S., Fender, R. P., Tzioumis, A. K., et al. 2000, A&A, 359, 251 Corbet, R. H. D., Cheung, C. C., Kerr, M., & Ray, P. S. 2012, in Fourth International Fermi Symposium Proceedings, ed. T. J. Brandt, N. Omodei, & C. Wilson-Hodge, eConf C121028, 21 Cordes, J. M. & Lazio, T. J. W. 2002, arXiv:astro-ph/0207156 Coti Zelati, F., Baglio, M. C., Campana, S., et al. 2014, MNRAS, 444, 1783 Cram, L. E. & Mullan, D. J. 1985, ApJ, 294, 626 Cranmer, S. R. & Saar, S. H. 2011, ApJ, 741, 54 Cromartie, H. T., Camilo, F., Kerr, M., et al. 2016, ApJ, 819, 34 Cromartie, H. T., Fonseca, E., Ransom, S. M., et al. 2020, Nature Astronomy, 4, 72 Cumming, A. 2003, in Astronomical Society of the Pacific Conference Series, Vol. 302, Radio Pulsars, ed. M. Bailes, D. J. Nice, & S. E. Thorsett, 277 Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, 2MASS All Sky Catalog of point sources. D’Amico, N., Possenti, A., Manchester, R. N., et al. 2001a, ApJ, 561, L89 D’Amico, N., Possenti, A., Manchester, R. N., et al. 2001b, ApJ, 561, L89 Damour, T. 2009, Astrophysics and Space Science Library, Vol. 359, Binary Systems as Test-Beds of Gravity Theories, ed. M. Colpi, P. Casella, V. Gorini, U. Moschella, & A. Possenti, 1 Davey, S. & Smith, R. C. 1992, MNRAS, 257, 476 de Jager, O. C. & Büsching, I. 2010, A&A, 517, L9 de Jager, O. C., Raubenheimer, B. C., & Swanepoel, J. W. H. 1989, A&A, 221, 180 148 de Martino, D., Falanga, M., Bonnet-Bidaud, J. M., et al. 2010, A&A, 515, A25 de Martino, D., Casares, J., Mason, E., et al. 2014a, MNRAS, 444, 3004 de Martino, D., Casares, J., Mason, E., et al. 2014b, MNRAS, 444, 3004 Deller, A. T., Moldon, J., Miller-Jones, J. C. A., et al. 2015, ApJ, 809, 13 Demorest, P. B., Pennucci, T., Ransom, S. M., Roberts, M. S. E., & Hessels, J. W. T. 2010, Nature, 467, 1081 Deneva, J. S., Ray, P. S., Camilo, F., et al. 2016, ApJ, 823, 105 Desai, S., Armstrong, R., Mohr, J. J., et al. 2012, ApJ, 757, 83 Dickey, J. M. & Lockman, F. J. 1990, ARA&A, 28, 215 Draghis, P., Romani, R. W., Filippenko, A. V., et al. 2019, ApJ, 883, 108 Drake, A. J. 2006, AJ, 131, 1044 Drake, A. J., Djorgovski, S. G., Mahabal, A., et al. 2009, ApJ, 696, 870 Drake, A. J., Graham, M. J., Djorgovski, S. G., et al. 2014, ApJS, 213, 9 Eastman, J., Siverd, R., & Gaudi, B. S. 2010, PASP, 122, 935 Eggleton, P. P. 1983, ApJ, 268, 368 Fruchter, A. S., Gunn, J. E., Lauer, T. R., & Dressler, A. 1988, Nature, 334, 686 Fruscione, A., McDowell, J. C., Allen, G. E., et al. 2006, in Proc. SPIE, Vol. 6270, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 62701V Fukugita, M., Ichikawa, T., Gunn, J. E., et al. 1996, AJ, 111, 1748 Gaensler, B. M., Madsen, G. J., Chatterjee, S., & Mao, S. A. 2008, Publications of the Astronomical Society of Australia, 25, 184 Gallo, E., Degenaar, N., & van den Eijnden, J. 2018, MNRAS, 478, L132 Gallo, E., Fender, R. P., & Pooley, G. G. 2003, MNRAS, 344, 60 Gallo, E., Miller-Jones, J. C. A., Russell, D. M., et al. 2014, MNRAS, 445, 290 Gandhi, P., Rao, A., Johnson, M. A. C., Paice, J. A., & Maccarone, T. J. 2018, submitted to ApJ, arXiv:1804.11349 [astro-ph.HE] Gentile, P. A., Roberts, M. S. E., McLaughlin, M. A., et al. 2014, ApJ, 783, 69 Ginzburg, S. & Quataert, E. 2020, MNRAS, 495, 3656 Goldreich, P. & Julian, W. H. 1969, ApJ, 157, 869 149 Grant, C. E., Bautz, M. W., & Virani, S. N. 2002, in Astronomical Society of the Pacific Conference Series, Vol. 262, The High Energy Universe at Sharp Focus: Chandra Science, ed. E. M. Schlegel & S. D. Vrtilek, 401 Gray, R. O. & Corbally, C. J. 2014, AJ, 147, 80 Greene, J., Bailyn, C. D., & Orosz, J. A. 2001, ApJ, 554, 1290 Greisen, E. W. 2003, in Astrophysics and Space Science Library, Vol. 285, Information Handling in Astronomy - Historical Vistas, ed. A. Heck, 109 Halpern, J. P., Gaidos, E., Sheffield, A., Price-Whelan, A. M., & Bogdanov, S. 2013, The As- tronomer’s Telegram, 5514 Halpern, J. P., Strader, J., & Li, M. 2017, ApJ, 844, 150 Heinke, C. O., Jonker, P. G., Wijnands, R., Deloye, C. J., & Taam, R. E. 2009, ApJ, 691, 1035 Hessels, J. W. T., Ransom, S. M., Stairs, I. H., et al. 2006, Science, 311, 1901 Hewett, P. C. & Wild, V. 2010, MNRAS, 405, 2302 HI4PI Collaboration, Ben Bekhti, N., Flöer, L., et al. 2016, A&A, 594, A116 Hill, A. B., Szostek, A., Corbel, S., et al. 2011, MNRAS, 415, 235 Hogg, D. W. & Foreman-Mackey, D. 2018, ApJS, 236, 11 Hulse, R. A. & Taylor, J. H. 1975, ApJ, 195, L51 Hynes, R. I., Zurita, C., Haswell, C. A., et al. 2002, MNRAS, 330, 1009 Illarionov, A. F. & Sunyaev, R. A. 1975, A&A, 39, 185 Jennings, R. J., Kaplan, D. L., Chatterjee, S., Cordes, J. M., & Deller, A. T. 2018, ApJ, 864, 26 Johansson, J., Thomas, D., & Maraston, C. 2010, MNRAS, 406, 165 Johnson, T. J., Ray, P. S., Roy, J., et al. 2015, ApJ, 806, 91 Kalberla, P. M. W., Burton, W. B., Hartmann, D., et al. 2005, A&A, 440, 775 Kandel, D., Romani, R. W., & An, H. 2019, ApJ, 879, 73 Kaplan, D. L., Bhalerao, V. B., van Kerkwijk, M. H., et al. 2013, ApJ, 765, 158 Kaspi, V. M. & Beloborodov, A. M. 2017, ARA&A, 55, 261 Keane, E. F., Barr, E. D., Jameson, A., et al. 2018, MNRAS, 473, 116 Kerr, M. 2011, ApJ, 732, 38 Kluzniak, W., Ruderman, M., Shaham, J., & Tavani, M. 1988, Nature, 334, 225 150 Kong, A. K. H., Huang, R. H. H., Cheng, K. S., et al. 2012, ApJ, 747, L3 Kramer, M., Lyne, A. G., Hobbs, G., et al. 2003, ApJ, 593, L31 Kreidberg, L., Bailyn, C. D., Farr, W. M., & Kalogera, V. 2012, ApJ, 757, 36 Landolt, A. U. 1992, AJ, 104, 340 Lattimer, J. M. & Prakash, M. 2001, ApJ, 550, 426 Lattimer, J. M. & Steiner, A. W. 2014, ApJ, 784, 123 Li, K.-L., Kong, A. K. H., Hou, X., et al. 2016, ApJ, 833, 143 Li, K.-L., Hou, X., Strader, J., et al. 2018, ApJ, 863, 194 Linares, M. 2014, ApJ, 795, 72 Linares, M. 2018, MNRAS, 473, L50 Linares, M., Shahbaz, T., & Casares, J. 2018, ApJ, 859, 54 Lipunov, V. M. 1987, Ap&SS, 132, 1 Lorimer, D. R. & Kramer, M. 2004, Handbook of Pulsar Astronomy, Vol. 4 Lorimer, D. R. & Kramer, M. 2012, Handbook of Pulsar Astronomy Lorimer, D. R., Yates, J. A., Lyne, A. G., & Gould, D. M. 1995, MNRAS, 273, 411 Margalit, B. & Metzger, B. D. 2017, ApJ, 850, L19 Marigo, P., Girardi, L., Bressan, A., et al. 2008, A&A, 482, 883 Maron, O., Kijak, J., Kramer, M., & Wielebinski, R. 2000, A&AS, 147, 195 Marsh, F. M., Prince, T. A., Mahabal, A. A., et al. 2017, MNRAS, 465, 4678 Matsuoka, M., Kawasaki, K., Ueno, S., et al. 2009, PASJ, 61, 999 Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1996, ApJ, 461, 396 McConnell, O., Callanan, P. J., Kennedy, M., et al. 2015, MNRAS, 451, 3468 McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127 Migliari, S. & Fender, R. P. 2006, MNRAS, 366, 79 Migliari, S., Fender, R. P., Rupen, M., et al. 2003, MNRAS, 342, L67 Migliazzo, J. M., Gaensler, B. M., Backer, D. C., et al. 2002, ApJ, 567, L141 151 Monet, D. G., Levine, S. E., Canzian, B., et al. 2003, AJ, 125, 984 Morin, J. 2012, in EAS Publications Series, Vol. 57, EAS Publications Series, ed. C. Reylé, C. Charbonnel, & M. Schultheis, 165 Mucciarelli, A., Salaris, M., Lanzoni, B., et al. 2013, ApJ, 772, L27 Munari, U. & Zwitter, T. 1997, A&A, 318, 269 Nolan, P. L., Abdo, A. A., Ackermann, M., et al. 2012a, ApJS, 199, 31 Nolan, P. L., Abdo, A. A., Ackermann, M., et al. 2012b, ApJS, 199, 31 Orosz, J. A. & Hauschildt, P. H. 2000a, A&A, 364, 265 Orosz, J. A. & Hauschildt, P. H. 2000b, A&A, 364, 265 Orosz, J. A. & van Kerkwijk, M. H. 2003, A&A, 397, 237 Özel, F. & Freire, P. 2016, ARA&A, 54, 401 Papitto, A. & Torres, D. F. 2015, ApJ, 807, 33 Papitto, A., Ferrigno, C., Bozzo, E., et al. 2013, Nature, 501, 517 Pfahl, E., Rappaport, S., Podsiadlowski, P., & Spruit, H. 2002, ApJ, 574, 364 Pickles, A. J. 1998, Publications of the Astronomical Society of the Pacific, 110, 863 Pitkin, M. 2018, JOSS, 3, 538 Plotkin, R. M., Miller-Jones, J. C. A., Gallo, E., et al. 2017a, ApJ, 834, 104 Plotkin, R. M., Bright, J., Miller-Jones, J. C. A., et al. 2017b, ApJ, 848, 92 Podsiadlowski, P., Rappaport, S., & Pfahl, E. D. 2002, ApJ, 565, 1107 Pons, J. A., Link, B., Miralles, J. A., & Geppert, U. 2007, Phys. Rev. Lett., 98, 071101 Price-Whelan, A. M., Hogg, D. W., Foreman-Mackey, D., & Rix, H.-W. 2017, ApJ, 837, 20 Radhakrishnan, V. & Srinivasan, G. 1982, Current Science, 51, 1096 Ransom, S. M. 2013, in IAU Symposium, Vol. 291, Neutron Stars and Pulsars: Challenges and Opportunities after 80 years, ed. J. van Leeuwen, 3 Ransom, S. M., Ray, P. S., Camilo, F., et al. 2011, ApJ, 727, L16 Reichart, D., Nysewander, M., Moran, J., et al. 2005, Nuovo Cimento C Geophysics Space Physics C, 28, 767 Rivera Sandoval, L. E., Hernandez Santisteban, J. V., Degenaar, N., et al. 2017, ArXiv e-prints, arXiv:1708.07041 [astro-ph.HE] 152 Rivera Sandoval, L. E., Hernández Santisteban, J. V., Degenaar, N., et al. 2018, MNRAS, 476, 1086 Roberts, M. S. E. 2011, in American Institute of Physics Conference Series, ed. M. Burgay, N. D’Amico, P. Esposito, A. Pellizzoni, & A. Possenti, Vol. 1357, 127 Roberts, M. S. E., Mclaughlin, M. A., Gentile, P., et al. 2014, Astronomische Nachrichten, 335, 313 Roberts, M. S. E., McLaughlin, M. A., Gentile, P. A., et al. 2015, 2014 Fermi Symposium - eConf C14102.1, arXiv:1502.07208 [astro-ph.HE] Romani, R. W., Filippenko, A. V., & Cenko, S. B. 2015a, ApJ, 804, 115 Romani, R. W., Graham, M. L., Filippenko, A. V., & Kerr, M. 2015b, ApJ, 809, L10 Romani, R. W. & Sanchez, N. 2016, ApJ, 828, 7 Roy, J., Ray, P. S., Bhattacharyya, B., et al. 2015, ApJ, 800, L12 Sabbi, E., Gratton, R., Ferraro, F. R., et al. 2003, ApJ, 589, L41 Sanchez, N. & Romani, R. W. 2017, ApJ, 845, 42 Sault, R. J., Teuben, P. J., & Wright, M. C. H. 1995, in Astronomical Society of the Pacific Conference Series, Vol. 77, Astronomical Data Analysis Software and Systems IV, ed. R. A. Shaw, H. E. Payne, & J. J. E. Hayes, 433 Savonije, G. J. 1987, Nature, 325, 416 Scargle, J. D. 1982, ApJ, 263, 835 Schlafly, E. F. & Finkbeiner, D. P. 2011a, ApJ, 737, 103 Schlafly, E. F. & Finkbeiner, D. P. 2011b, ApJ, 737, 103 Schroeder, J. & Halpern, J. 2014, ApJ, 793, 78 Shahbaz, T., Linares, M., & Breton, R. P. 2017, MNRAS, 472, 4287 Shajn, G. & Struve, O. 1929, MNRAS, 89, 222 Smith, D. A., Guillemot, L., Kerr, M., Ng, C., & Barr, E. 2017, arXiv e-prints, arXiv:1706.03592 Spitkovsky, A. 2011, Astrophysics and Space Science Proceedings, 21, 139 Stairs, I. H. 2003, Living Reviews in Relativity, 6, 5 Stappers, B. W., van Kerkwijk, M. H., Bell, J. F., & Kulkarni, S. R. 2001, ApJ, 548, L183 Steiner, A. W., Lattimer, J. M., & Brown, E. F. 2013, ApJ, 765, L5 Stinebring, D. R. & Condon, J. J. 1990, ApJ, 352, 207 153 Stovall, K., Lynch, R. S., Ransom, S. M., et al. 2014, ApJ, 791, 67 Strader, J., Chomiuk, L., Sonbas, E., et al. 2014, ApJ, 788, L27 Strader, J., Li, K.-L., Chomiuk, L., et al. 2016, ApJ, 831, 89 Strader, J., Smith, G. H., Larsen, S., Brodie, J. P., & Huchra, J. P. 2009, AJ, 138, 547 Strader, J., Chomiuk, L., Cheung, C. C., et al. 2015, ApJ, 804, L12 Strader, J., Swihart, S., Chomiuk, L., et al. 2019, ApJ, 872, 42 Stroh, M. C. & Falcone, A. D. 2013, ApJS, 207, 28 Swihart, S. J., Strader, J., Chomiuk, L., & Shishkovsky, L. 2019, ApJ, 876, 8 Swihart, S. J., Strader, J., Johnson, T. J., et al. 2017, ApJ, 851, 31 Swihart, S. J., Strader, J., Shishkovsky, L., et al. 2018, ApJ, 866, 83 Swihart, S. J., Strader, J., Urquhart, R., et al. 2020, ApJ, 892, 21 Tang, S., Kaplan, D. L., Phinney, E. S., et al. 2014, ApJ, 791, L5 Tauris, T. M. 2012, Science, 335, 561 Tauris, T. M. & Savonije, G. J. 1999, A&A, 350, 928 Tauris, T. M. & van den Heuvel, E. P. J. 2006, Formation and evolution of compact stellar X-ray sources, 623 Tetarenko, A. J., Bahramian, A., Sivakoff, G. R., et al. 2016, MNRAS, 460, 345 The Fermi-LAT collaboration. 2019, arXiv e-prints, arXiv:1902.10045 The LIGO Scientific Collaboration, the Virgo Collaboration, Abbott, R., et al. 2020, arXiv e-prints, arXiv:2006.12611 Tody, D. 1986, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 627, Instrumentation in astronomy VI, ed. D. L. Crawford, 733 Tudor, V., Miller-Jones, J. C. A., Patruno, A., et al. 2017, MNRAS, 470, 324 Tudose, V., Fender, R. P., Linares, M., Maitra, D., & van der Klis, M. 2009, MNRAS, 400, 2111 Urpin, V. & Konenkov, D. 1997, MNRAS, 284, 741 van den Heuvel, E. P. J. 2019, in IAU Symposium, Vol. 346, IAU Symposium, ed. L. M. Oskinova, E. Bozzo, T. Bulik, & D. R. Gies, 1 van Staden, A. D. & Antoniadis, J. 2016, ApJ, 833, L12 Venter, C., Kopp, A., Harding, A. K., Gonthier, P. L., & Büsching, I. 2015, ApJ, 807, 130 154 Viganò, D., Pons, J. A., & Miralles, J. A. 2012, Computer Physics Communications, 183, 2042 Viganò, D., Rea, N., Pons, J. A., et al. 2013, MNRAS, 434, 123 Vogt, S. S., Allen, S. L., Bigelow, B. C., et al. 1994, in Proc. SPIE, Vol. 2198, Instrumentation in Astronomy VIII, ed. D. L. Crawford & E. R. Craine, 362 Wadiasingh, Z., Harding, A. K., Venter, C., Böttcher, M., & Baring, M. G. 2017, ApJ, 839, 80 Walter, F. M., Battisti, A., Towers, S. E., Bond, H. E., & Stringfellow, G. S. 2012, PASP, 124, 1057 Wang, Z., Archibald, A. M., Thorstensen, J. R., et al. 2009, ApJ, 703, 2017 Werner, G. R., Uzdensky, D. A., Cerutti, B., Nalewajko, K., & Begelman, M. C. 2016, ApJ, 816, L8 Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914 Wu, J., Orosz, J. A., McClintock, J. E., et al. 2016, ApJ, 825, 46 Yao, J. M., Manchester, R. N., & Wang, N. 2017, ApJ, 835, 29 Zahn, J. P. 1977, A&A, 500, 121 Zharikov, S., Kirichenko, A., Zyuzin, D., Shibanov, Y., & Deneva, J. S. 2019, MNRAS, 489, 5547 155