Astrophysical jets from boosted compact objects

Ramiro Cayuso, 2, ∗ Federico Carrasco, 3, † Barbara Sbarato, ‡ and Oscar Reula § Facultad de Matemática, Astronomı́a, F́ısica y Computación, Universidad Nacional de Córdoba. Instituto de F́ısica Enrique Gaviola, CONICET. Ciudad Universitaria (5000), Córdoba, Argentina. Perimeter Institute for Theoretical Physics, Waterloo, ON, N2L2Y5, Canada. Departament de F́ısica & IAC3, Universitat de les Illes Balears and Institut d’Estudis Espacials de Catalunya, Palma de Mallorca, Baleares E-07122, Spain. (Dated: May 2, 2019)


I. INTRODUCTION
Enormous amounts of energy, in the form of Poynting winds or highly collimated relativistic jets, are often observed in various astrophysical scenarios. Such energetic phenomena are believed to be powered by compact objects like black holes (BH) and neutron stars (NS), from the interactions with strong and large-scale magnetic fields on their surrounding magnetospheres. In the seminal works of Goldreich & Julian [1] and Blandford & Znajek [2] (describing pulsars and active galactic nuclei, respectively), it was first demonstrated that the vicinity of these spinning compact objects would be filled with a tenuous plasma. In such rarefied environments, the electromagnetic force dominates over particle inertia and leads to a great simplification in the problem, allowing to capture the basic mechanism that taps rotational energy by means of the electromagnetic field. While pulsars admits a classical interpretation as Faraday disks [3], in the black hole scenario the energy is instead extracted in a form of generalized Penrose process (see e.g. [4]) known as the Blandford-Znajek mechanism. This low-inertia limit of relativistic magnetohydrodynamics, referred as forcefree electrodynamics (FFE), has been -since then-widely used to study global properties of neutron stars and black holes magnetospheres, like for instance Refs. [5][6][7][8][9][10][11].
In the force-free approximation, when there is a perturbation of an otherwise constant magnetic field, the dynamics makes these perturbation travel preponderantly along the magnetic field lines, thus carrying energy with them along this direction. In this work, we simulate a couple of astrophysically relevant situations where this happens, which consist on a black hole or a neutron star moving through a plasma-filled region of constant magnetic field. Galactic mergers could provide a likely scenario for the black hole case [12,13], since the resulting circumbinary disk of the merged galaxy will anchor magnetic field lines, some of which traverse the central region where the binary -and eventually the final supermassive black hole-moves. Another example could be a BH-NS binary, in which the black hole would move through the magnetic field of a neutron star [14,15]. In such cases, we expect the black hole to loose some kinetic energy, transforming it by enlarging its mass but also into electromagnetic energy propagated away by the jets. There has been a number of previous numerical studies on this scenario, [10,16,17], which we use as starting point for the present work. All of them analyze the problem from the point of view of the stationary magnetic field, namely, in their numerical grid the black hole moves and creates the jets which carry the energy. The advantage is that they can readily compute an approximate -since their time direction is a not Killing direction for the background geometry-energy flux. It is precisely this absence of a timelike Killing vector field and the correspondingly lack of a conserved positive-definite energy, which permits to have energy transport via jets in this approximation where the background is fixed. The disadvantage is that they can not model high speed black holes for they move too quickly outside the grid. In our case we choose to describe the problem from the black hole static geometry. The black hole sees a boosted magnetic field and the corresponding electric field, the interaction of its geometry with that electromagnetic configuration generates a stationary solution which takes energy away through jets. In our case we do have a background Killing vector field, and so conservation of the energy it defines, but this is not the energy an observer for which the (uniform) magnetic field is at rest would see. Thus, we also have to define approximate energy fluxes corresponding to these -for our description-moving observers. The energies so defined are transported away, as expected.
The other situation we model is that of a neutron star, also moving on a region of uniformly magnetized plasma. This could happen if a neutron star orbits near an active supermassive black hole, where both strong magnetic fields and force-free plasma are expected around the central region. It could also be relevant in the context of electromagnetic precursor signals from neutron stars mergers [18][19][20], the likely progenitors of gammaray burst. We consider here an idealized setting where the neutron star is represented by a perfectly conducting spherical surface and there is no field generated at the stellar interior. This might be regarded as the limiting case in which the exterior magnetic field is much stronger than the one associated to the star, so that the later can be neglected. We defer the inclusion of the star's own magnetic field and rotation to a more detailed analysis on a future work. A similar behavior to the boosted (nonspinning) black hole case is found, although the details of the operating mechanisms are not the same. Here, kinetic energy from the motion is transformed into Alfvén waves, sourced by the boundary conditions at the conductor. One nice aspect of this problem is that it allows to take the flat spacetime limit, in which the boosted time direction is also a Killing direction. This means there is no ambiguity on defining the energy fluxes used for the description; and thus, it might help gaining some insight into the previous -more delicate-black hole scenario.
In section II, we describe the setup for both cases: the numerical scheme, geometry and evolution equations; the initial data, the boundary conditions and, finally, the energy fluxes definitions. With all of these information one should be able to reproduce our results unambiguously. Except for the boundary conditions and energy fluxes, the setting is very similar to the one in [16,17]. In section III we present the results of our simulations, where different aspects of the problem were explored. Conclusions and perspectives are drawn on Section IV. Throughout, we adopt geometrized units in which c = G = 1 and Lorentz-Heaviside units for the electromagnetic field.

