Magnetospheres of black hole-neutron star binaries

We perform force-free simulations for a neutron star orbiting a black hole, aiming at clarifying the main magnetosphere properties of such binaries towards their innermost stable circular orbits. Several configurations are explored, varying the orbital separation, the individual spins and misalignment angle among the magnetic and orbital axes. We find significant electromagnetic luminosities, $L\sim 10^{42-46} \, [B_{\rm pole}/ 10^{12}{\rm G}]^2 \, {\rm erg/s}$ (depending on the specific setting), primarily powered by the orbital kinetic energy, being about one order of magnitude higher than those expected from unipolar induction. The systems typically develop current sheets that extend to long distances following a spiral arm structure. The intense curvature of the black hole produces extreme bending on a particular set of magnetic field lines as it moves along the orbit, leading to magnetic reconnections in the vicinity of the horizon. For the most symmetric scenario (aligned cases), these reconnection events can release large-scale plasmoids that carry the majority of the Poynting fluxes. On the other hand, for misaligned cases, a larger fraction of the luminosity is instead carried outwards by large-amplitude Alfv{\'e}n waves disturbances. We estimate possible precursor electromagnetic emissions based on our numerical solutions, finding radio signals as the most promising candidates to be detectable within distances of $\lesssim 200$\,Mpc by forthcoming facilities like the Square Kilometer Array.


I. INTRODUCTION
The spectacular combined detection of gravitational waves (GW) and electromagnetic (EM) radiation from the binary neutron star (BNS) merger GW170817 [1][2][3] has initiated the new era of multi-messenger astronomy. This single multi-messenger observation led to important breakthroughs in our understanding of astrophysics and even fundamental physics. It confirmed that BNS mergers give rise to short gamma-ray bursts and the production of heavy elements via r-process; as well as it has provided constraints on proprieties of matter at nuclear densities, an estimation of the Hubble constant, and stringent tests of general relativity, among others (e.g., [4]).
The presence of matter and its support for strong magnetic fields make binary systems involving at least one neutron star (NS) the most likely sources for such simultaneous detection of GW and EM signals. The interaction of the NS magnetosphere with its companion during the inspiral phase is expected to generate EM emissions. As they would be produced within the relatively cleaner environment preceding the merger, these precursor EM counterparts -not yet detected-could provide crucial information about the sky localization of the source and other physical parameters of the system which cannot be accurately obtained from the GW observation alone. For black hole-neutron star (BHNS) binaries these precursor EM signals could be particularly relevant, since in most cases these systems are not expected to produce strong post-merger EM emissions: it is known for high enough binary mass-ratios ( 5) and for a slowly spinning black hole (BH), the NS is likely to be swallowed by the BH without tidal disruption (e.g., [5]), and thus, not favoring the ejection of material from the system nor the formation of an accretion disk to support further EM counterparts. In such cases, the existence of precursor EM signals could be the only way to distinguish a BHNS from a binary BH system of the same masses.
A NS in a compact binary system at late stages prior to merger is likely to be surrounded by a tenuous plasma, well described by the force-free (FF) approximation [6]. The arguments around this assumption date back to the pioneering work of Goldreich & Julian in the context of pulsars [7], recently adapted to binaries scenarios in [8]. Such plasma environment enclosing the binary would allow to tap kinetic energy from the orbital motion (or from individual spins) and, eventually, a part of this EM energy will be re-processed within the magnetosphere leading to emissions on different bands of the EM spectrum. These central ideas have been well established for isolated pulsars (e.g., [7,[9][10][11][12]) and spinning BHs in active galactic nuclei (e.g., [13,14]).
A few theoretical models have been invoked to describe the energy extraction mechanisms in binary scenarios. One of them focuses on the EM energy-loss produced by the orbital motion of the NS magnetic dipole moment [15]. The other focuses on the interaction of the companion moving across the magnetized medium surrounding the NS, typically modeled by the unipolar induction (UI) effect. This idea was introduced long time ago for moving conductors such as satellites around a planet [7,16] and later applied to compact binaries [17][18][19][20][21][22], with the BH acting as a battery in a DC circuit with the NS for the case of a BHNS binary.
Several numerical studies have shown how, indeed, a fraction of the available kinetic energy of compact binary systems can be transferred into the surrounding FF plasma. General relativistic (GR) simulations matching magnetohydrodynamics stellar interiors with an exterior FF magnetosphere were performed for BNS [23][24][25] and, more recently, for BHNS binaries [26]. A similar GR FF approach, but assuming a helical Killing vector field, has been carried for the BHNS scenario as well [27]. On the other hand, special relativistic simulations have focused on the energy extraction associated with magnetic dipole radiation by the NS orbital motion [28] and the flaring events associated with the twisted magnetic fluxtube arising from the relative motion of two NSs [29]. This paper aims at further clarifying the main magnetospheric properties of BHNS binaries in their last orbits up-to the innermost stable circular orbit of the system. We assume that the NS is endowed with a dipolar magnetic field and moving around the BH through an FF plasma environment. Instead of performing full GR simulations, we approximate the problem by a fixed Kerr spacetime. Thus, while we are neglecting the curvature of the NS, we are prescribing its trajectory over the BH background. The former assumption might be acceptable acknowledging the dominant dynamical role of the boundary conditions, i.e., anchoring the magnetic field into the NS crust. While the latter is based on the fact that binary trajectories would be essentially determined by the energy-loss due to GWs, which are energetically predominant over any realistic EM interaction [15,20]. This simplified setting allows us to conduct rather inexpensive, very accurate, numerical simulations for a detailed description of the magnetospheric response and to explore different relevant configurations and parameters.
In this way, we extend our previous special relativistic study [28] to investigate how the intense curvature of the BH affects the resulting BHNS magnetosphere. In this paper, we focus only on the binary in circular orbits at several separations close to the innermost stable circular orbit. We first consider the case in which the magnetic moment is aligned with respect to the orbital angular momentum and both compact objects are non-spinning. Then, we independently incorporate their spins (aligned with the orbital axis), and later we consider various inclinations of the dipole magnetic moment. Besides providing estimations for the amount of kinetic energy that can be extracted from the binary as Poynting luminosity (and/or dissipation) on each setting, we assess how the EM energy is distributed within the magnetosphere and pay particular attention to, e.g., the formation of current sheets (CSs). Finally, relying on recent progresses in pulsar theory (e.g., [30][31][32][33][34]), we infer possible EM precursor emissions from the attained numerical solutions.
The article is organized as follows. In Sec. II we setup the problem and describe our numerical implementation. Results are presented in Sec. III, first focusing on irrotational and aligned BHNS binaries magnetospheres, and then incorporating the effects of spin and misalignment. Possible observational implications of our results are also discussed. We summarize and conclude in Sec. IV. Throughout this paper, unless otherwise stated, we use the geometrical units of G = 1 = c, where G is the gravitational constant and c is the speed of light. The Latin indexes i, j, and k denote spatial components while other Latin indexes (a, b, c, ..) the spacetime ones.
II. SETUP

A. General Setting
As discussed above, we want to consider a star that is following certain trajectory over a Kerr background. The Kerr metric is parametrized by the mass M and spin a, and can be written in the Kerr-Schild form as g ab = η ab + H a b , where η ab is the flat metric and a is a null co-vector with respect to both η ab and g ab . We shall describe the trajectory of the orbiting NS in the Cartesian coordinates {x a } = {t, x, y, z} associated with the flat part of the metric 1 . In these coordinates, the metric function H takes the form and the co-vector a reads a = 1, rx + ay r 2 + a 2 , ry − ax r 2 + a 2 , z r .
(4) a is often defined by η ab b , and with this definition, we have a a = 0 and k k = 1. Note that g ab is written as In the 3 + 1 formulation, the line element is written as with representing the lapse function, shift vector, and spatial metric, respectively. The trajectory is defined by the orbital separation r o (t) and the phase ϕ o (t) (being Ω o (t) ≡φ o (t) the associated angular velocity). However, our numerical domain will be centered on the NS instead (as in [28]), thus describing the dynamics for an adapted foliation with coordinates {x a } = {t,x,ŷ,ẑ}. The coordinates transformation into this "co-moving" foliation is defined by: y =ŷ + r o (t) sin(ϕ o (t)), and thus, Therefore, the line element (5) in the new coordinates is written simply as where a new piece of the shift vector due to the orbital motion is In the above, to simplify the notation, we have dropped all time dependencies and defined v o := r o Ω o . "(x(x))" on the metric components above implies the same functional form of the previous ones, but evaluated at the "displaced" spatial points as seen in the new coordinate system.

B. Evolution Equations
Force-free electrodynamics (FFE) is governed by Maxwell's equations, together with the FF condition, Here, F ab is the EM tensor, F * ab its dual, and ∇ b the covariant derivative with respect to g ab . In this regime, remarkably, the EM field can be evolved autonomously without keeping track of the plasma degrees of freedom. This is achieved by finding an expression for j a in terms of the EM field.
Following [5], we can write the corresponding evolution equations for this system in a pseudo-conservative form: where E i := √ γE i and B i := √ γB i are the densitized variables with E a := F ab n b and B a := −F * ab n b (n a is the unit time-like vector normal to the spatial hypersurfaces), and their associated fluxes are given by: with abc = n d dabc being the induced volume element on the spatial hypersurfaces. The FF current reads, where Q := √ γρ e = ∂ k E k is the electric charge density 2 .
However, this formulation of the theory is only weakly hyperbolic [36]. Therefore, one has to use either a symmetrized version of the theory like the one we developed in [37] or a simplified version of the equations without the terms parallel to B i in the current (22), i.e., while numerically enforcing E ·B = 0 (like in [12,38,39]). We choose this simplified version of FFE since it is much straightforward to implement and has shown to produce very similar results to those with the more elaborated system [37] (this was tested on the similar setting of [28]).

C. Initial and Boundary Data
The initial configuration is set by a magnetic dipolar field and vanishing electric field. We approximate this EM field from its flat spacetime vector potential solution, centered in our adapted coordinates {x a }, where µ is the magnetic dipole-moment. We focus primarily on the cases in which the magnetic axis is aligned with the orbital angular momentum, but we also consider other scenarios in which these two axes are not aligned. In our setup, theẑ axis is always perpendicular to the orbital plane, and thus, the misalignment is attained by just tilting the initial magnetic moment by an angle χ along thex-ẑ plane. We will consider relevant (fixed) values of the orbital separation r o , setting the orbital angular velocity accordingly as (e.g., chapter 12 of [40]) Here we incorporate the effect of the BH spin on the orbital motion phenomenologically mimicking the angular velocity for a particle orbiting the spinning BH [40]. The NS is gradually set in motion (at this fixed r o ) until it reaches the corresponding orbital frequency (26) and, from then on, the system follows its prescribed circular orbit trajectory.
The boundary condition at the NS surface is derived as in [28] by assuming the perfectly conducting condition, which can be easily generalized to incorporate the NS spin at frequency Ω * by, whereφ a ≡ (∂φ) a . Thus, the resulting condition on the electric field measured by a fiducial observer in this adapted foliation (i.e.Ê a := F abn a ) can be written aŝ

D. Numerical Implementation
We evolve the simplified version of FFE described above, with the coordinate system centered on the NS and with the inclusion of a divergence cleaning field φ, i.e.: Our numerical scheme to solve these equations is based on the multi-block approach [41][42][43][44], relying on higherorder finite difference operators with an adapted Kreiss-Oliger dissipation [45] and Runge-Kutta method for integration in time. The numerical approach, including the treatment given to boundary conditions and CSs, is exactly the same as in [28]; hence, we refer the readers to [28] (and to the references therein) for further details.
In the present version, the FF condition E · B = 0 is enforced through the projection, Eî → Eî − E · B B 2 Bî, at each Runge-Kutta sub step.
Since the BH interior is contained within our numerical domain, we need to regularize the metric functions inside the outer horizon located in the Boyer-Lindquist coordinates at r + = M BH + M 2 BH − a 2 . We shall do it by first rewriting H and i as where r i := {x, y, z} and v i := {y, −x, (a/r)z}. Then, with λ being less than unity (typically set to λ = 0.6), and p(r/r c ) a large polynomial in r/r c built to guarantee smoothness across the matching surface r = r c .
Of course, such modifications to the metric inside the BH horizon do not affect the physics outside, since these two regions are causally disconnected. However, we further prevent the production and amplification of highfrequency numerical noises inside the BH by damping the electric and magnetic fields through the following replacements: with damping parameter κ = 60M −1 BH in our simulations. We solve the system in a region between an interior sphere at radiusr = R * , which represents the NS surface (i.e., R * denotes the NS radius), and an exterior spherical surface located atr ≈ 128R * . The domain is represented by a total of 6 × 7 sub-domains, with 6 patches to cover the angular directions and 7 being the number of spherical shells expanding in radius. These spherical shells do not overlap each other, and have more resolution in the inner regions: from layer to layer, the radial resolution is decreased by a factor 2. Typically, we adopt a resolution with total grid numbers of Nθ × Nφ × Nr with Nφ = 2Nθ = 320, while Nr = 560 to span the whole computational domain. To give an idea, for a binary of mass ratio q(= M BH /M NS ) = 3 with orbital separation r o = 10M BH at this resolution, the BH is covered by 20 × 20 × 35 grid-points.

E. Analysis Quantities
We monitor the EM energy and its associated fluxes. In FFE the four-momentum, p a = −T ab EM t b , is conserved (i.e. ∇ a p a = 0) 3 in stationary spacetimes where a is a Killing field. Hence, we measure the corresponding energy and fluxes by where the Poynting luminosity L is being integrated on spherical surfaces of radius r around the BH and, We also monitor the charge distribution and electric currents present during the dynamics. Thus, we look at the FF current density along the magnetic field, as seen by an observer moving in the direction normal to the spatial slices n a ,

III. RESULTS
We aim at understanding in detail the magnetospheric properties of the steady states of BHNS binary systems at given circular orbits near their innermost stable circular orbit. We thus extend our previous studies [28], paying a particular attention to the effect that the intense curvature in the vicinity of the BH has on the magnetosphere: e.g., in the magnetic field topology, the induced electric currents and charge distributions, and in the obtained Poynting luminosity. Finally, we elaborate on the implications of these results to potential EM observations.

A. Circular Orbits: general features
We first consider circular orbits at fixed separation r o = 10M BH , for which both compact objects are nonspinning and the magnetic moment of the NS is aligned to the orbital angular momentum. As mentioned in Sec. II C, in the simulations, the NS is gradually set into motion until it attains its final orbital angular velocity given by Eq. (26). The system settles into a quasistationary regime after about 2 orbits. In Fig. 1 we illustrate the main features of the solution at 3 orbital periods, showing both compact objects NS/BH as grey/black surfaces together with carefully chosen magnetic field lines. The contour-plots in the figure represent regions of very large charge density (exceeding the Goldreich-Julian (GJ) values [46]) signaling the presence of strong CSs. One may distinguish, in principle, two kinds of CSs in the picture. There is a CS that follows a spiral pattern over the orbital plane, associated with discontinuities of the toroidal magnetic field across the plane, being the sign of the charge density (depicted by the red/blue contours in the plot) indicative of the orientation of this "jump". And there is another, out-of-the-plane CS, produced by a discontinuity of the poloidal magnetic field at an outer side of the BH horizon. It can be seen that the dynamical effect of the curvature is quite strong for a set of magnetic field lines, which bend towards the BH, suffering extreme twisting of about ∼ π and leading to X-point magnetic reconnections. Such continuous reconnections taking place near the BH produce closed magnetic loops (i.e., plasmoids) that are released at relativistic speeds, transporting electromagnetic energy away from the system 4 . The large-scale plasmoids can be clearly appreciated at the right side of the plots. For the larger loops, produced closer to the BH, the field strengths are higher (as also noticed in [26]). On the other hand, the smaller ones, involving weaker magnetic fields, enclose the blue contour component of the equatorial CS (see Fig. 1).
The Poynting flux distribution concentrates along the same spiral arm of the red component on the equatorial CS, as can be seen in Fig. 2. The reason for this is that this region of the magnetic loops is where most of the twisting concentrates, and thus, the strength of the magnetic field, and hence the EM energy, are largely enhanced. The bottom panel of Fig. 2 presents a skymap distribution of such EM fluxes on a distant sphere around the BH. The intensity peaks on a narrow band in the azimuthal direction and within ∼ 60°from the orbital plane. This would suggest a lighthouse effect being established with the setup orbital frequency, consistent with previous results (e.g., [27,28]).
We compute the luminosity for our late-time numerical solutions by integrating (38) at some distant radius around the BH in the wave-zone (typically, we chose r = 40M BH ). We normalize them with the EM luminosity given by the MD radiation formula [15], although it should be taken only as a reference value, since there is no particular reason to chose this scaling in the present GR modeling. Convergence of the luminosity is assessed by considering the results with three different numerical resolutions. Since CSs extend to a region very close to the BH horizon, one does not expect formal convergence of the numerical solutions in the wave-zone, as they would generally depend on the dissipation treatment given to these CSs. However, by keeping a fixed dissipation rate (controlled by a single parameter in our case), we achieved reasonable numerical convergence for the luminosity. The results are summarized in Fig. 3, in which the luminosity is shown as a function of the integration radius for the three resolutions. We notice that the lowest resolution with N θ = 160 (the one used throughout this work) is already practically convergent as it differs in less than 5% with respect to the higher resolution runs. The luminosity decreases slightly with the increase of the extraction radius, because a part of the Poynting flux is dissipated on the CSs. Figure 4 shows the luminosity as a function of the orbital separation. A reference scaling is presented as well, in order to compare with the results found in the absence of the BH in [28]. First, we notice that the luminosity found here is 15 to 30 times larger than those from the single orbiting NS scenario, and that the scaling is here steeper towards the closest binary orbits. This indicates that the new effect contributed by the BH curvature is dominant, and thus, the scaling associated with it certainly differs from the reference, Eq. (40) (even if accounting for the relativistic corrections found in [28]). It is important to remark that the luminosity predicted in the unipolar induction (UI) model posses the same scaling as L MD for a non-spinning BH, and only differs by a factor ∼ 5.6 (larger) for a binary of mass-ratio q = 3. Thus, as also found in [26], the resulting luminosity is significantly higher than the theoretical expectations from the UI mechanism, especially near the last stable orbits.
An estimated electric dipole moment in the near zone (approximated as r<c/Ωo Q x i d 3 x, with respect to the CoM) yields higher values by a factor of 5-10 than the magnetic dipole moment of the NS. This indicates that the induced charges and currents within the magnetosphere are likely to be the origin of the high EM lumi-nosity obtained. As magnetic reconnections occur near the BH, the reconnection site is closer to the NS than in the scenario without BH studied in [28] where the CS begin at ∼ c/Ω o . This implies that the magnetic field is here much stronger, and may explain the enhanced luminosity.
We also note the normalized luminosity depends only very weakly on the mass-ratio, q, when re-scaling the problem by the BH mass (at a fixed value of M NS ). Therefore, these results are quite generic in that respect, and can be extended directly to binaries with higher mass-ratios.
As shown by Eq. (40) Thus, near the innermost stable circular orbits, the luminosity L is likely to be enhanced up to ∼ 10 42 erg/s for the non-spinning BHs. For rapidly spinning BHs with the spin axis aligned with the orbital angular momentum, the smaller orbital radius can be taken. The luminosity is computed at r = 40MBH and normalized by LMD (see Eq. (40)). A reference scaling (dotted lines) from previous results on a single orbiting NS (i.e., without the BH curvature) [28], re-scaled here by a factor of 20, is included.
allowed (see below, Sec. III B). We now continue our detailed description of the magnetosphere, by observing the late time numerical solutions at a particular plane intersecting both compact objects. Figure 5 shows charge density, magnetic field strength, and parallel electric currents (Eq. (39)), together with some magnetic field lines projected onto this co-plane. The CSs are clearly found in the images, both the spiral CS on the orbital plane and the out-of-theplane component (seen at the left side of the BH). The figure also presents a clean view of the plasmoids (i.e., their projection onto the co-plane): see, for instance, all the closed magnetic loops to the left of the BH in the left panel. As seen before, both charge and current densities are very intense (far beyond the GJ values) along the CSs. There is also significant current flowing along those field lines that connect the BH to the NS. Such "electric circuit", depicted more clearly on the right panel of Fig. 5, is highly reminiscent of the UI mechanism (e.g., [17,18,20,21]).
One may then ask whether the simplified description of the UI model is able to capture the main dynamics and to explain the luminosity which we measure in the wave zone. The main question is perhaps how the energy provided by the DC circuit could be transported to large distances, since Alfvén waves would -in principle-remain trapped within the flux-tube connecting the two compact objects. In this respect, the magnetic reconnections which we observe in the vicinity of the BH (not taken into account in the UI model) may play a fundamental role in allowing EM energy to escape. Such extreme bending FIG. 5. Electric charge and current on the co-plane (x-z plane) for a BHNS binary in a circular orbit at an orbital separation ro = 10MBH. The numerical solution on the co-plane after 2 orbits is displayed. Charge density distribution (top-left panel) and parallel electric currents (bottom-left panel) are normalized by ΩoB/2πc and ΩoB/2π, respectively. Top-right plot shows the magnetic field strength, normalized by its value at the poles while bottom-right panel is a zoom-in version of the bottom-left one. The circle with its center at x = z = 0 shows the NS surface, and the thick black circle (located at x < 0) shows the BH horizon; The green arrows in the bottom-right are plotted to emphasize the direction of the currents, which are suggestive of the electronic circuit predicted by the UI model. The open circle around the origin shows the NS and the circle of its center at x = −10MBH is the BH horizon.
of magnetic field lines near the BH equator has been noticed before in the context of boosted BHs moving across a magnetized plasma [47,48], and identified as a pure gravitational effect in [49]. Here, this phenomenon could circumvent the theoretical expectation that asserts the azimuthal twist of the flux-tube in the DC circuit for a BHNS binary is 1 [20].
Nevertheless, we do find an approximate realization of the UI mechanism in these strong currents flowing inside the flux-tube, which can be regarded as an additional contribution to the energy budget (i.e., comple-mentary to the luminosity measured at larger distances). This contribution is expected to propitiate the generation of gamma-rays from synchrotron/curvature radiation, as well as thermal X-ray emissions from the accelerated particles that hit the NS surface forming hot-spots (see, e.g., [19]). tures due to the impact of magnetic reconnections (especially those taking place near the BH). There are two spots at each hemisphere, linked to the ingoing/outgoing currents of the 'instantaneous' electric circuit. We here estimate the local temperature of such hot-spots as expected from bombardment of relativistic particles, following a simple model from [50] in which the kinetic energy deposition is estimated by assuming that a fixed fraction κ of the total current density | J| is carried by relativistic particles with averaged Lorentz factorγ traveling towards the stellar surface (see also [51] for details). The rate at which this energy is deposited is finally equated to the power radiated as a black body, yielding: where m e , σ, and e denote the electron mass, Stefan-Boltzmann constant and elementary charge, respectively. Plugging the numbers from our simulations, we get  (43) for a binary of mass-ratio q = 3 at r o = 10M BH assuming that the current density is 20 times larger than the GJ value. We note that the temperature could be higher if the electron-positron pair creation occurs near the polar region efficiently.

B. Spin effects
In this section, we pay attention to the cases in which there is BH or NS spin aligned to the orbital angular momentum. We first study the effect of the NS spin, which is expected to be less relevant from the astrophysical perspective, since in the late inspiral phase tidal locking is unlikely to be achieved and the orbital frequency is usually much higher than the NS spin frequency [52]. Here, we only look at the most extreme case in which the spin of the NS is as high as the orbital frequency (i.e., Ω * = Ω o ) in a circular orbit at fixed orbital separation r o = 10M BH . For this case, the magnetosphere co-rotates with the orbital motion, and hence, the magnetic field lines penetrate the BH horizon in a stationary manner. This situation is highly different from that considered in the previous section, and thus, the properties of the magnetosphere are also expected to be modified significantly.
We find an EM luminosity of L ∼ 56L MD (at r = 40M BH ), which is about 50% higher than the one obtained for the same setting with non-spinning NS. This luminosity is exactly twice the spin-down luminosity of an isolated pulsar of the same period (∼ 4 ms for a binary of mass-ratio q = 3), estimated as L pulsar ≈ µ 2 Ω 4 * [12]. Very similar enhancement with respect to the pulsar spindown luminosity was already observed in [28] for analogous configuration but without the BH. This suggests that the dynamical role played by the BH is here subdominant in contrast to the case of irrotational BHNS binaries discussed in the previous section. Indeed, the co-rotation of the NS with the orbit diminishes the relative velocity among the plasma and the BH, dramatically reducing its impact on the magnetic-field structure. The magnetosphere, thus, looks overall rather similar to that of an aligned pulsar but with the spiral arm structure in its Poynting flux resulting from the orbital motion. An equatorial CS develops near the BH (similar to the irrotational case) and close to the LC at the opposite side of the NS, with all the magnetic field lines opening-up beyond that point (see Fig. 7).
There is, however, some additional enhancement of the magnetic field strength for those field lines threading the BH, associated with the current-carrying flux-tube between the two compact objects. The polar cap regions at the NS surface, being dominated in this case by the NS spin, are slightly modified by these flux-tube currents.
A sky-map distribution of Poynting flux is displayed in Fig. 8 for such tidally locked BHNS binary system. It looks very similar to that of the irrotational case (i.e., bottom panel of Fig. 2) but with its intensity being even FIG. 7. Magnetosphere of a BHNS binary in a circular orbit at an orbital separation ro = 10MBH with synchronized NS spin, Ω * = Ωo, after 3 orbits. Top panel: Charge density distribution (left) and parallel electric current (right) on the co-plane, normalized by ΩoB/2πc and ΩoB/2π, respectively. The circle with its center at x = z = 0 shows the NS surface, and the thick black circle (located at x < 0) shows the BH horizon. Bottom panel: 3D view of the magnetosphere (left) with contours of the charge density signaling the equatorial CS. Induced polar cap on the stellar surface (right), indicated by a color map of the magnetospheric currents reaching the NS. more concentrated towards the orbital plane.
In the following, we switch-on the BH spin in the direction perpendicular to the orbit and explore its impact on the magnetosphere. The dimensionless spin a/M BH is varied from −0.99 to 0.99 for a circular orbit at fixed orbital separation r o = 10M BH . Generally speaking, we find that the effect of the orbital motion rather than that of the BH spin determines the magnetosphere and EM luminosity, even for the near-extreme Kerr BHs. In particular, if only the BH spin is present (i.e., without orbital motion) we find a luminosity consistent with the Blandford-Znajek mechanism within the flux-tube connecting both objects; but no Poynting fluxes outside the near-zone. Thus, the orbital motion, with its non-trivial plasma effects close to the BH, provides the mechanism allowing the EM energy to escape to large distances.
A prograde BH spin -with respect to the orbit-tends to increase the dynamical effects produced by the orbital motion, while a retrograde spin acts in the opposite way. It is found, from Fig. 9, that the prograde spin is capable of enhancing the luminosity (up to a factor of ∼ 2), whereas a retrogade BH spin slightly (by ∼ 20%) diminishes its value. This is again in contrast with the expectations from the UI model (see, e.g., [19]), in which the trend of L UI /L MD with BH spin is very different from the one found here.
In order to understand this behavior, we plot in Fig. 10 the toroidal magnetic field -as defined from the two Killing vectors ∂ t and ∂ φ in Kerr spacetime-, We find that the extra gravitational frame-dragging effect contributes to further twisting of the magnetic field lines near the BH with the prograde spin, leading to stronger toroidal discontinuities that reconnect in a region closer to the BH horizon. On the other hand, the toroidal magnetic field is reduced in the retrograde case and the reconnection site is pushed farther away from the BH, even outside of the ergoregion (see bottom panel of Fig. 10).
For a high prograde spin of a 0.9M BH , the stable circular orbits are allowed down to 2.3M BH . For such orbits, L MD reaches the order of 10 42 erg/s for plausible values of µ, M BH , and M NS , as we showed in Eq. (41). Thus, the EM luminosity L could exceed 10 44 erg/s for plausible parameters of BHNSs with a high BH spin. For the nearly extreme spin case of a = 0.99M BH , L is likely to reach the order of 10 46 erg/s at the innermost stable circular orbit.

C. Misalignment effects
We now consider situations in which the magnetic moment of the NS is not aligned with the orbital angular momentum, and both compact objects are non-spinning. For such irrotational misaligned binaries, the magnetosphere acquires a phase dependency along the orbit, since the position of the BH relative to the inclined dipole of the NS makes an important difference in the structure of the magnetosphere. Thus, in contrast to the quasistationary aligned configurations, the misaligned solutions are instead quasi-periodical.
The phase dependency is illustrated by the sequence of snapshots presented in the top panel of Fig. 11 for a χ = 15°misaligned magnetosphere, showing parallel currents and projected magnetic field lines on the co-plane. Analogous plots (bottom panel) are also presented for the solutions at 3 orbits for other inclination angles of χ = {30°, 60°, 90°}. Although these solutions are intrinsically three-dimensional and very complex in structure, the co-plane view of Fig. 11 allows us to appreciate their most relevant aspects: (i) there are strong CSs wriggling about the dipole's equators, similar to those seen in [28] (cf. their figure 11) and also reminiscent to those of misaligned pulsars; (ii) magnetic reconnections take place at these CSs; (iii) there are very strong electric currents threading the BH horizon.
The lack of symmetry across the orbital plane makes the previous X-point reconnections less likely to occur and, instead, a larger fraction of the Poynting flux is carried by large amplitude Alfvén waves propagating along magnetic field lines that are anchored to the NS. There are, thus, fewer and smaller plasmoids being produced. Figure 12 illustrates such disturbances propagating on a χ = 15°misaligned magnetosphere at a particular instant along the binary trajectory 5 . As it can be seen, they resemble the poloidal discontinuities discussed in Sec. III A, which are produced in the vicinity of the BH due to a combination of its intense gravitational field and the orbital motion.
This mechanism would suggest a potentially different EM signal as compared to the aligned settings. For instance, analogous large amplitude Alfvén wave-packages were linked to fast radio burst originating in magnetars, modeled as coherent radio emission that arises from the decay of such disturbances as they move outwards [53].
The EM luminosity, measured on a fixed distant surface from the BH, is displayed in Fig. 13 as a function of time for several misalignment angles χ. It shows again the phase-dependent character of the misaligned solutions, here regarding the extraction of orbital energy by means of the surrounding plasma. It is likely that the enhancement of the luminosity (at specific moments along the orbit) relative to the aligned case arises due to higher strength of the magnetic field lines threading the horizon when the inclination of the dipole points towards the BH. At this orbital separation, we find an increase (at least 6 ) by a factor of ∼ 2 for the orthogonal case (i.e., χ = 90°).

D. Implication to observations
The magnetosphere generated around the BH and NS can be a source for the EM emission in a wide variety of wavelengths (e.g., [8]). We here consider possible EM emissions assuming that the magnetic field strength at the NS pole is B p = 10 12 G, paying primary attention to radio waves.
For radio waves to be emitted from the source, the angular frequency has to be higher than the plasma cutoff frequency [54] 6 We note that there is a gradual increase on the peak values with time, especially for larger misalignment angles. This could be a spurious numerical effect (associated with, e.g., lack of resolution) or it could be indicating a longer relaxation timescale for these configurations.
Our results show that the charge number density is where B is the local magnetic-field strength and f e is a numerical factor of ∼ 1-20: cf. Fig. 5. Thus the value for the cut-off frequency, ν p = ω p /(2π), is MHz.(46) Figure 5 shows that f e is larger than unity for the regions of dense plasma, and hence, radio waves with the frequency ν 1 GHz are not likely to be efficiently emitted from a region of strong magnetic fields with B 10 −3 B p (= 10 9 G for B p = 10 12 G), which is found in the light cylinder of r cΩ −1 o . By contrast, outside the light cylinder, the magneticfield strength becomes weaker than 10 −3 B p . Even in such region, a CS of the spiral pattern illustrated in Fig. 1 is generated. Such region can be a source for the strong radio emission with the frequency less than 1 GHz. Assuming that a fraction eff of the total luminosity goes into the radio emission, the spectral flux density at ν = 1 GHz is written as where we supposed the case in which the BH has a small spin and the binary is in a close orbit, that can result in the luminosity of 10 42 erg/s. However, as already mentioned, the luminosity can be enhanced by two orders of magnitude if the BH spin is high 0.9, and hence, a close circular orbit is allowed. Here, we note that the typical radio emission efficiency of pulsars is 10 −5 -10 −4 (see, e.g., Fig. 10 of [8]). With the Square Kilometer Array, whose sensitivity is ∼ 1 mJy for 0.35-1.05 GHz (SKA1-MID band1) and 0.95-1.76 GHz (SKA1-MID band2) [55,56], radio waves from BHNSs in close orbits with a distance of 200 Mpc may be detected in the future, in particular for the binaries of rapidly spinning BHs. As indicated by the recent observational results [57,58], strong hard X-rays and gamma-rays are likely to be emitted from the magnetosphere of BHNSs in close orbits with a high efficiency. However, for the extragalactic sources with the distance of D 100 Mpc, the detection for them is not very optimistic even if the EM luminosity is as high as 10 44 erg/s (e.g., see the estimation in [8]). Even if the BH spin is moderately high, e.g., ∼ 0.9, for which the total maximum EM luminosity can be ∼ 10 44 (B p /10 12 G) 2 erg/s, the detection by the telescopes currently working such as Swift Alert Telescope and Fermi Gamma-ray burst Monitor is possible only for the case that the magnetic field strength at the NS pole is very high 10 13 G and the emission efficiency in X-rays and gamma-rays is close to unity.

IV. CONCLUSIONS
In this paper, we have numerically studied the common magnetosphere of an NS in circular orbits around a BH. We constructed the magnetosphere assuming the NS to be well approximated by a perfectly conducting surface endowed with a dipole magnetic field, moving through an FF plasma on a fixed Kerr background. Although the spacetime geometry is only partially described this way, we expect it to be a good approximation that captures the main phenomena linked to the curvature near the BH. Such simplified approach has allowed us to analyze the magnetospheric response of these binary systems in great detail and to explore several relevant configurations and parameters.
We found a few generic features of these magnetospheres, present in all of the configurations explored in this work. The most noticeable one is the development of strong CSs that initiate near the BH horizon and extend to large distances, typically following an spiral arm pattern. There is also an electric circuit forming along the flux-tube that connects both compact objects, which can be recognized as an approximate realization of the UI effect. Additionally, we discovered that the relative motion of the BH respect to the surrounding magnetized plasma produces large twists for a set of magnetic field lines threading the horizon, leading to magnetic reconnections in its vicinity. The details of all these features, of course, depend on the specific setting and, especially, on the orientation of the magnetic and orbital axes. For instance, in the aligned irrotational BHNS binaries studied in Sec. III-A, X-point reconnections take place immediately to one side of the BH horizon, releasing largescale plasmoids that carry away most of the EM energy by the Poynting fluxes. By contrast, the lack of symmetry on misaligned configurations make reconnections less likely and a larger fraction of the EM energy is instead transported outwards by large-amplitude Alfvén disturbances. These two mechanisms, not accounted within the UI model, are responsible for allowing the EM energy to escape from the near-region and be transported to longer distances (where the luminosity is being computed).
All the above properties can be linked to various EM emissions, as it is the case in the context of isolated pulsar magnetospheres. The electric currents flowing in the flux-tube has been associated with high-energy emissions, either thermal X-rays from Joule heating (e.g., [24]) or nonthermal gamma-rays from synchrotron/curvature radiation [19]; also, this region can generate thermal emissions from relativistic particles bombarding and heating of the NS surface (e.g., [19,50]); even fast radio bursts from coherent curvature radiation were proposed to originate on this region in the late inspiral phase of BNS systems [59]. On the other hand, both large-scale plasmoids of aligned settings and the sharp Alfvén disturbances found in misaligned cases, which are associated with the binary motion, can be potential sources of radio emissions (e.g., [26,29,53]). And finally, the spiral CSs are arguably the most promising source of detectable precursor EM signals for BHNS binaries in close orbits.
As estimated in Sec. III D, for the most favorable scenarios with high BH spins 0.9, with which the EM luminosity can reach L ∼ 10 44−46 [B p /10 12 G] 2 erg/s, radio waves could be detected by forthcoming facilities like the Square Kilometer Array at distances of 200 Mpc, assuming typical radio emission efficiencies of pulsars. On the other hand, in the high-energy band, albeit strong, the signals are less likely to be detectable by the current instruments at distances 100 Mpc, unless the magnetic field strength is as high as B p = 10 13 G.