# A Joint X-ray and Optical Study of the Massive Redback Pulsar J2215+5135

Andrew G. Sullivan Kavli Institute for Particle Astrophysics and Cosmology,
Department of Physics, Stanford University, Stanford, CA 94305, USA  Roger W.
Romani Kavli Institute for Particle Astrophysics and Cosmology, Department of
Physics, Stanford University, Stanford, CA 94305, USA

###### Abstract

PSR J2215+5135 (J2215) is a ‘redback’ spider pulsar, where the intrabinary
shock (IBS) wraps around the pulsar rather than the stellar-mass companion.
Spider orbital light curves are modulated, dominated by their binary companion
thermal emission in the optical band and by IBS synchrotron emission in the
X-rays. We report on new XMM-Newton X-ray and U-band observations of J2215. We
produce orbital light curves and use them to model the system properties. From
our optical light curve modeling, we find a neutron star mass $M_{NS}=2.15\pm
0.10M_{\odot}$, somewhat lower than previously reported. From the X-ray
analysis, we find that the IBS still wraps around the pulsar, but with a
wind/companion wind momenta ratio unusually close to unity, implying a flatter
IBS geometry than seen in other spiders. Estimating the companion wind
momentum and speed from the X-ray light curve, we find a companion mass-loss
rate of ${\dot{M}}_{c}\gtrsim 10^{-10}$ M⊙ yr-1, so that J2215 may become an
isolated millisecond pulsar in $\sim 1$ Gyr. Our X-ray analyses support models
of magnetic reconnection and particle acceleration in the highly magnetized
relativistic IBS.

Pulsars (1306) – Binary pulsars (153)

##  1 Introduction

Spider pulsars found in compact binaries with a low-mass companion star have
orbital periods $P_{b}\lesssim 1$ day. The companion mass determines whether
the spider falls into the black widow ($M_{c}<0.1$ M⊙) or redback
($M_{c}\approx 0.1-0.4$ M⊙) category. In these systems, pulsar gamma-ray
emission and the relativistic particles from the pulsar wind irradiate the
companion, consequently driving off a massive stellar wind (Kluzniak et al.,
1988; van Paradijs et al., 1988; Roberts, 2013; Hui & Li, 2019). The pulsar
wind and companion wind collide to form an intrabinary shock (IBS). In
redbacks which typically have higher stellar wind momentum fluxes, the
companion wind dominates the pulsar wind so the IBS wraps around the pulsar,
while in black widows with lower stellar wind momentum fluxes, the IBS wraps
around the companion, as the pulsar wind dominates the stellar wind (Romani &
Sanchez, 2016; Wadiasingh et al., 2017; Kandel et al., 2019).

X-ray observations may probe the IBS while optical observations typically
reveal the companion heating. Particles in the shocked pulsar wind accelerate
and emit synchrotron X-rays (Kandel et al., 2019, 2021). The light curves are
characterized by two caustic peaks per orbital phase, associated with the
beamed emission from relativistic particles traveling tangent to the observing
line of sight. Thermal emission from the companion is generally boosted by re-
processed pulsar gamma-rays on the ‘day’ side of the companion and sometimes
by precipitating IBS particles reprocessed in surface hot spots. Consequently,
X-ray light curve and spectral analyses reveal the IBS structure and particle
acceleration mechanisms, while optical analyses examine heating of the
companion surface.

PSR J2215+5135 (J2215 hereafter) is a redback millisecond pulsar (MSP) with
spin period $P_{s}=2.61$ s and spin-down power $\dot{E}=5\times 10^{34}$
$I_{45}$ erg s-1 in a $P_{b}=4.14$ hr orbit with a $M\approx 0.3$ M⊙ companion
(Hessels et al., 2011). The source has dispersion measure of 69.2 pc cm-3
(Hessels et al., 2011), which corresponds to an estimated distance $\sim 3$
kpc (Yao et al., 2017, e.g.). First detected as an unidentified gamma-ray
source by the Fermi Large Area Telescope (Fermi-LAT) (Atwood et al., 2009) and
discovered by 350 MHz Green Bank Telescope followup observations (Ray et al.,
2012; Hessels et al., 2011), J2215 has been studied optically (Breton et al.,
2013; Schroeder & Halpern, 2014; Romani et al., 2015; Linares et al., 2018),
in radio (Hessels et al., 2011; Broderick et al., 2016), in gamma-rays (Romani
et al., 2015), and briefly in X-rays (Gentile et al., 2014; Linares, 2014)
prior to this work.

Optical light curve modeling suggests that J2215 contains a particularly
massive neutron star $\gtrsim 2$ M⊙ (Linares et al., 2018; Kandel & Romani,
2020). Kandel & Romani (2020) estimate a neutron star mass of $M_{NS}=2.24\pm
0.09$ M⊙ using a companion surface optical model that contains a magnetic pole
hot spot. Such a massive neutron star constrains properties of the dense
matter equation of state (Lattimer & Prakash, 2007; Steiner et al., 2013;
Brandes et al., 2023). Previous modeling is sensitive to the poorly
constrained companion hot spot, adding uncertainty to the binary inclination
$i$ and, by extension, $M_{NS}$ estimates. Also, previous X-ray observations,
which suggested widely separated light curve peaks and a flat IBS, relied on
an XMM-Newton exposure with strong background flaring and only two binary
orbits. These factors motivate further optical and X-ray study.

In this paper, we present the results of light curve analyses of new XMM-
Newton J2215 X-ray and U band observations. In sec. 2, we summarize the new
data, presenting the X-ray and U light curves. In sec. 3, we discuss our
analysis of the collected X-ray spectra. In sec. 4, we show the results of
optical companion heating and X-ray IBS light curve analyses. We discuss our
results and their implications in sec. 5.

##  2 XMM-Newton Observations

The XMM-Newton Observatory (Jansen et al., 2001) performed three new
observations of J2215, the first two on 2022 June 8-11 (ObsIDs 0900770101 and
0900770201) and the third on 2022 December 1 (ObsID 0900770301) lasting a
total of 120.7 ks. These observations provided both X-ray and U band data,
which are presented here for the first time. We supplement these new
observations with an archival 54.9 ks XMM-Newton X-ray observation (ObsID
0783530301). We perform data reduction using tools of the XMM-Newton
Scientific Analysis System (SAS) (Jansen et al., 2001).

###  2.1 X-rays

We process the EPIC-PN (PN) and EPIC-MOS (MOS) X-ray data using the SAS tools
epproc and emproc. We perform standard barycenter corrections using the SAS
barycen tool and filter out background flaring events. After this processing,
the total exposure times for PN and the MOS cameras are 100.9 ks and 128.5 ks
across 14 orbits. We extract the source from a circular region with radius
$\sim 30"$ and select circular background regions on the same detector chip as
the source with radii 4-5 times larger than the source aperture. We obtain
background subtracted X-ray light curves for the three detectors from the SAS
tool epiclccorr.

For the orbital parameters, we use the Fermi-LAT third gamma-ray pulsar
catalog ephemeris valid from 2008 August 4 until 2019 December 22 for all
observations (Smith et al., 2023), as this is the most recent published
ephemeris for J2215. The times of ascending node (TASC) deviate from the mean
by less than 0.003 in phase for the entire ephemeris range, so we expect no
significant shift in the TASC at our epoch (although 1.5-2 yr outside the
ephemeris range). The observed phase of the Optical Monitor (OM) U light
curve, to be discussed in the following subsection, is in good agreement with
the ephemeris prediction.

We convert the PN and MOS light curves from count-rate and count-rate error to
flux in each detector. We combine the data and distribute it into 26 phase
bins. We compute the average flux weighted by the inverse square of the errors
in each bin. Flux errors are estimated by the standard error of a variance-
weighted sample. Our combined 0.5-10 keV X-ray light curve is shown in the top
panel of Fig. 1.

Figure 1: Binned X-ray and optical orbital light curves of J2215. Less-finely
binned light curves are shown by faded markers. (Top) The combined PN and MOS
0.5-10 keV X-ray light curve. We identify four distinct regions for which we
separately fit spectra: two distinct peaks, P1 and P2 corresponding to the two
phase intervals in which the flow of the IBS shock is tangent to our line of
sight, a bridge (B) between the two peaks, and an off region around the
companion inferior conjunction where the X-ray emission comes principally from
the IBS nose. (Bottom) The OM U filter light curve for J2215 from the new XMM-
Newton observations.

###  2.2 Optical

Contemporaneous OM U filter data was taken in both image and fast mode. The
fast mode 10”$\times$10” images were obtained from a series of exposures over
each observation. We use the SAS tool omfchain to obtain U filter light
curves. For each exposure, we extract the source from a circular aperture with
radius of 3 detector pixels in the fast mode image and determine the
background from the average count rate of the corresponding image mode image
with all point sources removed. After processing, the OM exposure time is 91
ks. We convert the OM U data to flux using the published OM U calibration for
A-type stars. We place the data into flux bins by computing the variance
weighted sample mean and standard error. We also use the Fermi-LAT third
gamma-ray pulsar catalog ephemeris for the orbital parameters. The OM U light
curve is shown in the bottom panel of Fig. 1.

##  3 Spectral Analysis

We divide our X-ray light curve into four phase intervals as labeled in Fig. 1
and perform spectral analysis of these phase bins. We identify Peak 1 (P1;
$0.44<\Phi<0.58$) and Peak 2 (P2; $0.87<\Phi<1.07$) as the two characteristic
IBS X-ray peaks (Romani & Sanchez, 2016; Kandel et al., 2019). We also
demarcate a bridge region (B; $0.58<\Phi<0.87$) between P1 and P2 and an off
region ($0.07<\Phi<0.44$) corresponding to the companion at inferior
conjunction. Here non-IBS fluxes (e.g. from the magnetosphere or a pulsar wind
nebula) may dominate, as only weakly beamed X-rays from the IBS nose are
expected. We use the CIAO tool Sherpa (Freeman et al., 2001; Doe et al., 2007)
to fit non-thermal power-law spectral models to the XMM data with spectral
indices $\Gamma_{P1}$, $\Gamma_{P2}$, and $\Gamma_{B}$. We include an
additional phase-independent power-law component throughout the orbit with
power-law index $\Gamma_{0}$. The spectral model in counts s-1 keV-1 in each
phase region is

$f_{r}(E)=e^{-N_{H}\sigma(E)}\left(K_{r}E^{-\Gamma_{r}}+K_{0}E^{-\Gamma_{0}}\right),$
(1)

where $N_{H}$ is the equivalent hydrogen column, $\sigma(E)$ is the photo-
electric cross section, $K$ is the power-law normalization factor, and
$r=\\{P1,P2,B\\}$ demarcates the orbital phase region. For our absorption
model, we employ the photoelectric cross sections of Balucinska-Church &
McCammon (1992) and the elemental abundances of Anders & Grevesse (1989). In
the off phase, $K_{r}=0$. We include the underlying power-law to account for
phase-independent emission. The results of our fit are shown in Table 1.
Because the spectral indices of P1, P2, and B are not significantly different,
we also give the average spectral index across the IBS $\Gamma_{IBS}$. Fig. 2
shows the spectra across the IBS as well as in the off phases as observed by
PN.

