# Supplementary Information for Random-access quantum memory using chirped
pulse phase encoding

James O’Sullivan  These authors have contributed equally to this work  London
Centre for Nanotechnology, UCL, 17-19 Gordon Street, London, WC1H 0AH, UK
Oscar W. Kennedy  These authors have contributed equally to this work  London
Centre for Nanotechnology, UCL, 17-19 Gordon Street, London, WC1H 0AH, UK
Kamanasish Debnath  Department of Physics and Astronomy, Aarhus University,
DK-8000 Aarhus C, Denmark  Joseph Alexander  London Centre for Nanotechnology,
UCL, 17-19 Gordon Street, London, WC1H 0AH, UK  Christoph W. Zollitsch  London
Centre for Nanotechnology, UCL, 17-19 Gordon Street, London, WC1H 0AH, UK
Mantas Šimėnas  London Centre for Nanotechnology, UCL, 17-19 Gordon Street,
London, WC1H 0AH, UK  Akel Hashim  Lawrence Berkeley National Laboratory,
Berkeley, CA 94720, USA  Christopher N. Thomas  Cavendish Laboratory,
University of Cambridge, JJ Thomson Ave, Cambridge CB3 0HE, UK  Stafford
Withington  Cavendish Laboratory, University of Cambridge, JJ Thomson Ave,
Cambridge CB3 0HE, UK  Irfan Siddiqi  Lawrence Berkeley National Laboratory,
Berkeley, CA 94720, USA  Klaus Mølmer  Department of Physics and Astronomy,
Aarhus University, DK-8000 Aarhus C, Denmark  John J. L. Morton  London Centre
for Nanotechnology, UCL, 17-19 Gordon Street, London, WC1H 0AH, UK  Department
of Electrical and Electronic Engineering, UCL, Malet Place, London, WC1E 7JE,
UK

##  I Methods

###  I.1 Device fabrication

A float-zone silicon wafer of natural isotopic abundance is ion implanted with
bismuth at a chain of energies targeting a density of $10^{17}$ cm-3 bismuth
in the top 1 $\upmu$m of the sample. This is annealed at 900∘C for 5 minutes
to incorporate the bismuth into the silicon matrix forming spin-active donors.
A planar niobium microresonator is patterned on the top of the sample by
liftoff. The Nb film is 100 nm thick and has a field dependent frequency (see
Ref. o2020spin) of 7.093 GHz at a magnetic field of 46 mT. The resonator is a
lumped element design comprising a pair of parallel capacitive plates and a
narrow doubled-back inductor wire, designed to enhance the magnetic field
generated by the resonator close to the substrate surface. The sample is
placed inside a copper cavity, as shown in Fig. S1(d). Two antennae extend
into the box; one short stub antenna (high insertion loss $\gtrsim 30$ dB, low
coupling) is connected to a 30 dB attenuated microwave in-line and used to
apply microwave excitations. Another long antenna (large coupling) is used to
collect microwave signals and is connected to the amplification chain. The
asymmetry results in a large collection efficiency of spin-echo signals.

A schematic of the device is shown in Fig. S1(a), the simulated bismuth
implantation profile in Fig. S1(b), the simulated microwave magnetic fields
caused by the zero point fluctuations in the resonator in Fig. S1(c) and a
schematic of the 3D cavity in Fig. S1(d).

Figure S1:  The prototype memory. (a) Illustration of the silicon wafer with
near-surface implanted Bi layer and Nb resonator patterned on the surface. (b)
Simulated Bi implantation profile as a function of depth below the substrate
surface. (c) Finite element simulations of the magnetic field generated in the
vicinity of the resonator inductor (d) Illustration of the copper sample box
with asymmetric antennae.  Figure S2: Schematic of the spectrometer and
dilution fridge setup. Microwave pulse signals are generated using a vector
signal generator (VSG) and an arbitrary waveform generator (AWG). A microwave
generator (MWG) is used to drive the Josephson parametric amplifier (JPA). The
return signals are down-converted and detected with a digitizer. A signal at
half the up-conversion local oscillator (LO) frequency provided by the VSG is
doubled with a frequency doubler and used for down-conversion. Coaxial cable
types inside the fridge have been colour coded.

###  I.2 Microwave measurement

The measurement setup is shown in Fig. S2. The copper sample box is mounted at
the base plate of a dilution refrigerator at 100 mK. The output signal is
routed through two circulators to a Josephson parametric amplifier (JPA), a
quantum limited amplifier, which amplifies the signal in reflection. The
amplifier is driven via a directional coupler by a dedicated microwave source.
Calibration of the JPA is given in the next section. The signal was further
amplifed at 4 K by a high electron mobility transistor (HEMT) and again at
room temperature before being mixed down in frequency by an IQ mixer and
detected using a digitiser.

Pulsed microwave signals are generated by an arbitrary waveform generator
(AWG), connected to a vector signal generator (VSG). These pulses are routed
through either a high gain path (solid state amplifier +35 dB 3W max output)
or an attenuating path (10 dB attenuation) both signal routes are actively
gated using fast microwave switches. Signals are routed through the fridge
before being mixed down at room temperature.

Continuous wave measurements are performed using a vector network analyser
(VNA). We measured $S_{21}$ transmission between the antennae which is
modulated by the microresonator. To characterize the resonator quality factor
we fit the modulation (in linear magnitude) with a Breit-Wigner-Fano function:

$S_{\rm 21Lin}(f)=K\frac{q\kappa/2+f-f_{0}}{\kappa^{2}/4+(f-f_{0})^{2}}+mx+c$
(1)

where $K$ determines the size of the modulation, $q$ is an asymmetry
parameter, $f_{0}$ is the central frequency of the resonator, $\kappa$ is the
FWHM of the cavity and the $mx+c$ term is an approximation to the background
transmission. This fits our resonance notch well and an example when the
resonator is off resonance with the spins is shown in Fig. S3(c)

Figure S3: Effect of spins on (a) resonator frequency and (b) resonator
linewidth $\kappa$ as a function of magnetic field as the spin line passes
through the resonator. Data are fit simultaneously to equations in Ref.
Eisuke2011electron $g_{\rm ens}\sim$120 kHz and $\gamma\sim 2.4$ MHz. (c)
Measurement of transmission magnitude through the copper cavity using a VNA at
46mT (off resonance with spin line). The superconducting resonator gives a
characteristic Fano resonance response which we fit to extract a centre
frequency for this resonator of 7.09395 GHz and a Q factor of 18400.

###  I.3 System Calibration

Figure S4: (a) $T_{1}$ by inversion recovery. A WURST pulse is used to invert
the ensemble before a wait (indicated on the x axis) and a detection sequence
at constant high power and $\tau$. (b) $T_{2}$ measurement of the sample at
the same weak excitation strength as used in Fig. 2 main text. The pulse
sequence (inset) is a weak excitation followed by two 100 $\upmu$s WURST
pulses. Increasing $\tau$ changes the echo time allowing $T_{2}$ to be
measured. Fitting stretched or single exponentials gives different $T_{2}$
giving a combined uncertainty of $T_{2}=0.7(2)$ ms (c) Rabi oscillations using
a Gaussian $\theta$-pulse of duration 8 $\upmu$s and FWHM 4 $\upmu$s and a
two-WURST silenced echo detection sequence. The oscillations are heavily
damped due to a large inhomogeneity in Rabi frequency across the ensemble.

We measure $T_{1}$, $T_{2}$ and the Rabi frequency of the hybrid system and
show these measurements in Fig. S4. $T_{1}$ is measured by inversion recovery
where the WURST pulse inverts the spin ensemble and we measure the time taken
for the inverted ensemble to return to the ground state and find $T_{1}\sim$14
s. $T_{2}$ is measured using two WURST pulses to refocus a weak excitation
(same strength as in Fig. 2). Varying $\tau$ allows us to measure the effect
of dephasing on the echo amplitude. We fit a single quadrature of the echo
decay returning a coherence time of 0.7(2) ms where most of the uncertainty
arises due to different values returned fitting a stretched or single
exponential to the data.

To measure the Rabi frequency, $\Omega$, we appply a $\theta$ rotation pulse
followed 10 ms later by a detection sequence with VSG output power of -20 dBm
sent through the high gain path. The $\theta$-pulse is chosen to match the
duration and shape of the excitations used in the memory sequence in Fig. 3.
The heavily damped envelope to sinusoidal Rabi oscillations (shown in Fig. S4)
is due to a large inhomogeneity in single spin coupling, $g_{0}$, across the
ensemble and introduces substantial ($\sim 50\%$) uncertainty to this
measurement. We fit a decaying cosine function to the echo amplitude, allowing
us to extract an approximate $\pi$ pulse amplitude of 0.435 (where 1 is the
maximum pulse amplitude at the chosen VSG output power) giving $\Omega=125$
kHz.

