# Multiphonon quantum dynamics in cavity optomechanical systems

Mihai A. Macovei<EMAIL_ADDRESS>Institute of Applied Physics, Academiei
str. 5, MD-2028 Chişinău, Moldova  Adriana Pálffy<EMAIL_ADDRESS>Department of Physics, Friedrich-Alexander-Universität Erlangen-Nürnberg,
91058 Erlangen, Germany

###### Abstract

The multiphonon quantum dynamics in laser-pumped cavity optomechanical samples
containing a vibrating mirror is investigated. Especially, we focus on
dispersive interaction regimes where the externally applied coherent field
frequency detuning from the optical resonator frequency is not equal to the
mirror’s oscillating frequency or to its multiples. As a result, for
moderately strong couplings among the involved subsystems, the quantum
dynamics of this complex system is described by multiphonon absorption or
emission processes, respectively. Particularly, we demonstrate efficient ways
to monitor the phonon quantum dynamics via photon detection. The possibility
to extract the relevant sample parameters, for instance, the coupling strength
between the mechanical mirror and the electromagnetic field, is also
discussed.

##  I Introduction

Cavity optomechanical systems that couple electromagnetic field radiation with
nanomechanical or micromechanical motion of a vibrating component have proven
their potential in numerous different applications, e.g., optical networks,
quantum memories, quantum metrology, quantum amplifiers, quantum sensing or
gravimetry walls_milb ; optm1 ; qsens ; brag ; optm ; natc ; grav . Many of
these applications, whether advancing fundamental quantum physics or
technological in nature, rely on the fact that the photon-phonon coupling
renders possible the quantum cooling of quantized motion cool ; cool1 ; cool2
; cool3 ; cool4 ; cool5 . Cavity optomechanical systems also allow the
creation and control of quantum macroscopic Schrödinger cat states in cavities
with a moving mirror schr1 ; schr2 . Furthermore, the output of an externally
pumped cavity optomechanical system shows a clear evidence of
electromagnetically induced transparency phenomenon, called in this case
optomechanically induced transparency eit1 ; eit2 ; eit3 . Based on this
effect, a photon switch effect was demonstrated eits . Transferred to a
different frequency regime by a proper optomechanical interface,
optomagnetically induced transparency may help to control light-matter
interaction of x-rays via optical photons eitx .

Generally speaking, optomechanical systems may be viewed as pumped optical
resonators containing Kerr-like nonlinear elements ker1 . As a consequence,
the cavity optomechanical samples exhibit bi- or multi-stability, multiple
photon blockades, and various types of entanglement or squeezing phenomena
ker1 ; ker2 ; marq ; mma ; mblock ; ker3 ; ker4 ; ker5 . The photon blockade
effect in optomechanical systems, i.e. preventing multiple photons from
entering the cavity at the same time due to strong photon-photon interactions,
was considered from theory side in Ref. phblk1 . Cross-Kerr interactions
between the optical and phonon modes lead to generation of mechanical cat
states mecat , for instance. Since photons typically do not mutually interact,
enhanced photon-photon interactions in these systems would be of great
interest in quantum computation and quantum information processing luk1 ; luk2
; chuang . Furthermore, in the single-photon strong-coupling regime and good-
cavity limit, the cavity response shows several resolved resonances at
multiples of the mechanical frequency, respectively phblk2 . This effect can
be used for instance to measure the mechanical frequency if other involved
parameters are known. Also, strong laser driving in an optomechanical setup
was analyzed in stlas , where a transition from sub-Poissonian to super-
Poissonian photon statistics was determined.

Most of the above mentioned results were obtained under the condition that the
externally applied coherent field frequency detuning from the optical
resonator frequency is equal or close to the mirror’s oscillating frequency.
In contrast, here we shall focus on the different regime when the above
mentioned resonance condition is not met. The laser frequency detuning from
the cavity one is considered unequal to the frequency of the mirror’s
vibrations or to its multiples. In addition, we consider the mechanical
frequency is larger than the detuning and than the corresponding decay rates
in the sample, and is commensurable with but larger than the coupling strength
among the interacting subsystems. We demonstrate that under these
circumstances, the quantum dynamics of the oscillating mirror has a multi-
phonon nature in the sense that the quantities describing its evolution in
steady-state are characterized by absorption or emission of many mechanical
oscillation quanta. We show that in the proposed setup the corresponding
phonon dynamics, proper to the mechanical part, follows that of the photon
one. The calculation of the mean-phonon as well as mean-photon numbers
demonstrated that detecting the leaking photons from the optical cavity one
can monitor the vibrations of the moving subsystem. Correspondingly, the
quantum nature of these processes is established through the second-order
phonon-phonon or photon-photon correlation functions. Furthermore, our results
show that the sample’s parameters like the phonon-photon coupling strengths
can be extracted from the multi-peak structure of the mean-photon number
quantum dynamics. Finally, the parameter range needed to observe this behavior
is close to those for photon blockade effect in optomechanical systems phblk1
; phblk2 , and within reach of experiments exp1 ; exp2 ; exp3 .

This paper is organized as follows. In Sec. II we describe the analytical
approach and the system of interest, while in Sec. III we present and analyse
the obtained numerical results. The article concludes with a summary and
outlook in Sec. IV.

##  II Analytical approach

We describe our sample using the master equation approach under Born-Markov
approximations where the coupling of the relevant degrees of freedom to their
environmental counterparts is weak, whereas the photon and phonon memory
effects are negligible walls_milb ; reww . Appropriate unitary transformations
performed further will allow us to follow the multi-phonon quantum dynamics of
the mechanical part or the corresponding photon dynamics as a function of the
ratio of the coupling strength $g$ over mechanical oscillation frequency
$\omega$. This way one can distinguish the corresponding quantum dynamics of
the photon-phonon subsystems, respectively, when single- or many-phonons are
involved.

The master equation describing a standard laser pumped cavity optomechanical
setup, in the Born-Markov approximations and in a frame rotating at the
external laser field frequency $\omega_{L}$ is given by

$\displaystyle\frac{d}{dt}\rho+\frac{i}{\hbar}[H,\rho]$ $\displaystyle=$
$\displaystyle-\frac{\kappa_{a}}{2}[a^{\dagger},a\rho]-\frac{\kappa_{b}}{2}(1+\bar{n})[b^{\dagger},b\rho]$
(1) $\displaystyle-$
$\displaystyle\frac{\kappa_{b}}{2}\bar{n}[b,b^{\dagger}\rho]+H.c.,$

where $\rho$ is the density matrix and the Hamiltonian $H$ is given by the
expression

$\displaystyle H=\hbar\Delta a^{\dagger}a+\hbar\omega
b^{\dagger}b+\hbar\epsilon(a^{\dagger}+a)+\hbar ga^{\dagger}a(b+b^{\dagger}).$
(2)

The coherent evolution of the examined system is described by the second term
of the left-side part of the master equation (1). The damping effects of the
involved photon and phonon subsystems are characterized by the right-side part
of this equation, with $\kappa_{a}$ and $\kappa_{b}$ being the corresponding
photon or phonon damping rates, respectively. Here,
$\bar{n}=(\exp[\hbar\omega/(k_{B}T)]-1)^{-1}$ is the mean-phonon number due to
the thermal bath environment at temperature $T$ and at the vibration frequency
$\omega$ of the mechanical resonator, while $k_{B}$ is the Bolzmann’s
constant. $a^{\dagger}(b^{\dagger})$ is the creation operator of a photon
(phonon), whereas $a(b)$ is the corresponding photon (phonon) annihilation
operator, respectively, satisfying the standard bosonic commutation relations
$[a,a^{\dagger}]=1$, $[b,b^{\dagger}]=1$, and
$[b^{\dagger},b^{\dagger}]=[b,b]=0$, $[a^{\dagger},a^{\dagger}]=[a,a]=0$. The
first and the second components from the Hamiltonian (2) account for the free
energies of the optical and mechanical resonators, respectively, with
$\Delta=\omega_{c}-\omega_{L}$ being the detuning of the optical cavity
frequency $\omega_{c}$ from the laser one. The third term in (2) describes the
laser pumping effects of the optical resonator’s mode with $\epsilon$ being
the corresponding amplitude. The last term accounts for the interaction among
the optical and mechanical motion degrees of freedom characterized by the
coupling strength $g$.

In the following, we consider that $\omega>g>\epsilon$, and do not yet impose
any conditions on $\kappa_{a}$ and $\kappa_{b}$. Note that when the
optomechanical coupling is comparable to or larger than the optical decay rate
and the mechanical frequency, the steady state of the mechanical oscillator
can develop a nonclassical strongly negative Wigner density nwign . So, we
proceed to perform a unitary transformation

$\displaystyle U=e^{\chi a^{\dagger}a(b-b^{\dagger})},$ (3)

in the master equation (1) leading to the following new bosonic operators

$\displaystyle\bar{b}^{\dagger}=\chi\bar{a}^{\dagger}\bar{a}+b^{\dagger},~{}~{}~{}\bar{b}=\chi\bar{a}^{\dagger}\bar{a}+b,$
(4)

and

$\displaystyle\bar{a}^{\dagger}=a^{\dagger}e^{\chi(\bar{b}-\bar{b}^{\dagger})},~{}~{}~{}\bar{a}=ae^{-\chi(\bar{b}-\bar{b}^{\dagger})}.$
(5)

Upon the unitary transformation (3), the master equation (1) takes the form:

$\displaystyle\frac{d}{dt}\bar{\rho}$ $\displaystyle+$
$\displaystyle\frac{i}{\hbar}[\bar{H},\bar{\rho}]=-\frac{\kappa_{a}}{2}[\bar{a}^{\dagger}e^{-\chi(\bar{b}-\bar{b}^{\dagger})},\bar{a}e^{\chi(\bar{b}-\bar{b}^{\dagger})}\bar{\rho}]$
(6) $\displaystyle-$
$\displaystyle\frac{\kappa_{b}}{2}(1+\bar{n})[(\bar{b}^{\dagger}-\chi\bar{a}^{\dagger}\bar{a}),(\bar{b}-\chi\bar{a}^{\dagger}\bar{a})\bar{\rho}]$
$\displaystyle-$
$\displaystyle\frac{\kappa_{b}}{2}\bar{n}[(\bar{b}-\chi\bar{a}^{\dagger}\bar{a}),(\bar{b}^{\dagger}-\chi\bar{a}^{\dagger}\bar{a})\bar{\rho}]+H.c.,$

where

$\displaystyle\bar{H}$ $\displaystyle=$
$\displaystyle\hbar\Delta\bar{a}^{\dagger}\bar{a}+\hbar\omega\bar{b}^{\dagger}\bar{b}+\hbar\epsilon\bigl{(}\bar{a}e^{\chi(\bar{b}-\bar{b}^{\dagger})}+\bar{a}^{\dagger}e^{-\chi(\bar{b}-\bar{b}^{\dagger})}\bigr{)}$
(7) $\displaystyle-$
$\displaystyle\omega\chi^{2}(\bar{a}^{\dagger}\bar{a})^{2},$

with

$\displaystyle\chi=\frac{g}{\omega}.$ (8)