Figure 2: The observed PN X-ray spectra with the equivalent hydrogen column $N_{H}$ fixed at $0.09\times 10^{22}$ cm-2. The plotted points are the spectra in the IBS phase region ($0.44<\Phi<1.07$; including P1, P2, and B) as well as in the off region ($0.07<\Phi<0.44$). The solid line shows the full model in the IBS phases $f_{tot}$ and the two dashed lines show the absorbed power-law models for the IBS component $f_{IBS}$ with power-law index $\Gamma_{IBS}$ and the phase-independent component $f_{0}$ with power-law index $\Gamma_{0}$. Note that in the off region, the only model component is the absorbed power-law $f_{0}$. | Fixed $N_{H}$ | Free $N_{H}$
---|---|---
$N_{H}$ ($10^{22}$ cm-2) | 0.09 | $0.10\pm 0.02$
$\Gamma_{P1}$ | $0.99\pm 0.11$ | $1.00\pm 0.11$
$\Gamma_{P2}$ | $0.85\pm 0.12$ | $0.86\pm 0.12$
$\Gamma_{B}$ | $0.85\pm 0.15$ | $0.86\pm 0.12$
$\Gamma_{IBS}$ | $0.94\pm 0.08$ | $0.96\pm 0.08$
$\Gamma_{0}$ | $1.43\pm 0.08$ | $1.48\pm 0.12$
$F_{tot}$ ($10^{-14}$ erg cm-2 $s^{-1}$) | 15.1 | 15.2
$F_{IBS}$ ($10^{-14}$ erg cm-2 $s^{-1}$) | 9.2 | 10.0
$F_{0}$ ($10^{-14}$ erg cm-2 $s^{-1}$) | 5.3 | 5.2
$\chi^{2}$/DoF | 1.03 | 1.03
Table 1: X-ray spectral model fits. We fit a phase-independent power-law with
index $\Gamma_{0}$ throughout the orbit and add additional power-laws with
separate spectral indices in the P1, P2, and B phase regions. The total
measured flux in the IBS phase region $0.44<\Phi<1.07$ is $F_{tot}$. The
measured flux of the specific IBS component is $F_{IBS}$, while the phase-
independent component is $F_{0}$. The equivalent hydrogen column $N_{H}$ is
fixed at $0.09\times 10^{22}$ cm-2, the value inferred from 3D dust maps
(Green et al., 2018), in the left column, while allowed to vary in the right.

##  4 Light Curve Modelling

###  4.1 Optical Models

Previous optical modeling of J2215 has been conducted by Linares et al. (2018)
and Kandel & Romani (2020) using g’r’i’ data from two observing nights at the
William Herschel Telescope (WHT). These light curve analyses show strong
heating asymmetries on the companion star, suggesting the presence of a hot
spot, which could plausibly be a magnetic pole (Sanchez & Romani, 2017).
Kandel & Romani (2020) include a hot spot at the coordinates
$(\theta_{HS},\phi_{HS})$ defined from the nose on the stellar surface. At
that location, the stellar temperature increases by a factor $(1+A_{HS})$ from
$T\approx 5600$ K to $T\approx 8000$ K with a Gaussian profile over radius
$r_{HS}$. Their best-fit model gives $i=68.9^{\circ}$ and consequently a very
high $M_{NS}=2.24$ M⊙. These results, however, are sensitive to the hot spot
parameters. Additionally, they were obtained using a simplified limb-darkening
model, which does not account for the limb-darkening of the reprocessed pulsar
heating luminosity (Claret & Bloemen, 2011; Kandel & Romani, 2023).

Figure 3: The best-fit companion optical light curve model plotted against the
g′r′i′ data from WHT and the OM U data shifted by the estimated calibration
offsets along with the residuals. Phase $0\leq\Phi\leq 1$ shows WHT night 1
while $1\leq\Phi\leq 2$ shows WHT night 2.

We refit the WHT g’r’i’ data along with the new OM U data using the ICARUS
light curve modeling code (Breton et al., 2013) including the extensions
described in Kandel & Romani (2020) and Kandel & Romani (2023). We fit the
light curves in each band simultaneously, using the same Fermi-LAT ephemeris
since it is valid for the WHT data taken 2014 August 11 and 2014 September 1.
We perform Monte Carlo sampling of the parameters using the library
PyMultinest (Feroz & Hobson, 2008; Feroz et al., 2009; Buchner et al., 2014;
Feroz et al., 2019), which uses nested sampling to search the parameter space.

The model is very similar to the best-fit hot spot model of Kandel & Romani
(2020); we augment a direct pulsar heating model with a surface hot spot, as
well as add a flat spectrum veiling flux $f_{2}$ to the WHT night 2 data. This
veiling flux was required by Kandel & Romani (2020), and is retained in our
fits. Our fits prefer calibration offsets $\delta m_{r}\approx 0.03$ to both
nights in the WHT r’ band and $\delta m_{U}\approx 0.66$ to the OM U data. We
have attempted to find the origin of this large U offset. As noted in Sec.
2.2, we obtain the OM U AB magnitude by converting from count rate to flux
with the SAS-recommended U conversion factor of $1.7\times 10^{-16}$ erg cm-2
Å-1 count-1 for a spectral type A stars. This spectral type is appropriate for
the companion near optical maximum. Using archival XMM OM U exposures in
fields covered by the Sloan Digital Sky Survey (SDSS), where catalog $u$
magnitudes are available, we have extracted stellar count rates following the
same procedure used for J2215. From these data we find a conversion $2.7\pm
0.2\times 10^{-16}$ erg cm-2 Å-1 count-1, which would suggest that the
companion U fluxes are brighter by a factor of $\sim 1.6$ (or offset by $\sim
0.5$ mag) with respect to the our XMM calibrated flux, close to the offset
needed.

We have also examined the calibrated 2014 Keck LRIS spectra of Romani et al.
(2015), whose gri colors are consistent with the 2014 WHT photometry. We find
the Keck fluxes at wavelengths below $4000$ Å consistent with BT-Settl model
atmosphere spectra (Allard, 2014) at the day side temperature of the star. A
change to the thermal companion flux (unlike the IBS non-thermal flux) seems
unlikely and would require a concomitant change in the heating power. Our
X-ray flux is consistent with past values and there is no significant change
to the (dominant) Fermi-LAT gamma-ray heating. Consequently, we consider a U
calibration offset the most likely explanation. While fitting for $\delta
m_{U}$ slightly decreases the predictive power of our model, the OM U light
curve shape nevertheless constrains the hot spot geometry and allows a good
model fit. Of course, precisely calibrated fluxes and simultaneous UV/optical
photometry would provide more secure constraints.

The fit-estimated parameters are shown in Table 2, while the best-fit model
light curves are shown in Fig. 3. The $\chi^{2}/$DoF=1.22, while not perfect,
is lower than 1.28 as reported by Kandel & Romani (2020). While the overall
convergence of the fit has improved, our constraints on $i$ are weaker. This
is a consequence of the different limb-darkening prescription used. Allowing
the heated luminosity to be reprocessed in addition to the internal stellar
flux tends to increase the flux on the night side of the star, thus requiring
larger $i$ values to reproduce the night-day flux difference. The constraints
on the hot spot parameters, however, have improved by a factor of 2, a result
of both the different limb-darkening prescription and the additional U band
data.

The larger value of $i$ implies a slightly lower neutron star mass than found
by (Kandel & Romani, 2020). To determine the neutron star mass (also shown in
Table 2) we use the Kandel & Romani (2020) radial velocity measurement
$K_{c}=429.8\pm 3.9$ km/s. Our new mass estimate $M_{NS}=2.15\pm 0.10$ M⊙ is
within a standard deviation of $M_{NS}=2.24\pm 0.09$ M⊙. Consequently, J2215
likely remains $>2$ M⊙ and is thus useful for constraining the neutron star
equation of state.

Parameters | Results | KR20 Results
---|---|---
$i$ (deg) | 71.8 ${}^{+2.5}_{-2.0}$ | 68.9${}^{+1.9}_{-1.8}$
$M_{NS}$ (M⊙) | 2.15${}^{+0.10}_{-0.10}$ | $2.24^{+0.09}_{-0.09}$
$f_{c}$ | 0.930${}^{+0.009}_{-0.009}$ | 0.94${}^{+0.01}_{-0.01}$
$L_{H}/10^{34}$ (erg s${}^{-1})$ | 2.50${}^{+0.09}_{-0.08}$ | 2.6${}^{+0.1}_{-0.1}$
$T_{N}$ (K) | 5592${}^{+18}_{-16}$ | 5682${}^{+14}_{-15}$
$d$ (kpc) | $3.32^{+0.04}_{-0.04}$ | $3.30^{+0.04}_{-0.04}$
$f_{2}$ ($\mu$Jy) | $1.7^{+0.2}_{-0.2}$ | 1.2${}^{+0.2}_{-0.2}$
$\delta m_{r}$ | $0.027^{+0.001}_{-0.001}$ | $0.031^{+0.001}_{-0.001}$
$\delta m_{U}^{\dagger}$ | $0.664^{+0.035}_{-0.036}$ | –
$\theta_{HS}$ (deg) | $340.5^{+6.0}_{-5.0}$ | $324.0^{+10.1}_{-8.7}$
$\phi_{HS}$ (deg) | $71.6^{+5.0}_{-10.0}$ | $73.2^{+5.0}_{-8.4}$
$A_{HS}$ | $0.35^{+0.16}_{-0.09}$ | $0.6^{+0.4}_{-0.2}$
$r_{HS}$ (deg) | $19.1^{+3.3}_{-6.4}$ | $15.5^{+3.9}_{-6.8}$
$\chi^{2}$/DoF | 311/255 | 297/232
Table 2: Parameter results for optical modeling of WHT+OM J2215 data including
the 68% confidence intervals. $f_{c}$ is the Roche lobe filling factor,
$L_{H}$ is the pulsar heating luminosity assuming an equatorially concentrated
heating flux, and $d$ is the source distance. The hot spot coordinates
$(\theta_{HS},\phi_{HS})$ are defined so that $(0,0)$ refers to the sub-pulsar
point at the companion nose. The dust extinction value $A_{V}=0.40$ (Green et
al., 2018) is used. The right column shows the fit results of Kandel & Romani
(2020) (KR20) which use a different gravity darkening model and do not include
U data. †$\delta m_{U}$ is an offset from the standard SAS flux calibration;
however, we see evidence for a systematic calibration shift of at least 0.5
mag.

###  4.2 X-ray Models

Figure 4: The best-fit IBS X-ray light curves. One model fit fixes $i$ at the best-fit value from the optical modeling, while the other model fit leaves $i$ as a free parameter. Parameters | Fixed $i$ | Free $i$
---|---|---
$i$ (deg) | 71 | 51.8 ${}^{+16.9}_{-8.5}$
$\beta$ | 2.12${}^{+0.252}_{-0.23}$ | 1.52${}^{+0.51}_{-0.32}$
$f_{v}$ | 5.7${}^{+5.4}_{-2.5}$ | 3.4${}^{+3.2}_{-1.3}$
$k$ | 0.092${}^{+0.031}_{-0.026}$ | 0.208${}^{+0.150}_{-0.105}$
$\delta\Phi$ | 0.027${}^{+0.009}_{-0.008}$ | 0.032${}^{+0.010}_{-0.010}$
$\dot{N}_{0}$ ($10^{31}$ s-1) | 3.8${}^{+0.3}_{-0.5}$ | 3.2${}^{+0.3}_{-0.4}$
$F_{b}$ ($10^{-14}$ erg cm-2 s-1) | 4.7${}^{+1.7}_{-1.4}$ | 5.0${}^{+1.4}_{-1.4}$
$\dot{M}_{w}$ ($10^{-10}$ $I_{45}$ M⊙ yr-1) | 1.8${}^{+1.4}_{-0.9}$ | 2.0${}^{+1.1}_{-0.9}$
$\chi^{2}$/DoF | 29/20 | 22/19
Table 3: Parameter results from IBS modeling of the X-ray light curve.
$\dot{M}_{w}$ is computed assuming $\dot{E}=5.2\times 10^{34}$ $I_{45}$ erg
s-1 (Smith et al., 2023).