Using $T_{1}$ and $\Omega$ we calibrate the photon number in the resonator.
The Purcell effect limits $1/T_{1}=4g_{0}^{2}/\kappa$, allowing a measure of
average spin-resonator coupling $g_{0}$. The different $\kappa$ dependence
from Ref. bienfait2016controlling is due to the definition of $\kappa$ being
the HWHM in this work and FWHM in Ref. bienfait2016controlling. At the centre
of the line $\kappa\sim 400$ kHz, $T_{1}=14.7$ s giving
$g_{0}\sim\leavevmode\nobreak\ 80$ Hz ($g_{0}$ varies across the ensemble and
this number is indicative). The Rabi frequency can be written as
$\Omega=2g_{0}\sqrt{n}$ which allows us to calibrate the photon number $n\sim
5.7\times 10^{5}$ for the $\pi$ pulse above. We can rescale this to the photon
number to the memory sequence in Fig. 3 (main text) based on an amplitude of
0.2 at a power $\sim$20 dB lower than the Rabi measurement and find $\langle
n\rangle\sim$ 1200 photons for the memory protocol. Due to the uncertainty in
Rabi measurements these photon numbers should be treated as indicative powers
accurate to approximately $50\%$.

We note that in Fig. 2(d) main text, that repeated WURST pulses before a pulse
sequence can cause a reduction in echo amplitude which implies that WURST
pulses cause some saturation of the ensemble. However, comparing the $T_{2}$
decay shown here and the $T_{M}$ (Fig. 2(e) main text) acquired when using
repeated WURST pulses, we see that repeated WURST pulses actually increase the
final echo amplitude due to dynamical decoupling. These two observations
appear initially at odds. However, we understand this discrepancy using a toy
model where we divide spins into three zones, (i) weakly coupled spins, far
from the resonator and unaffected by WURST pulses, (ii) spins a moderate
distance from the resonator, in turn moderately coupled and undergo imperfect
inversion under WURST pulses (iii) spins close and strongly coupled to the
resonator which are inverted by WURST pulses with high fidelity. Repeated
WURST pulses will saturate zone (ii) which would contribute some of the echo
strength. When operating this memory with a long string of WURST pulses, zone
(ii) would always remain saturated, and therefore not contribute to the
protocol. This means that in future devices, sufficient cooperativity for high
fidelity memory transfer must be achieved with only spins in zone (iii).

The use of a JPA in measurements was essential to perform experiments at low
photon numbers. Such amplifiers are prone to saturation and nonlinearity with
large input signals. We calibrate the JPA to determine the onset and extent of
the nonlinear regime, and any additional distortions that may be present. For
comparison, the same measurements were repeated with the JPA turned off, using
only a HEMT at 4 K for amplification. We confirm that for Fig. 3 (main
text)the echo signals are well into the linear regime but the input
excitations are in the nonlinear regime of the JPA and are prone to
distortion. The HEMT remains linear and un-distorted at all powers, as
expected. JPA calibrations are shown in Fig. S5. The inhomogeneous linewidth
of the ensemble is $\sim$2.5 MHz (see Ref. o2020spin for further details).

Figure S5: (a) Measured pulse amplitude as a function of input pulse
amplitude, measured using the JPA (blue) and the HEMT without JPA (orange).
(b) Low power regime of the same experiment.

To further increase signal we chose an operating temperature to maximise the
population of the ground state of our chosen bismuth transition ($\sim$10 %).
Pulsed experiments are run at long a shot-repetition rate of 160 s to minimise
the saturation of spins. These considerations are to improve the cooperativity
between the spin ensemble and resonator which contributes to the read/write
efficiency of the memory and thus improve signal - approximately 17 % one way
efficiency as shown in the main text. This means that when measuring $\sim$200
photon input pulses, the output pulse is actually $0.17^{2}\times 200\sim 6$
photon output pulses - without considering any $T_{2}$ decay. Despite our
optimized antenna configuration, we may fail to capture all of this output
pulse, and also find that the JPA is working slightly sub-optimally (given the
elevated temperature). The combination of these effects means that to observe
a low-noise echo at low powers we typically required hundreds of averages.
Recent work using single microwave photon detectors albertinale2021detecting
is a promising route to improve sensitivity. In tandem it is necessary to
consider routes to improve memory efficiency in order to achieve useful
memories. Promising approaches to this include the use of nuclear-spin 1/2
species (allowing the ensemble to be 100 % thermally polarised at mK), optimal
resonator design, increasing filling factor whilst maintaining large single-
spin coupling necessary for WURST inversion pulses and the use of 28Si
substrates which would improve coherence times and also narrow inhomogeneous
spin linewidths.

###  I.4 Cooperativity

Cooperativity in a hybrid micro-resonator spin system can be measured by the
dispersive shift to the resonator frequency ($f$), and the increased half-
width half max of the resonator frequency response ($\kappa$)
Eisuke2011electron:

$f=f_{0}-\frac{g_{\rm ens}^{2}\Delta}{\Delta^{2}+\gamma^{2}/4}$ (2)
$\kappa=\kappa_{0}+\frac{g_{\rm ens}^{2}\gamma/2}{\Delta^{2}+\gamma^{/}4}$ (3)

where $\kappa_{0}$ and $f_{0}$ are respectively the half width and frequency
of the resonator in the absence of spins, and the detuning
$\Delta=(B_{0}-B_{R})\times\partial f/\partial B_{0}$, where $B_{R}$ is the
magnetic field at which spins and resonator are resonant. Fitting equations 2
and 3 to our data gives $C\sim 0.067$ as shown in Fig. S3.

We also extract the cooperativity from $T_{2}$ measurements with CPMG
sequences in Fig. 2 (main text) by modelling the echo intensity as a function
of echo number and considering the reduction of energy in the echo field with
repeated emission. In A[BB]nA sequences the echo is emitted after the final A
pulse. As (i) WURST pulses efficiently refocus spins and (ii) only one echo
forms, we attribute any change in echo intensity to dephasing of the spins. We
fit the echo intensity $A_{\rm sil}(t)$ to extract the dephasing time $T_{2}$
using

$A_{\rm sil}(t)=\sqrt{A_{0}^{2}\exp(-2t/T_{2})+K^{2}},$ (4)

where $A_{0}$ is the echo amplitude at time $t=0$, $T_{2}$ is the dephasing
time and $K$ is the background giving a dephasing time of $T_{2}=2.0\pm 0.2$
ms for a refocusing rate of 7.1 kHz. We model the echo amplitude in A[AA]nA
sequences by assuming that every time an echo forms, a constant fraction of
the energy stored in the echo field is lost and that in between echo emission
the echo field dephases as in A[BB]nA sequences. We model the echo amplitude
by

$A(t)=\sqrt{A_{0}^{2}\exp(-2t/T_{2})\times(1-\eta_{\rm em})^{N}+K^{2}}$ (5)

where $\eta_{\rm em}$ is the efficiency with which energy is emitted from the
echo field to the cavity, $N$ is the number of echoes which have occurred
before time $t$. This is fit with one free parameter, $\eta_{\rm em}$, and is
shown shown in Fig. 2.

Using Equation (17) from Ref. afzelius2013proposal we relate the one-way
efficiency to the cooperativity by

$\eta_{\rm em}=\frac{4C}{(1+C)^{2}}$ (6)

giving $C=0.047$ in close agreement to the value extracted in the more typical
method by measuring effects on the resonator.

###  I.5 WURST Pulses

WURST pulses are a type of adiabatic fast passage which have been used to
control spin ensembles o2013wurst. They have advantages in the context of
semi-classical control as, in the adiabatic limit, they impart faithful $\pi$
pulses across a wide bandwidth of spins irrespective of the coupling between
the spin and the drive field which may be inhomogeneous across the ensemble.
By chirping the frequency of the control pulse, the effective magnetic field
observed from the rotating frame of the un-driven spin flips from pointing
along one pole of the Bloch sphere to pointing along the other pole following
a path along the surface of the Bloch sphere. In this reference frame the spin
precesses around this effective magnetic field, and if that precession rate is
sufficiently high relative to the rate of the field reversal (i.e. the pulse
is sufficiently strong relative to the chirp rate) then the spin will
adiabatically track this magnetic field and undergo a $\pi$ inversion pulse.
This is described well in section 2.2 of Ref. o2013wurst. In this work we use
WURST pulses of order 20.