II. SETUP
We are interested on modeling the magnetosphere of a compact object (BH or NS) that travels across a uniform magnetic field by solving the equations of force-free electrodynamics. The code used here was first described in [21] for black holes and later extended in [22], by developing appropriate boundary conditions for the perfectly conducting surface of a neutron star. Since we adopt the reference frame of the central object, its motion relative to the uniform magnetic field will be accomplished through suitable boundary conditions at the external surface of the domain. We shall look for stationary solutions obtained by evolving the fields until they do not change appreciably. The resulting state is determined only by boundary conditions, the background geometry, and to some extent on the handling of the electric field growth on the current sheets that the dynamics generates.
Although a detailed description of our numerical implementation can be found on previous works [21,22], we shall start this section by briefly summarizing its basic features along with information about the metric and the set of evolution equations employed. Then, we shall describe initial data and boundary conditions for the two scenarios we want to study. And finally, we shall discuss the energy fluxes definitions used to analyze the results.

A. Numerical Implementation
We evolve a particular version of force-free electrodynamics derived at [23], which has some improved properties in terms of well posedness and involves the full forcefree current density 1 . More concretely, we shall consider the evolution system given by Eqs. (8)-(9)-(10) in [22]. The numerical scheme to solve these equations is based on the multi-block approach [26][27][28][29], in which the numerical domain is built from several non-overlapping grids where only grid-points at their boundaries are sheared. The equations are discretized at each individual subdomain by using difference operators constructed to satisfy summation by parts (SBP). In particular, we employ difference operators which are sixth-order accurate on the interior and third-order at the boundaries. Numerical dissipation is incorporated through the use of adapted Kreiss-Oliger operators. These compatible difference and dissipation operators were both taken from Ref. [30]. Penalty terms [27][28][29] are added to the evolution equations at boundary points. These terms penalize possible mismatches between the different values the characteristic fields take at the interfaces, providing a consistent way of communicate information between the different blocks: essentially, the outgoing characteristic modes of one grid are matched onto the ingoing modes of its neighboring grids. At each subdomain, it is possible to find a semidiscrete energy defined by both a symmetrizer of the system at the continuum and a discrete scalar product (with respect to which SBP holds). The summation by parts property of the operators allows one to obtain an energy estimate, up to outer boundary and interface terms left after SBP. The penalties are constructed so that they make a contribution to the energy estimate which cancels inconvenient interface terms, thus providing an energy estimate which covers the whole integration region across grids. Such semi-discrete energy estimates -provided an appropriate time integrator is chosen-guarantee the stability of the numerical scheme [27-29, 31, 32]. A classical fourth order Runge-Kutta method is used for time integration in our code.
We use a particular multiple patch infrastructure that has been equipped with the Kerr metric, as in Ref. [26]. This provides a numerical domain that is perfectly adapted to the geometry of the problems, having two global inner/outer boundaries with spherical topology 2 . 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 . For visual representation, throughout this article, we will present our results in the Cartesian coordinates {t, x, y, z} associated with the flat part of the metric 3 . In these coordinates, the metric function H takes the form H = 2M r r 2 + a 2 z 2 /r 2 (1) and the co-vector a reads a = 1, Typically we solve in a region between an interior sphere whose radio is either inside the black hole or represents a perfectly conducting boundary, and an exterior spherical surface at r = 162M . This region is covered by a total of 9 × 6 grids, being 9 the number of layers. The typical resolution used on each of these grids is of 41 × 41 grid-points in the angular directions and N r = 101 points for the radial one. The grids layers do not cover regions of identical radial extension, having more resolution near the inner boundary than in the asymptotic region: from layer to layer we decrease the effective radial resolution by a factor 1.3. In some cases, we have increased the resolution of the individual grids to 61 × 61 and N r = 151.
In order to handle current sheets, we use a rather standard approach in which electric field is effectively dissipated to maintain the condition that the plasma is magnetically dominated (i.e. B 2 − E 2 > 0), discussed in [21]. At the current sheets the magnetic field presents a jump discontinuity in the vertical direction in the y component while the electric field has a spike in the x component. When high order finite difference operators are used in such discontinuous regions, the fields behave in an unsatisfactory manner. Indeed, high order operators tend to give a noisier results. To overcome this issue, the precision of the finite difference operators is reduced from 6th to 2nd order for those grids covering the region where current sheets form. Thus, providing a substantial improvement in the quality of the numerical approximation.
B. Initial Data