We model the X-ray emission with the IBS synchrotron emission prescription of
Kandel et al. (2019) implemented in the ICARUS IBS code. The code discretizes
the IBS (assumed to be thin) into triangular tiles of constant angular size as
viewed from the pulsar, representing different zones from which synchrotron
radiation is emitted. An electron and positron (hereafter $e^{\pm}$)
population is injected into the IBS with energy spectrum in the flow frame

$\dot{N}(\gamma_{e})d\gamma_{e}=\dot{N}_{0}\gamma_{e}^{-p}d\gamma_{e},$ (2)

where $\gamma_{e}$ is the particle Lorentz factor in range
$\gamma_{min}<\gamma_{e}<\gamma_{max}$, $\dot{N}_{0}$ (in e s-1) is a global
normalization coefficient related to the particle injection rate and $p$
depends on the particle acceleration mechanism. After injection at a
particular tile, the $e^{\pm}$ population is advected to adjacent tiles with
bulk Lorentz factor $\Gamma_{bulk}$. As the particles are advected,
synchrotron losses cool the population. The residence time of the particle
population in a particular zone is computed in the bulk flow frame as in
Sullivan & Romani (2023). The synchrotron spectrum at each tile is calculated
from the time-averaged particle spectrum in that tile. See Kandel et al.
(2019) for further details on the emission prescription.

We assume an equatorially concentrated pulsar wind momentum flux
$\propto\sin^{2}\theta_{\star}$ in accordance with striped wind models (e.g.
Michel, 1973). $\theta_{\star}$ denotes the angle between the pulsar spin axis
(assumed to be aligned with the orbital angular momentum axis since spiders
pulsars are strongly recycled) and a position in the wind. The pulsar wind
collides with a spherical companion wind (Kandel et al., 2019). The stellar
wind to pulsar wind momentum ratio $\beta=\dot{M_{w}}v_{w}c/\dot{E}_{PSR}$
(Romani & Sanchez, 2016; Kandel et al., 2019) governs the overall shock
geometry and sets the X-ray peak separation. The relative heights of the peaks
are influenced by $i$ as well as the ratio of the stellar wind speed to the
orbital speed $f_{v}=v_{w}/v_{orb}$. Low values of $f_{v}$ cause the IBS shock
structure to be swept back by the orbital motion, with the geometry tracing
out an Archimedean spiral (see Romani & Sanchez (2016) for details on
implementation in ICARUS), rather than a symmetric conical structure. The
prescription for the bulk Lorentz factor in the post-shock flow is

$\Gamma_{Bulk}(s)=\Gamma_{nose}\left(1+k\frac{s}{r_{0}}\right),$ (3)

where $s$ is the arclength from the nose to a given point on the IBS, $r_{0}$
is the nose-standoff distance from the pulsar, $\Gamma_{nose}$ is the bulk
Lorentz factor at the nose, and $k$ is a scaling parameter that controls the
flow speed increase. The bulk flow direction is assumed to be parallel to the
contact discontinuity.

We fit for $\beta$, $f_{v}$, and $k$. We perform fits with fixed
$i=71^{\circ}$, the global minimum value from the optical fits, and with $i$
left as a free parameter. We fix $\Gamma_{nose}=1.1$, consistent with
numerical models (Bogovalov et al., 2008; Dubus et al., 2015; Kandel et al.,
2019). At the IBS nose, we set
$B_{IBS}=3B_{LC}\left(R_{LC}/r_{0}\right)\approx 70$ G (where $B_{LC}$ is the
magnetic field at the light cylinder radius $R_{LC}$; the factor of 3 is due
to relativistic shock compression; Sullivan & Romani, 2023), based on the
inferred orbital separation $a\approx 1.5\times 10^{11}$ cm, $P_{s}=2.61$ and
$\dot{P}_{s}=2.38\times 10^{-20}$ (Smith et al., 2023). We set $p=0.76$ so
that our IBS model spectrum has $\Gamma_{IBS}=0.94$ in the $0.5-10$ keV band,
consistent with the actual spectral results shown in Table 1. Our fits prefer
a forward phase shift $\delta\Phi\sim 0.03$ from the radio/gamma-ray
determined ephemeris and from our optical phase. We leave $\delta\Phi$ as a
fit parameter. We also fit the injected particle spectrum normalization
$\dot{N}_{0}$ and a phase independent flux $F_{b}$ added as an offset. We add
the offset to account for additional X-ray flux present in the source,
possibly due to a pulsar wind nebula (PWN) or the magnetosphere; this
component is also present in the spectral fits, as noted above.

Our best-fit light curves are shown in Fig. 4 and the parameter estimates are
shown in Table 3. Our model fits the data moderately well, with
$\chi^{2}/$DoF$=1.45$ when fixing $i$ and $\chi^{2}/\text{DoF}=1.16$ when
leaving $i$ free. This corresponds to an Akaike information criterion (AIC)
value of 41 for the fixed $i$ model and 36 for the free $i$ model. Departures
from a simple light curve in the phase range $\Phi=0.6-0.7$ contribute a large
factor to the $\chi^{2}$. In examining the orbit-to-orbit variations, we find
that this phase range has somewhat higher variability than the rest of the
light curve, but we are unable to attribute the structure in this phase range
to IBS flares (as have been seen in other spider pulsars; e.g. An et al.,
2017). Thus, the flux dip at $\Phi\approx 0.6$ might represent an absorption
component, or, alternatively, the flux increase at $\Phi\approx 0.65$ might
represent an extra emission component, possibly structure in a thick post-
shock flow not captured by our thin-shock geometry. We have left these points
in the fit to be conservative; excising these points barely changes the fixed
$i$ best-fit parameters. We also note the low bin at $\Phi=0.25$. While at the
appropriate phase for a companion eclipse, its significance is low.
Furthermore, we would not expect significant occultation of the extended IBS
at our moderate $\sim 70^{\circ}$ inclination.

The model fitting gives values $\beta\approx 2$ and $f_{v}\approx 3-6$,
signaling a rather flat IBS and slow companion wind. If left free, $i$
decreases to $52^{+17}_{-9}$ deg to accommodate the bridge flux (including the
$\Phi\approx 0.65$ component). Such low $i$ implies $M_{NS}=3.8^{+3.9}_{-1.4}$
M⊙; the central value is unphysical, although less than $2\sigma$ different
from our optical $i$ mass determination. Indeed, $i$ itself is only in $\sim
1\sigma$ disagreement with the optical value. We thus consider the fixed $i$
IBS fit more reliable.

##  5 Discussion and Conclusion

Our optical and X-ray model fits, while imperfect, provide useful estimates of
the parameters of J2215. Having used the radial velocity $K_{c}$ from Kandel &
Romani (2020), we believe that these new values should supersede the results
of that paper and those in the less realistic modeling of Linares et al.
(2018). Our revised value for the mass is just about $1\sigma$ less than the
Kandel & Romani (2020) estimate, but, at $M_{NS}=2.15\pm 0.10$ M⊙, remains
more than $1\sigma$ larger than 2 M⊙.

With well-detected light curve peaks, our X-ray results are much more robust
than those of Romani & Sanchez (2016). Our modeling suggests $\beta\lesssim 2$
and $f_{v}<10$. From these results and $\dot{E}=5.2\times 10^{34}$ $I_{45}$
erg s-1 (Smith et al., 2023), we estimate the companion mass-loss rate

$\dot{M}_{w}=\frac{\beta\dot{E}}{f_{v}v_{orb}c},$ (4)

showing the results in Table 3. $\dot{M}_{w}\approx 2\times 10^{-10}$ $I_{45}$
M⊙ yr-1 suggests that the companion should evaporate in $\sim 1.5/I_{45}$ Gyr,
a modest fraction of the age of the universe, and comparable to the
$\tau_{c}=1.2$ Gyr characteristic spindown age. This is short enough for J2215
to evolve into an isolated MSP (Kluzniak et al., 1988). However, with $\beta$
unusually close to unity, the IBS location and thus mass-loss rate are
sensitive to the fit details; these include the small unexplained phase shift.

The estimated value of the phase-independent flux from the light curve fits
$F_{b}\approx 5\times 10^{-14}$ erg cm-2 s-1 is in good agreement with the
phase-independent flux estimated from the spectral fitting $F_{0}$. These
fluxes correspond to a luminosity of $7\times 10^{31}$ erg s-1 at $d=3.3$ kpc
and could originate from magnetospheric emission or, more likely, from a
faint, compact PWN ( e.g. Li et al., 2008). The spectral index
$\Gamma_{0}\approx 1.5$ is quite consistent with a PWN source, while the
luminosity of $10^{-3}{\dot{E}}$ is similar to other observed PWNe
efficiencies. The slightly larger value of $F_{0}$ would include a small
amount of emission from the wings of the IBS peaks in the off-phase.

It is also useful to compare our measured energy fluxes with the pulsar spin-
down power budget ${\dot{E}}=5\times 10^{34}$ $I_{45}$ ${\rm erg\,s^{-1}}$.
Several neutron star masses have been securely measured above $2M_{\odot}$,
pointing towards a quite stiff equation of state. We can estimate
$I_{45}\approx[0.8-1.2](M/M_{\odot})^{3/2}$ (Lattimer & Schutz, 2005), i.e.
$2.5-3.8$ for $2.15$ $M_{\odot}$. The ${\rm sin}^{2}\theta_{\ast}$-distributed
heating luminosity $L_{H}\approx 2.5\times 10^{34}$ ${\rm erg\,s^{-1}}$ is in
good accord with the isotropic GeV gamma-ray luminosity measured by the Fermi
LAT, $L_{\gamma}=3.3\times 10^{34}$ ${\rm erg\,s^{-1}}$, suggesting that our
‘direct heating’ is dominated by the GeV photons. This requires that the
radiation be directed toward the companion at the pulsar spin equator, as
expected from the spin-orbit aligned recycled pulsar outer magnetosphere/wind
emission model (Philippov & Spitkovsky, 2018). This gamma-ray direct heating
power is a substantial fraction of the spin-down energy budget, but can be
comfortably accommodated for $I_{45}=2-3$. This also reduces our estimated
evaporation timescale below 1 Gyr.

Our IBS fit also estimates the $e^{\pm}/B$ pulsar wind power $\dot{E}_{pw}$,
which is processed by the shock into a power-law electron distribution,

$\dot{E}_{pw}=\frac{\dot{N}_{0}}{2-p}m_{e}c^{2}(\gamma^{2-p}_{\rm
max}-\gamma_{\rm min}^{2-p}).$ (5)

The observed 0.5-10 keV X-ray spectrum comes from $10^{4}<\gamma<10^{5}$
$e^{\pm}$ particles. Our fit $\dot{N}_{0}\sim 3-4\times 10^{31}$ s-1 thus
gives a minimum $\dot{E}_{pw}\gtrsim 5\times 10^{31}$ erg s-1 for $p\approx
0.8$. This is well below the available power ${\dot{E}}-L_{H}\approx 12\times
10^{34}$ ${\rm erg\,s^{-1}}$ for $I_{45}\approx 3$. To avoid saturating the
residual spin-down power, there is an upper limit of $\gamma_{\rm max}\sim
5.5\times 10^{7}$. Higher energy X-ray and gamma-ray measurements can further
probe the high $\gamma$ $e^{\pm}$ population (which cools rapidly). If
magnetic reconnection is the primary acceleration mechanism (as suggested by
the very hard observed spectra), the maximum electron energy should be
$\gamma_{\rm max}\sim\sigma\gamma_{\rm min}$, where
$\sigma=B^{2}/(4\pi\gamma_{\rm min}n_{0}m_{e}c^{2})$ is the pulsar wind
magnetization and $n_{0}$ is the e± number density (Sironi & Spitkovsky, 2011,
2014). As expected, the J2215 pulsar wind should be very strongly magnetized
with $10\lesssim\sigma\lesssim 6\times 10^{7}$. Finally, the companion hot
spot indicates extra localized heating. Sanchez & Romani (2017) attribute
companion hot spots to IBS-energized particles precipitating to the companion
magnetic poles. In our fits, the small $\beta$ allows the relatively flat IBS
shock to capture $\sim 40-50\%$ of the $\dot{E}_{pw}$ flux, while the thermal
emission of the hot spot represents $1\times 10^{32}$ ${\rm erg\,s^{-1}}$.
Only a small fraction $\sim 2\times 10^{-3}$ of the maximum available IBS-
processed pulsar wind power needs to reach the surface to heat the companion
magnetic pole. Simultaneous optical/UV observations would be useful to re-
check the U band normalization and further constrain the hot spot
contribution.