Figure S6: WURST-20 pulses with (a) the pulse amplitude and frequency shift
for a WURST pulse with bandwidth 2 MHz and duration 100 $\upmu$s. (b) The I
and Q quadratures of the same WURST pulse. (c) The X and Y field quadratures
in the cavity once the WURST pulse is filtered by a $\kappa=$200 kHz cavity.

####  I.5.1 Input Output Theory

Wideband uniform rate smooth truncation (WURST) pulses can be characterized by
frequency and amplitude modulation (FM and AM respectively)

$\displaystyle f_{\rm W}(t)=-\frac{\Gamma_{\rm W}}{2}+\frac{\Gamma_{\rm
W}}{T_{\rm W}}t,$ (7) $\displaystyle A_{\rm
W}(t)=1-\absolutevalue{\sin(\pi\frac{t-T_{\rm W}/2}{T_{\rm W}})}^{N},$

where $\Gamma_{\rm W}$ is the bandwidth of the WURST pulse, $T_{\rm W}$ is the
WURST pulse duration and $N$ is the WURST pulse index, in this work we use
$N=20$ pulses.

The FM results in a time varying phase of the output pulse

$\displaystyle\phi(t)=\int_{0}^{t}f_{\rm W}(t^{\prime})dt^{\prime}+\phi_{0},$
(8)

where $\phi_{0}$ is the phase of the WURST pulse. The I and Q pulse
quadratures applied by the AWG are

$\displaystyle I(t)=A_{\rm W}(t)\sin\phi(t),$ $\displaystyle Q(t)=A_{\rm
W}(t)\cos\phi(t),$ (9)

The WURST pulse is modulated by the cavity described by input/output theory
gardiner1985input; ranjan2020pulsed where

$\displaystyle\dot{X}(t)=\sqrt{\kappa_{C}}I(t)-\frac{\kappa}{2}X(t),$
$\displaystyle\dot{Y}(t)=\sqrt{\kappa_{C}}Q(t)-\frac{\kappa}{2}Y(t),$ (10)

where $\kappa_{C}$ is the coupling linewidth of the cavity and $\kappa$ is the
loaded linewidth of the cavity. Increasing the bandwidth of the WURST pulse
above the bandwidth of the cavity shortens the effective WURST pulse duration.
We show an example WURST pulse in FM/AM (Fig. S6a), consequent IQ quadratures
(Fig. S6b) and after modulation by a cavity (Fig. S6c)

####  I.5.2 Limits on WURST Pulses

In Fig. 4, a region of suitable WURST pulses — WURST pulses which inverts a
large fraction of the spin ensemble — is mapped out based on both theory and
experiment described in this section. Using pulse sequences such as that in
Fig. 1(a) (excitation followed by two identical WURST pulses) we have measured
the refocusing efficiency. We measure the amplitude of the first (nominally
silenced) and second echo and present them as a function of the WURST pulse
bandwidth and amplitude in Fig. S7. A ‘good’ WURST pulse results in no echo
after the first WURST, and a loud echo after the second.

By thresholding the lines we determine regions where WURST pulses
adiabatically refocus a large fraction of spins. The constraint that the
second echo must have a suitable magnitude (85 mV) gives the experimental
boundary to suitable WURST pulses shown in Fig. 4 (main text) based on
adiabaticity.

This is supported by a theoretical treatment doll2013adiabatic.

$Q_{\rm min}=2\pi\nu^{2}/R\gg 1,$ (11)

where $Q_{\rm min}$ is an adiabaticity factor, $R$ is the chirp rate and $\nu$
is the nutation frequency (how quickly the spin precesses around the effective
field). $\nu$ is minimal (and equal to the Rabi frequency) when the WURST
frequency is the same as the spin frequency. We choose a minimum $g_{0}$ for
spins to be refocused and solve this equation for $Q_{\rm min}=1$ finding a
lower limit on pulse strength.

In the homodyne scheme used in this work there is a limit placed by the finite
frequency resolution we can imprint onto WURST pulses. The effective duration
of the WURST pulse, $T_{\rm W,eff}=\kappa T_{\rm W}/\Gamma_{\rm W}$, is the
time the WURST pulse of duration $T_{\rm W}$ and total bandwidth $\Gamma_{\rm
W}$ takes to chirp across the cavity frequency $\kappa$. The inverse of this
time gives an indication of the minimum frequency that can be resolved by
modulating the pulse. This frequency must be less than the cavity bandwidth
placing a limit on the maximum $\Gamma_{\rm W}$ which can be used,

$\Gamma_{\rm W}\ll\kappa^{2}T_{\rm W}=37\leavevmode\nobreak\ \mathrm{MHz}.$
(12)

The maximum spectrometer power limits the power of WURSTs we can apply and is
shown as the top boundary in Fig. 4. This line cannot be increased to
arbitrarily high powers using higher power amplifiers as the superconducting
resonators will limit the maximum power. There is also a minimum WURST
bandwidth imposed by the cavity bandwidth.

####  I.5.3 WURST pulse distinctiveness

Figure S7: Top (bottom) panel – amplitude of the first (second), nominally
silenced (unsilenced), echo as a function of the pulse bandwidth for 200us
WURST pulses. Inset to the bottom panel is the unsilenced data for low
bandwidths, showing that there is a minimum bandwidth, similar to the cavity
width, required to maximise the echo amplitude by refocusing all the spins in
the cavity. The adiabaticity limit is determined by thresholding the
unsilenced echo intensity relative to 85 mV shown as a black dashed line. Non-
monotocity of echo intensity and uncertainty in echo amplitude give
uncertainty in the thresolded value shown in Fig. 4.

We partition the space of suitable WURSTs into ‘distinct’ WURST pulses to
estimate the memory capacity. We partition this region based upon two-WURST
echo sequences where the first and second WURST pulse have different
properties. Inset to Fig. 4 (main text) is the echo amplitude as the amplitude
and chirp rate of the B pulse is varied. This sequence results in strong
echoes along a line with positive gradient. In the main text we explain that
this line of equivalent pulses is due to the acquisition of a dynamic phase.
Together with other two dimensional and one dimensional maps varying a B-pulse
parameter we interpolate the gradient and width of the line of equivalent
pulses across the space of WURST pulses giving good inversion and use this
interpolation to count WURST pulses. This interpolation is shown in Fig. S8
with the results in Fig. 4.

Figure S8: (a) The line of strong echoes from 2D maps of AB-echo sequences. We
fit a line to these data and extract their gradient. (b) The gradient of the
lines fit in (a) as a function of their extrapolated bandwidth at amplitude
$A$=1 (black circles). The gradient determined when B-pulse has been swept in
1D in both amplitude and bandwidth are inferred and plot on the same axes
(grey triangles) showing good agreement. We fit these data with $\partial
A/\partial R=C/R$ with $C$ a fit parameter and return the red dashed line and
$C=0.93$. (c) Width of the line of equivalent WURST pulses as a function of
WURST pulse parameters. We fit a plane to this data to allow interpolation of
the width as a function of bandwidth and amplitude.

####  I.5.4 Theoretical description of WURST pulses

The cavity is driven by a chirped field which results in a temporal intensity
profile of the field inside the cavity, and the coupled equations for the
spins and the cavity field must in general be solved numerically. Here we
shall make some observations for a simplified model, that will qualitatively
explain most features seen in the experiments. The spins have a central
frequency $\omega_{c}$ and each spin is characterized by its detuning
$\delta=\omega-\omega_{c}$. We characterize a WURST pulse by its chirp rate
$R$ and pulse centre $t_{0}$, yielding its frequency $\omega(t)=R(t-t_{0})$
and phase $\phi(t)=\frac{1}{2}R(t-t_{0})^{2}$ in the frame rotating at
$\omega_{c}$. In that frame the Hamiltonian for each given spin is
$\delta\sigma_{z}/2$, and the time dependent complex Rabi frequency is
$\Omega(t)=\Omega_{0}\exp(-i\phi(t))$.

The spin transitions are resonant when $R(t-t_{0})=\delta$, i.e., at the time
$t_{\delta}\equiv t_{0}+\delta/R$. In the interaction picture with respect to
the spin excitation energy $\delta\sigma_{z}/2$, the Hamiltonian of a single
spin reads,

$H_{I}(t)=\frac{\Omega_{0}}{2}(e^{-i\phi(t)}e^{i\delta
t}\sigma^{+}+e^{i\phi(t)}e^{-i\delta t}\sigma^{-})$ (13)

where $\phi(t)-\delta t=\frac{1}{2}R(t-t_{\delta})^{2}+\delta
t_{0}-\delta^{2}/2R$.