Boosted Black Hole
We consider a black hole moving with velocity v o respect to a reference frame in which the magnetic field is asymptotically uniform and along the z-axis, while the electric field vanishes. The boost direction is not necessarily orthogonal to this magnetic field, so we shall study the evolution of the system for different alignments between the black hole velocity and the asymptotic magnetic field. Since we adopt the reference frame in which the black hole is at rest, the electric and magnetic fields should arise from a Lorentz transformation of the electromagnetic field from the frame in which the magnetic field is uniform and the electric field vanishes. That is, where γ = 1 √ 1−v 2 is the Lorentz factor. In Kerr-Schild Cartesian coordinates, the uniform magnetic field would read while the velocity can be generically written as being χ the angle among the velocity and the y-axis.
We emphasize the fact that the steady state solutions would only depend on the boundary conditions, namely on the asymptotic boosted fields, and not on the particular way we have chosen to set the initial configuration in the interior.

Boosted perfectly conducting sphere
We shall consider a perfectly conducting surface as an idealized model of a neutron star. In the present work, we have assumed that the interior magnetic field of the star is several orders of magnitude weaker than the exterior uniform magnetic field. Thus, we propose an initial data that comes from the configuration of an asymptotically uniform magnetic field around a superconducting sphere, which ensures a vanishing normal component of the magnetic field at the stellar surface. Concretely, the initial condition is given by where R is the radius of the star. Hence the field is tangent to the stellar surface and asymptotically matches a uniform magnetic field along the z-axis, as we wanted. Notice, however, that it does not satisfies the magnetically dominated plasma condition at the poles where the magnetic field vanishes. Since the force-free equations would break down at these points, we simply chose a grid that does not contain them, which has proven to be enough for achieving well-behaved solutions.

C. Boundary Conditions
As mentioned before, our numerical domain is bounded by two inner/outer spherical surfaces, where boundary conditions needs to be specified. The treatment given in the present article to these boundaries was previously described and employed in Refs. [21,22]. However, we find important to briefly summarize here the main aspects 4 and explain how these boundary conditions are applied in this new astrophysical context. Generally speaking, physical conditions are imposed by fixing appropriately all the incoming characteristic (physical) modes via the penalty method. Whereas for the constraints, in this case the divergence-free condition ∇ · B = 0, we adopt a method presented in [21] (see also [34]) which restricts possible incoming violations at both boundary surfaces.
Our implementation of the outer boundary conditions consist on setting the incoming physical modes according to a fixed source U ext = ( E ext , B ext ) that we control. This idea appear motivated on the interfaces treatment and was already employed in Ref. [21], where U ext represented a uniform magnetic field sourced by a distant accretion disk. We shall use this strategy again here, with U ext now being the boosted electromagnetic configurations that threatens the compact object magnetosphere. Thus, for the black hole case, U ext is given by the boosted uniform magnetic field (see eqns. (5)-(8)). While for the neutron star, the source is given by the boost of its initial configuration (9) evaluated at the boundary. In this later case, we shall introduce the boost smoothly to its final boost velocity v 0 (at time t f ) by using a time-dependent The treatment of the inner boundary is, on the other hand, very different between the black hole and neutron star scenarios. In the first case, the inner edge of the domain is simply placed inside the black hole horizon where all characteristic modes points inward (i.e. they are all outgoing from our numerical domain perspective), and hence, there are no incoming modes to be prescribed. In the neutron star case, the inner edge of our domain is placed at the stellar surface which is assumed to behave as an idealized perfect conductor. In the present work, we have further assumed that the interior magnetic field of the star is negligibly small with respect to the one of the external magnetized plasma and also that the star is not rotating. Thus, the boundary conditions reduces to with α, β i and h ij being the lapse function, shift vector and the intrinsic metric on the spatial slices, respectively. The normal magnetic field is keep fixed to zero by enforcing it at each Runge-Kutta substep, as done in Ref. [35]. While the electric field components are imposed through the penalty method, by fixing the incoming physical modes to a rather involved combination of outgoing modes. We refer the interest reader to [22] (in particular Sec. II-C and the Appendix) for further details on this.