J2215 remains one of the most massive neutron stars known. Its high inferred
mass-loss rate makes this system a plausible isolated MSP progenitor, while
its unusually flat IBS may reprocess a large fraction of the spin-down power.
The hard X-ray spectrum and sharp peaks provide an important probe of
relativistic shock dynamics and particle acceleration. Higher energy
observations (e.g. with NuSTAR) should further constrain the IBS $e^{\pm}$
population. Since the prominent IBS peaks indicate beamed emission and thus
moderate bulk Lorentz factors in the IBS flow, the various orbital phases
probe different portions of the shock with varying obliquity. Thus, IBS-
dominated spider pulsars, like J2215, show emission from a range of
relativistic shock geometries. Future phase-resolved spectrum and polarization
measurements of these natural laboratories can be compared with numerical
simulations to advance our understanding of oblique, strongly magnetized
relativistic shocks.

## Acknowledgements

The authors thank Dinesh Kandel for help using ICARUS and for clarifications
on the previous optical modeling as well as the anonymous referee for helpful
feedback. This work was supported in part by NASA grant 80NSSC22K1506. A.S.
acknowledges the support of the Stanford University Physics Department
Fellowship and the National Science Foundation Graduate Research Fellowship.

## References

  * Allard (2014) Allard, F. 2014, in Exploring the Formation and Evolution of Planetary Systems, ed. M. Booth, B. C. Matthews, & J. R. Graham, Vol. 299, 271–272, doi: 10.1017/S1743921313008545
  * An et al. (2017) An, H., Romani, R. W., Johnson, T., Kerr, M., & Clark, C. J. 2017, ApJ, 850, 100, doi: 10.3847/1538-4357/aa947f
  * Anders & Grevesse (1989) Anders, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197, doi: 10.1016/0016-7037(89)90286-X
  * Atwood et al. (2009) Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071, doi: 10.1088/0004-637X/697/2/1071
  * Balucinska-Church & McCammon (1992) Balucinska-Church, M., & McCammon, D. 1992, ApJ, 400, 699, doi: 10.1086/172032
  * Bogovalov et al. (2008) Bogovalov, S. V., Khangulyan, D. V., Koldoba, A. V., Ustyugova, G. V., & Aharonian, F. A. 2008, MNRAS, 387, 63, doi: 10.1111/j.1365-2966.2008.13226.x
  * Brandes et al. (2023) Brandes, L., Weise, W., & Kaiser, N. 2023, arXiv e-prints, arXiv:2306.06218, doi: 10.48550/arXiv.2306.06218
  * Breton et al. (2013) Breton, R. P., van Kerkwijk, M. H., Roberts, M. S. E., et al. 2013, ApJ, 769, 108, doi: 10.1088/0004-637X/769/2/108
  * Broderick et al. (2016) Broderick, J. W., Fender, R. P., Breton, R. P., et al. 2016, MNRAS, 459, 2681, doi: 10.1093/mnras/stw794
  * Buchner et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125, doi: 10.1051/0004-6361/201322971
  * Claret & Bloemen (2011) Claret, A., & Bloemen, S. 2011, A&A, 529, A75, doi: 10.1051/0004-6361/201116451
  * Doe et al. (2007) Doe, S., Nguyen, D., Stawarz, C., et al. 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, 543
  * Dubus et al. (2015) Dubus, G., Lamberts, A., & Fromang, S. 2015, A&A, 581, A27, doi: 10.1051/0004-6361/201425394
  * Feroz & Hobson (2008) Feroz, F., & Hobson, M. P. 2008, MNRAS, 384, 449, doi: 10.1111/j.1365-2966.2007.12353.x
  * Feroz et al. (2009) Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601, doi: 10.1111/j.1365-2966.2009.14548.x
  * Feroz et al. (2019) Feroz, F., Hobson, M. P., Cameron, E., & Pettitt, A. N. 2019, The Open Journal of Astrophysics, 2, 10, doi: 10.21105/astro.1306.2144
  * Freeman et al. (2001) Freeman, P., Doe, S., & Siemiginowska, A. 2001, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4477, Astronomical Data Analysis, ed. J.-L. Starck & F. D. Murtagh, 76–87, doi: 10.1117/12.447161
  * Gentile et al. (2014) Gentile, P. A., Roberts, M. S. E., McLaughlin, M. A., et al. 2014, ApJ, 783, 69, doi: 10.1088/0004-637X/783/2/69
  * Green et al. (2018) Green, G. M., Schlafly, E. F., Finkbeiner, D., et al. 2018, MNRAS, 478, 651, doi: 10.1093/mnras/sty1008
  * Hessels et al. (2011) Hessels, J. W. T., Roberts, M. S. E., McLaughlin, M. A., et al. 2011, in American Institute of Physics Conference Series, Vol. 1357, Radio Pulsars: An Astrophysical Key to Unlock the Secrets of the Universe, ed. M. Burgay, N. D’Amico, P. Esposito, A. Pellizzoni, & A. Possenti, 40–43, doi: 10.1063/1.3615072
  * Hui & Li (2019) Hui, C. Y., & Li, K. L. 2019, Galaxies, 7, 93, doi: 10.3390/galaxies7040093
  * Jansen et al. (2001) Jansen, F., Lumb, D., Altieri, B., et al. 2001, A&A, 365, L1, doi: 10.1051/0004-6361:20000036
  * Kandel & Romani (2020) Kandel, D., & Romani, R. W. 2020, ApJ, 892, 101, doi: 10.3847/1538-4357/ab7b62
  * Kandel & Romani (2023) —. 2023, ApJ, 942, 6, doi: 10.3847/1538-4357/aca524
  * Kandel et al. (2019) Kandel, D., Romani, R. W., & An, H. 2019, ApJ, 879, 73, doi: 10.3847/1538-4357/ab24d9
  * Kandel et al. (2021) —. 2021, ApJ, 917, L13, doi: 10.3847/2041-8213/ac15f7
  * Kluzniak et al. (1988) Kluzniak, W., Ruderman, M., Shaham, J., & Tavani, M. 1988, Nature, 334, 225, doi: 10.1038/334225a0
  * Lattimer & Prakash (2007) Lattimer, J. M., & Prakash, M. 2007, Phys. Rep., 442, 109, doi: 10.1016/j.physrep.2007.02.003
  * Lattimer & Schutz (2005) Lattimer, J. M., & Schutz, B. F. 2005, ApJ, 629, 979, doi: 10.1086/431543
  * Li et al. (2008) Li, X.-H., Lu, F.-J., & Li, Z. 2008, ApJ, 682, 1166, doi: 10.1086/589495
  * Linares (2014) Linares, M. 2014, ApJ, 795, 72, doi: 10.1088/0004-637X/795/1/72
  * Linares et al. (2018) Linares, M., Shahbaz, T., & Casares, J. 2018, ApJ, 859, 54, doi: 10.3847/1538-4357/aabde6
  * Michel (1973) Michel, F. C. 1973, ApJ, 180, L133, doi: 10.1086/181169
  * Philippov & Spitkovsky (2018) Philippov, A. A., & Spitkovsky, A. 2018, ApJ, 855, 94, doi: 10.3847/1538-4357/aaabbc
  * Ray et al. (2012) Ray, P. S., Abdo, A. A., Parent, D., et al. 2012, arXiv e-prints, arXiv:1205.3089, doi: 10.48550/arXiv.1205.3089
  * Roberts (2013) Roberts, M. S. E. 2013, in Neutron Stars and Pulsars: Challenges and Opportunities after 80 years, ed. J. van Leeuwen, Vol. 291, 127–132, doi: 10.1017/S174392131202337X
  * Romani et al. (2015) Romani, R. W., Graham, M. L., Filippenko, A. V., & Kerr, M. 2015, ApJ, 809, L10, doi: 10.1088/2041-8205/809/1/L10
  * Romani & Sanchez (2016) Romani, R. W., & Sanchez, N. 2016, ApJ, 828, 7, doi: 10.3847/0004-637X/828/1/7
  * Sanchez & Romani (2017) Sanchez, N., & Romani, R. W. 2017, ApJ, 845, 42, doi: 10.3847/1538-4357/aa7a02
  * Schroeder & Halpern (2014) Schroeder, J., & Halpern, J. 2014, ApJ, 793, 78, doi: 10.1088/0004-637X/793/2/78
  * Sironi & Spitkovsky (2011) Sironi, L., & Spitkovsky, A. 2011, ApJ, 741, 39, doi: 10.1088/0004-637X/741/1/39
  * Sironi & Spitkovsky (2014) —. 2014, ApJ, 783, L21, doi: 10.1088/2041-8205/783/1/L21
  * Smith et al. (2023) Smith, D. A., Abdollahi, S., Ajello, M., et al. 2023, ApJ, 958, 191, doi: 10.3847/1538-4357/acee67
  * Steiner et al. (2013) Steiner, A. W., Lattimer, J. M., & Brown, E. F. 2013, ApJ, 765, L5, doi: 10.1088/2041-8205/765/1/L5
  * Sullivan & Romani (2023) Sullivan, A. G., & Romani, R. W. 2023, ApJ, 959, 81, doi: 10.3847/1538-4357/ad09ae
  * van Paradijs et al. (1988) van Paradijs, J., Allington-Smith, J., Callanan, P., et al. 1988, Nature, 334, 684, doi: 10.1038/334684a0
  * Wadiasingh et al. (2017) Wadiasingh, Z., Harding, A. K., Venter, C., Böttcher, M., & Baring, M. G. 2017, ApJ, 839, 80, doi: 10.3847/1538-4357/aa69bf
  * Yao et al. (2017) Yao, J. M., Manchester, R. N., & Wang, N. 2017, ApJ, 835, 29, doi: 10.3847/1538-4357/835/1/29


# On the Control of Flying Qubits††thanks: The material in this paper was
partially presented at the 20th IFAC World Congress, July 13-17, 2020, Berlin,
Germany. RBW acknowledges supports by National Key Research and Development
Program of China (Grant No. 2017YFA0304300) and NSFC grants (Nos. 61833010 and
61773232). GFZ acknowledges supports in part by Hong Kong Research Grant
council (RGC) grants (No. 15208418, No. 15203619) and Shenzhen Fundamental
Research Fund, China under Grant No. JCYJ20190813165207290.

Wen-Long Li Department of Automation, Tsinghua University, Beijing, 100084,
China Guofeng Zhang Department of Applied Mathematics, The Hong Kong
Polytechnic University, Hong Kong Re-Bing Wu Department of Automation,
Tsinghua University, Beijing, 100084, China, and Beijing National Research
Center for Information Science and Technology, Beijing, 100084, China

###### Abstract