Assuming a near perfect adiabatic chirp, the unitary evolution operator of the
full chirp process is given in the $\sigma_{z}$ eigenstate basis:

$U=\left(\begin{array}[]{cc}0&-ie^{-i\theta_{\delta}}\\\
-ie^{i\theta_{\delta}}&0\\\ \end{array}\right)$ (14)

where $\theta_{\delta}=\phi_{W}-\delta t_{0}$. $\phi_{W}$ contains a term
$\delta^{2}/2R$ and a term that does not depend on $\delta$ but will in
general depend on the chirp rate $R$ and on the coupling strength $\Omega_{0}$
(see below).

For initial spin states with a (small) excitation amplitude of $\epsilon$ at
$t=0$ , the $\delta$-dependence of the mean value of the $\sigma^{+}$ operator
is $\propto\epsilon e^{2i\delta t_{0}}$ (disregarding the contribution from
$\phi_{W}$), and in the Schrödinger picture, multiplying by $e^{-i\delta t}$
we obtain $\sigma^{+}\propto\epsilon e^{i\delta(2t_{0}-t)}$. If it was not for
the detuning dependence of $\phi_{W}$, these terms would rephase at the time
$t=2t_{0}$ as in the conventional $\pi-$ pulse echo. Now, instead, the spins
may have very different excitation phases and the integral over $\delta$
vanishes and the echo is silenced.

All WURST pulses are described by the same form (14), and the action of two
subsequent pulses (in the interaction picture) at times $t_{0}$ and $t_{1}$ is
readily found,

$U=U_{1}U_{0}=\left(\begin{array}[]{cc}-a^{*}b&0\\\ 0&-ab^{*}\\\
\end{array}\right)$ (15)

With $a=e^{-i\theta_{\delta}}$, with parameters $R_{0},\Omega_{0}$ and
$b=e^{-i\theta_{\delta}}$, with parameters $R_{1},\Omega_{1}$, we get
$ab^{*}=e^{i(\phi_{W_{0}}-\phi_{W_{1}})+i\delta(t_{1}-t_{0})}$ and if the
pulse parameters are identical, this simplifies to
$ab^{*}=e^{i\delta(t_{0}-t_{1})}$. Thus, after two identical WURST pulses, a
weak spin excitation returns with a global minus sign and a detuning dependent
phase factor. In the Schrödinger picture, the linear detuning dependence of
the phase disappears at $t=2(t_{1}-t_{0})$ and we obtain an echo.

Assuming perfect inversion, a sequence with an even number of WURST pulses
yields a final phase which is the sum of all WURST phases $\theta_{\delta}$
with $+(-)$ signs if they are raising (lowering) pulses, cf., (13). Two pulses
with different chirp rates or strengths will not cause an echo due to the
nontrivial dependence of the phase difference $\phi_{W_{0}}-\phi_{W_{1}}$. The
quadratic variation $\delta^{2}(1/2R_{0}-1/2R_{1})$ due to the phase chirp may
thus cause destructive interference over the distribution of $\delta$ and the
dynamical phase associated with the energies of the adiabatic eigenstates
during the chirped adiabatic process may cause destructive interference over a
range of values for the coupling strength.

As we observe in the experiments, variation of the chirp rate and strength
permit addressing of separate storage modes. The field amplitude $\Omega_{0}$,
and the chirp rate $R$ both appear in the dressed state energies and hence in
the dynamical phases, and the two control parameters turn out not lead to
exploration of a two dimensional space of storage modes . An analytical theory
for the phase textures indicated in Fig.1(b) would require solution of the
driven two-level dynamics by a chirped frequency field with a time and
frequency dependent Rabi frequency $\Omega_{0}(t)$ imposed by the finite
bandwidth cavity mode. While this is not possible, we can solve the dynamics
numerically and do so to elucidate two effects shown in Fig.
LABEL:fig:theory_efficiency and Fig. LABEL:F2 respectively.

###  I.6 High power memory sequences

Figure S9: The memory sequence shown in Fig.3 was performed at higher input
power ($\sim 15000$ photons per excitation pulse); the raw transient trace of
the sequence is shown in (e). In Figs.(a-d) the same experiment was repeated
with only one of each of the 4 input excitation pulses, A,B,C and D turned on.
This allows us to unambiguously determine which echo originates from which
excitation pulse. We see that a single echo appears in the expected read slot
in each case. We also see that no other echoes appear when only one excitation
is turned on. In this power regime, as well as amplifying (non-linearly), the
JPA also imparts poorly characterized phase shifts seen particularly when
comparing (c,d,e).

To further confirm the memory protocol functions as intended, we ran identical
memory protocols cycling different excitations on and off so that in each run
only one excitation was used. To speed up this measurement we increased the
power of the input pulses and the results are shown in Fig. S9 where we also
show a replication of the full memory protocol at higher power. In each of the
sequences where there is only one excitation stored, we only retrieve one
echo, and the echo occurs when we would expect the echo to form further
confirming that the memory protocol works.


# Searching for Triple Systems Unbound by Supernovae

Karina Barboza,1 Christopher S. Kochanek,1,2
1Department of Astronomy, The Ohio State University, 140 West 18th Avenue,
Columbus OH 43210, USA
2Center for Cosmology and AstroParticle Physics, The Ohio State University,
191 W. Woodruff Avenue, Columbus OH 43210, USA E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

A large fraction of massive stars are found in higher order systems where the
presence of a tertiary may significantly modify the system’s evolution. In
particular, it can lead to increased numbers of compact object binaries and
accelerate their mergers with important implications for gravitational wave
observations. Using Gaia, we constrain the number of Galactic supernovae that
produce unbound triples. We do this by searching 8 supernova remnants for
stars with consistent Gaia parallaxes and paths intersecting near the center
of the supernova remnant at a time consistent with the age of the remnant. We
find no candidates for unbound triple systems. Combined with prior work, less
than 11.4% of supernovae leave behind unbound triples at a 90% confidence
limit. The absence of such systems limits their role in the evolution of
massive stars and the formation of merging compact objects.

###### keywords:

supernova remnants – triples – compact objects

††pubyear: 2023††pagerange: Searching for Triple Systems Unbound by
Supernovae–Searching for Triple Systems Unbound by Supernovae

##  1 Introduction

Most massive stars have at least one stellar companion (e.g., Evans et al.,
2006; Kobulnicky & Fryer, 2007; Sana et al., 2014; Moe & Stefano, 2017) and
are expected to interact at some point during their evolution (Sana et al.,
2012). For massive stars, triple systems are likely more common than binaries
(Moe & Stefano, 2017). While these surveys have provided valuable insights
into the multiplicity of massive stars, the path to mergers or interacting
compact objects is long and intricate. To understand the expected properties
of supernovae (SN) and compact objects, we must understand the evolution of
such systems.

Binary population synthesis (BPS) models are used to try to understand and
predict the evolution of binary populations (e.g., Belczynski et al., 2008;
Eldridge et al., 2017). Recent expansions in BPS studies now encompass the
tracking of triple or quadruple systems with their significantly richer
dynamical properties (e.g., Toonen et al., 2020; Hamers et al., 2021). The
introduction of a tertiary to a system’s evolution can lead to enhanced tidal
interactions (Fabrycky & Tremaine, 2007), mass transfer in eccentric orbits
(Thompson, 2011; Shappee & Thompson, 2013; Toonen et al., 2020), increase the
number of compact object binaries (Thompson, 2011; Shappee & Thompson, 2013)
and accelerate their mergers (Perets & Fabrycky, 2009). These accelerated
mergers can then produce transients such as gamma-ray bursts, Type Ia
supernovae (SNe), and gravitational wave mergers (Wen, 2003; Perets &
Fabrycky, 2009; Thompson, 2011; Katz & Dong, 2012; Hamers et al., 2013;
Kimpson et al., 2016; Zwart & van den Heuvel, 2016; Abbott et al., 2021).

Aside from the standard issues in stellar evolution, there are complexities in
BPS models such as mass transfer, common envelope evolution, neutron star and
black hole natal kicks (see, e.g., the review of BPS models by Han et al.,
2020), and understanding the connections between a core collapse and its
outcomes (supernovae, failed supernovae, neutron stars, or black holes see,
e.g., Pejcha & Thompson, 2015; Kochanek, 2014; Sukhbold et al., 2016). These
physical effects have significant uncertainties, but they are crucial to
predictions, especially for the population of merging compact objects. One way
to reduce the uncertainties is to determine the statistics of binary and
higher order systems at the time of stellar death rather than birth. This can
be done through studies of Galactic SNRs.