The last term in the Hamiltonian (7) describes the vibration-induced Kerr-like
nonlinearity effects. Furthermore, the master equation (6) exhibits resonance
conditions if $\Delta=\pm k\omega$, $\\{k=1,2,\cdots\\}$. In what follows, we
shall focus on the case when $\Delta\not=\pm k\omega$, but rather
$\omega>\Delta$.

For the parameter regime of interest, we can proceed by expanding the
exponents in Eqs. (6-7) in the Taylor series using the small parameter
$\chi\ll 1$,

$\displaystyle
e^{\pm\chi(\bar{b}-\bar{b}^{\dagger})}=\sum^{\infty}_{n=0}\frac{(\pm\chi)^{n}}{n!}(\bar{b}-\bar{b}^{\dagger})^{n}.$
(9)

Further, using the bosonic operator identity

$\displaystyle(A+B)^{n}=\sum^{n}_{k}\frac{n!}{k!(\frac{n-k}{2})!}\bigl{(}-C/2\bigr{)}^{\frac{n-k}{2}}\sum^{k}_{r=0}\binom{k}{r}A^{r}B^{k-r},$
(10)

where $[A,B]=C$ and $[A,C]=[B,C]=0$, whereas $k$ is odd for an odd $n$ and
even for an even $n$, respectively, one can simplify the expression (9)
depending on the assumed approximations. Performing a unitary transformation

$\displaystyle
V=\exp[i(\Delta\bar{a}^{\dagger}\bar{a}+\omega\bar{b}^{\dagger}\bar{b})t],$
(11)

in the master equation (6) and avoiding any resonances in the system, i.e.
$\Delta\pm k\omega\not=0$, $\\{k=1,2,\cdots\\}$, one can neglect then all
time-dependent terms $\propto e^{\pm ik\omega t}$ in the master equation. This
approximation additionally requires that
$\omega\gg\\{\kappa_{a},\kappa_{b}\\}$. We keep those terms oscillating at
$e^{\pm i\Delta t}$, which accounts to further assuming that $\Delta/\omega\ll
1$. As a result, the exponent expression (9), which enters in the Hamiltonian
(7), for instance, takes the following form in this case:

$\displaystyle
e^{\pm\chi(\bar{b}-\bar{b}^{\dagger})}=\sum^{\infty}_{n=0}\sum^{n}_{m=0}\frac{(-1)^{m}\chi^{2n}}{(m!)^{2}(n-m)!}\frac{\bar{b}^{m}\bar{b}^{\dagger
m}}{2^{n-m}}.$ (12)

Respectively, the exponent expression entering in the optical resonator’s
damping in Eq. (6), i.e.

$\displaystyle
e^{\chi(\bar{b}-\bar{b}^{\dagger})}\bar{a}\bar{\rho}\bar{a}^{\dagger}e^{-\chi(\bar{b}-\bar{b}^{\dagger})}=$
$\displaystyle\sum^{\infty}_{n_{1}n_{2}=0}\frac{(-\chi)^{n_{1}}\chi^{n_{2}}}{n_{1}!n_{2}!}\bigl{(}\bar{b}^{\dagger}-\bar{b}\bigr{)}^{n_{1}}\bar{a}\bar{\rho}\bar{a}^{\dagger}\bigl{(}\bar{b}^{\dagger}-\bar{b}\bigr{)}^{n_{2}},$

acquires the form:

$\displaystyle
e^{\chi(\bar{b}-\bar{b}^{\dagger})}\bar{a}\bar{\rho}\bar{a}^{\dagger}e^{-\chi(\bar{b}-\bar{b}^{\dagger})}$
$\displaystyle=$
$\displaystyle\sum^{\infty}_{n_{1}n_{2}=0}\sum^{n_{1}}_{k_{1}}\sum^{n_{2}}_{k_{2}}\frac{(-\chi)^{n_{1}}\chi^{n_{2}}}{k_{1}!k_{2}!}\frac{(1/2)^{(n_{1}-k_{1})/2}(1/2)^{(n_{2}-k_{2})/2}}{[(n_{1}-k_{1})/2]![(n_{2}-k_{2})/2]!}\sum^{k_{1}}_{r_{1}=0}\sum^{k_{2}}_{r_{2}=0}(-1)^{r_{1}+r_{2}}\binom{k_{1}}{r_{1}}\binom{k_{2}}{r_{2}}$
(13) $\displaystyle\times$ $\displaystyle\bar{b}^{r_{1}}\bar{b}^{\dagger
k_{1}-r_{1}}\bar{a}\bar{\rho}\bar{a}^{\dagger}\bar{b}^{r_{2}}\bar{b}^{\dagger
k_{2}-r_{2}}e^{i(k_{1}-2r_{1})\omega t}e^{i(k_{2}-2r_{2})\omega t},$

where, again, $k_{i}$, $\\{i=1,2\\}$, is odd for an odd $n_{i}$ and even for
an even $n_{i}$, respectively. It is easy to observe that the above expression
is time-independent if

$k_{1}-2r_{1}+k_{2}-2r_{2}=0.$ (14)

Expressions (12-14) have to be introduced in the master equation (6) and one
can already recognize the multiphonon nature of the cavity optomechanical
dynamics in the off-resonance situation considered here.

Once we have arrived at a time-independent master equation, we can obtain the
corresponding equation for photon-phonon distribution function, namely,
$P_{n_{1}n_{2},m_{1}m_{2}}$=$\langle
n_{1},m_{1}|\bar{\rho}|m_{2},n_{2}\rangle$, where indices $n(m)$,
$\\{n,m=0,1,2,\cdots\\}$, refers to photon(phonon) subsystem, respectively,
with $|n(m)\rangle$ being the corresponding Fock state. Actually, that
equation is diagonal with respect to phonon degrees of freedom [see Eq. (6)
with expressions (12)-(14) and (20)]. Therefore, the final distribution
function will be represented as $P_{n_{1}n_{2},mm}=\langle
n_{1},m|\bar{\rho}|m,n_{2}\rangle$. Furthermore, the $N$-phonon processes are
described by terms proportional to $\chi^{2N}$, see Appendix A. In the
presence of corresponding damping effects, we can calculate the photon-phonon
distribution function $P_{n_{1}n_{2},mm}$ numerically in steady-state.

From the master equation (6), modified based on relations (12-14), one obtains
an infinite number of equations for the photon-phonon distribution function
$P_{n_{1}n_{2},mm}$. In order to solve the infinite system of equations for
$P_{n_{1}n_{2},mm}$, we truncate it at a certain maximum value
$\\{n=n_{max},m=m_{max}\\}$ so that a further increase of its value, i.e.
$\\{n_{max},m_{max}\\}$, does not modify the obtained results if other
involved parameters are being fixed. Thus, using the operator relations (4-5)
and keeping the time-independent terms only, the optical resonator’s steady-
state mean quanta number can be expressed as

$\displaystyle\langle
a^{\dagger}a\rangle=\sum^{n_{max}}_{n=0}\sum^{m_{max}}_{m=0}nP_{nm},$ (15)

while the mechanical vibrational mean-phonon number is

$\displaystyle\langle
b^{\dagger}b\rangle=\sum^{n_{max}}_{n=0}\sum^{m_{max}}_{m=0}(m+\chi^{2}n^{2})P_{nm}.$
(16)

Here $P_{nm}\equiv P_{nn,mm}=\langle n,m|\bar{\rho}|m,n\rangle$ with

$\displaystyle\sum^{n_{max}}_{n=0}\sum^{m_{max}}_{m=0}P_{nm}=1.$ (17)

The second-order photon-photon correlation function glauber  is defined as

$\displaystyle g^{(2)}_{a}(0)=\bigl{(}1/\langle
a^{\dagger}a\rangle^{2}\bigr{)}\sum^{n_{max}}_{n=0}\sum^{m_{max}}_{m=0}n(n-1)P_{nm},$
(18)

whereas the second-order phonon-phonon correlation function is calculated via
the expression:

$\displaystyle g^{(2)}_{b}(0)$ $\displaystyle=$
$\displaystyle\bigl{(}1/\langle
b^{\dagger}b\rangle^{2}\bigr{)}\sum^{n_{max}}_{n=0}\sum^{m_{max}}_{m=0}\bigl{(}m(m-1)+4\chi^{2}mn^{2}$
(19) $\displaystyle+$ $\displaystyle\chi^{4}n^{4}\bigr{)}P_{nm}.$

We note here that in a cavity optomechanical setup, second-order phonon-phonon
correlation functions are measured experimentally, as demonstrated e.g. in
Ref. exp4 .

(a)(b)

Figure 1:  (a) The steady-state optical cavity mean photon number $\langle
n\rangle\equiv\langle a^{\dagger}a\rangle$ as well as (b) the corresponding
mechanical resonator mean phonon number $\langle m\rangle\equiv\langle
b^{\dagger}b\rangle$ as a function of $\Delta/\omega$. Here $\chi=0.1$,
$\epsilon/\omega=0.02$, $\kappa_{a}/\omega=2\cdot 10^{-3}$ and
$\kappa_{b}/\omega=2\cdot 10^{-5}$, and only single-phonon processes were
included. Also, $\bar{n}=1$ for solid lines, while $\bar{n}=0$ for short-
dashed lines, respectively.

(a)(b)

Figure 2:  The steady-state second-order correlation function for (a) photons,
$g^{(2)}_{a}(0)$, and (b) for phonons, $g^{(2)}_{b}(0)$, as a function of
$\Delta/\omega$. The solid lines are for $\bar{n}=1$, while the short-dashed
one for $\bar{n}=0$, respectively. Other parameters are as in Fig. 1.

##  III Numerical Results and Discussion

In Fig. 1 we present the numerical results of the steady-state mean photon and
phonon numbers for a generic system with parameters within the considered
regime with $\chi=0.1$, $\epsilon/\omega=0.02$, $\kappa_{a}/\omega=2\cdot
10^{-3}$ and $\kappa_{b}/\omega=2\cdot 10^{-5}$. These parameters are
accessible experimentally. For instance, a similar set with $\omega\sim
1\rm{MHz}$, $g\sim 0.3\rm{MHz}$, $\kappa_{a}\sim 0.2\rm{MHz}$ and
$\kappa_{b}\sim 150\rm{Hz}$ was reported in Ref. exp1 . For our case, also a
good-cavity regime with $\omega\gg\kappa_{a}$ would be of interest, since a
too strong photon loss washes out the predicted features. Furthermore, we
considered here only one-phonon process, which corresponds to keeping terms up
to $\chi^{2}$ in Eq. (6), respectively. Few interesting features can be
observed: (i) In these parameter ranges, the steady-state phonon dynamics is
quite sensitive on external temperature variations, (ii) the photon dynamics
follows that of phonon one (see Appendix B), and vice versa, though with a
different magnitude, and (iii) the peak frequency intervals equal the Kerr-
like nonlinearity $\omega\chi^{2}$, see Hamiltonian (7). The latter result can
be intuitively understood from the free-photon Hamiltonian plus the Kerr-like
contribution, i.e.
$\bar{H}_{0}=\hbar\bar{a}^{\dagger}\bar{a}(\Delta-\omega\chi^{2}\bar{a}^{\dagger}\bar{a})$,
from where it follows that the effective detuning vanishes at
$\Delta_{k}=k\omega\chi^{2}$ or $|\Delta_{k}-\Delta_{k\pm
1}|=\omega\chi^{2}\equiv g^{2}/\omega$. This reveals that the induced Kerr
nonlinearity shifts the frequency of the corresponding photon Fock states,
which can be detected by scanning the frequency of the cavity photons. Thus,
the mean-phonon number dynamics can be extracted via detection of the mean-
photon number’s quantum dynamics, while the coupling constant $g$ can be
estimated if the mechanical motion frequency $\omega$ is known (or vice
versa).