The control of flying quantum bits (qubits) carried by traveling quantum
fields is crucial for coherent information transmission in quantum networks.
In this paper, we develop a general framework for modeling the generation,
catching and transformation processes of flying qubits. We introduce the
quantum stochastic differential equation (QSDE) to describe the flying-qubit
input-output relations actuated by a standing quantum system. Under the
continuous time-ordered photon-number basis, the infinite-dimensional QSDE is
reduced to a low-dimensional deterministic non-unitary differential equation
for the state evolution of the standing system, and the outgoing flying-qubit
states can be calculated via randomly occurring quantum jumps. This makes it
possible, as demonstrated by examples of flying-qubit generation and
transformation, to analyze general cases when the number of excitations is not
reserved. The proposed framework lays the foundation for the design of flying-
qubit control systems from a control theoretic point of view, within which
advanced control techniques can be incorporated for practical applications.
Keywords: quantum control; flying qubits; quantum stochastic differential
equation

##  1 Introduction

Quantum computation has been evidently demonstrated to be surpassing classical
supercomputers in the near future [37]. Towards large-scale and distributed
quantum information processing (QIP) [9], the efficient control of networked
quantum systems is highly demanding [5]. In the past decades, the control of
“standing” components (namely the nodes, e.g., atoms or resonators [21]) for
on-site QIP had been extensively studied [4, 6, 45, 25, 15, 42]. However, the
control of “flying” components for information transmission, as are often
called flying qubits, has received much less attention [31].

In quantum networks, the flying qubits are encoded by the number of photons or
the poarization of single photons in travelling electromagnetic fields [14, 3,
9] ranging from optical regime (e.g., from quantum dots or NV-centers [28, 27,
10]) down to the microwave regime (e.g., from superconducting artificial atoms
[22]). Roughly speaking, the flying-qubit control problems can be categorized
into the following three classes (see Fig. 1 for illustration).

Figure 1: Schematics of flying-qubit control problems actuated by a standing
quantum system: (a) the generation of flying qubits; (b) the catching of
flying qubits; (c) the transformation of flying qubits.

The first class is the generation of flying qubits for the purpose of coding
and sending quantum information, the generated flying qubits are required to
be on-demand (in the sense that single photons can be triggered at any
appointed time) and be flexibly tunable in the frequency [35] or the photon
shapes [26, 34, 12, 2, 1] so as to match remote network nodes [34, 2, 12]).
The second class is the catching of incoming flying qubits for the purpose of
receiving (and storing) quantum information, which is usually done by exciting
the standing system from its ground state [43, 33]. The third class is the
transformation of flying qubits (e.g., manipulating the photon number of the
shape of incoming flying qubits) for general purposes of information
processing [41, 29].

The control of flying qubits need to be actuated by a standing quantum system
(e.g., an atom or a resonator [46]) of which the incoming and outgoing flying
qubits form quantum inputs or outputs. In the generation of flying qubits, the
standing system receives vaccum quantum input and flying qubits are yielded as
its quantum output. As an inverse process of the generation, the standing
system for catching flying qubits receives flying qubits as its quantum input
and yields vacuum quantum output. The transformation process, which is also
called scattering processes in quantum physics [44], can be taken as a
generalization of the generation and catching processes, in which both the
quantum input and output are non-empty.

To control the flying-qubit input-output processes, one also need classical
controls, i.e., tunable time-dependent parameters with some classical devices,
e.g., electromagnetic fields that coherently drive the standing system [34,
11] and tunable couplers that incoherently vary the interaction of the
standing system to the waveguide [36]. The control can also be indirectly done
by nonlocal spectral filtering using ancilla qubits [1] or coherent quantum
feedback loops [7, 48]. To model the underlying control dynamics, the Quantum
Stochastic Differential Equation (QSDE) [23, 13] is a powerful tool for
analyzing and filtering the input-output responses [49, 47, 16, 40, 8, 48],
based on which flying-qubit generaters can be designed with linear or
nonlinear finite-level systems [17, 16, 19], as well as the catching of flying
qubits via two-level systems with tunable waveguide coupling [32]. However,
most existing results [46, 39] are heavily dependent on the conservation of
energy (or the excitation number) that substantially reduces the effective
dimensionality of system. This symmetry will be broken in the presence of
time-dependent coherent driving controls as they usually lead to infinite
excitations.

In a preliminary version of this paper [30], we found that, under an infinite
series expansion, the QSDE approach can be extended to flying-qubit generation
systems whose number of excitations is not preserved. Here, we will show that
this proposed approach can be further applied to catching and transformation
problems with non-vacuum quantum inputs, forming a unified framework for the
analysis and design of flying-qubit control systems. The overall contributions
we made are as follows. First, we introduce the time-ordered photon-number
basis to represent flying-qubit states and their joint state with the standing
system. The resulting expansion can be taken as a generalization of the
representation proposed in Refs. [19, 32], but it is applicable to arbitrary
states whose number of excitations are not conservative. Second, as summarized
in Theorem 3.1 and Corollary 3.1 in Sec. 3, we reduce the original infinite-
dimensional QSDE to a low-dimensional deterministic non-unitary evolution
equation (see Eq. (11) in Sec. 3) associated with the standing system. This
equation can be efficiently solved to calculate the state of outgoing flying
qubits via random quantum jumps (see Eq. (19) in Sec. 3), which clearly
interpret the physical generation process of flying qubits. The resulting
model is fully consistent with the results obtained by scattering theory in
Ref. [11] for flying-qubit generation, but the derivation is much more concise
and easier to generalize to complex quantum networks via $(S,L,H)$ formulation
[18]. This is quite beneficial to quantum control engineering owing to the
prevalent use of the state-space models in control system design. Third, we
show that the flying-qubit catching and transformation problems can be unified
as generation problems via an auxiliary quantum system that virtually produces
the incoming flying qubits, thus forming a systematic QSDE-based framework for
treating flying-qubit control problems.

The remainder of this paper is organized as follows. In Sec. 2, necessary
preliminaries will be provided for the state representation of flying qubits.
In Sec. 3, the model based on a non-unitary evolution equation is derived for
calculating the steady-state of the outgoing flying qubits with vacuum or non-
vacuum quantum inputs. In Sec. 4, we demonstrate the model with applications
to flying-qubit generation and transformation problems, respectively, with a
two-level standing quantum system. Finally, concluding remarks are made in
Sec. 5.

##  2 The Quantum State Representation of Flying Qubits

Consider flying qubits carried by single photons in quantized electromagnetic
fields traveling in waveguides. The field generally contains continuous
bosonic modes whose annihilation and creation operators satisfy the singular
commutation relations
$[b(\omega),b^{{\dagger}}(\omega^{\prime})]=\delta(\omega-\omega^{\prime})$,
where $\omega\in\mathbb{R}$ is the mode frequency. In practice, we often
represent the field with the inverse Fourier transform [13]

$b_{\tau}=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}d\omega\
e^{-i\omega\tau}b(\omega)$ (1)

$b(\omega)$, which describes the annihilation of a photon at time $\tau$. The
temporal quantum field. $b_{\tau}$ satisfies the singular commutation relation
$[b_{\tau},b^{{\dagger}}_{\tau^{\prime}}]=\delta(\tau-\tau^{\prime})$, and can
be further used to construct the Ito-type quantum Wiener increment

$dB_{\tau}=\int_{\tau}^{\tau+d\tau}b_{s}ds,\quad d\tau>0,$ (2)

that obeys Ito’s rule

$\begin{split}&dB_{\tau}dB^{{\dagger}}_{s}=\delta_{\tau,s}d\tau,\\\
&dB^{{\dagger}}_{\tau}dB_{s}=dB_{\tau}dB_{s}=dB^{{\dagger}}_{\tau}dB^{\dagger}_{s}=0.\end{split}$
(3)

Using the Ito-type increment $dB^{{\dagger}}_{\tau}$, we can define the state
of an arbitrary single flying-qubit as follows:

$\ket{1_{\xi}}_{F}=\int_{-\infty}^{\infty}\xi^{\tau}dB^{\dagger}_{\tau}\ket{\Omega}_{F},$
(4)

where $\ket{\Omega}_{F}$ is the vacuum state of the quantum field and the
time-dependent shape function $\xi^{\tau}d\tau$ indicates the probability
amplitude of observing the photon between $\tau$ and $\tau+d\tau$. Similarly,
we can define the $\ell$-photon state [11],

$\ket{\ell_{\xi}}_{F}=\int_{-\infty}^{\infty}\int_{-\infty}^{\tau_{\ell}}\cdots\int_{-\infty}^{\tau_{2}}\xi^{\tau_{1},\ldots,\tau_{\ell}}dB^{\dagger}_{\tau_{1}}\cdots
dB^{\dagger}_{\tau_{\ell}}\ket{\Omega}_{F},$ (5)

and more generally, any flying-qubit state can be expanded as the
superposition of multi-photon states:

$\ket{\xi}_{F}=\sum_{\ell=0}^{\infty}\alpha_{\ell}\ket{\ell_{\xi}}=\sum_{\ell=0}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\tau_{\ell}}\cdots\int_{-\infty}^{\tau_{2}}\xi^{\tau_{1},\ldots,\tau_{\ell}}dB^{\dagger}_{\tau_{1}}\cdots
dB^{\dagger}_{\tau_{\ell}}\ket{\Omega}_{F},$ (6)

in which $\xi$ is a scalar when $\ell=0$ and the coefficients $\alpha_{\ell}$
has been absorbed by $\xi^{\tau_{1},\cdots,\tau_{\ell}}$. Here, for
simplicity, we abuse the symbolic use of $\xi$ to represent all components in
the superposition of multi-photon wavepackets
$\\{\xi,\xi^{\tau_{1}},\xi^{\tau_{1},\tau_{2}},...\\}$. Under this
representation, the probability of observing $\ell$ flying qubits in
$\ket{\xi}$ can be calculated by the integral

$P_{\ell}=\int_{-\infty}^{\infty}\int_{-\infty}^{\tau_{\ell}}\cdots\int_{-\infty}^{\tau_{2}}|\xi^{\tau_{1},\cdots,\tau_{\ell}}|^{2}d\tau_{\ell}\cdots
d\tau_{1}.$ (7)

The summation $\sum_{\ell=0}^{\infty}P_{\ell}$ of these probabilities must be
equal to 1, and this can be easily verified by Ito’s rule from the
normalization condition ${}_{F}\langle\xi|\xi\rangle_{F}=1$.

##  3 Flying-Qubit Dynamics via Classical Controls

Throughout this paper, we always assume that the classical controls on the
standing system are imposed from the initial time $t_{0}=0$, and the standing
system is coupled to a chiral waveguide in which all flying qubits propagate
rightwards [35]. To analyze the time evolution of the joint system, we expand
the time-dependent joint state of the system and the flying qubits as follows:

$\ket{\Psi(t)}=\sum_{\ell=0}^{\infty}\int_{0}^{t}\int_{0}^{\tau_{\ell}}\cdots\int_{0}^{\tau_{2}}|\psi^{\tau_{1},\cdots,\tau_{\ell}}(t)\rangle\otimes
dB^{\dagger}_{\tau_{1}}\cdots dB^{\dagger}_{\tau_{\ell}}\ket{\Omega}_{F},$ (8)

where $|\psi^{\tau_{1},\cdots,\tau_{\ell}}(t)\rangle$ is the standing system’s
correlated state with the field when $\ell$ photons are generated at
$\tau_{1},\cdots,\tau_{\ell}$. The modified lower and upper bounds of the
first integral indicate that, physically, flying qubits can only be generated
between time $0$ and the present time $t$. The expansion contains terms that
possess different numbers of excitations, which can be taken as the
generalization of the state representation used in Ref. [32] that involves one
and only one excitation.

In the following, we will introduce the Quantum Stochastic Differential
Equations (QSDEs) to model the flying-qubit control processes with vacuum or
non-vacuum quantum inputs, respectively.

###  3.1 Standing system with vacuum quantum inputs