For SNRs with known compact objects, finding bound systems is relatively
straight forward, as the companion and the remnant are co-located and, in the
three known cases, interacting (Hillwig & Gies, 2008, SS433; HESS J0632$+$057
in MWC 148, Hinton et al. 2009; and 1FGL J1018.6$-$5856 in 2MASS
J10185560$-$5856489, Corbet et al. 2011). Kochanek et al. (2019) searched for
surviving but non-interacting binaries and found none in 23 SNRs.

Detecting unbound binaries or higher order systems in supernova remnants
(SNRs) is challenging because the remnant and the former companion are not co-
located. This can become even more challenging in the case where there is no
identified remnant, as illustrated by the many studies searching for a former
stellar companion to the Type Ia SN Tycho (e.g., Ruiz-Lapuente, 2004;
Kerzendorf et al., 2018; Ihara et al., 2007). Typically, the former companion
should remain close to the center of the SNR. However, the distortion and
asymmetry of many SNRs makes it challenging to accurately determine the search
region from the structure of the remnant. Boubert et al. (2017) and Lux et al.
(2021) have attempted such searches for core-collapse SNRs, but this approach
likely has a high false positive rate (see Kochanek, 2021).

In the case where a remnant is identified and the SNR is very young, the
process becomes straightforward because a surviving companion must be nearly
co-located with the remnant, as in the cases of the Crab, Cas A, and SN 1987A
(Kochanek, 2018). If the remnant has a measured proper motion (and preferably
a parallax), this can be combined with stellar proper motions and parallaxes
to search for companions whose back projected trajectories intersect with that
of the remnant near the center of the SNR at a reasonable age. This was done
by Dinçel et al. (2015), and expanded by Fraser & Boubert (2019), Kochanek
(2021), and Kochanek (2023). To date, there is only one good candidate for an
unbound binary companion, HD 37424 in S147 (Dinçel et al. 2015). The present
statistics are that 77% (58-90% at 90% confidence) are not binaries at death,
9% (4-19%) are bound binaries, 6% (2-13%) are interacting binaries, $<9\%$ are
non-interacting binaries and 12% (3-31%) are unbound binaries (Kochanek 2023).

Finding unbound triples is actually easier than finding unbound binaries
because it is not dependent on even identifying the remnant, let alone having
its proper motions or parallaxes. One can simply look for a pair of stars with
consistent parallaxes whose back projected trajectories intersect with the
remnant near the center of the SNR. This was first done by Kochanek (2021),
who studied 10 SNRs and found no candidates, which implied that fewer than
18.9% core-collapse SN lead to unbound triples at 90% confidence. Here we
expand this search for unbound triples to 8 additional SNRs using Gaia DR3. We
outline the method in §2, discuss the 8 SNRs and the search results in §3 and
summarize them in §4.

##  2 Identifying Stellar Companions

For each SNR, we search Gaia DR3 (Gaia Collaboration et al., 2023, 2016) for
stars within a radius $R/R_{SNR}<1$ where $R_{SNR}$ is the radius of the SNR
and with a magnitude limit $G_{lim}$ chosen to reach stars with $M\gtrsim
M_{\odot}$ given the distance and extinction to the SNR. We discuss the
adopted distances and extinctions in §3. We identify candidate triples using
the methods from Kochanek (2021), outlined in this section and summarized in
Table 3. While all the SNRs contain neutron stars, we primarily search for
possible unbound pairs of stars.

First we restrict our candidates to those that can intersect within the
neutron star’s path in the time frame $-t_{m}<t<0$ for a neutron star velocity
of $|v_{ns}|<v_{max}$. The positional difference between the neutron star and
the star can be written as
$\Delta\overrightarrow{\text{x}}=d\Delta\overrightarrow{\alpha}$ for a
distance $d$ where
$\Delta\overrightarrow{\alpha}=\overrightarrow{\alpha}_{ns}-\overrightarrow{\alpha}_{*}$,
the velocity of the star is
$\overrightarrow{v}_{*}=d\overrightarrow{\mu}_{*}$, so the neutron star
velocity ($v_{ns}$) needed to intercept the compact object at time t is

$\centering
v^{2}_{ns}=v^{2}_{*}+2t^{-1}\overrightarrow{v}_{*}\cdot\Delta\overrightarrow{x}+\Delta
x^{2}t^{-2}.\@add@centering$ (1)

The angular coordinate differences incorporate the necessary cos$\delta$ term
for the separation in RA and the angular separations are small enough to just
work in the tangent plane. We keep the stars which can intercept the remnant
at $R<R_{SNR}/4$ with $|v_{ns}|\leq v_{max}=400$ km/s and $-10^{5}\leq t\leq
0$ years.

Next we search for pairs of these stars with consistent parallax estimates. If
the parallaxes are $\varpi_{i}\pm\sigma_{\varpi,i}$ where $i=1,2$, then their
joint parallax is

$\centering\varpi=\left[\sum_{i}\frac{\varpi_{i}}{\sigma^{2}_{\varpi,i}}\right]\left[\sum_{i}\frac{1}{\sigma^{2}_{\varpi,i}}\right]^{-1}\@add@centering$
(2)

with a goodness of fit

$\centering\chi^{2}_{\varpi}=\sum_{i}\left(\frac{\varpi_{i}-\varpi}{\sigma^{2}_{\varpi,i}}\right)^{2}.\@add@centering$
(3)

Pairs meeting the requirement $\chi^{2}_{\varpi}<9$ are kept as potential
candidates for unbound triple systems.

The paths of a pair of unbound stars need not precisely intersect at a
particular time because the disrupted binaries had a finite size. We introduce
a maximum binary separation $a_{max}$ which corresponds to an angle of

$\centering\sigma_{a}=10\aas@@fstack{\prime\prime}0\left(\frac{a_{max}}{10^{4}\text{AU}}\right)\left(\frac{\varpi}{\text{mas}}\right)\@add@centering$
(4)

We must also account for the uncertainties in the proper motions. We do this
by treating the true proper motions as an intermediate variable. For example,
for the motion in RA we can let the current separation of two objects be
$\Delta\alpha_{12}=\alpha_{1}-\alpha_{2}$, their measured current motions be
$\mu_{\alpha,i}$, and their true proper motions be $\mu^{T}_{\alpha,i}$. The
cosine declination terms on $\Delta\alpha$ are again implicit. The measured
proper motions impose constraints on the true proper motion, leading to the
fit statistic

$\displaystyle\chi^{2}_{\alpha,12}=\frac{(\mu_{\alpha,1}-\mu^{T}_{\alpha,1})^{2}}{\sigma^{2}_{\alpha,1}}+\frac{(\mu_{\alpha,2}-\mu^{T}_{\alpha,2})^{2}}{\sigma^{2}_{\alpha,2}}+$
(5)
$\displaystyle\frac{\left(\Delta\alpha_{12}+t\left(\mu^{T}_{\alpha,1}-\mu^{T}_{\alpha,2}\right)\right)^{2}}{\sigma^{2}_{\alpha}}.$

Here the first two terms use the measured proper motions to constrain the true
proper motions, while the last term constrains the closest approach distance
relative to the maximum binary size based on the true proper motions. We can
then optimize the fit statistic for the unknown true proper motions to obtain

$\displaystyle\chi^{2}_{\alpha,12}=\frac{(\Delta\alpha_{12}/t+\Delta\mu_{\alpha,12})^{2}}{\sigma^{2}_{\alpha,12}}$
(6)

where $\Delta\mu_{\alpha,12}=\mu_{\alpha,1}-\mu_{\alpha,2}$, and

$\displaystyle\sigma^{2}_{\alpha,12}=\sigma^{2}_{\alpha,1}+\sigma^{2}_{\alpha,2}+\sigma^{2}_{\alpha}/t^{2}.$
(7)

We now define $\Delta_{\delta,12}=\Delta\delta_{12}/t+\Delta\mu_{\delta,12}$
to obtain

$\displaystyle\chi^{2}_{\mu}=\frac{\Delta^{2}_{\alpha,12}}{\sigma^{2}_{\alpha,12}}+\frac{\Delta^{2}_{\delta,12}}{\sigma^{2}_{\delta,12}},$
(8)

with the declination term. A value of $\chi^{2}_{\mu}\simeq 3$ roughly
corresponds to the stars passing within a distance $a_{max}$ of each other
with proper motions that are each within $1\sigma$ of their measured values.