D. Fluxes
In this section we shall discuss the relevance of different quantities needed to describe the jets and the energy extraction process. Before computing any quantity it is important to recall that the electromagnetic field does not by itself define any four-momentum, a fourvector and the energy-momentum tensor are necessary to build this quantity. The different choices of that fourvector define different four-momenta that can be thought to be related to different observers. Thus, the result obtained by computing the electromagnetic flux (as an spatial component of the four-momentum) in the BH's frame will certainly be different from the one obtained in the plasma frame of reference. The electromagnetic flux expression for the BH's frame comes naturally, it is associated to the time-like Killing vector field of the background spacetime geometry that allows to construct a conserved four-momentum current. To obtain fluxes this vector is contracted to the normal to some spacelike hypersurface, usually the boundary of a sphere at some radius r = R in Kerr-Schild coordinates. We will refer to this quantity as Poynting flux, defined as follows, where p a is the four-momentum defined by, where T ab is the Energy-momentum tensor and k a is the Killing vector field related to stationarity (see for instance Appendix B of [21]). The factor √ −g appears from the normalization at the surface r = R. Now the task is to find a quantity that is representative of the electromagnetic energy flux in the reference frame in which the BH is not static, namely the frame where the asymptotic magnetic field is constant and at rest. To represent an observer that moves relative to the BH with velocity v a , we can take its four-velocity n a to be, One can thus define the four-momentum p a for this observer as, where n a is the normal vector to the equal-time hypersurface Σ t of the space-time foliation.
Since the background geometry describes a curved spacetime, this boosted quantity is rather arbitrary: it does not have the same normalization as the Killing vector field nor gives a conserved current. Nevertheless it acquires meaning as an asymptotic quantity, far away from the BH we can use the fact that the geometry is approximately flat there and so the boosted time direction would approach a boosted Killing vector of the underlying asymptotic geometry. Thus, for large distances we can define the four-momentum of an observer that is in the plasma rest frame (i.e. with velocity v a = −v a ) as in (15). Now using this four-momentum we can define a quantity that is representative of the electromagnetic flux in this frame, we can do so by contracting p a with a vector N a that is both of norm unity and orthogonal to n a . By proposing an Ansatz of the form N a = (N a + ζv a + ξk a ) (where k a is the Killing vector field associated with stationarity and N a is the radial unit vector) and using the conditions: we can obtain the value of the constants ξ and ζ as, where q = (−α 2 + β 2 ), v β = v i β i and α & β are respectively the Lapse and Shift of the spacetime foliation. It can be easily checked that for large values of r (i.e. where α → 1 and β → 0) the vector N a is simply the Lorentz transformation of the N a vector, and so the surface determined by the vectors N a is asymptotically a boosted sphere. Now we can finally define the electromagnetic energy flux for the boosted frame as where h is the determinant of induced metric to the surface in which this flux is computed.
It is very important to stress that, even though the expression (19) can be evaluated in the whole numerical domain, it's value will only be representative of the electromagnetic flux far away from the BH, where n a approaches an asymptotic Killing vector field and consequently the p a is an asymptotically conserved fourmomentum. This downside of not being able to properly define the electromagnetic flux (for the plasma's reference frame) in the whole numerical domain is not a consequence of our choice of the BH rest reference frame. In previous works, [16,17,36] a measurement of the electromagnetic flux is also given accurately only far away from the Black Hole(s) system since the Killing field needed to construct the conserved four momentum is only an asymptotic concept.
Another aspect that has to be taken into account is the fact that actually this system is not isolated, since uniform magnetic and electric fields are present as a background in the whole numerical domain and they would remain so as a consequence of the boundary conditions imposed at the outer boundary. Special care has to be taken in order to distinguish the electromagnetic flow generated by the BH's interaction with the fields, from the ubiquitous flow generated by the background. In order to subtract adequately this background radiation we take the same approach as in [17,36], that is, we subtract the value of the field's initial condition to the stationary values of the electromagnetic field (final configuration) before computing the value of the electromagnetic energy flux. Another interesting approach to this problem is to compute asymptotic fluxes using the plane wave structure of force-free electrodynamics, that is, its fast and Alfvén propagation modes. Mensurable quantities can be constructed for each plane wave mode and study them separately. For instance, we shall look at the radial fluxes Φ ± A (Alfvén modes) and Φ ± T (fast magnetosonic) of the final stationary solutions. Following appendix A of Ref. [22] and assuming our solutions reasonably satisfy the constraint, i.e. E · B ≈ 0, we get: where σ ±