We start from the simple case that the quantum input of the standing system is
empty (e.g., in its vacuum state), i.e., the joint system is initially at
$\ket{\Psi(0)}=\ket{\psi_{0}}\otimes\ket{\Omega}_{F}$ with $|\psi_{0}\rangle$
as the initial state of the standing system. Under this circumstance, the
evolution of the joint state is governed by the following QSDE of the joint
state $|\Psi(t)\rangle$ [20, 32]:

$d|\Psi(t)\rangle=\left[-iH_{\rm
eff}(t)dt+L(t)dB^{{\dagger}}_{t}\right]|\Psi(t)\rangle,$ (9)

where $L(t)$ represents the coupling operator of the standing system to the
flying qubits, and incoherent classical controls can be realized by varying
$L(t)$ in time. The effective non-Hermitian Hamiltonian

$H_{\rm eff}(t)=H(t)-\frac{i}{2}L^{\dagger}(t)L(t),$ (10)

where $H(t)$ includes the standing system’s internal Hamiltonian and
interaction Hamiltonian associated with the driving fields that serve as
coherent classical control functions.

Now we present the main Theorem for the analysis of flying-qubit control
dynamics.

######  Theorem 3.1.

Let $V(t)$ be the propagator of the differential equation

$\dot{V}(t)=-iH_{\rm eff}(t)V(t),$ (11)

with $V(0)=\mathbb{I}$. Then, the correlated states

$|\psi^{\tau_{1},\cdots,\tau_{\ell}}(t)\rangle=V(t)\tilde{L}(\tau_{\ell})\cdots\tilde{L}(\tau_{1})|\psi_{0}\rangle,$
(12)

for $\ell=0,1,\ldots$, where
$\tilde{L}(\tau_{k})=V^{-1}(\tau_{k})L(\tau_{k})V(\tau_{k})$.

_Proof._ Differentiate both sides of Eq. (8), we have

$\begin{split}d|\Psi(t)\rangle=&~{}\sum_{\ell=0}^{\infty}\int_{0}^{t}\int_{0}^{\tau_{\ell}}\cdots\int_{0}^{\tau_{2}}|\dot{\psi}^{\tau_{1},\cdots,\tau_{\ell}}(t)\rangle
dt\otimes dB^{{\dagger}}_{\tau_{1}}\cdots
dB^{{\dagger}}_{\tau_{\ell}}\ket{\Omega}_{F}\\\
&+dB^{{\dagger}}_{t}\sum_{\ell=1}^{\infty}\int_{0}^{t}\int_{0}^{\tau_{\ell-1}}\cdots\int_{0}^{\tau_{2}}|\psi^{\tau_{1},\cdots,\tau_{\ell-1},t}(t)\rangle\otimes
dB^{{\dagger}}_{\tau_{1}}\cdots
dB^{{\dagger}}_{\tau_{\ell-1}}\ket{\Omega}_{F},\end{split}$ (13)

where the dot “$\cdot$” represents the partial derivative with respect to time
variable $t$. Replacing Eqs. (8) and (13) into Eq. (9), we can obtain

$|\dot{\psi}^{\tau_{1},\cdots,\tau_{\ell}}(t)\rangle=-iH_{\rm
eff}(t)|\psi^{\tau_{1},\cdots,\tau_{\ell}}(t)\rangle$ (14)

by comparing the terms with $dt$, where $\ell\geq 0$. The comparison of terms
with $dB^{{\dagger}}_{t}$ gives the boundary condition

$|\psi^{\tau_{1},\cdots,\tau_{\ell-1},t}(t)\rangle=L(t)|\psi^{\tau_{1},\cdots,\tau_{\ell-1}}(t)\rangle.$
(15)

This condition indicates the emission of a photon (or a flying qubit)
accompanied with the lowering of the standing system’s energy levels by
operator $L(t)$, which is also called quantum jumps in open quantum system
theory [45]. The transition of wavepacket functions from
$|\psi^{\tau_{1},\cdots,\tau_{\ell-1}}(t)\rangle$ to
$|\psi^{\tau_{1},\cdots,\tau_{\ell-1},t}(t)\rangle$ implies that the number of
photons in the field is increased by 1, where the additional photon is emitted
at time $t$ from the standing system.

Observing that all vector functions
$|\psi^{\tau_{1},\cdots,\tau_{\ell}}(t)\rangle$ share the same differential
equation (14) with respect to time $t$, we denote by $V(t)$ their common non-
unnitary propagator steered by the effective Hamiltonian $H_{\rm eff}(t)$, as
in Eq. (11). Consequently, let

$G(s,s^{\prime})=V(s)V^{-1}(s^{\prime})$ (16)

be the transition operator between any $s,s^{\prime}\leq t$. From Eqs. (11)
and (14) we have $|\psi(t)\rangle=G(t,0)|\psi_{0}\rangle$ and

$|\psi^{\tau_{1},\cdots,\tau_{\ell}}(t)\rangle=G(t,\tau)|\psi^{\tau_{1},\cdots,\tau_{\ell}}(\tau)\rangle,$
(17)

for any $0\leq\tau_{1}\leq\cdots\leq\tau_{\ell}\leq\tau\leq t$. Repeatedly
using Eq. (17) and Eq. (15), we have

$\displaystyle|\psi^{\tau_{1},\cdots,\tau_{\ell}}(t)\rangle$ $\displaystyle=$
$\displaystyle
G(t,\tau_{\ell})|\psi^{\tau_{1},\cdots,\tau_{\ell}}(\tau_{\ell})\rangle$ (18)
$\displaystyle=$ $\displaystyle
G(t,\tau_{\ell})L(\tau_{\ell})|\psi^{\tau_{1},\cdots,\tau_{\ell-1}}(\tau_{\ell})\rangle$
$\displaystyle\vdots$ $\displaystyle=$ $\displaystyle
G(t,\tau_{\ell})L(\tau_{\ell})G(\tau_{\ell},\tau_{\ell-1})L(\tau_{\ell-1})G(\tau_{\ell-1},\tau_{\ell-2})\cdots
L(\tau_{1})G(\tau_{1},0)|\psi_{0}\rangle.$

The compact form (12) is obtained from Eq. (18) by Eq. (16). $\Box$

######  Remark 3.1.

Theorem 3.1 characterizes quantum jumps that is well-known in the theory of
open quantum systems [45], along which single photons are sequentially
released at random time instants $\tau_{1},\ldots,\tau_{\ell}$. These quantum
jumps intervene in the non-unitary evolution (14) of
$|\psi^{\tau_{1},\cdots,\tau_{\ell}}(t)\rangle$ by
${L}(\tau_{1}),\ldots,{L}(\tau_{\ell})$, forming stochastic and discontinuous
evolution trajectories that can be numerically calculated to analyze and
design the flying-qubit generation processes.

Theorem 3.1 provides the entire solution for the transient process of
generating flying qubits, which is reduced to Eq. (11) on the space of the
standing quantum system. To calculate the steady-state of the outgoing flying
qubits, we have the following corollary.

######  Corollary 3.1.

If the standing system asymptotically decays to its ground state $|0\rangle$,
then the outgoing flying-qubit state can be decomposed into the superposition
of multi-photon states represented by Eq. (6), where

$\xi^{\tau_{1},\cdots,\tau_{\ell}}=\langle
0|V(\infty)\tilde{L}(\tau_{\ell})\cdots\tilde{L}(\tau_{1})|\psi_{0}\rangle$
(19)

for $\ell=0,1,2,\cdots$.

_Proof._ The steady-state of the outgoing flying qubits can be obtained from
the asymptotic limit $t\rightarrow\infty$ of the joint state
$|\Psi(t)\rangle$. When the atom decays to its ground state $|0\rangle$, only
ground-state components in all $\ket{\psi^{\tau_{1},\cdots,\tau_{\ell}}(t)}$
may survive. Hence, the asymptotic state of $|\Psi(t)\rangle$ can be
decomposed as

$|\Psi(\infty)\rangle=|0\rangle\otimes\sum_{\ell=0}^{\infty}\int_{0}^{\infty}\int_{0}^{\tau_{\ell}}\cdots\int_{0}^{\tau_{2}}\xi^{\tau_{1},\cdots,\tau_{\ell}}dB^{{\dagger}}_{\tau_{1}}\cdots
dB^{{\dagger}}_{\tau_{\ell}}\ket{\Omega}_{F},$ (20)

where $\xi^{\tau_{1},\cdots,\tau_{\ell}}={}\langle
0|\psi^{\tau_{1},\cdots,\tau_{\ell}}(\infty)\rangle$. The outgoing flying-
qubit state can be directly observed from the flying-qubit part of (19).
$\Box$

######  Remark 3.2.

The assumption made in Corollary 3.1, i.e., the standing system decays to its
ground state, always holds when the system is a two-level atom that is
constantly coupled to the waveguide and is driven by a finite-duration driving
pulse. Under this circumstance, the standing system will spontaneously decay
to its ground state after the control pulse is turned off, and eventually is
disentangled from the flying qubits. It is possible that, under elaborately
designed coherent and incoherent controls, or with multi-level atoms that have
non-unique stable states (e.g., a $\Lambda$-type atom), the standing system
remains entangled with the flying qubits. This has many intriguing
applications in quantum networks, e.g., distributing quantum entanglements
between remote nodes in the quantum network.

###  3.2 Standing system with non-vacuum quantum inputs

The above analysis is well suited for the generation of flying qubits, but it
is not directly applicable to the catching or transformation problems because
the quantum input of incoming flying qubits is not initially empty. As
schematically shown in Fig. 2, we can introduce an auxiliary system $A$ that
is cascaded to the standing system $B$, where the auxiliary system is
engineered to generate a flying qubit at the same state of the incoming flying
qubit [19]. In this way, the analysis can be equivalently done with the joint
$AB$ system with vacuum quantum inputs, to which Theorem 3.1 and Corollary 3.1
become applicable.

Figure 2: Schematics of a flying-qubit control system $B$ with quantum inputs
generated from an auxiliary system $A$.

Let $H_{A}(t)$ and $H_{B}(t)$ be the Hamiltonians of the $A$-system and
$B$-system, and $L_{A}(t)$ and $L_{B}(t)$ be their coupling operators to the
waveguide, respectively. According to the $(S,L,H)$ formalism presented in
[18], the series product of systems $A$ and $B$ has the following equivalent
Hamiltonian

$H(t)=H_{A}(t)\otimes\mathbb{I}_{B}+\mathbb{I}_{A}\otimes
H_{B}(t)+\frac{1}{2i}\left[L_{A}(t)\otimes
L_{B}^{\dagger}(t)-L_{A}^{\dagger}(t)\otimes L_{B}(t)\right],$ (21)

and the coupling operator

$L(t)=L_{A}(t)\otimes\mathbb{I}_{B}+\mathbb{I}_{A}\otimes L_{B}(t).$ (22)

The series product representation makes it possible to apply Theorem 3.1 to
analyze flying-qubit catching or transformation problems, because the joint
$AB$ system receives vacuum quantum inputs via the subsystem $A$. Note that
the auxiliary system $A$ needs to be carefully designed so as to yield the
incoming flying qubits to be caught or transformed.

##  4 Examples

This section will illustrate how the proposed model is applied to the control
of flying qubits. For simplicity, the standing system is chosen as a two-level
qubit, whose controlled Hamiltonian is

$\displaystyle H(t)$ $\displaystyle=$
$\displaystyle\epsilon(t)\sigma_{+}\sigma_{-}+\frac{u(t)}{2}\sigma_{+}+\frac{u^{*}(t)}{2}\sigma_{-},$
(23) $\displaystyle L(t)$ $\displaystyle=$
$\displaystyle\sqrt{2\gamma(t)}\sigma_{-},$ (24)