Because of equation (7), the optimization of $\chi^{2}_{\mu}$ is non-linear in
$t^{-1}$. To solve this problem, we keep the time in equation (7) fixed to
determine a new optimal time, use this new time to update the equation, and
repeat the process to convergence. The positional uncertainties are excluded
because they are unimportant, but they could be incorporated as an error term
such as $\sigma_{\alpha}$.

With the uncertainties fixed, the closest approach time is

$t_{m}=-\frac{\Delta\alpha_{12}^{2}\sigma_{\delta,12}^{2}+\Delta\delta_{12}^{2}\sigma_{\alpha,12}^{2}}{\Delta\alpha_{12}\Delta\mu_{\alpha,12}\sigma_{\delta,12}^{2}+\Delta\delta_{12}\Delta\mu_{\delta,12}\sigma_{\alpha,12}^{2}}$
(9)

with a goodness of fit of

$\displaystyle\chi_{12}^{2}=\frac{(\Delta\delta_{12}\Delta\mu_{\alpha,12}-\Delta\alpha_{12}\Delta\mu_{\delta,12})^{2}}{\Delta\alpha_{12}^{2}\sigma_{\delta,12}^{2}+\Delta\delta_{12}^{2}\sigma_{\alpha,12}^{2}}.$
(10)

We require $\chi^{2}_{\mu}<9$, an intercept time of $-10^{5}\leq t_{m}\leq 0$
years and an intersection point whose radius relative to the SNR center
satisfies $R/R_{SNR}<1/4$. We can then check the consistency of the candidates
with the estimated age and distance of the SNR. Figure 1 illustrates a model
of a good candidate. The dashed circle is the SNR search radius $R_{SNR}/4$,
the red and green stars are the current positions of the stars, and the orange
circle is the intersection point. This is a model candidate because we find no
pairs meeting all the selection criteria.

We employed PARSEC isochrones from Pastorelli et al. (2020) to categorize the
stars into different mass ranges, based on their positions on an extinction-
corrected $M_{G}$ versus $B_{P}-R_{P}$ color-magnitude diagram. This is the
only aspect of the calculation involving distances, and employing the inverse
parallax is sufficiently significantly accurate because of the steepness of
stellar mass-luminosity relations. We grouped stars by applying limits of 1,
2, 3, 4 and 5$M_{\odot}$, with the color-dependent absolute magnitude limits
given in Table 2. For example a star with a $B_{P}-R_{P}<0$ and $M_{G}<-0.4$
should have a mass $>5M_{\odot}$.

Figure 1: An example of a pair which could be a good candidate based on the
geometry of the SNR G$065.7+01.2$. The current positions of the stars are
given by the red and green stars and the position at the optimal intersection
time is given by the orange circle. The blue cross mark is the estimated
center of the SNR, the magenta plus sign is the location of the remnant, and
the dotted circle is the $R_{SNR}/4$ region we require for an acceptable
candidate. Unfortunately, we found no candidates which met all the
requirements.

##  3 Sample and Results

From the 23 SNRs in Table 2 of Kochanek et al. (2019), we examined eight:
G065.7$+$01.2, G078.2$+$02.1, G106.3$+$02.7, G114.3$+$00.3, G119.5$+$10.2,
G189.1$+$03.0, G266.2$-$01.2, G291.0$-$00.1. The positions and radii of these
SNRs are shown in Table 1. Of the remaining 15 remnants, 12 were previously
studied in Kochanek (2021) and Kochanek (2017), and the remaining 3
(G034.7$-$00.4, G320.4$-$01.2, G332.4$-$00.4) required using such faint stars
that the searches were dominated by false positives. We adopt the SNR centers
and diameter from Green (2019). For the stellar extinctions we used mwdust
(Bovy et al., 2016), which is based on the which combine the Drimmel et al.
(2003), Marshall et al. (2006), and Green (2019) dust models. These models
generally agree well with X-ray estimates of the column density N(H) for the
distance of the SNR.

The selection process is outlined in Table 3. For the first step, we queried
GAIA DR3 Gaia Collaboration et al. (2023) to find all the stars within the
radius of the SNR where the limiting magnitude $G_{lim}$ is chosen based on
the distance and extinction of the remnant. From these stars, in step 2 we
selected stars whose paths intersect the neutron star with
$R/R_{SNR}<1/4,|v_{ns}|<400$ km/s and $-10^{5}\leq t\leq 0$ years (equation
1). As shown in Table 3, this is where we lost most of the initial sample.
Next in step 3, we determined which pairs of stars had consistent parallaxes
$\chi_{\varpi}^{2}<9$ (equations 2 and 3). For step 4, we determined which of
these pairs of stars could intercept with $\chi_{\mu}^{2}<9$ (equation 8).
Finally, our last constraint was that star pairs intersect within $R_{SNR}/4$
and at a time between $-10^{5}<t_{m}<0$ years, which none of our candidates
satisfied.

As we tested the search, we discovered two pairs which had present day
separations of $2.6\times 10^{4}$ AU and $4.9\times 10^{5}$ AU with
statistically consistent parallaxes and proper motions. These are almost
certainly wide binaries. We then searched all the pairs for cases which have
present day separations $<10^{6}$ AU and statistically consistent parallaxes
and proper motions. We defined statistically consistent as
$\chi^{2}_{b}=\chi^{2}_{\varpi}+\chi^{2}_{\mu,\alpha}+\chi^{2}_{\mu,\delta}<9$
where $\chi^{2}_{\varpi}$ is simply Eqn 3 and $\chi^{2}_{\mu,\alpha}$ and
$\chi^{2}_{\mu,\delta}$ are the analogous fit statistics for the two
components of the proper motions. Of all pairs found in step 3, 11 have a
present day separation $<10^{6}$ AU with $\chi^{2}_{b}<9$.

SNR G065.7+01.2 contains the resolved pulsar wind nebula (PWN) DA 495 and
central object J1952.2+2925 (Arzoumanian et al., 2004). Kothes et al. (2008)
calculate the age of the remnant to be $\sim 20,000$ yrs. The first star pair
candidate is located at 3.3 kpc, a distance consistent with that of the
remnant. The distance to the second star pair candidate is 6.5 kpc which is
not consistent with the remnant. While these two candidates meet the two
$\chi^{2}$ requirements, neither candidate intersects on a reasonable time
scale or close to the remnant’s center. Both of these star pair candidates are
almost certainly wide binaries since they have present day separations of
$3.1\times 10^{5}$ AU and $4.9\times 10^{5}$ AU respectively, and
statistically consistent parallaxes and proper motions.

The SNR G078.2+02.1 is associated with the $\gamma$-ray pulsar PSR J2021+4026
(Abdo et al., 2009). Trepl et al. (2010) estimated an age of 6.6 kyr which is
in agreement with the age from Mavromatakis (2003) and an order of magnitude
less than the $\tau\sim 77$ kyr spin-down age of the pulsar (Trepl et al.,
2010). Of the 21 candidates found in step 2, none met the two $\chi^{2}$
requirements.

SNR G106.3+02.7 contains the radio pulsar PSR J2229+6114. Halpern et al.
(2001b) determine the characteristic age of the pulsar as 10,460 years. None
of the 11 candidates from step 2 met the two $\chi^{2}$ requirements.

The remnant G114.3+00.3 is associated with radio pulsar PSR B2334+61 (Dewey et
al., 1985).Cioffi et al. (1988) estimate the age of the remnant to be 7700
years. None of the stars from our initial search made it past step 2.

G119.5+10.2 (CTA 1) contains the X-ray source RX J0007.0+7303 (Halpern et al.,
2004), later found to be a $\gamma$-ray (Abdo et al., 2009) and X-ray pulsar
(Lin et al., 2010). The pulsar has a spin-down age of $\tau_{c}\sim 14,000$
yrs which is consistent with the age of the remnant (10,000–20,000 yrs, Lin et
al., 2010). None of the 3 candidates from step 2 met the two $\chi^{2}$
requirements.

The SNR G189.1+03.0 (IC 443) contains the neutron star CXO J061705.3+222127
and the PWN G189.22+2.90 (Keohane et al., 1997). Chevalier (1999) estimate an
age of $\sim 30,000$ years for the remnant. None of the 4 candidates from step
2 met the two $\chi^{2}$ requirements.

SNR G266.2-01.2 (Vela Jr.) contains the X-ray source CXOU J085201.4-461753
(Pavlov et al., 2001). The characteristic age of the remnant is 140 kyr
according to Kramer et al. (2003) while Allen et al. (2014) estimate an age of
2.4-5.1 kyr. The one candidate left after step 2 did not meet the two
$\chi^{2}$ requirements.