A. Boosted black hole
In this section we shall explore different parameters of the problem such as boost velocity v, black hole spin a and inclination angle χ among the boost direction and the y-axis (i.e. χ represents the departure of the direction of the motion from the case in which it is orthogonal to the asymptotic magnetic field). Lets start first by the case where the asymptotic magnetic field is orthogonal to the boost velocity v a (e.i. χ = 0). We shall study some aspects of the solution for different magnitudes of this velocity, such as the topology of the electric and magnetic field configurations, the power of the electromagnetic flux and the development of a current sheet. Later, we shall study the dependence of the luminosity on the misalignment between the boost direction and the asymptotic magnetic field. And finally, we will analyze the effects that including rotation has on the electromagnetic flux.

Orthogonal velocity
We present a late time configuration for the case of a Schwarzschild black hole, boosted with velocity v = 0.5 orthogonal to the asymptotic magnetic field. The general structure of the solution is depicted on Fig. 1, where it can be seen that magnetic field lines are disturbed by the passage of the black hole, leaving a trail behind it. Similarly, electric field lines, shown at the z = 0 plane, are also dragged by the black hole as it moves along the y-axis, but from the opposite side. The way in which the magnetic field is pulled towards the black hole, result in a discontinuity of the y-component at the equatorial plane, thus producing a current sheet with strong electric fields. In Fig. 2 we have plotted the quantity B 2 −E 2 B 2 , which when close to zero signals the location this current sheet. In such regions, the numerical mechanism that effectively dissipates electric field is actively operating to avoid violations of the magnetic domination condition. We see that for the present case, the current sheet extends behind the black hole's motion, up to approximately 6M . Figure 3 displays the radial electromagnetic energy flux density on the x = 0 plane, as measured in the two relevant reference frames of the problem, namely: the one where the observer is at rest with respect to the asymptotic magnetic field (plasma frame, top image) and the one in which the black hole is at rest (BH frame, bottom image). Both the flux in the plasma frame, i.e. p a N a , and the flux measured in the BH frame, p r , exhibit a pair of highly collimated jets emerging from the black hole. These jets form an angle with the z-axis given by, θ jet = tan −1 (γv), in the co-moving frame; and equivalently, θ jet = tan −1 (v), in the plasma frame. Such misalignment between the collimated energy flux and the original magnetic field orientation is expected and has been reported previously in [17]. At a first glance, we see the main difference between the two fluxes is at their noncollimated components, especially in front of the black hole where it gives negatives values in the BH frame. Figure 4 presents again the radial electromagnetic energy flux density as measured in the BH's frame, but now the initial background fields has not been subtracted from the final stationary configuration. The same pair of jets as in Fig. 3 can be seen, except that they are somewhat hidden now by a mainly dipolar flux density distribution arising from the background electromagnetic field being boosted against the black hole. It is worth emphasizing that integrating this flux around the BH horizon gives a rather small but negative value, thus showing that there is no energy extraction from the black hole. This is consistent with the fact that there is no ergoregion here and, hence, the Blandford-Znajek mechanism is not possible. The net positive flow of electromagnetic energy in the plasma frame, on the other hand, must arise from the available energy due to the relative motion between the magnetized plasma and the black hole. We turn next to a more quantitative analysis and consider, first, how does the emitted jet power changes with the distance from the black hole horizon. Thus, we measure the integrated flux in the plasma frame, Φ , as a function of radius. The integration is performed on a surface determined by the normal vector N a (that in this frame represents spheres), within a cylindrical region of 60M diameter enclosing the collimated jet. Following the notation of Refs. [16,17], we shall compute this quantity in physical units, respect to a representative system in which a black hole of mass M = 10 8 M is immersed on a magnetic field of strength B o = 10 4 G. That is, the results will be expressed proportional to (M 8 B 4 ) 2 ≡ ( M 10 8 M ) 2 ( B 10 4 G ) 2 , allowing for an easy translation to any pair of physical values M and B. Figure 5 presents the behavior of Φ in the range r = 70 − 150M , for a black hole moving at speed v = 0.5. It can be seen that the emitted power drops approximately 20% between r = 70M and r = 150M . A function of the form Φ ∼ Φ ∞ (1 + σr −1 ), also shown in the plot, fits the numerical data very well. The asymptotic value for the collimated flux is Φ ∞ = 3.44 × 10 44 erg/s, and  The results obtained for several boost velocities for both reference frames are summarized on Fig.6. For the plasma frame (top figure), as expected from previous numerical experiments in the regime v ≤ 0.2 [17], further supported on theoretical arguments [37,38] later, the emitted power for non-relativistic speeds goes as ∝ v 2 (see red dashed curve). However, we find that for larger boost velocities the correct dependence is instead given by ∝ γv 2 (solid black line), which fits the numerical values very well for the whole range of velocities explored. To the best of our knowledge, no one has pointed out this behavior before, which may have important observational consequences on astrophysical scenarios were such relativistic speeds are plausible. Meanwhile, from the BH's frame, the behavior with the boost velocity is instead given by ∝ γ 4 v 2 .