(a)(b)

Figure 3:  The steady-state behaviors of (a) the mean-phonon number $\langle
m\rangle$ as well as (b) the second-order phonon-phonon correlation function
$g^{(2)}_{b}(0)$ as a function of $\chi=g/\omega$. The dashed blue line is
plotted for single-phonon whereas the green solid one for three-phonon
processes, respectively. Here $\Delta/\omega=0.09$ and $\bar{n}=1$ while other
parameters are as in Fig. (1).

Figure 2 shows the corresponding second-order correlation functions. The
photon statistics changes from sub-Poissonian ($g^{(2)}_{a}(0)<1$) for
$\Delta/\omega\lesssim 0.05$, to super-Poissonian ($g^{(2)}_{a}(0)\geq 2$) in
the interval 0.06-0.08 and to quasi-coherent features
($g^{(2)}_{a}(0)\gtrapprox 1$ ) for $\Delta/\omega>0.08$. On the other hand,
the phonon statistics exhibits quasi-thermal to super-Poissonian properties.
Also in this case, the curves are characterized by a multi-peak structure. For
smaller and negative $\Delta/\omega$, the lines in Figs. (1) and (2) remain
flat. Notice that in the absence of vibrations, i.e. when $g=0$, the cavity
mean-photon number and its correlation function are given by the following
expressions in the steady-state:

$\displaystyle\langle a^{\dagger}a\rangle$ $\displaystyle=$
$\displaystyle\epsilon^{2}/\bigl{(}\Delta^{2}+(\kappa_{a}/2)^{2}\bigr{)},~{}~{}{\rm
and}$ $\displaystyle\langle a^{\dagger 2}a^{2}\rangle$ $\displaystyle=$
$\displaystyle\epsilon^{4}/\bigl{(}\Delta^{2}+(\kappa_{a}/2)^{2}\bigr{)}^{2},$

leading to $g^{(2)}_{a}(0)\equiv\langle a^{\dagger 2}a^{2}\rangle/\langle
a^{\dagger}a\rangle^{2}=1$. For an independent boson mode in thermal
equilibrium with its surrounding thermal bath we obtain

$\displaystyle\langle b^{\dagger}b\rangle=\bar{n},~{}~{}{\rm and}~{}~{}\langle
b^{\dagger 2}b^{2}\rangle=2\bar{n}^{2},$

with $g^{(2)}_{b}(0)\equiv\langle b^{\dagger 2}b^{2}\rangle/\langle
b^{\dagger}b\rangle^{2}=2$.

In Fig. 3 we investigate the mean phonon number and its statistics as a
function of $\chi=g/\omega$ for a fixed ratio $\Delta/\omega=0.09$. As can be
seen in Fig. 1, for the particular coupling strength $g$ chosen ($\chi=0.1$)
the mean-quanta numbers are quite small at this detuning. Figure 3 compares
the case of single-phonon processes, i.e. keeping terms up to $\chi^{2}$ in
Eq. (6), with the case of three-phonon processes, for which we consider all
terms up to $\chi^{6}$ in Eq. (6). As a first observation, we notice that the
mean-phonon number increases with the ratio $g/\omega$. The peaks are slightly
higher for three-phonon processes, but the differences are barely noticeable
in Fig. 3a. Furthermore, the maxima occur at
$\chi_{k}$=$\sqrt{\Delta/(k\omega)}$, which for $\Delta/\omega=0.09$, give
k=6, 5, 4 and 3. Also, the steady-state mean phonon number is direct
proportional to the mean thermal phonon number $\bar{n}$ due to the
environmental reservoir, see Appendix B. The second-order phonon-phonon
correlation function changes accordingly and exhibits quasi-thermal features
for $\bar{n}=1$. In this case, single-phonon and three-phonon processes become
distinguished at higher ratios of $g/\omega$, see Fig. 3(b). Thus, multiphonon
processes occur for higher ratios of $g/\omega$. It can be shown that also in
this case the mean photon number dynamics is following the mean phonon number
dynamics, albeit with a smaller scaling. Finally, we have also calculated
higher-order phonon correlation functions multix ; multiy , i.e.
$g^{(k)}_{b}(0)=\langle b^{\dagger k}b^{k}\rangle/\langle
b^{\dagger}b\rangle^{k}$ as a function of $\Delta/\omega$ with
$k\in\\{3,4\\}$. Our numerical results show that these curves closely follow
the ones of $g^{(2)}_{b}(0)$, however with a larger magnitude. For instance,
the peak at $\Delta/\omega=0.07$ appears for the third-order (fourth order)
correlation function at the same position, but with a value of approx. 7.6
(39).

Our results confirm that the presence of a vibrating mirror changes
significantly both the photon’s as well as the phonon quantum dynamics of an
off-resonant pumped cavity optomechanical system. Moreover, the parameter
range required to observe this behavior is within experimental reach exp1 ;
exp2 ; exp3  and is close to those for photon blockade effect predicted in
optomechanical systems phblk1 ; phblk2 . Finally, the analytical approaches
proper to cavity optomecanics with a moving mirror equally apply to other
related samples, like e.g., hybrid metal-dielectric cavities hdcav , plasmon-
excitonic polaritons pep , superconducting qubits and quantum circuits qcirc
or other types of nanomechanical resonators nmr , respectively, rendering our
developed analytical approach relevant also for these systems.

##  IV Conclusions

We have investigated a cavity optomechanical setup where the detuning of the
external coherent electromagnetic field frequency from the optical resonator’s
one is different from the frequency (or its multiples) of one of the vibrating
mirror. As a result, the quantum dynamics of this complex system is
accompanied by multiple phonon absorption and emission processes. We have
computed the mean-quanta numbers and the corresponding second-order
correlation functions and described their properties. Particularly, we have
found that the inter-peak frequency intervals observed in the quantum dynamics
of the photon and phonon subsystem as a function of detuning equal the
vibration-induced Kerr-like non-linearity. The photon mean-number dynamics
follows that of the mean-phonon one which is convenient for monitoring the
phonon quantum dynamics by photon detection. The corresponding second-order
correlation functions also exhibit a multi-peak structure and are completely
different from the fixed-mirror case. The photon-photon correlation function
may exhibit sub-Poissonian to super-Poissonian photon statistics. The
corresponding correlation function for phonons lies within quasi-thermal to
thermal phonon statistics, characterized by super-Poissonian features.
Furthermore, the second order phonon correlation function can be used to
distinguish among the one- and few-phonon processes for stronger coupling
strengths among the involved interacting subsystems. For a more detailed
investigation of non-classical signatures, future work could focus on the
Wigner function for representative parameter sets, or on so-called bundle
correlation functions, $g^{(2)}_{2}(\tau)=\langle b^{\dagger 2}(0)b^{\dagger
2}(\tau)b^{2}(\tau)b^{2}(0)\rangle/\langle(b^{\dagger
2}b^{2})(0)\rangle\langle(b^{\dagger 2}b^{2})(\tau)\rangle$ multix ; multiy ,
which could reveal highly correlated phonon behaviour.

###### Acknowledgements.

MM acknowledges the financial support by the Moldavian National Agency for
Research and Development, grant No. 20.80009.5007.07. AP gratefully
acknowledges support from the Heisenberg Program of the Deutsche
Forschungsgemeinschaft (DFG).

##  Appendix A Expanding the exponents up to $\chi^{4}$

In order to demonstrate the multi-phonon nature of the combined quantum
dynamics, let us expand the exponents in the expression (13), i.e.
$e^{\chi(\bar{b}-\bar{b}^{\dagger})}\bar{a}\bar{\rho}\bar{a}^{\dagger}e^{-\chi(\bar{b}-\bar{b}^{\dagger})}$,
up to $\chi^{4}$ and keep the time-independent terms only, that is,

$\displaystyle
e^{\chi(\bar{b}-\bar{b}^{\dagger})}\bar{a}\bar{\rho}\bar{a}^{\dagger}e^{-\chi(\bar{b}-\bar{b}^{\dagger})}=\bar{a}\bar{\rho}\bar{a}^{\dagger}+\biggl{(}\bar{a}\bar{\rho}\bar{a}^{\dagger}\bigl{\\{}\frac{\chi^{2}}{2!}(1-2\bar{b}\bar{b}^{\dagger})+\frac{\chi^{4}}{4!}(6\bar{b}^{2}\bar{b}^{\dagger
2}-12\bar{b}\bar{b}^{\dagger}+3)\bigr{\\}}+H.c.\biggr{)}+\chi^{2}(1+\chi^{2})\biggl{(}\bar{b}\bar{a}\bar{\rho}\bar{a}^{\dagger}\bar{b}^{\dagger}$
$\displaystyle+\bar{b}^{\dagger}\bar{a}\bar{\rho}\bar{a}^{\dagger}\bar{b}\biggr{)}-\frac{\chi^{4}}{2}\biggl{(}\bar{b}\bar{a}\bar{\rho}\bar{a}^{\dagger}\bar{b}\bar{b}^{\dagger
2}+\bar{b}^{\dagger}\bar{a}\bar{\rho}\bar{a}^{\dagger}\bar{b}^{2}\bar{b}^{\dagger}+H.c.\biggr{)}+\frac{\chi^{4}}{4}\biggl{(}\bar{b}^{2}\bar{a}\bar{\rho}\bar{a}^{\dagger}\bar{b}^{\dagger
2}+(1-2\bar{b}\bar{b}^{\dagger})\bar{a}\bar{\rho}\bar{a}^{\dagger}(1-2\bar{b}\bar{b}^{\dagger})+\bar{b}^{\dagger
2}\bar{a}\bar{\rho}\bar{a}^{\dagger}\bar{b}^{2}\biggr{)}.$ (20)

When considering the photon-phonon distribution function $P_{n_{1}n_{2},mm}$,
then from (20) one can observe that terms proportional at least to $\chi^{2}$
contribute to single-phonon processes, whereas those proportional to
$\chi^{4}$ \- to two-phonon processes, respectively. For instance, the last
term from the first line of expression (20) accounts for single-phonon
processes where the phonon number changes by $\pm 1$, i.e. $m\pm 1$. On the
other side, the last term from (20) describes two-phonon processes where the
phonon number modifies by $\pm 2$, $m\pm 2$. Thus, the expression (20),
expanded up to $\chi^{4}$ terms, describes simultaneously single- and two-
phonon processes, respectively. Also, single-phonon processes are influenced
by the second-order ones, i.e. there is a contribution to single-phonon
processes coming from $\chi^{4}$ terms.

Generalizing now, one can state that $N$-phonon processes are described by
terms proportional to $\chi^{2N}$. There is a contribution to $N$-phonon
processes which arises from higher order terms proportional to
$\chi^{2(N+1)}$. This way, depending on a certain power of $\chi$, one can
investigate the multiphonon quantum dynamics of the cavity optomechanical
system.