G291.0-00.1 is associated with the X-ray pulsar CXOU J111148.6-603926 (Slane
et al., 2012). The age of this remnant is estimated to be 1300 years (Harrus
et al., 1998). None of the stars from our initial search made it past step 2.

##  4 Conclusions

We searched for unbound triple systems in the 8 SNRs listed in Table 1 and
found none. Combined with the 10 systems considered by Kochanek (2021), we
have now searched M = 18 SNRs and found no good candidates for unbound
triples. This limits the fraction $f_{T}$ of triple systems becoming unbound
at stellar death to the binomial probability

$\centering P(<f_{T})=1-(1-f_{T})^{1+M}.\@add@centering$ (11)

Thus, less than $f_{T}<11.4$% of SNe leave behind unbound massive triples at a
90% confidence limit. This limit along with the measurements of bound and
unbound binary fractions from Kochanek et al. (2019), Kochanek (2021), and
Kochanek (2023) provide constraints on multiplicity that should inform and
test BPS models. Recently, Burdge et al. (2024) found that the black hole low
mass X-ray binary V404 Cygni was likely a triple system, and discuss the
challenges of keeping such a system bound. We note that these limits on binary
systems effectively apply to the abundance of surviving triple systems as
well.

## Acknowledgements

CSK is supported by the NSF grants AST-1908570 and AST-2307385. KB is
supported by the Heising-Simons Foundation through grant number 2022-3533. We
are thankful to be able to think and work on the ancestral and contemporary
territory of the Shawnee, Potawatomi, Delaware, Miami, Peoria, Seneca,
Wyandotte, Ojibwe, and many other Indigenous peoples.

## Data Availability

All data used in this paper are publically available.

Table 1: SNR candidates and their properties: The column densities N(H) are,
in order, from Arzoumanian et al. (2008), Hui et al. (2015), Halpern et al.
(2001a), McGowan et al. (2006), Slane et al. (1997), Gaensler et al. (2006),
Pavlov et al. (2001), Slane et al. (2012). Distance estimates are from Karpova
et al. (2016), Landecker et al. (1980), Halpern et al. (2001b), McGowan et al.
(2006), Pineault et al. (1993), Fesen (1984), Liseau et al. (1992), Slane et
al. (2012).

SNR Center $R_{SNR}$ d N(H) $G_{lim}$ Remnant (’) (kpc) ($cm^{-2}$) (mag)
G$065.7+01.2$ 19:52:10 $+$29:26:00 $11$ $3.0\pm 2.0$ $(2-7)\times 10^{21}$
$15.5$ PSR J$1952.2+2925$ G$078.2+02.1$ 20:20:50 $+$40:26:00 $30$ $1.5\pm 0.5$
$6.4^{+0.8}_{-1.8}\times 10^{21}$ $15$ PSR J$2021+4026$ G$106.3+02.7$ 22:27:30
$+$60:50:00 $19$ 3.0 $6.3\pm 1.3\times 10^{21}$ $16$ PSR J$2229+6114$
G$114.3+00.3$ 23:37:00 $+$61:55:00 $35$ $3.2\pm 1.7$ $(0.26-1.0)\times
10^{22}$ $17$ PSR B$2334+61$ G$119.5+10.2$ 00:06:40 $+$72:45:00 $45$ $1.4\pm
0.3$ $2.8^{+0.6}_{-0.5}\times 10^{21}$ $12$ RX J$0007.0+7303$ G$189.1+03.0$
06:17:00 $+$22:34:00 $23$ $1.7\pm 0.3$ $(7.2\pm 0.6)\times 10^{21}$ $14.5$ PWN
G$189.22+2.90$ G$266.2-01.2$ 08:52:00 $-$46:20:00 $60$ $0.7\pm 0.2$ $(3.0\pm
1.0)\times 10^{21}$ $11$ CXOU J$085201.4-461753$ G$291.0-00.1$ 11:11:54
$-$60:38:00 $7$ $3.5-5.0$ $(6.7\pm 0.7)\times 10^{21}$ $16.5$ CXOU
J$111148.6-603926$

Table 2: CMD Selection Limits: Mag gives the upper limit on $M_{G}$ for stars above the listed mass. Mass | Color | Mag | Color | Mag | Color | Mag | Color | Mag
---|---|---|---|---|---|---|---|---
$>1M_{\odot}$ | $<1.5$ | $\hphantom{-}5.5$ |  |  |  $1.5$-$2.0$ | $-1-3(B_{P}-R_{P}-1)$ | $>2$ | $-4$
$>2M_{\odot}$ | $<0.7$ | $\hphantom{-}2.0$ |  $0.7$-$1.1$ | $\hphantom{-}3.0$ |  $1.1$-$2.0$ | $-1-3(B_{P}-R_{P}-1)$ | $>2$ | $-4$
$>3M_{\odot}$ | $<0.7$ | $\hphantom{-}1.0$ |  $0.7$-$1.1$ | $\hphantom{-}2.0$ |  $1.1$-$2.0$ | $-1-3(B_{P}-R_{P}-1)$ | $>2$ | $-4$
$>4M_{\odot}$ | $<0.0$ | $\hphantom{-}0.5$ |  $0.0$-$1.3$ | $-1.9$ |  $1.3$-$2.0$ | $-1-3(B_{P}-R_{P}-1)$ | $>2$ | $-4$
$>5M_{\odot}$ | $<0.0$ | $-0.4$ |  $0.0$-$1.3$ | $-1.9$ |  $1.3$-$2.0$ | $-1-3(B_{P}-R_{P}-1)$ | $>2$ | $-4$
Mag gives the upper limit on $M_{G}$ for each $B_{P}-R_{p}$ color range.
Table 3: Star pair candidates during the selection process.

Steps G$065.7+01.2$ G$078.2+02.1$ G$106.3+02.7$ G$114.3+00.3$ G$119.5+10.2$
G$189.1+03.0$ G$266.2-01.2$ G$291.0+00.1$ 1 initial stars stars 700 554 392
6336 132 252 88 1090 2 intersect NS stars 92 21 11 0 3 4 1 0 3 consistent
$\varpi$ pairs 62 0 0 0 0 0 0 0 4 intersect pairs 2 0 0 0 0 0 0 0 5 time/
radius pairs 0 0 0 0 0 0 0 0