Misaligned case
Now, we consider situations in which the exterior magnetic field and the black hole velocity are not orthogonal. In order to do that, we take the boost velocity to lay on the y − z plane and parametrize different orientations by the angle χ it forms with the y-axis (see expression (8)). For this particular scenario, we have focused in the case of a black hole moving at speed v = 0.5.
A representative late time configuration of the magnetic field is shown in Fig. 7, corresponding to a black hole traveling with inclination angle χ = −π/4. Notice that the field topology is similar to the one of the orthogonal case (shown in Fig. 1), but now the asymptotic field is rotated an angle Θ = tan −1 ( B y B z ) = within the y − z plane, as seen from the black hole's frame. Such rotation is induced by the Lorentz transformation and can be straightforwardly computed from equation (5). Figure 8 shows the EM energy flux density p a N a in the x = 0 plane, corresponding to a late time solution of a black hole moving with a direction determined by χ = π/4. In contrast to Fig. 3, we see that -as expectedthe solution has lost the reflection symmetry respect to the z = 0 plane, exhibiting now a pair of asymmetric jets. To study the dependence of the jet power on inclination, we vary the angle χ from 0 to π/2. It can be seen, in Fig. 9, that the power for each individual upper/lower jet highly depends on the inclination angle: they can be up to ≈ 17% higher than in the orthogonal case (χ = 0) and vanishing for χ = π/2. Figure 9 also shows the net collimated power, i.e. the sum of the contributions from both jets, along with a fitting ∝ cos 2 (χ), which fits very well with the numerical data. The same cos 2 (χ) behaviour was observed in [21] for the stationary Kerr BH case, with the exception that χ represents the inclination angle between the asymptotic magnetic field and the BH rotation axis. This behaviour is expected, since the electric field in the BH frame is ∝ cos(χ). By performing simulations for negative inclination angles we observed that, as expected, there is a symmetry between the lower and upper jets, in which the upper jet power for a given angle χ o equals the power of the lower jet at the opposite angle, i.e. −χ o . The simulations performed for different angles χ also show that the direction of each jet is shifted respect to the χ = 0 case, this displacement is shown in Fig.10, which also presents a fitting of the numerical data with functions ∝ cos(χ − δ) (for the up-  per jet) and ∝ cos(χ + δ) (for the lower one). For the boost velocity employed here (i.e. v = 0.5), we find that δ = 0.17 ± 0.01.