where $\sigma_{\pm}$ are the standard Pauli raising and lowering operators.
The functions $\epsilon(t)$ and $u(t)$ are the detuning frequency and the
envelope of the coherent driving field on the qubit. The tunable coupling
strength $\gamma(t)$ plays the role of incoherent classical controls. In the
following, we discuss the flying-qubit generation and transformation problems,
which correspond to the cases with vacuum and non-vacuum quantum inputs,
respectively.

###  4.1 The generation of flying qubits

We will first apply Theorem 3.1 to reproduce a well-known result of generating
arbitrary single-photon wavepackets via the incoherent control of tunable
waveguide coupling function [17, 20, 19], with a slight generalization in that
a time-variant phase is allowed in the wavepacket function.

######  Proposition 4.1.

Given a differentiable single-photon wavepacket
$\xi^{\tau}=\nu(\tau)e^{i\phi(\tau)}$, where $\nu(\tau)$ and $\phi(\tau)$ are
the amplitude and phase functions of $\xi^{\tau}$. Then, a flying qubit at
state $\ket{1_{\xi}}$ can be generated without coherent driving $u(t)$ when
the standing qubit is initially prepared at the excited state $|1\rangle$ and

$\gamma(\tau)=\frac{\nu^{2}(\tau)}{2\int_{\tau}^{\infty}\nu^{2}(s)ds},\quad\epsilon(\tau)=-\frac{d\phi(\tau)}{d\tau}.$
(25)

_Proof._ According to Theorem 3.1 and Corollary 3.1, it is straightforward to
calculate the wavepacket of outgoing single-photon state as follows:

$\xi^{\tau}=\sqrt{2\gamma(\tau)}\exp\left[-\int_{0}^{\tau}\gamma(s)ds\right]\exp\left[-i\int_{0}^{\tau}\epsilon(s)ds\right].$
(26)

Comparing with the shape function $\xi^{\tau}=\nu(\tau)e^{i\phi(\tau)}$, we
have $-\int_{0}^{t}\epsilon(\tau)d\tau=\phi(t)$ and

$\nu(\tau)=\sqrt{2\gamma(\tau)}\exp\left[-\int_{0}^{\tau}\gamma(s)ds\right]$
(27)

for all $0<\tau<t$. Combined with the normalization condition of $\nu(\tau)$,
the amplitude and phase conditions in the proposition can be derived. $\Box$

The proof of Proposition 4.1 is relatively simple because the controls
$\gamma(t)$ and $\epsilon(t)$ preserve the number of excitations in the joint
system. Next, we consider the case with a coherent driving pulse that does not
preserve the number of excitations. Assume that the driving control posed
along the $x$-axis is resonant with the qubit (i.e., $\epsilon(t)\equiv 0$)
and the coupling strength $\gamma(t)\equiv\gamma$. Under a rectangular driving
control pulse on $0\leq t\leq T<\infty$, the effective Hamiltonian can be
written as

$H_{\rm
eff}(t)=\begin{cases}\frac{\Omega}{2}\sigma_{x}-i\frac{\gamma}{2}\sigma_{+}\sigma_{-},&0<t<T\\\
-i\frac{\gamma}{2}\sigma_{+}\sigma_{-},&t>T\end{cases}$ (28)

where $\sigma_{x}=\sigma_{+}+\sigma_{-}$ and $\Omega$ is the power of the
driving pulse.

Depending on the magnitude of the driving power $\Omega$, we discuss the
outgoing flying-qubit states in three regimes, namely the strong-driving
regime ($\Omega>\gamma$), the balanced-driving regime ($\Omega=\gamma$) and
the weak-driving regime ($\Omega<\gamma$). According to Theorem 3.1 and
Corollary 3.1, we can obtain, for the example of strong-driving case, that

$\xi=\langle 0|V(\infty)|0\rangle=e^{-\frac{\gamma
T}{2}}\cos\left(\frac{\omega T}{2}-\phi\right)$ (29)

and the single-photon packet