## References

  * Abbott et al. (2021) Abbott R., et al., 2021, The Astrophysical Journal Letters, 915, L5
  * Abdo et al. (2009) Abdo A. A., et al., 2009, ApJS, 183, 46
  * Allen et al. (2014) Allen G. E., Chow K., DeLaney T., Filipović M. D., Houck J. C., Pannuti T. G., Stage M. D., 2014, The Astrophysical Journal, 798, 82
  * Arzoumanian et al. (2004) Arzoumanian Z., Safi-Harb S., Landecker T. L., Kothes R., 2004, ApJ, 610, L101
  * Arzoumanian et al. (2008) Arzoumanian Z., Safi-Harb S., Landecker T. L., Kothes R., Camilo F., 2008, The Astrophysical Journal, 687, 505
  * Belczynski et al. (2008) Belczynski K., Kalogera V., Rasio F. A., Taam R. E., Zezas A., Bulik T., Maccarone T. J., Ivanova N., 2008, The Astrophysical Journal Supplement Series, 174, 223
  * Boubert et al. (2017) Boubert D., Fraser M., Evans N. W., Green D. A., Izzard R. G., 2017, A&A, 606, A14
  * Bovy et al. (2016) Bovy J., Rix H.-W., Green G. M., Schlafly E. F., Finkbeiner D. P., 2016, The Astrophysical Journal, 818, 130
  * Burdge et al. (2024) Burdge K. B., et al., 2024, arXiv e-prints, p. arXiv:2404.03719
  * Chevalier (1999) Chevalier R. A., 1999, ApJ, 511, 798
  * Cioffi et al. (1988) Cioffi D. F., McKee C. F., Bertschinger E., 1988, ApJ, 334, 252
  * Corbet et al. (2011) Corbet R. H. D., et al., 2011, The Astronomer’s Telegram, 3221, 1
  * Dewey et al. (1985) Dewey R. J., Taylor J. H., Weisberg J. M., Stokes G. H., 1985, ApJ, 294, L25
  * Dinçel et al. (2015) Dinçel B., Neuhäuser R., Yerli S. K., Ankay A., Tetzlaff N., Torres G., Mugrauer M., 2015, Monthly Notices of the Royal Astronomical Society, 448, 3196
  * Drimmel et al. (2003) Drimmel R., Cabrera-Lavers A., López-Corredoira M., 2003, A&A, 409, 205
  * Eldridge et al. (2017) Eldridge J. J., Stanway E. R., Xiao L., McClelland L. A. S., Taylor G., Ng M., Greis S. M. L., Bray J. C., 2017, Publications of the Astronomical Society of Australia, 34
  * Evans et al. (2006) Evans C. J., Lennon D. J., Smartt S. J., Trundle C., 2006, Astronomy &; Astrophysics, 456, 623
  * Fabrycky & Tremaine (2007) Fabrycky D., Tremaine S., 2007, The Astrophysical Journal, 669, 1298
  * Fesen (1984) Fesen R. A., 1984, ApJ, 281, 658
  * Fraser & Boubert (2019) Fraser M., Boubert D., 2019, ApJ, 871, 92
  * Gaensler et al. (2006) Gaensler B. M., Chatterjee S., Slane P. O., van der Swaluw E., Camilo F., Hughes J. P., 2006, The Astrophysical Journal, 648, 1037–1042
  * Gaia Collaboration et al. (2016) Gaia Collaboration et al., 2016, A&A, 595, A1
  * Gaia Collaboration et al. (2023) Gaia Collaboration et al., 2023, A&A, 674, A1
  * Green (2019) Green D. A., 2019, Journal of Astrophysics and Astronomy, 40, 36
  * Halpern et al. (2001a) Halpern J. P., Gotthelf E. V., Leighly K. M., Helfand D. J., 2001a, ApJ, 547, 323
  * Halpern et al. (2001b) Halpern J. P., Camilo F., Gotthelf E. V., Helfand D. J., Kramer M., Lyne A. G., Leighly K. M., Eracleous M., 2001b, The Astrophysical Journal, 552, L125
  * Halpern et al. (2004) Halpern J. P., Gotthelf E. V., Camilo F., Helfand D. J., Ransom S. M., 2004, The Astrophysical Journal, 612, 398
  * Hamers et al. (2013) Hamers A. S., Pols O. R., Claeys J. S. W., Nelemans G., 2013, Monthly Notices of the Royal Astronomical Society, 430, 2262
  * Hamers et al. (2021) Hamers A. S., Rantala A., Neunteufel P., Preece H., Vynatheya P., 2021, Monthly Notices of the Royal Astronomical Society, 502, 4479
  * Han et al. (2020) Han Z.-W., Ge H.-W., Chen X.-F., Chen H.-L., 2020, Research in Astronomy and Astrophysics, 20, 161
  * Harrus et al. (1998) Harrus I. M., Hughes J. P., Slane P. O., 1998, The Astrophysical Journal, 499, 273
  * Hillwig & Gies (2008) Hillwig T. C., Gies D. R., 2008, The Astrophysical Journal, 676, L37
  * Hinton et al. (2009) Hinton J. A., et al., 2009, ApJ, 690, L101
  * Hui et al. (2015) Hui C. Y., et al., 2015, The Astrophysical Journal, 799, 76
  * Ihara et al. (2007) Ihara Y., Ozaki J., Doi M., Shigeyama T., Kashikawa N., Komiyama Y., Hattori T., 2007, Publications of the Astronomical Society of Japan, 59, 811–826
  * Karpova et al. (2016) Karpova A., Zyuzin D., Danilenko A., Shibanov Y., 2016, Journal of Physics: Conference Series, 769, 012002
  * Katz & Dong (2012) Katz B., Dong S., 2012, The rate of WD-WD head-on collisions may be as high as the SNe Ia rate (arXiv:1211.4584)
  * Keohane et al. (1997) Keohane J. W., Petre R., Gotthelf E. V., Ozaki M., Koyama K., 1997, The Astrophysical Journal, 484, 350
  * Kerzendorf et al. (2018) Kerzendorf W. E., Long K. S., Winkler P. F., Do T., 2018, Monthly Notices of the Royal Astronomical Society, 479, 5696
  * Kimpson et al. (2016) Kimpson T. O., Spera M., Mapelli M., Ziosi B. M., 2016, Monthly Notices of the Royal Astronomical Society, 463, 2443
  * Kobulnicky & Fryer (2007) Kobulnicky H. A., Fryer C. L., 2007, The Astrophysical Journal, 670, 747
  * Kochanek (2014) Kochanek C. S., 2014, The Astrophysical Journal, 785, 28
  * Kochanek (2017) Kochanek C. S., 2017, Monthly Notices of the Royal Astronomical Society, 473, 1633
  * Kochanek (2018) Kochanek C. S., 2018, MNRAS, 473, 1633
  * Kochanek (2021) Kochanek C. S., 2021, Monthly Notices of the Royal Astronomical Society, 507, 5832
  * Kochanek (2023) Kochanek C. S., 2023, MNRAS, 519, 3865
  * Kochanek et al. (2019) Kochanek C. S., Auchettl K., Belczynski K., 2019, Monthly Notices of the Royal Astronomical Society, 485, 5394
  * Kothes et al. (2008) Kothes R., Landecker T. L., Reich W., Safi-Harb S., Arzoumanian Z., 2008, The Astrophysical Journal, 687, 516
  * Kramer et al. (2003) Kramer M., et al., 2003, MNRAS, 342, 1299
  * Landecker et al. (1980) Landecker T. L., Roger R. S., Higgs L. A., 1980, A&AS, 39, 133
  * Lin et al. (2010) Lin L. C. C., Huang R. H. H., Takata J., Hwang C. Y., Kong A. K. H., Hui C. Y., 2010, The Astrophysical Journal, 725, L1
  * Liseau et al. (1992) Liseau R., Lorenzetti D., Nisini B., Spinoglio L., Moneti A., 1992, A&A, 265, 577
  * Lux et al. (2021) Lux O., Neuhäuser R., Mugrauer M., Bischoff R., 2021, Astronomische Nachrichten, 342, 553
  * Marshall et al. (2006) Marshall D. J., Robin A. C., Reylé C., Schultheis M., Picaud S., 2006, Astronomy & Astrophysics, 453, 635–651
  * Mavromatakis (2003) Mavromatakis F., 2003, A&A, 408, 237
  * McGowan et al. (2006) McGowan K. E., Zane S., Cropper M., Vestrand W. T., Ho C., 2006, The Astrophysical Journal, 639, 377
  * Moe & Stefano (2017) Moe M., Stefano R. D., 2017, The Astrophysical Journal Supplement Series, 230, 15
  * Pastorelli et al. (2020) Pastorelli G., et al., 2020, MNRAS, 498, 3283
  * Pavlov et al. (2001) Pavlov G. G., Sanwal D., Kızıltan B., Garmire G. P., 2001, The Astrophysical Journal, 559, L131
  * Pejcha & Thompson (2015) Pejcha O., Thompson T. A., 2015, The Astrophysical Journal, 801, 90
  * Perets & Fabrycky (2009) Perets H. B., Fabrycky D. C., 2009, The Astrophysical Journal, 697, 1048
  * Pineault et al. (1993) Pineault S., Landecker T. L., Madore B., Gaumont-Guay S., 1993, AJ, 105, 1060
  * Ruiz-Lapuente (2004) Ruiz-Lapuente P., 2004, The Astrophysical Journal, 612, 357–363
  * Sana et al. (2012) Sana H., et al., 2012, Science, 337, 444
  * Sana et al. (2014) Sana H., et al., 2014, The Astrophysical Journal Supplement Series, 215, 15
  * Shappee & Thompson (2013) Shappee B. J., Thompson T. A., 2013, The Astrophysical Journal, 766, 64
  * Slane et al. (1997) Slane P., Seward F. D., Bandiera R., Torii K., Tsunemi H., 1997, The Astrophysical Journal, 485, 221
  * Slane et al. (2012) Slane P., et al., 2012, The Astrophysical Journal, 749, 131
  * Sukhbold et al. (2016) Sukhbold T., Ertl T., Woosley S. E., Brown J. M., Janka H.-T., 2016, The Astrophysical Journal, 821, 38
  * Thompson (2011) Thompson T. A., 2011, The Astrophysical Journal, 741, 82
  * Toonen et al. (2020) Toonen S., Zwart S. P., Hamers A. S., Bandopadhyay D., 2020, Astronomy &; Astrophysics, 640, A16
  * Trepl et al. (2010) Trepl L., Hui C. Y., Cheng K. S., Takata J., Wang Y., Liu Z. Y., Wang N., 2010, MNRAS, 405, 1339
  * Wen (2003) Wen L., 2003, The Astrophysical Journal, 598, 419
  * Zwart & van den Heuvel (2016) Zwart S. F. P., van den Heuvel E. P. J., 2016, Monthly Notices of the Royal Astronomical Society, 456, 3401