Spinning black hole
In this section, we present the results of the simulations performed on a Kerr background, focusing on the effects the black hole rotation has in the emitted power. In Ref. [37], analytic vacuum Maxwell solutions were found for the field configurations in the vicinity of black hole which is both moving and spinning. The estimated luminosities from these solutions have shown that the effect of rotation would be subdominant respect to the one associated with the translation motion. This scenario has also been studied numerically in Ref. [17], within the forcefree approximation. Considering boost velocities up to v = 0.2, the authors of [17] have proposed a decomposition of the total luminosity as, where Φ spin and Φ boost represent the contributions of spin and linear motion, respectively. Thus, suggesting that the two mechanisms acts separately, with Φ spin being independent of the velocity v (only depending on a) and Φ boost being a constant which does not depend on spin. In Fig. 11 we show the luminosities at two different spin values: a = 0 (i.e. non spinning case) and a = 0.6, for velocities up to v = 0.2. The plot reproduces very well the results of Ref. [17] (specifically, their figure 2), which illustrates -through the fitting curves-the above mentioned behavior of the two separate contributions. This serves two purposes: it confirm their results by an independent approach to the problem, and on the other hand, it further validates our numerical implementation. Now, we shall explore what happens if one goes to larger speeds. To that end, we present in Fig. 12 the resulting luminosities at two spin values, a = 0 and a = 0.9, for velocities that ranges from v = 0.1 to v = 0.7. Surprisingly, we find that the curves that represent the non-spinning and the highly-spinning black holes tend to overlap for highly-relativistic boost velocities, v 0.5. It means the decomposition (22) made above no longer holds for such speeds, where Φ spin is seen to decrease (see Fig. 13). This can be explained, at least in part, by noticing that the power on the BZ mechanism diminishes with the inclination angle among the rotation axis and the asymptotic magnetic field. Thus, the angle θ jet (v) produced by the motion would tend to reduce the spin contribution from the total luminosity. Indeed, we have confirmed numerically that by aligning the spin axis to the upper jet (since it can not be aligned simultaneously to both jets) one gets a larger spin contribution Φ spin than compared to the case in which the rotation axis is perpendicular to the motion and, moreover, that this difference increments with boost velocity.

B. Boosted perfectly conducting sphere
We turn now to the idealized setup of a neutron star that is moving across a uniformly magnetized plasma in the orthogonal direction. The star has been modeled by a perfectly conducting spherical surface on a Schwarzschild spacetime, and it was further assumed to have no magnetic field on its own 5 . The star is smoothly brought to relative motion by gradually boosting the initial electromagnetic configuration at the outer boundary. After an initial dynamical transient, the numerical solutions reach a steady state showing collimated electromagnetic jets. We will analyze these solutions and how their jet power vary with the boost velocity v and stellar compactness C ≡ M/R. Of particular interest is the flat spacetime limit, M = 0, for that in such case the existence of the Killing vectors makes the notion of the boost and fluxes well defined everywhere, and thus allows for an interesting comparison with the previous black hole scenario.
In a force-free environment, conductors are shown to act as sources by imposing boundary conditions on the surrounding fields [39]. A similar setting to our star embedded in flat spacetime has been studied in the context of satellites (see e.g. [40]), where the problem is fairly well understood. The motion (ŷ) across a uniform magnetic field (ẑ) induces charge separation along the transverse (x) direction, which is conducted away through the plasma in the form of Alfvén waves. An stationary electric circuit, as the one depicted on Fig. 15, is then established. Such configurations gives rise to a significant damping on the motion of the object, as mechanical energy is converted to Alfvén radiation [40]. In Fig.14 we show the Poynting fluxes produced by the neutron star moving at speed v = 0.5, for the case discussed so far in flat spacetime (top panel), and also when including curvature effects setting a stellar compactness C = 0.2 (bottom panel). The overall qualitative picture is similar among these two cases, and also when comparing with the black hole scenario: there are positive electromagnetic fluxes along the (same) jet directions, where plasma currents are sustained by two counter-oriented twisted bundles of magnetic field lines (see right panel of Fig. 15). Even though the distortions on the magnetic field by the strong curvature of the BH are similar to the ones produced by the NS in flat spacetime, the underlying mechanism operating is quite different. As one might expect, the currents at the black hole horizon does not look like those at the conducting surface of the star. Moreover, there is no current sheet forming behind the star, as the one shown in Fig. 2 for the BH, when there is no curvature (i.e. M = 0). But when the mass of the neutron star is tuned-on (i.e. M = 0), an analogous current sheet emerges and the emitted power gets enhanced. Thus, indicating that a composition of the two effects is acting; namely, the one associated with the perfect conductor condition and the one due to spacetime curvature.
Quantitatively, the dependence of the collimated jet power on the boost speed is again ∝ γv 2 as in the black hole scenario (see Fig. 16). Even though there is no unambiguous way to compare luminosities among a black hole and a perfectly conducting sphere in flat spacetime, we choose to relate the stellar radius R with the BH mass M by setting R = 2M in geo-metric units. This way, we will be comparing the plasma frame luminosity produced by the black hole with the one of a NS whose surface is placed at the Schwarzschild radius of the BH, exploring different stellar compactness at this fixed radius (Fig. 17). We find that the emitted power is now larger by a factor between 1.4 and 3, depending on the compactness. A similar enhanced luminosity for the NS scenario was found for rotating compact objects [41], with the explanation there being traced to distinct effective resistances of their electronic circuits. This argument, reminiscent from the membrane paradigm, indicates that the perfect conductor (placed at the BH horizon) in flat spacetime would produce larger values (∼ 40% in this case) simply because the black hole behaves effectively as a poorer conductor. We observe that increasing the stellar compactness then leads to an interesting interplay between the effects related with the conducting surface and those of gravity. The luminosity quickly rises (up to ∼ 3L BH by C = 0.1) and then smoothly begins to drop, presumably approaching the value L BH at compactness C ≈ 0.5.