##  Appendix B Relationship among the optical and phonon modes

In what follows, we shall give details on how the mean cavity photon and mean
phonon numbers are interconnected in the dispersive interaction regime
investigated here. The master equation (6), containing terms up to $\chi^{2}$
for simplicity and in the considered approximations, reads

$\displaystyle\frac{d}{dt}\bar{\rho}(t)$ $\displaystyle+$
$\displaystyle\frac{i}{\hbar}[\bar{H},\bar{\rho}]=$ $\displaystyle-$
$\displaystyle\frac{\kappa_{a}}{2}\biggl{(}\bar{a}^{\dagger}\bar{a}\bar{\rho}-e^{\chi(\bar{b}-\bar{b}^{\dagger})}\bar{a}\bar{\rho}\bar{a}^{\dagger}e^{-\chi(\bar{b}-\bar{b}^{\dagger})}\biggr{)}$
$\displaystyle-$
$\displaystyle\frac{\kappa_{b}}{2}(1+\bar{n})\biggl{(}[\bar{b}^{\dagger},\bar{b}\bar{\rho}]+\chi^{2}[\bar{a}^{\dagger}\bar{a},\bar{a}^{\dagger}\bar{a}\bar{\rho}]\biggr{)}$
$\displaystyle-$
$\displaystyle\frac{\kappa_{b}}{2}\bar{n}\biggl{(}[\bar{b},\bar{b}^{\dagger}\bar{\rho}]+\chi^{2}[\bar{a}^{\dagger}\bar{a},\bar{a}^{\dagger}\bar{a}\bar{\rho}]\biggr{)}+H.c.,$

where

$\displaystyle\bar{H}$ $\displaystyle=$
$\displaystyle\hbar\Delta\bar{a}^{\dagger}\bar{a}+\hbar\omega\bar{b}^{\dagger}\bar{b}+\hbar\epsilon(\bar{a}+\bar{a}^{\dagger})\bigl{\\{}1-\chi^{2}(\bar{b}^{\dagger}\bar{b}+\bar{b}\bar{b}^{\dagger})/2\bigr{\\}}$
(22) $\displaystyle-$
$\displaystyle\omega\chi^{2}(\bar{a}^{\dagger}\bar{a})^{2}.$

The exponent term from the first line of the damping part of the above master
equation can be obtained from (20) while setting $\chi^{4}\to 0$. Then, in the
steady-state, from the master equation (LABEL:app2) one can easily show that

$\displaystyle\langle\bar{b}^{\dagger}\bar{b}\rangle=\bar{n}+\chi^{2}\frac{\kappa_{a}}{\kappa_{b}}\langle\bar{a}^{\dagger}\bar{a}\rangle,$
(23)

while from Eqs. (4), one finally arrives at the relationship among the steady-
state mean phonon and photon numbers, namely,

$\displaystyle\langle b^{\dagger}b\rangle$ $\displaystyle=$
$\displaystyle\langle\bar{b}^{\dagger}\bar{b}\rangle+\chi^{2}\langle\bar{a}^{\dagger
2}\bar{a}^{2}\rangle$ (24) $\displaystyle=$
$\displaystyle\bar{n}+\chi^{2}\frac{\kappa_{a}}{\kappa_{b}}\langle\bar{a}^{\dagger}\bar{a}\rangle+\chi^{2}\langle\bar{a}^{\dagger
2}\bar{a}^{2}\rangle.$

One can observe here that there is an almost linear dependence between mean
phonon and cavity photon numbers, respectively, since the last term may give
little contribution as long as $\chi\ll 1$. In addition, we have assumed that
$\kappa_{a}/\kappa_{b}\gg 1$ . Thus, detecting the cavity photons one can
estimate the mean phonon number and vice versa.

## References

  * (1) D. F. Walls, and G. J. Milburn, Quantum Optics (Springer, Berlin, Heidelberg, 2008).
  * (2) M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Cavity optomechanics, Rev. Mod. Phys. 86, 1391 (2014).
  * (3) C. L. Degen, F. Reinhard, and P. Cappellaro, Quantum sensing, Rev. Mod. Phys. 89, 035002 (2017).
  * (4) V. B. Braginskii, Classical and quantum restrictions on the detection of weak disturbances of a macroscopic oscillator, Sov. Phys. JETP 26, 831 (1968).
  * (5) I. Pikovski, M. R. Vanner, M. Aspelmeyer, M. S. Kim, and C. Brukner, Probing Planck-scale physics with quantum optics, Nature Physics 8, 393 (2012).
  * (6) G. A. Brawley, M. R. Vanner, P. E. Larsen, S. Schmid, A. Boisen, and W.P. Bowen, Nonlinear optomechanical measurement of mechanical motion, Nature Communications 7:10988, 1 (2016).
  * (7) S. Qvarfort, A. Serafini, P. F. Barker, and S. Bose, Gravimetry through non-linear optomechanics, Nature Communications 9:3690, 1 (2018).
  * (8) S. Gigan, R. Böhm, M. Paternostro, F. Blaser, G. Langer, J. B. Hertzberg, K. C. Schwab, D. Bäuerle, M. Aspelmeyer, and A. Zeilinger, Self-cooling of a micromirror by radiation pressure, Nature (London) 444, 67 (2006).
  * (9) O. Arcizet, P. F. Cohadon, T. Briant, M. Pinard, and A. Heidmann, Radiation-pressure cooling and optomechanical instability of a micromirror, Nature (London) 444, 71 (2006).
  * (10) I. Wilson-Rae, N. Nooshi, W. Zwerger, and T. J. Kippenberg, Theory of Ground State Cooling of a Mechanical Oscillator Using Dynamical Backaction, Phys. Rev. Lett. 99, 093901 (2007).
  * (11) F. Marquardt, J. P. Chen, A. A. Clerk, and S. M. Girvin, Quantum Theory of Cavity-Assisted Sideband Cooling of Mechanical Motion, Phys. Rev. Lett. 99, 093902 (2007).
  * (12) J. D. Teufel, T. Donner, D. Li, J. W. Harlow, M. S. Allman, K. Cicak, A. J. Sirois, J. D. Whittaker, K. W. Lehnert, and R. W. Simmonds, Sideband cooling of micromechanical motion to the quantum ground state, Nature (London) 475, 359 (2011).
  * (13) J. Chan, T. P. Mayer Alegre, A. H. Safavi-Naeini, J. T. Hill, A. Krause, S. Gröblacher, M. Aspelmeyer, and O. Painter, Laser cooling of a nanomechanical oscillator into its quantum ground state, Nature (London) 478, 89 (2011).
  * (14) S. Mancini, V. I. Man’ko, and P. Tombesi, Ponderomotive control of quantum macroscopic coherence, Phys. Rev. A 55, 3042 (1997).
  * (15) S. Bose, K. Jacobs, and P. L. Knight, Preparation of nonclassical states in cavities with a moving mirror, Phys. Rev. A 56, 4175 (1997).
  * (16) G. S. Agarwal, and S. Huang, Electromagnetically induced transparency in mechanical effects of light, Phys. Rev. A 81, 041803(R) (2010).
  * (17) S. Weis, R. Riviere, S. Deleglise, E. Gavartin, O. Arcizet, A. Schliesser, and T. J. Kippenberg, Optomechanically Induced Transparency, Science 330, 1520 (2010).
  * (18) H. Xionga, and Y. Wu, Fundamentals and applications of optomechanically induced transparency, Applied Phys. Rev. 5, 031305 (2018).
  * (19) G. S. Agarwal, and S. Huang, Optomechanical systems as single-photon routers, Phys. Rev. A 85, 021801(R) (2012).
  * (20) W.-T. Liao, and A. Pálffy, Optomechanically induced transparency of x-rays via optical control, Scientific Reports 7: 321, 1 (2017).
  * (21) S. Aldana, Ch. Bruder, and A. Nunnenkamp, Equivalence between an optomechanical system and a Kerr medium, Phys. Rev. A 88, 043826 (2013).
  * (22) P. D. Drummond, and D. F. Walls, Quantum theory of optical bistability. I. Nonlinear polarisability model, J. Phys. A 13, 725 (1980).
  * (23) F. Marquardt, J. G. E. Harris, and S. M. Girvin, Dynamical Multistability Induced by Radiation Pressure in High-Finesse Micromechanical Optical Cavities, Phys. Rev. Lett. 96, 103901 (2006).
  * (24) M. A. Macovei, Measuring photon-photon interactions via photon detection, Phys. Rev. A 82, 063815 (2010).
  * (25) A. Miranowicz, M. Paprzycka, Y.-x. Liu, J. Bajer, and F. Nori, Two-photon and three-photon blockades in driven nonlinear systems, Phys. Rev. A 87, 023809 (2013).
  * (26) D. Vitali, S. Gigan, A. Ferreira, H. R. Böhm, P. Tombesi, A. Guerreiro, V. Vedral, A. Zeilinger, and M. Aspelmeyer, Optomechanical Entanglement between a Movable Mirror and a Cavity Field, Phys. Rev. Lett. 98, 030405 (2007).
  * (27) L. Tian, Robust Photon Entanglement via Quantum Interference in Optomechanical Interfaces, Phys. Rev. Lett. 110, 233602 (2013).
  * (28) Z. R. Gong, H. Ian, Y.-x. Liu, C. P. Sun, and F. Nori, Effective Hamiltonian approach to the Kerr nonlinearity in an optomechanical system, Phys. Rev. A 80, 065801 (2009).
  * (29) P. Rabl, Photon Blockade Effect in Optomechanical Systems, Phys. Rev. Lett. 107, 063601 (2011).
  * (30) F. Zou, L.-B. Fan, J.-F. Huang, and J.-Q. Liao, Enhancement of few-photon optomechanical effects with cross-Kerr nonlinearity, Phys. Rev. A 99, 043837 (2019).
  * (31) K. Stannigel, P. Rabl, A. S. Sorensen, P. Zoller, and M. D. Lukin, Optomechanical Transducers for Long-Distance Quantum Communication, Phys. Rev. Lett. 105, 220501 (2010).
  * (32) D. E. Chang, V. Vuletic, and M. D. Lukin, Quantum nonlinear optics — photon by photon, Nature Photonics 8, 685 (2014).
  * (33) M. A. Nielsen, and I. L. Chuang, Quantum Computation and Quantum Information (Cambridge University Press, Cambridge, 2000).
  * (34) A. Nunnenkamp, K. Borkje, and S. M. Girvin, Single-Photon Optomechanics, Phys. Rev. Lett. 107, 063602 (2011).
  * (35) A. Kronwald, M. Ludwig, and F. Marquardt, Full photon statistics of a light beam transmitted through an optomechanical system, Phys. Rev. A 87, 013847 (2013).
  * (36) S. Gröblacher, K. Hammerer, M. R. Vanner, and M. Aspelmeyer, Observation of strong coupling between a micromechanical resonator and an optical cavity field, Nature 460, 724 (2009).
  * (37) R. Leijssen, and E. Verhagen, Strong optomechanical interactions in a sliced photonic crystal nanobeam, Scientific Reports 5:15974, 1 (2015).
  * (38) G. A. Peterson, S. Kotler, F. Lecocq, K. Cicak, X. Y. Jin, R.W. Simmonds, J. Aumentado, and J. D. Teufel, Ultrastrong Parametric Coupling between a Superconducting Cavity and a Mechanical Resonator, Phys. Rev. Lett. 123, 247701 (2019).
  * (39) M. Kiffner, M. Macovei, J. Evers, and C. H. Keitel, Vacuum-induced processes in multilevel atoms, Progress in Optics 55, 85 (2010).
  * (40) J. Qian, A. A. Clerk, K. Hammerer, and F. Marquardt, Quantum Signatures of the Optomechanical Instability, Phys. Rev. Lett. 109, 253601 (2012).
  * (41) R. J. Glauber, The Quantum Theory of Optical Coherence, Phys. Rev. 130, 2529 (1963).
  * (42) J. D. Cohen, S. M. Meenehan, G. S. MacCabe, S. Gröblacher, A. H. Safavi-Naeini, F. Marsili, M. D. Shaw, and O. Painter, Phonon counting and intensity interferometry of a nanomechanical resonator, Nature 520, 522 (2015).
  * (43) C. Sánchez Munoz, E. del Valle, A. González Tudela, K. Müller, S. Lichtmannecker, M. Kaniber, C. Tejedor, J. J. Finley, and F. P. Laussy, Emitters of N-photon bundles, Nat. Photonics 8, 550 (2014).
  * (44) Q. Bin, X.-Y. Lü, F. P. Laussy, F. Nori, and Y. Wu, N-Phonon Bundle Emission via the Stokes Process, Phys. Rev. Lett. 124, 053601 (2020).
  * (45) M. K. Dezfouli, R. Gordon, and S. Hughes, Molecular Optomechanics in the Anharmonic Cavity-QED Regime Using Hybrid Metal-Dielectric Cavity Modes, ACS Photonics 6, 1400 (2019).
  * (46) T. Neuman, J. Aizpurua, and R. Esteban, Quantum theory of surface-enhanced resonant Raman scattering (SERRS) of molecules in strongly coupled plasmon-exciton systems, Nanophotonics 9(2), 295 (2020).
  * (47) S. N. Shevchenko, and A. N. Omelyanchouk, Multiphoton transitions in Josephson-junction qubits, Low Temperature Physics 38, 283 (2012).
  * (48) Y. Greenberg, Y. Pashkin, and E. Il’ichev, Nanomechanical resonators, Physics-Uspekhi, 182, 407 (2012).