$\xi^{\tau}=\left\\{\begin{array}[]{ll}-\frac{i\sqrt{2\gamma}e^{-\frac{\gamma
T}{2}}}{\cos^{2}\varphi}\sin\frac{\omega{\tau}}{2}\cos\left(\frac{\omega(T-{\tau})}{2}-\varphi\right),&{\tau}<T\\\
-\frac{i\sqrt{2\gamma}e^{-\frac{\gamma T}{2}}}{\cos\varphi}\sin{\frac{\omega
T}{2}}e^{-\gamma({\tau}-T)},&{\tau}>T,\end{array}\right.$ (30)

where $\varphi=\arcsin{\frac{\gamma}{\Omega}}$ and
$\omega=\sqrt{\Omega^{2}-\gamma^{2}}$. The expressions for multi-photon
emission, as well as those in balanced-driving and weak-driving regimes, can
also be analytically solved under retangular driving pulses, but they will not
be listed here due to the limit of space.

Based on the above calculation, we display in Fig. 3(a) the shape
$|\xi^{\tau}|^{2}$ of the single-photon component generated by soft
$\pi$-pulses (i.e., $\Omega T=\pi$) in the strong, balanced and weak driving
regimes, respectively. All these wavepackets gradually rise from zero and the
steepness of the rising slope is proportional to the driving strength
$\Omega$. After the driving pulse is turned off, they spontaneously decay to
zero. We also calculate the probabilities of generating $0$, $1$, $2$ and more
photons according to Eq. (7) when the pulse power $\Omega$ increases from
weak- to strong-driving regimes. As is shown in Fig. 3(b), soft $\pi$-pulses
do not generate perfect single flying qubits due to the coexisting multi-
photon emission processes, unless in the hard-pulse limit that $\Omega$
approaches infinity.

Figure 3: (a) The shapes of single-photon wavepackets generated by
$\pi$-pulses in the strong, balanced and weak driving regimes; (b) the
probabilities of generating $0$, $1$, $2$ and more flying qubits under
$\pi$-pulses with driving strength varying from weak to strong coupling
regimes.

###  4.2 The transformation of flying qubits

Now let us see how a standing qubit responses to non-vacuum quantum inputs,
e.g., an incoming flying qubit at state $\ket{1_{\xi_{0}}}$, where
$\xi_{0}^{\tau}$ is its wavepacket. Due to the limit of space, we only discuss
the transformation problems here, as the catching problem can be taken as one
of its special case.

Suppose that the standing qubit is prepared at the ground state $|0\rangle$
and is controlled by the tunable coupling $\gamma(t)$ and the coherent driving
$u(t)$. We introduce the auxiliary two-lelvel system $A$ shown in Fig (2) to
generate the incoming flying qubit, which is initially excited and its
Hamiltonian $H_{0}(t)=\epsilon_{0}(t)\sigma_{+}\sigma_{-}$ and coupling
operator $L_{0}(t)=\sqrt{2\gamma_{0}(t)}\sigma_{-}$ are chosen according to
Proposition 4.1. From Eq. 21, the effective Hamiltonian of the $AB$ system is

$H_{\rm
eff}(t)=\frac{u(t)}{2}\sigma_{+}+\frac{u^{*}(t)}{2}\sigma_{-}-i\left[\gamma_{A}(t)\sigma_{+}\sigma_{-}\otimes\mathbb{I}_{2}+\mathbb{I}_{2}\otimes\gamma_{B}(t)\sigma_{+}\sigma_{-}+2\sqrt{\gamma_{0}(t)\gamma(t)}\sigma_{-}\otimes\sigma_{+}\right]$
(31)

and
$L(t)=\sqrt{2\gamma_{0}(t)}\sigma_{-}\otimes\mathbb{I}_{2}+\sqrt{2\gamma(t)}\mathbb{I}_{2}\otimes\sigma_{-}$.
The joint $AB$ system is initially prepared at state
$|\psi_{AB}(0)\rangle=|10\rangle$.

Consider the simple case that the coherent driving is a hard $\pi$ pulse
(i.e., $u(t)=\pi\delta(t)$) under which the system $B$ is instantaneously
flipped from $|0\rangle$ to $|1\rangle$, and hence $|\psi_{AB}(t)\rangle$ is
flipped from $|10\rangle$ to $|11\rangle$ at $t=0$. Hence, we only need to
calculate the joint system’s dynamics starting from the flipped state
$|\psi_{AB}(0)\rangle=|11\rangle$. It can be verified from Corollary 3.1 that
the wavepacket functions $\xi^{\tau_{1},\cdots,\tau_{\ell}}$ in the steady
state all vanish except when $\ell=2$, and the two-photon wavepacket can be
expressed as follows:

$\displaystyle\xi^{\tau_{1},\tau_{2}}=$ $\displaystyle
2\sqrt{{\gamma}_{0}(\tau_{2}){\gamma}(\tau_{1})}e^{-{\Gamma}_{A}(\tau_{2})-{\Gamma}_{B}(\tau_{1})}+2\sqrt{{\gamma}_{0}(\tau_{1}){\gamma}(\tau_{2})}e^{-{\Gamma}_{A}(\tau_{1})-{\Gamma}_{B}(\tau_{2})}$
$\displaystyle+2\sqrt{{\gamma}(\tau_{1}){\gamma}(\tau_{2})}e^{-{\Gamma}_{B}(\tau_{1})-{\Gamma}_{B}(\tau_{2})}\left[\Xi(\tau_{1})-\Xi(\tau_{2})\right],$
(32)

where

$\Gamma_{A,B}(t)=\int_{0}^{t}\gamma_{A,B}(s)ds,\quad\Xi(t)=2\int_{0}^{t}\sqrt{\gamma_{0}(s)\gamma(s)}e^{\Gamma_{B}(s)-\Gamma_{A}(s)}ds.$
(33)

Recall that the system $A$ is chosen to generate the desired single flying-
qubit wavepacket, i.e.,
$\xi_{0}^{\tau}=\sqrt{2{\gamma}_{0}(\tau)}e^{-{\Gamma}_{A}(\tau)}$. Let
$\xi_{1}^{\tau}=\sqrt{2\gamma(\tau)}e^{-\Gamma_{B}(\tau)}$ be the wavepacket
of a single flying-qubit spontaneously emitted from the system $B$ driven by
vacuum input. We can rewrite Eq. (4.2) in a more compact form:

$\xi^{\tau_{1},\tau_{2}}=\xi_{0}^{\tau_{1}}\xi_{1}^{\tau_{2}}+\xi_{0}^{\tau_{2}}\xi_{1}^{\tau_{1}}+\xi_{1}^{\tau_{1}}\xi_{1}^{\tau_{2}}\int_{\tau_{1}}^{\tau_{2}}\xi_{0}^{s}\xi_{1}^{s}\left[\int_{0}^{s}|\xi_{1}^{\alpha}|^{2}d\alpha\right]^{-1}.$
(34)

This form indicates that the incoming single photon knocks out an additional
photon from the excited qubit, as is also called stimulated emission [38], and
the two single photons are entangled in a complex manner.

Now we consider the soft driving control pulses, and assume that the incoming
single photon has exponentially decaying waveform
$\xi_{0}^{\tau}=\sqrt{2\gamma_{0}}e^{-\gamma_{0}\tau}$ (corresponding to
$\epsilon(t)=0$ and $\gamma_{0}(t)\equiv\gamma_{0}$ for the auxiliary system
$A$). The coherent driving control $u(t)=u_{x}(t)+iu_{y}(t)$ on the system $B$
is turned on over a fixed time interval to reshape the flying qubit, and the
joint system is initially prepared at state $\ket{\Psi(0)}=\ket{10}$. We apply
the genetic algorithm to optimize the waveform of $u(t)$ that maximizes the
single-photon emission probability $P_{1}$. Numerical simulations show that
the coherent control can improve the single-photon generation probability to
$95.22\%$ after 200 generations, yielding the optimized control pulse
displayed in Fig. 4(a). As is shown in Fig. 4(b), the shape of the outgoing
single-photon component is suppressed in the beginning and then revives,
showing the ability of the driving control on reshaping the incoming flying
qubit.

Figure 4: The transformation of a flying qubit with (a)-(b) optimized control
pulses and (c) the shapes of the incoming and outgoing flying qubits.

##  5 Conclusion

To conclude, we have developed a QSDE-based modeling method for the control of
flying qubits in networked quantum information processing systems. Under the
time-ordered photon-number basis, the proposed model can describe most general
cases when the standing quantum system is steered by coherent or incoherent
classical controls, and the quantum inputs can be either vacuum or non-vacuum.
As demonstrated by examples, the states of the outgoing flying qubits can be
efficiently computed by the derived low-dimensional non-unitary Schrödinger
equation and the associated quantum jumps. This QSDE-based approach can be
naturally applied to the analysis flying-qubit control dynamics in quantum
networks by combining the $(S,L,H)$ formulation.

The established framework also lays the foundation for practical design of
flying-qubit control systems, such as the precise shaping of flying qubits,
the suppression of undesired multi-photon emissions, or robust control of
flying qubits. As is exemplified in the numerical example, these design
problems can be formulated as properly defined optimal control problems, and
efficient algorithms need to be developed. These problems have been under
consideration in our future studies.

The proposed model can be naturally extended to cases when the standing
quantum system has multiple levels or is coupled to multiple waveguides.
However, although the analysis of flying-qubit catching and transformation
problems can still be handled by introducing multiple auxiliary systems
coupled to these waveguides, the resulting calculation will become extremely
hard as the equivalent system defined on the tensor of multiple Hilbert spaces
become large. Hence, a direct method that is free of auxiliary systems is
highly demanded. This has been shown to be possible from quantum scattering
theory [24, 44], but the QSDE-based approach, which is potentially
advantageous in the applications to complex networks, is still an open
question. This will be also explored in our future works.

## References

  * [1] V. Averchenko, D. Sych, C. Marquardt, and G. Leuchs.  Efficient generation of temporally shaped photons using nonlocal spectral filtering.  Physical Review A, 101(1):013808, 2020.
  * [2] V. Averchenko, D. Sych, G. Schunk, U. Vogl, C. Marquardt, and G. Leuchs.  Temporal shaping of single photons enabled by entanglement.  Physical Review A, 96:043822, 2017.
  * [3] J. I. Cirac, P. Zoller, H. J. Kimble, and H. Mabuchi.  Quantum state transfer and entanglement distribution among distant nodes in a quantum network.  Phys. Rev. Lett., 78:3221–3224, 1997.
  * [4] D. D’Alessandro.  Introduction to quantum control and dynamics.  Chapman & Hall/CRC, Taylor & Francis Group, 2008.
  * [5] D. DiVincenzo.  Quantum computation.  Science, 270:255–261, 10 1995.
  * [6] D. Dong and I. R. Petersen.  Quantum control theory and applications: a survey.  IET Control Theory Applications, 4(12):2651–2671, 2010.
  * [7] Z. Dong, L. Cui, G. Zhang, and H. Fu.  Wigner spectrum and coherent feedback control of continuous-mode single-photon fock states.  Journal of Physics A: Mathematical and Theoretical, 49(43):435301, 2016.
  * [8] Z. Dong, G. Zhang, and N. Amini.  Quantum filtering for a two-level atom driven by two counter-propagating photons.  Quantum Information Processing, 18(5):136, 2019.
  * [9] L.-M. Duan, M. D. Lukin, J. I. Cirac, and P. Zoller.  Long-distance quantum communication with atomic ensembles and linear optics.  Nature, 414(6862):413–418, 2001.
  * [10] M. D. Eisaman, J. Fan, A. Migdall, and S. V. Polyakov.  Single-photon sources and detectors.  Review of Scientific Instruments, 82(7):071101, 2011.
  * [11] K. Fischer, R. Trivedi, V. Ramasesh, I. Siddiqi, and J. Vučković.  Scattering into one-dimensional waveguides from a coherently-driven quantum-optical system.  Quantum, 2:69, 2018.
  * [12] P. Forn-Diaz, C. W. Warren, C. W. S. Chang, A. M. Vadiraj, and C. M. Wilson.  On-demand microwave generator of shaped single photons.  Physical Review Applied, 8(5):054015, 2017.
  * [13] C. W. Gardiner and M. J. Collett.  Input and output in damped quantum systems: Quantum stochastic differential equations and the master equation.  Physical Review A, 31:3761–3774, 1985.
  * [14] K. M. Gheri, K. Ellinger, T. Pellizzari, and P. Zoller.  Photon-wavepackets as flying quantum bits.  Fortschritte der Physik, 46:401–415, 1998.
  * [15] S. J. Glaser, U. Boscain, T. Calarco, C. P. Koch, W. Kockenberger, R. Kosloff, I. Kuprov, B. Luy, S. Schirmer, T. Schulte-Herbruggen, D. Sugny, and F. K. Wilhelm.  Training schrodinger’s cat: quantum optimal control strategic report on current status, visions and goals for research in europe.  European Physical Journal D, 69(12):279, 2015.
  * [16] J. Gough, M. James, and H. Nurdin.  Quantum trajectories for a class of continuous matrix product input states.  New Journal of Physics, 16(7):075008, 2014.
  * [17] J. Gough, M. James, H. Nurdin, and J. Combes.  Quantum filtering for systems driven by fields in single-photon states or superposition of coherent states.  Physical Review A, 86:043819, 2012.
  * [18] J. Gough and M. R. James.  The series product and its application to quantum feedforward and feedback networks.  IEEE Transactions on Automatic Control, 54(11):2530–2544, 2009\.
  * [19] J. Gough and G. Zhang.  Generating nonclassical quantum input field states with modulating filters.  EPJ Quantum Technology, 2(15), 2015.
  * [20] J. E. Gough, M. R. James, and H. I. Nurdin.  Quantum filtering for systems driven by fields in single photon states and superposition of coherent states using non-markovian embeddings.  Quantum Information Processing, 12(3):1469–1499, 2013.
  * [21] X. Gu, A. Kockum, A. Miranowicz, Y.-x. Liu, and F. Nori.  Microwave photonics with superconducting quantum circuits.  Physics Reports, 718:1–102, 2017.
  * [22] A. Houck, D. Schuster, J. Gambetta, J. Schreier, B. Johnson, J. Chow, L. Frunzio, J. Majer, M. Devoret, S. Girvin, and R. Schoelkopf.  Generating single microwave photons in a circuit.  Nature, 449:328–331, 10 2007.
  * [23] R. L. Hudson and K. R. Parthasarathy.  Quantum ito’s formula and stochastic evolutions.  Communications in Mathematical Physics, 93(3):301–323, 1984.
  * [24] D. Hurst and P. Kok.  Analytic few-photon scattering in waveguide QED.  Physical Review A, 97:043850, 2018.
  * [25] K. Jacobs.  Quantum measurement theory and its applications.  Cambridge University Press, 2014.
  * [26] M. Keller, B. Lange, K. Hayasaka, W. Lange, and H. Walther.  Continuous generation of single photons with controlled waveform in an ion-trap cavity system.  Nature, 431(7012):1075–1078, 2004.
  * [27] A. Kuhn, M. Hennrich, and G. Rempe.  Deterministic single-photon source for distributed quantum networking.  Physical Review Letters, 89:067901, 09 2002.
  * [28] C. Kurtsiefer, S. Mayer, P. Zarda, and H. Weinfurter.  Stable solid-state source of single photons.  Physical Review Letters, 85(2):290–293, 2002.
  * [29] V. Leong, M. Seidler, M. Steiner, A. Cere, and C. Kurtsiefer.  Time-resolved scattering of a single photon by a single atom.  Nature Communications, 7:13716, 2016.
  * [30] W.-L. Li, G. Zhang, and R.-B. Wu.  The dynamical model of flying-qubit control systems.  In 20th IFAC World Congress, volume 50, pages 11755 – 11759, 2020\.  20th IFAC World Congress.
  * [31] M. Lucamarini, G. Di Giuseppe, D. Vitali, and P. Tombesi.  Open-loop and closed-loop control of flying qubits.  Journal of Physics B, 44(15):154005, 2011.
  * [32] H. I. Nurdin, M. R. James, and N. Yamamoto.  Perfect single device absorber of arbitrary traveling single photon fields with a tunable coupling parameter: A qsde approach.  In 55th IEEE Conference on Decision and Control (CDC), pages 2513–2518, 2016.
  * [33] Y. Pan, G. Zhang, W. Cui, and M. R. James.  Single photon inverting pulse for an atom in a cavity.  In 54th IEEE Conference on Decision and Control (CDC), pages 6429–6433, 2015.
  * [34] M. Pechal, L. Huthmacher, C. Eichler, S. Zeytinoglu, A. Abdumalikov, S. Berger, A. Wallraff, and S. Filipp.  Microwave-controlled generation of shaped single photons in circuit quantum electrodynamics.  Physical Review X, 4(4):041010, 2013.
  * [35] Z. H. Peng, S. E. De Graaf, J. S. Tsai, and O. V. Astafiev.  Tuneable on-demand single-photon source.  Nature Communications, 7(12588), 2016.
  * [36] M. Pierre, I. Svensson, S. Sathyamoorthy, G. Johansson, and P. Delsing.  Storage and on-demand release of microwaves using superconducting resonators with tunable coupling.  Applied Physics Letters, 104(23):232604, 2014.
  * [37] J. Preskill.  Quantum Computing in the NISQ era and beyond.  Quantum, 2:79, 2018.
  * [38] E. Rephaeli and S. Fan.  Stimulated emission from a single excited atom in a waveguide.  Physical Review Letters, 108:143602, 2012.
  * [39] J. Shen and S. Fan.  Theory of single-photon transport in a single-mode waveguide: I. coupling to a cavity containing a two-level atom.  Physical Review A, 79(2):023837, 2009.
  * [40] H. T. Song, G. F. Zhang, and Z. R. Xi.  Continuous-mode multi-photon filtering.  SIAM Journal on Control and Optimization, 54(3):1602–1632, 2016\.
  * [41] B. Srivathsan, G. Gulati, A. Cere, B. Chng, and C. Kurtsiefer.  Reversing the temporal envelope of a heralded single photon using a cavity.  Physical Review Letters, 113(16):163601, 2014.
  * [42] D. Stefanatos and E. Paspalakis.  A shortcut tour of quantum control methods for modern quantum technologies.  EPL (Europhysics Letters), 132(6):60001 (7pp), 2020.
  * [43] M. Stobinska, G. Alber, and G. Leuchs.  Perfect excitation of a matter qubit by a single photon in free space.  Europhysics Letters, 86(1):14007, 2009.
  * [44] Rahul Trivedi, Kevin Fischer, Shanshan Xu, Shanhui Fan, and Jelena Vuckovic.  Few-photon scattering and emission from low-dimensional quantum systems.  Phys. Rev. B, 98:144112, Oct 2018.
  * [45] H. Wiseman and G. Milburn.  Quantum measurement and control.  Cambridge University Press, 2010.
  * [46] W. Yao, R.-B. Liu, and L. J. Sham.  Theory of control of the spin-photon interface for quantum networks.  Physical Review Letters, 95:030504, 2005.
  * [47] G. Zhang.  Analysis of quantum linear systems’s response to multi-photon states.  Automatica, 50(2):442 – 451, 2014.
  * [48] G. Zhang.  Single-photon coherent feedback control and filtering, pages 1–4.  Springer, London, 2020.
  * [49] G. Zhang and M.R. James.  On the response of quantum linear systems to single photon input fields.  IEEE Transactions on Automatic Control, 58(5):1221–1235, 2011.