C. Decomposing the solutions in physical modes
Here we want to analyze the radial fluxes associated with the physical propagation modes, as defined by equations (20) and (21). To that end, the fluxes are plotted in figure 18, for the three main cases under consideration. Namely, the black hole and the perfectly conducting sphere in flat/Schwazschild spacetime, moving at speed v = 0.5 along the y-axis. Fast magnetosonic modes move at the speed of light without reference to the background magnetic field, thus here they seems to carry the contribution from the relative motion between the magnetized plasma and the object. On the other hand, Alfvén modes show funnels of positive radiation along the jets. But they also show similarly collimated regions of nega- tive (i.e. incoming) flux at the opposite side of the y = 0 plane. This is a bit puzzling, considering there is no electric charge density nor currents at those regions, as opposed to what happens with the positive components (outgoing) along the jets. It might be related, however, with the fact that we are computing these fluxes on the frame of the moving objects and not in the plasma frame; in that frame the net sum of the different contributions to the energy should be approximately zero when integrating on (e.g. spherical) surfaces enclosing the object.

IV. DISCUSSION
In this paper, we have studied jets arising from black holes and neutron stars moving across a magnetized force-free plasma environment. Even though, as one may expect, there is no conserved (Killing) energy being extracted from the moving objects -like e.g. emerging from the BH horizon-, we find however a positive (outgoing) net flux of approximate energy at the plasma frame. Such energy represents the one measured by an observer at rest with respect to the (asymptotic) uniform magnetic field. We see that part of the available energy from the relative motion between the magnetized plasma and the object gets transferred to the electromagnetic field, producing two counter-oriented twisted bundles of magnetic fieldlines which induce stationary currents and support the highly collimated energy fluxes found; namely, the jets.
By working on the co-moving reference frame, we were able to explore a wide range of boost velocities, finding a dependence of the luminosity of the form, L ∝ γ v 2 , in both scenarios. Would there be an astrophysical situation where relative high velocities appear between compact objects and magnetic fields, the gamma factor would become preponderant. We have also explored in some detail the other relevant parameters of the problem, like the orientation of the motion respect to the asymptotic magnetic field or the inclusion of black hole spin. We looked not only to the total energy flux, but also the contributions from each of the different force-free modes, namely Alven and fast modes. This way we were able to better understand the character of the process.
Comparing a black hole with a perfectly conducting sphere on flat spacetime, we have concluded that the overall effects are quite similar, although there are subtle but important differences among the two mechanisms. Clearly, the horizon does not behave as a perfect conductor 6 and, moreover, there is a strong current sheet emerging and playing an important role in the black hole case. We find that the perfect conductor generates about 40% larger luminosity than a black hole -when placed at the BH horizon-, in agreement with the familiar arguments from the membrane paradigm which states a BH would posses an effective finite conductivity. Furthermore, we saw that when the mass of the neutron star is turnedon, a nontrivial superposition of these two mechanisms operates; interestingly, producing larger luminosities at intermediate values of the stellar compactness.
As mentioned before, we left the inclusion of the stellar magnetic field and rotation to a future work, where the interaction among this field and the external one can be explored in detail. This configurations could be used to mimic the presence of a NS companion orbiting in the context of a binary merger and to systematically study precursor electromagnetic signals along the lines of Refs. [18,20].

V. ACKNOWLEDGMENTS
We acknowledge financial support from CONICET, SeCyT-UNC and MinCyT-Argentina. F.C. acknowledge support from the Spanish Ministry of Economy and Competitiveness grants AYA2016-80289-P and AYA2017-82089-ERC (AEI/FEDER, UE). This work used computational resources from Mendieta Cluster (Centro de Computación de Alto Desempeño, Universidad Nacional de Córdoba), Pirayu Cluster (supported by the Agencia Santafesina de Ciencia, Tecnología e Innovación, Gobierno de la Provincia de Santa Fe, Proyecto AC-00010-18) and Centro de Cómputos de Alto Rendimiento (Ce-CAR). Which are all part of the Sistema Nacional de Computación de Alto Desempeño -MinCyT-Argentina.