# From Feast to Famine: A Systematic Study of Accretion onto Oblique Pulsars
with 3D GRMHD Simulations

Ariadna Murguia-Berthier NASA Hubble Fellow  Center for Interdisciplinary
Exploration & Research in Astrophysics (CIERA), Physics & Astronomy,
Northwestern University, Evanston, IL 60202, USA  Kyle Parfrey Princeton
Plasma Physics Laboratory, Princeton, NJ 08540, USA  Alexander Tchekhovskoy
Center for Interdisciplinary Exploration & Research in Astrophysics (CIERA),
Physics & Astronomy, Northwestern University, Evanston, IL 60202, USA  Jonatan
Jacquemin-Ide Center for Interdisciplinary Exploration & Research in
Astrophysics (CIERA), Physics & Astronomy, Northwestern University, Evanston,
IL 60202, USA

###### Abstract

Disk-fed accretion onto neutron stars can power a wide range of astrophysical
sources ranging from X-ray binaries, to accretion powered millisecond pulsars,
ultra-luminous X-ray sources, and gamma-ray bursts. A crucial parameter
controlling the gas–magnetosphere interaction is the strength of the stellar
dipole. In addition, coherent X-ray pulsations in many neutron star systems
indicate that the star’s dipole moment is oblique relative to its rotation
axis. Therefore, it is critical to systematically explore the 2D parameter
space of the star’s magnetic field strength and obliquity, which is what this
work does, for the first time, in the framework of 3D general-relativistic
magnetohydrodynamics. If the accretion disk carries its own vertical magnetic
field, this introduces an additional factor: the relative polarity of the disk
and stellar magnetic fields. We find that depending on the strength of the
stellar dipole and the star–disk relative polarity, the neutron star’s jet
power can either increase or decrease with increasing obliquity. For weak
dipole strength (equivalently, high accretion rate), the parallel polarity
results in a positive correlation between jet power and obliquity, whereas the
anti-parallel orientation displays the opposite trend. For stronger dipoles,
the relative-polarity effect disappears, and jet power always decreases with
increasing obliquity. The influence of the relative polarity gradually
disappears as obliquity increases. Highly oblique pulsars tend to have an
increased magnetospheric radius, a lower mass accretion rate, and enter the
propeller regime at lower magnetic moments than aligned stars.

accretion — pulsars — neutron stars — magnetohydrodynamics — relativistic jets
— general relativity

††software: harmpi (Tchekhovskoy, 2019), matplotlib (Hunter, 2007), numpy
(Harris et al., 2020), scipy (Virtanen et al., 2020), hdf5 (The HDF Group,
1997-NNNN), yt (Turk et al., 2011)

##  1 Introduction

Accretion onto neutron stars (NSs) can power a wide range of astrophysical
sources ranging from X-ray binaries (Walter et al., 2015), to accretion
powered millisecond pulsars (Archibald et al., 2009), ultra-luminous X-ray
sources (Bachetti et al., 2014; Israel et al., 2017), magnetar-powered
superluminous supernovae (Woosley, 2010) and gamma-ray bursts (Duncan &
Thompson, 1992; Metzger et al., 2011). The fact that NSs have a surface and
show periodic variability provides a unique opportunity to turn them into
astrophysical laboratories of accretion by directly measuring their spins and
spin derivatives: this enables direct measurements of the accretion flow and
magnetospheric torques experienced by the star, a luxury not afforded to us by
black holes (BHs) where spin measurements are scarce and often come with
significant uncertainties (Remillard & McClintock, 2006; Steiner et al., 2009;
Reynolds, 2014). However, these data-rich NS observations are _interpretation-
limited_ due to the highly nonlinear gas-magnetosphere interactions that
confound analytic models and the scarcity of first-principles 3D models that
describe the interaction between the NS magnetosphere and infalling gas.

Simulating the gas-magnetosphere interactions in the NS context is indeed a
challenging problem, because it involves enormous density and magnetization
contrasts: light relativistically magnetized magnetospheric plasma interacting
with heavy weakly magnetized infalling gas. The two regions require seemingly
mutually exclusive approaches: a force-free approach capable of handling the
high magnetization of the magnetosphere (but incapable of handling the gas)
and a general-relativistic magnetohydrodynamical (GRMHD) approach capable of
handling moderately magnetized fluid (but incapable of handling the extreme
magnetization of the magnetosphere).

Several groups have performed studies of gas-magnetosphere interactions in the
NS context (Parfrey et al., 2017; Parfrey & Tchekhovskoy, 2017; Das et al.,
2022; Parfrey & Tchekhovskoy, 2023; Das & Porth, 2023). Recently, Parfrey &
Tchekhovskoy (2023) self-consistently simulated in 3D a neutron star
surrounded by an accretion disk using a hybrid GRMHD approach, which combines
the advantages of the force-free and GRMHD techniques. They studied how the
relative orientation of the large-scale NS and disk magnetic field vectors
affected magnetized NS accretion, as a function of NS dipole strength. They
found that similar to 2D simulations (Parfrey & Tchekhovskoy, 2017), the
reconnection of disk and magnetospheric magnetic fields can aid in opening
closed stellar magnetic field lines and result in an increased NS spindown
power, launching a relativistic jet. However, in contrast to 2D, even when the
disk and magnetospheric fields are parallel to each other and cannot
reconnect, the field line opening can still happen due to the _interchange
slingshot_ , a version of the magnetic interchange instability. They also
found that this mechanism only activates when the NS magnetic fields are
strong enough to push the accretion flow away (i.e., in the propeller regime),
or close to doing so. Therefore, both the field strength and its orientation,
or the tilt of the NS magnetosphere, crucially affect the disk-magnetosphere
interactions.

Indeed, the magnetic and rotational axes of accreting NSs are typically tilted
relative to one other, with tilt angles sometimes reaching as high as
$60^{\circ}$ (e.g., XSS J12270-4859, de Martino et al., 2014). How does such a
large tilt affect the system? Does it decrease the magnetospheric radius at
constant mass accretion rates (Wang, 1997; Bozzo et al., 2018)? Does it
increase the pulsar’s spindown power similar to isolated pulsars (Spitkovsky,
2006; Tchekhovskoy et al., 2013, 2016)? If so, can this make it harder for gas
to accrete and, possibly, even bring about the propeller state, in which the
disk–magnetosphere interaction expels the accretion flow?

Due to the lack of symmetries and the non-linearity of gas–magnetosphere
interactions, direct multidimensional simulations provide an attractive line
of attack on this critical problem. Romanova et al. (2003, 2004, 2013, 2021)
explored the interaction between a non-relativistic oblique dipole and an
accretion flow via an $\alpha-$prescription, and Romanova et al. (2012) via
the magnetorotational instability (MRI, Balbus & Hawley 1991; Hawley & Balbus
1991; Balbus & Hawley 1998). They found that obliquity dramatically affects
the geometry of the incoming streams, mass accretion rate, neutron star
spindown power, and the radiative signatures from the hotspots produced by the
infalling streams on the stellar surface.

In this paper, we perform the first systematic exploration of magnetized
accretion onto a rapidly rotating oblique NS dipolar magnetosphere, for a wide
range of obliquities and stellar dipolar strengths, using the 3D hybrid-GRMHD
approach. This enables us to self-consistently handle relativistically
magnetized magnetospheric plasma alongside dense weakly magnetized non-
relativistic gas, and to resolve the MRI. In Sec. 2, we discuss the problem
setup and numerical method, in Sec. 3 we investigate the effects of the
obliquity angle (Sec. 3.1) and stellar magnetic dipole moment (Sec. 3.2), and
finish with the special case of a perpendicular rotator (Sec. 3.3). We
conclude in Sec. 4.

##  2 Numerical Method and Problem Setup

We solve the equations of GRMHD using a modified version of the harmpi code
(Tchekhovskoy, 2019), which is based on the harm code (Gammie et al., 2003;
Noble et al., 2006):

$\displaystyle\nabla_{\mu}(\rho u^{\mu})$ $\displaystyle=0,$ (1)
$\displaystyle\nabla_{\mu}T^{\mu}_{\nu}$ $\displaystyle=0,$ (2)
$\displaystyle\nabla_{\mu}{{}^{{}^{*}}\\!\\!F}^{\mu\nu}$ $\displaystyle=0,$
(3)

where $\rho$ is the fluid-frame mass density, $u^{\mu}$ is the 4-velocity, and
${{}^{{}^{*}}\\!\\!F}^{\mu\nu}=b^{\mu}u^{\nu}-b^{\nu}u^{\mu}$ is the dual of
the electromagnetic tensor. We set $G=M_{*}=c=1$ and work in Lorentz-Heaviside
units such that the fluid-frame magnetic pressure is $b^{2}/2$, where
$b^{\mu}={{}^{{}^{*}}\\!\\!F}^{\nu\mu}u_{\nu}$ is the fluid-frame magnetic
field 4-vector. The stress energy tensor takes the form:

$T^{\mu\nu}=\rho
hu^{\mu}u^{\nu}+Pg^{\mu\nu}+{{b^{2}}}u^{\mu}u^{\nu}+\frac{1}{2}{{b^{2}}}g^{\mu\nu}-b^{\mu}b^{\nu}\;,$
(4)

where $h=\left(1+\epsilon+P/\rho\right)$ is the specific enthalpy, $P$ is the
gas pressure, $\epsilon$ is the specific internal energy density, and
$g_{\mu\nu}$ is the metric tensor.

Figure 1: Density slices through the $\mu-\Omega$ (top row) and equatorial
(bottom row) planes of our simulations with a magnetic moment of $\mu=20$
reveal that the increase in magnetic obliquity brings the system closer to the
propeller state. We show the system at $t=14{,}200\,r_{\rm g}/c\approx
113\,P_{*}$, for both parallel and anti-parallel configurations. [Left 2
panels]: At low obliquity, $\chi=10^{\circ}$, the gas reaches the NS. [Right 2
panels]: However, at high obliquity, $\chi=60^{\circ}$, the magnetosphere
pushes most of the gas outside of the corotation radius, shown with the dashed
lines. As a result, the system enters the propeller regime, and the gas
struggles to reach the NS. Directed pink lines trace out the magnetic field in
the image plane. The solid cyan lines indicate the light cylinder. Figure 2:
Different quantities as a function of time for $\mu=20$. The dashed lines
correspond to simulations in which the magnetic field of the torus and the NS
are anti-parallel, and the solid line represents simulations in which they are
parallel. The magenta line represents the perpendicular simulation. From top
to bottom, the quantities are: (a) mass accretion rate at the NS surface; (b)
open magnetic flux at $r=R_{\rm LC}$ normalized to an isolated pulsar; (c) jet
power at $r=R_{\rm LC}$, similarly normalized; (d) total torque, similarly
normalized; (e) magnetospheric radius. See section 3 for a detailed
explanation of the quantities and how they were calculated. The horizontal
lines represent $\Psi_{\rm open}=\Psi_{\rm open,0}$, $L=L_{0}$, $N_{\rm
tot}=-|N_{0}|$, $r_{\rm mag}=r_{\rm co}$. Figure 3: Density snapshots of the
$\mu-\Omega$ plane (top panel) and the equatorial plane (bottom panel) for
simulations with $\chi=60^{\circ}$, in an anti-parallel configuration. We show
different magnetic moments at time $t=14{,}200\,r_{\rm g}/c$, or around 113 NS
spin periods. The pink lines represent the in-plane magnetic field component.
The solid cyan line is placed at the light cylinder, and the dashed cyan line
is the corotation radius. Figure 4: Density snapshots of a simulation with a
magnetic moment of $\mu=10$ at different times, where the star’s magnetic
moment is perpendicular to its rotation axis. Top row panels show different
planes, from left to right: (a) the plane perpendicular to $\mu-\Omega$, (b)
the $\mu-\Omega$ plane, and (c) a plane between those two. The bottom panel is
the equatorial plane. The pink lines represent the magnetic field of the in-
plane component. The solid cyan line is placed at the light cylinder, and the
dashed cyan line is the corotation radius.

We place the NS stellar surface at $r_{*}=4\,r_{\rm g}$, where $r_{\rm
g}=GM_{*}/c^{2}$ is the gravitational radius of the star. We adopt the angular
velocity of the NS to be $\Omega_{*}=0.05\,c/r_{\rm g}$, which places the
light cylinder at $R_{\rm LC}\equiv c/\Omega_{*}=20\,r_{\rm g}$ and sets the
NS rotational period, $P_{*}=2\pi/\Omega_{*}$. We use the Kerr metric, with a
dimensionless spin of $a=1/3$.

The Keplerian angular velocity, $\Omega_{\rm K}(r)=1/(r\sqrt{r}+a)$,
determines the corotation radius, $r_{\rm co}\approx 7.3\,r_{\rm g}$, at which
where the NS spin and Keplerian angular velocities are the same. We add a
dipolar magnetic field by setting the magnetic vector potential (Wasserman &
Shapiro, 1983):

$\displaystyle A_{\phi,\rm WS}(r,\theta,\phi)$
$\displaystyle=\Phi\sin{\theta}(\cos{\chi}\sin{\theta}-\sin{\chi}\cos{\theta}\cos{\phi}),$
(5) $\displaystyle A_{\theta,\rm WS}(r,\theta,\phi)$
$\displaystyle=-\Phi\sin{\phi}\sin{\chi},$ (6)

where

$\Phi=\frac{3\mu\sin^{2}\theta}{2M_{*}}\left[\lambda^{2}\ln(1-\lambda^{-1})+\lambda+\frac{1}{2}\right],$
(7)

$\lambda=r/2M_{*}$, and $\mu$ is the magnetic moment of the star in units of
$4\pi r_{\rm g}^{3}\sqrt{\rho_{\rm max}c^{2}}$.

For the rotating magnetized NS boundary conditions, we set the surface
electric fields in our constrained transport scheme such that at each time
step they result in the surface radial magnetic field, which is a
discretization of the initial dipolar magnetic flux rotated by an angle
$\Delta\phi=\Omega_{*}t$, where $t$ is the current simulation time
(Tchekhovskoy et al., 2013; Philippov et al., 2014).

We include an equilibrium hydrodynamic torus (Fishbone & Moncrief, 1976;
Chakrabarti, 1985), with the inner radius at $r_{\rm in}=40r_{\rm g}$ and the
radius of maximum pressure at $r_{\rm max}=60r_{\rm g}$. Inside the torus, we
add a single poloidal (i.e., in the $R$\- and $z$-directions) magnetic field
loop that follows the lines of constant density, described by the magnetic
vector potential, $A_{\phi,\rm disk}\propto\rho^{2}r^{5}$. We normalize the
disk magnetic field so that the ratio of the gas to magnetic pressure maxima
is $\max P/\max(b^{2}/2)=100$. Note that we have the freedom to orient the
disk magnetic field loop to be either parallel or anti-parallel to the dipolar
magnetic field of the NS, and in this work we will consider both magnetic
field orientations. We note that, unlike Parfrey & Tchekhovskoy (2023), in
these simulations we do not deform the magnetic field around the torus. This
results in accretion occurring earlier and a different profile of vertical
magnetic flux being dragged inward by the disk over time.

We adopt the approach of handling both MHD and force-free regions in a single
simulation in a self-consistent way (Parfrey & Tchekhovskoy, 2017, 2023).
Namely, we introduce a passive scalar, $\mathcal{F}$, that we evolve as

$\nabla_{\mu}(\mathcal{F}\rho u^{\mu})=0.$ (8)

The scalar enables the algorithm to distinguish the magnetospheric and gaseous
regions inside the light sphere $r\leq R_{\rm LC}$: $\mathcal{F}\approx 1$ in
the magnetospheric, or force-free, regions, and $\mathcal{F}\approx 0$ where
the accreting gas dominates. When $\mathcal{F}=0$ the GRMHD equations are used
without modification, while when $\mathcal{F}=1$ the non-force-free degrees of
freedom are driven towards force-free-like values over a characteristic
timescale; intermediate $\mathcal{F}$ values smoothly interpolate between
these extremes. The GRMHD evolution is unmodified in all cells for $r>R_{\rm
LC}$.

We performed a suite of simulations covering a wide range of dipolar magnetic
field strengths, $\mu$, and obliquity angles, $\chi$, where $\chi=0^{\circ}$
denotes an aligned rotator. Namely, we performed simulations for
$\mu=5,10,20,100$, each for several obliquity angles,
$\chi=0^{\circ},10^{\circ},60^{\circ},90^{\circ}$, and for both the anti-
parallel and parallel configurations of torus magnetic field.111Note that the
two torus magnetic field configurations are equivalent for $\chi=90^{\circ}$.
We also simulated the intermediate, $\chi=30^{\circ},45^{\circ}$, cases for
$\mu=20$ and anti-parallel torus magnetic field.

For comparison, we repeated all of the above models for an isolated pulsar
(i.e., without the initial torus). Our simulations have a resolution of
$N_{r}=192$, $N_{\theta}=128$, $N_{\phi}=64$. We evolve all simulations until
$t_{\rm F}=20{,}000\,r_{\rm g}/c\approx 150P_{*}$.

##  3 Results

###  3.1 Effect of obliquity angle

Figure 1 shows the vertical, $\mu-\Omega$, and equatorial, $x-y$, slices
through the density in our $\mu=20$ simulations for low, $\chi=10^{\circ}$,
and high, $\chi=60^{\circ}$, obliquities, at $t=14200r_{\rm g}/c\approx
113P_{*}$. Directed pink lines show the image plane magnetic field.

At the beginning of the simulation, magnetic stresses are induced in the torus
by its differential rotation, triggering the MRI. As angular momentum is
transported outwards, material from the torus is transported inwards, onto the
NS. When the obliquity angle is low, the pressure due to the star’s closed
magnetic field lines halts the accretion flow and redirects some of the
material to the poles. The radius at which this happens is known as the
magnetospheric radius. If the magnetic field the torus presents to the NS is
anti-parallel to the NS closed dipolar field, the two fields can reconnect:
this will open NS field lines that would be closed in an otherwise isolated
pulsar. The additional open field lines add to the pulsar wind, and create a
Poynting-flux-dominated relativistic jet. This can be seen on the panels for
$\chi=10^{\circ}$ (low obliquity), for the anti-parallel configuration, in the
$\mu-\Omega$ plane. We see ordered magnetic field lines forming what in 3D is
a helical structure. We can also see how material from the torus is accreted
onto the NS in a “ring”-like manner. On the other hand, if the magnetic field
lines of the torus and the dipolar field are parallel, reconnection is not
possible, but, for high enough magnetic moment (like $\mu=20$), magnetic field
lines can open due to the interchange slingshot instability (Parfrey &
Tchekhovskoy, 2023).

If we increase the obliquity angle, instead of a “ring”, material will be
accreted onto the NS via a single, thin stream onto each magnetic pole. This
is illustrated in the $\mu-\Omega$ panels of the $\chi=60^{\circ}$ simulation
in Figure 1. Although it is not seen in the shown snapshots, at several times,
and at several different planes, the thin stream will fall onto the NS and be
accreted. Similar results were found in non-relativistic MHD simulations by
Romanova et al. (2003, 2004, 2021). Even though the material from the torus
has an easier, natural path to the NS, the stream is thinner than the
“ring”-columns, and less material is accreted onto the NS.

Aside from the different accretion stream geometry, there is another
interesting effect. Since the equatorial region of the NS magnetic field is no
longer in the same plane as the material being accreted, magnetic reconnection
should occur less frequently. This results in more closed lines. Since there
is less additional open flux to strengthen the pulsar wind, this effect may
severely inhibit jet formation.

In Figure 2, we plot key quantities as a function of time. From top to bottom
panels, we show the mass accretion rate, open magnetic flux, jet power, NS
torque, and magnetospheric radius. We calculate the quantities as in Parfrey &
Tchekhovskoy (2023), and normalize some of them (open flux, jet power, and
torque) to those of an isolated pulsar at the same obliquity angle, obtained
from simulations without a torus and indicated with a “0” subscript.

The mass accretion rate is calculated by integrating the flux of mass at the
NS surface. The open magnetic flux, calculated at the light sphere, $r=R_{\rm
LC}$, is found by integrating over the force-free ($\sigma>1$) region,

$\Psi_{\rm
open}=\int_{0}^{2\pi}\int_{0}^{\pi}|B^{r}(r,\theta,\phi)|\sqrt{-g}\,{\rm
d}\theta\,{\rm d}\phi,$ (9)

where $\sigma=b^{2}/\rho h$ is the magnetization parameter and $B^{r}$ is a
measure of the radial magnetic field.

The jet power, also calculated at the light sphere and for the force-free
region, is measured as

$L=-\int_{0}^{2\pi}\int_{0}^{\pi}\left(b^{2}u^{r}u_{t}-b^{r}b_{t}\right)\sqrt{-g}\,{\rm
d}\theta\,{\rm d}\phi.$ (10)

The torque is calculated at the NS surface, and obtained by integrating the
hydrodynamic and electromagnetic contributions:

$N_{\rm
tot}=-\int_{0}^{2\pi}\int_{0}^{\pi}\left(b^{2}u^{r}u_{\phi}-b^{r}b_{\phi}+\rho
hu^{r}u_{\phi}\right)\sqrt{-g}\,{\rm d}\theta\,{\rm d}\phi.$ (11)

In this convention, a negative torque results in spindown of the star.

We define the magnetospheric radius as the location at which $\sigma=1$ in the
equatorial plane. Due to the geometry of the problem, we calculate two
different magnetospheric radii: $r_{\rm mag,1}$ is measured in the
$\mu-\Omega$ plane, while $r_{\rm mag,2}$ is defined in the plane
perpendicular to the $\mu-\Omega$ plane.

Figure 2 shows these quantities as a function of time for the simulations
performed with $\mu=20$. As the simulation starts, the MRI in the accretion
disk transports the angular momentum outward and results in the infall of the
torus material towards the NS. At $t\approx 7{,}500\,r_{\rm g}/c$, the
material reaches the surface of the NS, and at roughly $t=10{,}000\,r_{\rm
g}/c$ the quantities settle into an approximate steady state.

If we focus first on the magnetospheric radius (Figure 2e), it increases as we
increase obliquity, eventually becoming larger than the corotation radius.
When $r_{\rm mag}\gtrsim r_{\rm co}$, the pulsar enters the so-called
propeller regime, and little to no material is accreted as it would have to
overcome the centrifugal barrier.

High-obliquity pulsars tend to reach the propeller stage at lower magnetic
moments than their non-oblique counterparts. As seen in Figure 1, in the
$\mu-\Omega$ plane of the $\chi=60^{\circ}$ pulsars, the interactions with the
magnetosphere redirect the material from the torus onto each of the star’s
magnetic poles through a single, thin stream.

The magnetic field strength of a dipole is twice as high at its pole as at its
equator, and therefore tilting the pole increases the field strength at the
rotational equator in the $\mu$–$\Omega$ plane. This may partly explain why
the magnetospheric radius in the $\mu$–$\Omega$ plane ($r_{\rm mag,1}$) is
consistently larger than the measurement in the perpendicular plane ($r_{\rm
mag,2}$), where the field strength is unchanged from its aligned-rotator
value.

From Figure 2c, which shows the jet power, we can see that as we increase
obliquity the jet power decreases. At some angle (around $\chi=45^{\circ}$),
it resembles the isolated-pulsar value, and as we further increase the
obliquity angle the jet power becomes less than the isolated-pulsar spindown
power. The open magnetic flux in Figure 2b has a similar behaviour. Both of
these effects are partly due to the increase in the magnetospheric radius. At
high obliquity angles, the magnetospheric radius is larger than the corotation
radius the majority of the time. Material now has to break the centrifugal
barrier in order to accrete, thus, the mass accretion rate will be inhibited
(as can be seen in the mass accretion rate panel), and less magnetic flux will
be open.

We also note, when comparing different obliquity angles at the same magnetic
moment, that the isolated-pulsar spindown power is larger at higher obliquity.
Spitkovsky (2006) found a difference of two between aligned and orthogonal
rotators (the difference may be slightly higher; Tchekhovskoy et al. 2013).
Since in the isolated case $L\propto\mu^{2}$, the factor of two translates
into a difference in the effective magnetic moment of $\sqrt{2}$. In our
simulations, even though at $\mu=5$ all our runs are in the accretor (non-
propeller) regime, at $\mu=10$ and 20 we have the interesting case where low-
obliquity pulsars are in the accretor regime, whereas at larger obliquity
angles the pulsar is already in the propeller regime.

###  3.2 Effect of the magnetic moment and star–torus relative field
orientation

In Figure 3, we show simulations of a $\chi=60^{\circ}$, anti-parallel
magnetic field configuration with different magnetic moments. At lower
magnetic moments ($\mu=5$), we still see the formation of a single, thin
accretion stream in the $\mu-\Omega$ plane.

As we increase the magnetic moment, the pressure from the magnetosphere
increases, and the stream becomes thinner. There is very little dense material
inside the corotation radius, and thus the mass accretion rate drops. Notice
that at intermediate magnetic moment ($\mu=20$), there is still some
accretion, although it is severely inhibited. At even higher magnetic moment
($\mu=100$), the magnetospheric pressure is high enough that there is no
accretion onto the NS and the pulsar wind dominates. In our simulations at
$\mu=100$ the pulsars are fully in the propeller regime for all obliquity
angles.

As pointed out in Parfrey & Tchekhovskoy (2017, 2023), at low obliquity
angles, and particularly at low magnetic moments where there is no interchange
slingshot instability, there are major differences between the anti-parallel
and parallel configurations when the pulsar is far from the propeller state.
In the anti-parallel configuration, there is reconnection between the magnetic
field lines from the NS and the magnetic field lines of the torus with
opposite polarity. The reconnection causes more magnetic field lines to open
up, therefore increasing the power of the jet more than in the isolated
pulsar. Even though the difference between the two cases is more pronounced in
2D simulations, in 3D we also see this effect. As we increase obliquity,
though, the difference between the anti-parallel and parallel configurations
decreases, and at higher obliquity is not as influential as changing the
magnetic moment or the obliquity angle. Due to the geometry, in high-obliquity
pulsars reconnection between the magnetic field lines of the NS and the torus
is less effective. In that sense, high-obliquity pulsars are more similar to
the parallel configurations than the anti-parallel.

In addition, high obliquity angles are geometrically less favorable for flux
opening in general, as, beyond any given radius, less of the pulsar’s closed
flux passes through the region near the rotational equator, where it can
interact with the accretion flow. This would also suppress flux opening via
magnetic interchange.

These trends are clear in Figure 2, as the open magnetic fluxes — and
therefore jet powers — for the higher-obliquity cases are lower than for the
isolated pulsar case, just like for the small-$\mu$ parallel-orientation
configurations with low obliquity angles.

###  3.3 Perpendicular rotator

Figure 4 shows density snapshots of our perpendicular rotator,
$\chi=90^{\circ}$, for $\mu=10$. In the top panel, we show three different
planes. From left to right, we are showing (a) the plane perpendicular to
$\mu-\Omega$, (b) the $\mu-\Omega$ plane, and (c) a plane between these two.
We choose to show this particular magnetic moment, since at higher magnetic
moments, the pulsar is no longer in an accretor state. Similar to a
$\chi=60^{\circ}$ pulsar, accretion occurs in thin streams through the
magnetic poles. We can also note that accretion is not limited to the
$\mu-\Omega$ plane, or the plane perpendicular to the $\mu-\Omega$ plane. This
means that the hot-spot (where material falls onto the NS) is not a simple
ring, but has a rather complicated geometry. Indeed, unlike an aligned
rotator, the accretion does not occur in a “ring”-like manner.

##  4 Conclusions

In this work, we present the first systematic exploration, via 3D GRMHD
simulations, of gas interaction with rapidly rotating oblique NS
magnetospheres, spanning the full range of NS field strength (from the
accretor to propeller) and magnetic obliquity (from aligned to orthogonal).
For the first time for oblique NS accretion, we study the effect of disk
magnetic polarity: anti-parallel and parallel polarity relative to the NS
closed-zone magnetic field (able and unable to magnetically reconnect with the
NS closed-zone field, respectively).

###  4.1 Strong-field Accretor Regime

Let us first consider the case of relatively strong stellar dipolar magnetic
field, $\mu=20$, shown with green lines in Figure 5. Figure 5 shows that as
the inclination angle increases, the magnetosphere pushes the gas further away
(Figure 5e), which suppresses the accretion (Figure 5a), closes more magnetic
flux (Figure 5b), lowers jet luminosity (Figure 5c), and weakens spindown
torques (Figure 5d). Why does this happen? One possibility is that the
magnetic field strength of a dipole is twice as high at its pole as at its
equator: therefore tilting the pole increases the field strength at the
rotational equator in the $\mu$–$\Omega$ plane and might be able to push away
the gas more efficiently. Figure 5 (orange symbols) shows that the results are
similar for $\mu=10$.

While preparing this manuscript for publication, we became aware of a recent
preprint (Das & Porth, 2023) that also describes 3D GRMHD simulations of
accretion onto an oblique rotator. Their simulations also show that the jet
power decreases as the obliquity angle increases, therefore agreeing with our
results for $\mu\sim 10-20$.

Interestingly, there is little difference between the parallel and anti-
parallel torus magnetic field polarities for $10\lesssim\mu\lesssim 20$: for
instance, Figure 5b–e shows that for $\mu=20$, the green dashed and solid
lines, which show parallel and anti-parallel polarities, respectively, are
very close to each other.

A mechanism that leads to such similarity of the two magnetic polarities is
the _interchange slingshot_ instability (Parfrey & Tchekhovskoy, 2023), which
is a variation of the magnetic Rayleigh-Taylor instability that is activated
when the magnetospheric field is strong enough to push the gas close to, or
into, the propeller state. This instability, which is most vigorous when the
disk and magnetospheric magnetic field lines are parallel (so they cannot
reconnect), mixes the closed magnetospheric magnetic field lines into the
disk, causing them to open up and power the twin polar jets.

###  4.2 Weak-field Accretor Regime

With the help of the interchange slingshot instability, even the parallel
magnetic configuration can open magnetic field lines and launch jets for high
enough $\mu$ values. What happens at weaker stellar fields? We find that the
system can show more complex behavior. For example, at $\mu=5$, Figure 5 shows
that the two different polarities exhibit very different behavior. In the
anti-parallel configuration there is reconnection of the magnetic field lines
from the torus and the magnetosphere that increases the open magnetic flux
(Figure 5b), whereas in the parallel configuration there is no reconnection,
and no opening of the magnetic flux. This suppresses the jet power. However,
as $\chi$ increases and approaches the perpendicular case, the difference
between the two configurations becomes smaller.

Figure 5: Several quantities averaged in time, as a function of obliquity
angle. The quantities are the same as in Figure 2. The grey horizontal lines
represent $\Psi_{\rm open}=\Psi_{\rm open,0}$, $L=L_{0}$, $N_{\rm
tot}=-|N_{0}|$, $r_{\rm mag}=r_{\rm co}$. We averaged at times when the
quantities are in the steady state ($t>10{,}000\,r_{\rm g}/c$); the error
corresponds to the standard deviation over that time period. We show lines
that connect the values for simulations with $\mu=20$ (green), $\mu=5$ (blue)
for both parallel (solid line), and anti-parallel (dashed line)
configurations.

What causes the difference between the polarities to vanish at high obliquity?
Firstly, the two polarities are equivalent for an orthogonal rotator
($\chi=90^{\circ}$), and one would naturally expect the polarity-dependence to
disappear gradually as this limit is approached. In addition, Figure 5e shows
that high-obliquity pulsars approach the propeller regime, $r_{\rm mag}\gtrsim
r_{\rm co}$. For $\mu=5$, the transition to the propeller regime happens at
$\chi\simeq 60^{\circ}$, leading to the activation of the interchange
slingshot instability that opens up the magnetic flux and therefore increases
the jet power, as seen in Figure 5c. Remarkably, this leads to opposing jet
power trends with the obliquity angle for the two star–disk magnetic
orientations at $\mu=5$: whereas $L(\chi)$ decreases for the anti-parallel
polarity, $L(\chi)$ increases for the parallel polarity. We conclude that the
$\chi$-dependence of the various key quantities also significantly depends on
the magnetic moment and the relative magnetic polarity.

###  4.3 Propeller Regime

Our simulations show that high-obliquity pulsars enter the propeller regime
($r_{\rm mag}>r_{\rm co}$) at lower magnetic moments than low-obliquity or
aligned pulsars. If we focus on low magnetic moments ($\mu=5$), we can see
that in all cases, we have some non-zero mass accretion rate that decreases
with increasing $\chi$ (Figure 5a). High-obliquity pulsars accrete through a
thin stream, whereas low-obliquity pulsars accrete in a “ring”-like manner.
For $\mu=5$, $r_{\rm mag}\lesssim r_{\rm co}$ at all values of $\chi$ (Figure
5e), thus all pulsars we simulated for $\mu=5$ are in the accretor regime. If
we increase the magnetic moment while keeping the gas reservoir constant, the
difference between the low- and high-obliquity cases becomes more noticeable.
The mass accretion rate decreases, and the magnetospheric radius increases as
$\chi$ increases at constant $\mu$. In our simulations, the high-obliquity
pulsars reach the propeller ($r_{\rm mag}>r_{\rm co}$) regime between $\mu=5$
and $\mu=10$, whereas the aligned rotators change state between $\mu=20$ and
$\mu=100$. It is interesting that sometimes pulsars in the propeller regime
exhibit a non-zero mass accretion rate.

Through mapping out the parameter space, we explored the transitions between
the accretor and propeller regimes. For example, for $\mu=20$, Figure 5 shows
two transitions somewhere between $\chi\sim 30^{\circ}$ and $\chi\sim
45^{\circ}$: the magnetospheric radius passes through and exceeds the
corotation radius, and the jet power goes from being larger than the isolated-
pulsar spindown power (i.e. $L>L_{0}$) to being less than this value. These
intermediate states occur between the full accretor and full propeller states.

This research was supported in part through the computational resources and
staff contributions provided for the Quest high performance computing facility
at Northwestern University which is jointly supported by the Office of the
Provost, the Office for Research, and Northwestern University Information
Technology. A.M-B. is supported by NASA through the NASA Hubble Fellowship
grant HST-HF2-51487.001-A awarded by the Space Telescope Science Institute,
which is operated by the Association of Universities for Research in
Astronomy, Inc., for NASA, under contract NAS5-26555. KP, JJ, and AT
acknowledge support by NASA 80NSSC21K1746 grant. AT and JJ also acknowledge
support by NASA 80NSSC22K0938 and NSF AST-2009884 grants. AT was additionally
supported by NSF grants AST-2107839, AST-1815304, AST-1911080, OAC-2031997,
and AST-2206471. KP was supported in part by the Laboratory Directed Research
and Development Program at Princeton Plasma Physics Laboratory, a national
laboratory operated by Princeton University for the U.S. Department of Energy
under Prime Contract No. DE-AC02-09CH11466. Support for this work was also
provided by the National Aeronautics and Space Administration through Chandra
Award Number TM1-22005X issued by the Chandra X-ray Center, which is operated
by the Smithsonian Astrophysical Observatory for and on behalf of the National
Aeronautics Space Administration under contract NAS8-03060. This research was
also made possible by NSF PRAC award no. 1615281 at the Blue Waters sustained-
petascale computing project and supported in part under grant no. NSF
PHY-1125915. This research used resources of the Oak Ridge Leadership
Computing Facility, which is a DOE Office of Science User Facility supported
under Contract DE-AC05-00OR22725 via ALCC, INCITE, and Director Discretionary
allocations PHY129.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Archibald et al. (2009) Archibald, A. M., Stairs, I. H., Ransom, S. M., et al. 2009, Science, 324, 1411, doi: 10.1126/science.1172740
  * Bachetti et al. (2014) Bachetti, M., Harrison, F. A., Walton, D. J., et al. 2014, Nature, 514, 202, doi: 10.1038/nature13791
  * Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214, doi: 10.1086/170270
  * Balbus & Hawley (1998) —. 1998, Reviews of Modern Physics, 70, 1, doi: 10.1103/RevModPhys.70.1
  * Bozzo et al. (2018) Bozzo, E., Ascenzi, S., Ducci, L., et al. 2018, A&A, 617, A126, doi: 10.1051/0004-6361/201732004
  * Chakrabarti (1985) Chakrabarti, S. K. 1985, ApJ, 288, 1, doi: 10.1086/162755
  * Das & Porth (2023) Das, P., & Porth, O. 2023, arXiv e-prints, arXiv:2311.05301.  https://arxiv.org/abs/2311.05301
  * Das et al. (2022) Das, P., Porth, O., & Watts, A. L. 2022, MNRAS, 515, 3144, doi: 10.1093/mnras/stac1817
  * de Martino et al. (2014) de Martino, D., Casares, J., Mason, E., et al. 2014, MNRAS, 444, 3004, doi: 10.1093/mnras/stu1640
  * Duncan & Thompson (1992) Duncan, R. C., & Thompson, C. 1992, ApJ, 392, L9, doi: 10.1086/186413
  * Fishbone & Moncrief (1976) Fishbone, L. G., & Moncrief, V. 1976, ApJ, 207, 962, doi: 10.1086/154565
  * Gammie et al. (2003) Gammie, C. F., McKinney, J. C., & Tóth, G. 2003, ApJ, 589, 444, doi: 10.1086/374594
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357–362, doi: 10.1038/s41586-020-2649-2
  * Hawley & Balbus (1991) Hawley, J. F., & Balbus, S. A. 1991, ApJ, 376, 223, doi: 10.1086/170271
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Israel et al. (2017) Israel, G. L., Belfiore, A., Stella, L., et al. 2017, Science, 355, 817, doi: 10.1126/science.aai8635
  * Metzger et al. (2011) Metzger, B. D., Giannios, D., Thompson, T. A., Bucciantini, N., & Quataert, E. 2011, MNRAS, 413, 2031, doi: 10.1111/j.1365-2966.2011.18280.x
  * Noble et al. (2006) Noble, S. C., Gammie, C. F., McKinney, J. C., & Del Zanna, L. 2006, ApJ, 641, 626, doi: 10.1086/500349
  * Parfrey et al. (2017) Parfrey, K., Spitkovsky, A., & Beloborodov, A. M. 2017, MNRAS, 469, 3656, doi: 10.1093/mnras/stx950
  * Parfrey & Tchekhovskoy (2017) Parfrey, K., & Tchekhovskoy, A. 2017, ApJ, 851, L34, doi: 10.3847/2041-8213/aa9c85
  * Parfrey & Tchekhovskoy (2023) —. 2023, arXiv e-prints, arXiv:2311.04291.  https://arxiv.org/abs/2311.04291
  * Philippov et al. (2014) Philippov, A., Tchekhovskoy, A., & Li, J. G. 2014, MNRAS, 441, 1879, doi: 10.1093/mnras/stu591
  * Remillard & McClintock (2006) Remillard, R. A., & McClintock, J. E. 2006, ARA&A, 44, 49, doi: 10.1146/annurev.astro.44.051905.092532
  * Reynolds (2014) Reynolds, C. S. 2014, Space Sci. Rev., 183, 277, doi: 10.1007/s11214-013-0006-6
  * Romanova et al. (2021) Romanova, M. M., Koldoba, A. V., Ustyugova, G. V., et al. 2021, MNRAS, 506, 372, doi: 10.1093/mnras/stab1724
  * Romanova et al. (2004) Romanova, M. M., Ustyugova, G. V., Koldoba, A. V., & Lovelace, R. V. E. 2004, ApJ, 610, 920, doi: 10.1086/421867
  * Romanova et al. (2012) —. 2012, MNRAS, 421, 63, doi: 10.1111/j.1365-2966.2011.20055.x
  * Romanova et al. (2013) —. 2013, MNRAS, 430, 699, doi: 10.1093/mnras/sts670
  * Romanova et al. (2003) Romanova, M. M., Ustyugova, G. V., Koldoba, A. V., Wick, J. V., & Lovelace, R. V. E. 2003, ApJ, 595, 1009, doi: 10.1086/377514
  * Spitkovsky (2006) Spitkovsky, A. 2006, ApJ, 648, L51, doi: 10.1086/507518
  * Steiner et al. (2009) Steiner, J. F., McClintock, J. E., Remillard, R. A., Narayan, R., & Gou, L. 2009, ApJ, 701, L83, doi: 10.1088/0004-637X/701/2/L83
  * Tchekhovskoy (2019) Tchekhovskoy, A. 2019, HARMPI: 3D massively parallel general relativictic MHD code, Astrophysics Source Code Library, record ascl:1912.014.  http://ascl.net/1912.014
  * Tchekhovskoy et al. (2016) Tchekhovskoy, A., Philippov, A., & Spitkovsky, A. 2016, MNRAS, 457, 3384, doi: 10.1093/mnras/stv2869
  * Tchekhovskoy et al. (2013) Tchekhovskoy, A., Spitkovsky, A., & Li, J. G. 2013, MNRAS, 435, L1, doi: 10.1093/mnrasl/slt076
  * The HDF Group (1997-NNNN) The HDF Group. 1997-NNNN, Hierarchical Data Format, version 5
  * Turk et al. (2011) Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9, doi: 10.1088/0067-0049/192/1/9
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Walter et al. (2015) Walter, R., Lutovinov, A. A., Bozzo, E., & Tsygankov, S. S. 2015, A&A Rev., 23, 2, doi: 10.1007/s00159-015-0082-6
  * Wang (1997) Wang, Y. M. 1997, ApJ, 475, L135, doi: 10.1086/310481
  * Wasserman & Shapiro (1983) Wasserman, I., & Shapiro, S. L. 1983, ApJ, 265, 1036, doi: 10.1086/160745
  * Woosley (2010) Woosley, S. E. 2010, ApJ, 719, L204, doi: 10.1088/2041-8205/719/2/L204