Modeling Magnetically Channeled Winds in 3D: I. Isothermal Simulations of a Magnetic O Supergiant

Although fewer than 1% of stars end their lives in core-collapse supernovae, massive stars have had a profound influence on the chemical, ionization and star-formation history of the universe. The first very massive stars formed from the pristine H/He gas of the early universe, prior to the formation of the first galaxies. How massive stars evolve over cosmic time depends on understanding the role that binarity, mass transfer, mergers, rotation, metallicity, mass loss, and magnetic fields play on the pre-supernova evolution of individual systems (see Langer, 2012, for a review).

The last two decades have seen the discovery of an exciting new class of stars: magnetic massive stars. Strong, organized magnetic fields on O and early-B stars strongly influence angular momentum evolution, and impact their late-stage evolution prior to core-collapse. The occurrence rate of magnetic fields in O- and early-B stars is $7-10\%$ (Wade et al., 2016; Grunhut et al., 2017; Schöller et al., 2017). That said, the origin of the magnetic fields is not known. Among the massive stars with detectable magnetic fields, their magnetic geometries appear to be stable on timescales of many decades. Single magnetic massive stars rotate more slowly than non-magnetic stars, suggesting long-term magnetic braking (Wade et al., 2016).

In the presence of a strong, large-scale magnetic field the stellar wind is channeled toward the magnetic equator, the magnetosphere is closed at relatively small radii. The wind pulls open the field at high magnetic latitude, and at large radii. For weak fields, the field is unconfined and blown open by the wind. ud-Doula & Owocki (2002a) showed that the effect of the field on the stellar wind depends largely on just a single wind magnetic confinement parameter, $\eta_{\star}$, which characterizes the competition between field and wind outflow energy density.

where the equatorial field strength $B_{\mathrm{eq}}=\frac{1}{2}B_{\mathrm{p}}$ is half the polar field strength, $R_{\star}$ is the stellar radius, and $\dot{M}_{B=0}$ and $v_{\infty}$ are the mass-loss rate and terminal wind speed in the absence of a field.

In an extensive series of papers, 2D MHD simulations were used to understand the dynamics of magnetic wind channeling and the consequent hard X-ray emission seen with Chandra (ud-Doula & Owocki, 2002a; Owocki & ud-Doula, 2004; Gagné et al., 2005; ud-Doula et al., 2006, 2008, 2009, 2013, 2014).

Specifically, ud-Doula et al. (2008) showed how magnetic channeling depends on rotation and magnetic confinement. We hasten to point out that all this early modeling work was done within the context of aligned rotators, where the rotation axis and the dipole axis coincide. Indeed, that is the only case for which models can be built. Gagné et al. (2005) incorporated an energy equation to model the X-rays from $\theta^{1}$ Orionis C, the central star of the Orion Nebula. These early 2D models were able to explain many of the observational properties of the few magnetic O stars that were known at the time. More recently, Petit et al. (2013) compiled an exhaustive list of 64 confirmed magnetic OB stars with $T_{\rm eff}>16$ kK, along with their physical, rotational and magnetic properties.

Many of the observational properties of magnetic massive stars depend on the ratio of the Alfvén and the Kepler corotation radii (Petit et al., 2013). The Alfvén radius $R_{\mathrm{A}}$, which, at the magnetic equator, is the approximate boundary between the closed magnetosphere and the open field, can be expressed in terms of the magnetic confinement parameter (ud-Doula et al., 2008):

The torque of the magnetic field on the wind outflow maintains rigid body rotation within the magnetosphere. The rotation parameter $W$ is the ratio of the rotational and orbital velocities at the photosphere. For near-critical rotation, the star is oblate, but for $W\leq 0.5$ we use a spherical approximation (ud-Doula et al., 2008):

where $\Omega$ is the rotational angular frequency and $W=\Omega/\Omega_{\rm crit}$ is the critical rotation ratio. The outward centrifugal acceleration from rigid-body rotation will exactly balance the inward gravitational acceleration at the Kepler corotation radius.

Many magnetic O stars are slow rotators with approximately kG fields such that the Alfvén radius is smaller than the Kepler radius (Petit et al., 2013; Wade et al., 2016; ud-Doula & Nazé, 2016; Grunhut et al., 2017). In these so-called dynamical magnetospheres (DM), wind material fed into the magnetosphere falls back onto the photosphere within a dynamical time scale. But for the more rapidly rotating B stars, and, notably, the O7.5 III secondary in Plaskett’s star, $R_{\mathrm{K}}<R_{\mathrm{A}}$; wind mass fed into the magnetosphere builds up at the magnetic equator between the Kepler and Alfvén radii. These are the so-called centrifugal magnetospheres (CM) (Petit et al., 2013). Some mechanism must allow for mass and density to leak out of these centrifugal magnetospheres: slowly via diffusion, or sporadically, via centrifugal breakout events, wherein magnetic tension can no longer contain the built-up mass (Shultz et al., 2020). Owocki et al. (2020) show that the variable H$\alpha$ profiles of centrifugal magnetosphere stars strongly favor centrifugal breakout events.

Although 2D simulations have been used to successfully model aligned rotators, magnetically channelled stellar winds are intrinsically 3D. The magnetic dipole axis is nearly always tilted to the rotation axis, with typical obliquity angles between $30\arcdeg$ and $90\arcdeg$ (Wade et al., 2016). A schematic diagram of an oblique magnetic rotator is shown in Figure 1. Daley-Yates et al. (2019) use 3D isothermal MHD simulations performed with the PLUTO code to predict the radio emission from the winds of oblique magnetic rotators. Their simulations show the formation of a two-armed spiral structure.

In this work, our aim is to extend this modeling effort using the Riemann Geomesh MHD code (Balsara et al., 2019; Florinski et al., 2020). Rather than using a conventional spherical $(r,\theta,\phi)$ grid, Riemann Geomesh employs a recursive partitioning of the spherical icosahedron based triangular meshing of the sphere, and a non-linear radial grid to map the volume of the wind, typically out to $10-20R_{\star}$. The resultant mesh maps the surface of the sphere as uniformly as possible, which traditional Cartesian based meshing is unable to accomplish, especially near the rotational poles. The Riemann Geomesh code is therefore uniquely suited to simulate oblique magnetic rotators in a bias free fashion.

The classification of centrifugal magnetosphere (CM) and dynamical magnetosphere (DM) is based on an aligned magnetic rotator, as displayed in Figures 2a and 2c. This system of classifying magnetospheres is especially relevant to aligned magnetic rotators but including a tilt angle complicates this simple classification scheme, as can be seen in Figures 2b and 2d. Figure 2b shows a tilted “dynamical” magnetosphere. Figures 2b and 2d are just conjectured sketches intended to sensitize the reader that the inclusion of three dimensionality introduces a third important parameter into our discussion; i.e. the tilt between the rotation axis and the dipole axis. A major aim of this work is to examine the dynamics of these magnetospheres in 3D to understand the role of magnetic tilt.

In this paper, we present and analyze 3D simulations of magnetically channeled winds around oblique magnetic rotators. The first goal of this paper is to show that the geodesic mesh MHD code can perform accurate simulations with large tilt angles and rotation rates. The second goal of this paper is to examine the dependence of the quasi-steady state mass-loss rate with rotation, magnetic confinement, and magnetic tilt angle, and to test the prediction that the overall mass-loss rate should decrease with the magnetic field, and increase with the rotation rate (ud-Doula et al., 2009) . The third goal of this paper is to use the 3D simulations to look for and study the episodic breakout events which are predicted to occur in centrifugal magnetospheres (ud-Doula et al., 2006; Shultz et al., 2020). The corresponding angular momentum flux is also shown as a function of $\theta,\phi$ at the outer most boundary of the simulation. The fourth goal is to estimate spindown as the star loses mass and angular momentum over its lifetime. This is achieved with the help of mass flux and angular momentum flux obtained at the outer boundary of the simulation. ud-Doula et al. (2008) analyze angular momentum evolution and spindown for aligned magnetic rotators in 2D.

The outline of this paper is as follows: in section 2 we describe the geodesic mesh along with the value it provides for this work. In section 3 we describe the model, the boundary conditions, and the numerical simulations. In section 4 we present the quasi-steady state mass-loss rate for different magnetic tilt angles and rotation rates. In section 5 we show the episodic centrifugal breakout events in 3D, as well as the mass fallback close to the star. And in section 6 we plot angular momentum flux as a function of $\theta,\phi$, and show regions of high angular momentum flux, and estimate characteristic stellar spindown time. Lastly, in section 7 we present the summary and conclusion of this work.

The overarching problem consists of simulating the magnetically channeled, line driven winds around rotating massive stars. Due to the spherical aspect of the problem, it is best to carry out the simulation on a mesh that is uniquely suited for this problem.

One of the innovative aspects of this paper is that we have carried out the simulations using a geodesic mesh, and imposed a globally divergence-free formulation of the MHD equations. The other major advance in the simulation of line driven winds is that the magnetic field is split into two components: one the curl free, background magnetic field ${\bf B_{0}}$ and the other being the time evolving component ${\bf B_{1}}$. This would be reflective of the fact that the primary source of the magnetic field i.e., the star itself, lies outside the simulation domain. The variations in the magnetic field due to the material outflow is modeled with the help of the evolving field ${\bf B_{1}}$. The simulations are carried out in the rotating frame of reference of the star. Consequently, the simulations include the forces that arise due to rotation, which are centrifugal and coriolis forces. The complete description of the MHD equations employed in this manuscript is provided in Appendix A. In Subsection 2.1 we describe the structure of the geodesic mesh. In Subsection 2.2 we explain why it is particularly useful for this problem.

The schematic of the simulation setup is shown in Figure 1. The radius of the stellar surface $R_{\star}$ is taken as the inner boundary and the outer boundary $R_{\rm o}$ is chosen as $15R_{\star}$. The dipolar magnetic field is maintained throughout the simulation domain by means of the background magnetic field ${\bf B_{0}}$. The magnetic flux density is parametrized by $\eta_{\star}$ (see Eq. 1).

The 3D Riemann Geomesh simulations described here span a range of magnetic tilt angles $\zeta=0,45\arcdeg,75\arcdeg$, two rotation rates $W=0,0.5$, and three magnetic confinement parameters, $\eta_{\star}=0,10,50.$ The ten cases are listed in Table 1, with the corresponding Kepler and Alfvén radii, and the corresponding aligned-magnetosphere classification: CM for centrifugal magnetospheres, and DM for dynamical magnetospheres. The outgoing mass-loss rates in the last column are discussed below.

We first note that although we do not include our simulations at $W=0.25$, the behavior of the DMs can be seen in the $W=0$ simulations. Second, we did not consider $W>0.5$ because at very high rotation the photosphere becomes oblate (ud-Doula et al., 2008), and the spherical mesh is no longer appropriate. Third, we note that higher values of $\eta_{\star}$ lead to significantly longer compute times. High magnetic confinement simulations will be the subject of a future paper.

Figure 5 shows the integrated mass loss from the outer boundary of the simulation as a function of time. The figures show the mass-loss rate reaching a quasi-steady value after $100-120$ ksec. The mass-loss spikes seen in all simulations are the initial mass blow out, which occurs as the initially dense inner CAK wind exits the outer boundary.

In the first simulation with no magnetic field and no rotation, a steady state is reached quickly, and, reassuringly, the density, velocity and mass-loss profiles follow those predicted by CAK theory. These can be seen, respectively in Figures 6a, 6b and 6c as solid black lines (simulation) and black open circles (CAK theory). In the second simulation, rapid rotation is accounted for by means of centrifugal and Coriolis forces, as described in Appendix A. These forces act perpendicular to the $z$ rotation axis and tend to increase density and mass-loss rate, but decrease radial velocity, approaching the rotational equator in the $xy$-plane. Averaged over $\theta,\phi$, this can be seen as the yellow solid radial profiles in Figs. 6a-c. The $\sim 15\%$ increase in the mass loss rate, averaged over $\theta,\phi$ at the outer boundary, can be seen in the last column of Table 1, and is similar to the 2D simulation results obtained by ud-Doula et al. (2008). In the absence of a magnetic field, the increase in $\dot{M}$ is accompanied by an average decrease in $v_{\infty}$. This can be understood as a decrease in the escape velocity $v_{\rm esc}$ with increased rotation.

For the magnetic simulations with $\eta_{\star}=10,50$, $W=0,0.5$ and $\zeta=0\arcdeg,45\arcdeg,75\arcdeg$ (see Table 1) the initially dipolar field is stretched open near the magnetic poles, and remains closed near the magnetic equator, out to approximately the Alfvén radius, channeling wind material into the closed magnetosphere, thereby reducing the mass-loss rate at the outer boundary. The reduction in $\dot{M}$ is more pronounced for higher magnetic confinement, consistent with prior 2D results (ud-Doula et al., 2008). Table 1 also shows that, in all but one case ($\eta_{\star}=50$, $\zeta=75\arcdeg$), higher rotation leads to higher $\dot{M}$.

From Table 1 for $\eta_{\star}=10$, we see that as the tilt increases, the mass-loss rate is reduced. This behaviour suggests a dot-product $\cos\zeta$ variation involving the magnetic dipole moment ${\bf m_{B}}$ and the angular velocity ${\bf\Omega}$. That is, the moment arm of the mass outflow is reduced as tilt increases, thereby reducing the net mass-loss rate. The same is not true, however, at higher confinement, where the maximum $\dot{M}$ is achieved for $\zeta=45\arcdeg$.

In fact, the phenomenon depends upon the geometry of the tilt in addition to the competition between magnetic confinement and centrifugal forces. Therefore, a larger parameter study will be required to yield meaningful predictions. The complexity of the situation can be seen for $\eta_{\star}=50$, $W=0.5$ and $\zeta=45\arcdeg$, where the mass-loss rate is higher than both the $\zeta=0\arcdeg$ and $\zeta=75\arcdeg$ cases. Here we resort to the simulations to explain the situation. At high field $\eta_{\star}=50$, with no tilt $\zeta=0$, and rapid rotation $W=0.5$, much of the outflow close to the rotational equator is trapped in the magnetosphere. Because the centrifugal acceleration is greatest in the equatorial XY plane, $\dot{M}$ is not significantly increased by rapid rotation. This is evident in rows 7 and 8 of Table 1: the mass-loss rates at $W=0$ and $W=0.5$ are comparable. Examining Figure 9, with $W=0.5$, $\eta_{\star}=50$, and $\zeta=45\arcdeg$, most of the outflow near the rotational equator is along open field lines, thereby maximizing $\dot{M}$ and $\dot{J}$, the total rate of angular momentum loss. Figure 9 will be discussed in more detail in section 5. The total angular momentum loss is listed in Table 2 and discussed in section 6 below.

The Figure 6 shows the variation in log density, velocity and mass-loss rate as a function of radius at the end of each simulation. These are averaged over the $\theta,\phi$ to illustrate their radial variation. Outside the Alfvén radii, the mass-loss rate remains constant, while the density decreases rapidly and the velocity increases, as expected from CAK. In general, the presence of magnetic confinement $\eta_{\star}>1$ increases the overall radial velocity, while the density profile undergoes a reduction due to magnetic confinement. In the same way, in the presence of rotation, the density profile and the mass-loss rate increase, as the wind is flung out by rotational forces, consistent with prior 2D simulations (ud-Doula et al., 2008).

For the aligned rotator case, there is some well-developed theory for the mass loss case even when we have rotation and magnetic fields. We compare our simulations for aligned rotators with that theory in this paragraph. In figure 6d, the mass-loss rate for aligned magnetic field cases are presented. The mass-loss rate is normalized with respect to the non-magnetic and no-rotation case ($\dot{M}_{B=0}$). The green and purple circles in the plot are obtained from the analytic formula provided in equations 23 and 24 of ud-Doula et al. (2008) and repeated here as:

where, $R_{c}$ is the confinement radius of the magnetic closed loops. From our simulation results the confinement radius is observed to be of the form, $R_{c}\approx R_{*}+0.6(R_{A}-R_{*})$. We see from figure 6d that the agreement between our simulations and the theory is quite good. In the quasi-steady state, the simulations show the magnetic confinement of the wind and the episodic breakout events at the magnetic equator. This is detailed in section 5 below.

The mass outflow from the surface of the star is free to flow at the magnetic poles and it is confined near the equatorial region by the closed loops of the magnetic dipole field. Due to this, the wind coming from the either side of the equatorial region, guided by the magnetic field lines, meet and cause a field lines breaking at the magnetic equator. This overall behaviour of closed magnetic loops near the surface of the star and open field lines as we move farther can be seen in figures 7.

Figures 7a-f shows the density colorized with the magnetic field lines at different times of the simulation, for the cases of magnetic field strength $\eta_{\star}=10,50$ and no rotation. These figures show the expected dynamics of a magnetically channeled line driven wind. We see that each closed loop has two foot-points that connect to the star at similar moment arms. Owing to the radially outward line driving force, both foot-points have similar amounts of matter that is driven out from the surface of the star. Because Figs. 7a, b, c correspond to a smaller magnetic field than Figs 7d, e, f, we can clearly see that the former set of figures produce a smaller magnetosphere than the latter set of figures.

The matter driven out from the either side of the star meet at the equatorial region, forming a dense knot. The knot is at a much higher density than the rest of the wind, and hence experiences a greater gravitational attraction towards the surface of the star. However, the closed field lines just below the knot prevents the direct fallback of matter on to the star and hence the path of the density knot fallback is dictated by the magnetic field lines. This can be seen in figures 7c, d. These observations are inline with the 2D simulations of ud-Doula & Owocki (2002b). In addition, in this 3D simulation, it can be seen that the density knots at the either end of the magnetic equator fallback from the opposite sides of the field lines, thereby conserving the momentum of the star (figures 7c-f). Also, the overall outflow is guided by the magnetic field lines from the either sides of the pole, in an alternate fashion and hence the subsequent density fallbacks come from the either side of the magnetic equator. This is also inline with the observations made in ud-Doula & Owocki (2002b).

The next set of figures 8a-f, represent the simulations carried out with the high rotation and an inclined magnetic dipole. The magnetic field has the strength of $\eta_{\star}=10$ and the rotation rate is $W=0.5$, the rotation axis is the $z-$axis and the magnetic dipole makes a $\zeta=45\arcdeg$ inclination with the rotation axis. The dynamics are different and interesting in Figure 8, when comparing to that of Figure 7. Due to the tilt and high rotation, the streams of gas that come off from the two foot-points of a closed magnetic loop experiences unequal amount of centrifugal force. The magnetic foot-point that is closer to the rotational equator has more centrifugal force and is flung out faster. The other magnetic foot-point that is closer to the pole, does not have centrifugal assist when being flung outwards. Therefore the streamer that comes off from the foot-point that is closer to the rotational-equator dominates and hence that one dominant streamer is doing the flailing. This is in contrast to the case of aligned rotation, where, there are two streamers at either sides of the magnetic equator that go back and forth.

Similar lines of observations can be made for the figures 9a-f, that corresponds to $\eta_{\star}=50$, $W=0.5$, $\zeta=45\arcdeg$ case. In table 1, this case has $R_{A}=2.9$ and $R_{K}=1.59$ and it has been classified as a centrifugal magnetosphere (CM) based on the traditional explanations that have been developed so far. Upon observing the images at different times in Figure 9a-f, it can be seen that the magnetosphere successively filled and with clumps falling down on the star. This situation is similar to that of Figure 8a-f. Here, one would therefore be inclined to conclude that the simulation in Figure 9a-f is a dynamical magnetosphere (DM) even though it was classified as a CM. But a little reflection allows us to understand that $R_{A}$ should be modified by $cos\zeta$ so that $~{}R_{A}~{}cos\zeta\sim R_{K}$. Here, the tilt has brought the cosine-modified Alfven radius much closer to the Keplarian radius. This might explain two important features in Figure 9a-f. One, the clumps episodically keep filling the magnetosphere and falling back on to the star, similar to that of a dynamical magnetosphere case. Second, the majority of the magnetically streamlined mass outflow comes out from the magnetic foot point that is closer to the rotational-equator (xy-plane).

Figure 10 again shows one episode of the $\eta_{\star}=50$, $\zeta=75\arcdeg$ and $W=0.5$ simulation, as in Figure 9. We see one clump initiated in Figure 10a. Figures 10b-e show successive times when the clump is reined in by the magnetic field and in-falling material fills the magnetosphere. Figure 10f shows the precise moment when the clump impacts the star. Interestingly, at outer radii in Figure 10f we can see the start of another episode of clump initiation.

The episodes of clump formation and fallback can also be seen, when observing the mass-loss rate close to the stellar surface. Referring to Figure 11, the mass-loss rate experiences abrupt drops for the cases, i) $\eta=10,\zeta=0,W=0$ and ii) $\eta=10,\zeta=75,W=0.5$, where the density knot is formed by one side of the stream and that knot proceeds to fall onto star on the other side of the magnetic loop. The other cases do not show such a complete fallback and hence they do not exhibit these sharp drops mass-loss rate close to the stellar surface.

In summary, the tilted dipole simulations shown in Figs. 8, 9 and 10 suggest that the tilt should be factored in when deciding whether we have a DM or CM. Larger simulated data sets will eventually enable us to identify the boundary between DM and CM in cases where the magnetic dipole has a large tilt relative to the rotational axis.

All of the simulation results presented in this work have the rotation axis along the $z-$axis and the magnetic tilt oriented along the $xz-$plane. Owing to this, much of the interesting phenomena are best viewed in the $xz-$plane. In order to illustrate the 3D nature of the simulations, three cross sections are presented in Figure 12, where density and magnetic field lines are presented in the XY, XZ and YZ planes at an intermediate simulation time, for $\eta_{\star}=50$, $W=0.5$ and $\zeta=45\arcdeg$. In the next section, the angular momentum flux is described along with the calculation of mass loss rate and spin down time of the star.

Mass loss from rotating stars carries away angular momentum and leads to stellar spindown. In non-magnetic massive stars with line-driven winds, the mass-loss rate is large, but not enough to significantly increase their rotation periods during the main sequence because the angular momentum is only carried away by gas, and the main sequence lifetime is relatively short (e.g., Maeder & Meynet, 2000). In magnetic stars, most of the angular momentum is lost via the Poynting stress provided by the magnetic field, and not by the outflowing plasma itself (Weber & Davis, 1967; ud-Doula et al., 2009). Weber & Davis (1967) modeled the angular momentum loss of the Sun with the far-field approximation of a magnetic monopole. For stronger dipole fields, ud-Doula et al. (2009) showed that the angular momentum loss rate of an aligned magnetic rotator could be approximated as:

Similar to Eq. (11), the angular momentum loss rate due to magnetic braking is modeled with the help of the Maxwell stress tensor. The corresponding angular momentum flux is defined as,

With this preamble, the remainder of this section describes the angular momentum flux and angular momentum loss rate for the simulations listed in Table 1. Figure 13 displays the total angular momentum flux $dj_{\rm tot}/dt$ in the $xz-$plane for three tilt angles: (a) aligned rotation $\zeta=0\arcdeg$, (b) $\zeta=45\arcdeg$, and (c) $\zeta=75\arcdeg$. Figure 13a illustrates the total angular momentum flux for the case of aligned rotation. Far from the star, the angular momentum loss is expected to be maximal in the $xy-$plane (i.e. perpendicular to the axis of rotation) and minimal along the axis of rotation (i.e. the $z-$axis). This effect is shown in Figure 13a (i.e. in the left panel), where the large scale flow shows that $dj_{\rm tot}/dt$ is small along the rotation axis and substantially larger in the plane perpendicular to the rotation axis. However, in the in the zoomed right panel of Figure 13, we see that the angular momentum loss perpendicular to the plane of rotation is not maximal because closed magnetic loops near the stellar surface inhibit mass loss. Focusing on the right panel of Figure 13a, we see that the maximal angular momentum flux comes from the footpoints of the open magnetic field lines that connect to the equatorial outflow. Once the wind propagates past $R_{\rm A}$ the angular momentum flux is channeled towards the xy-plane.

Thus the angular momentum flux close to the surface of the star is maximal along the open field lines (${\bf B}_{\rm open}$) in the plane of rotation; i.e., where $|{\bf B}_{\rm open}\times{\bf\Omega}|$ is maximal. This result is especially evident for the case of of tilted-rotation with $\zeta=45\arcdeg$, shown in Figure 13b. In the zoomed panel, we see that the angular momentum loss close to the star is low along the magnetic equator, i.e. along the closed field lines, similar to the case of aligned rotation. Comparing the zoomed versions of Figs. 13a and 13b, in the $\zeta=45\arcdeg$ case the open field lines are present on either side of the magnetic equator. As a result of the magnetic tilt, the open field lines along the rotational equator, leading to have larger angular momentum flux close to the star. Conversely, the open field lines along the rotational $z-$axis have smaller angular momentum flux.

Similar observations can also be made for the $\zeta=75\arcdeg$ case (see Figure 13c), where the magnetic poles are nearly perpendicular to the axis of rotation. Hence, the total angular momentum loss close to the stellar surface is maximum near the magnetic poles, where the open field lines and the plane of rotation coincide.

Beyond the Alfvén radius, much of the angular momentum flux emerges along the magnetic equator because the angular momentum flux is channeled by the field lines. (see Figure 14). Figures 14a and 14b, show the total angular momentum flux at the outer boundary of the simulation, i.e. at a radius of $R\approx 15R_{\star}$. Figures 14a and 14b correspond to the case of aligned rotation. We see that most of the angular momentum flux is channeled by the magnetic field lines. Figures 14c-f, show the total angular momentum flux at the outer boundary of the simulation $R\approx 15R_{\star}$, for the case of tilted-rotation with $\zeta=45\arcdeg,75\arcdeg$. Here, similar to the aligned rotation, the angular momentum flux is influenced by the magnetic channeling of the wind and rotation. As we move farther from the star, the effect of rotation overpowers the strength of the magnetic field and the same can be observed from Figure 13b.

The presence of a dipolar magnetic field and its tilt with respect to the rotation axis can have notable influence on the mass outflow and the angular momentum outflow. The mass outflows from the O and B-type stars are very significant and thus influence the overall evolution of the star. Hence, it is inevitable that the presence of magnetic field and its orientation would play a significant role in the evolution of the O and B-type stars. The detailed simulation and analysis of the magnetic O and B-stars are discussed in the literature for the case of aligned rotation ud-Doula & Owocki (2002b); ud-Doula et al. (2008, 2009) with the help of 2D simulations.

For the simulation of spherical systems in full 3D, we have recently developed the Riemann Geomesh MHD code. The code is based on the spherical icosahedron-based meshing of the sphere. The surface of the sphere is mapped as uniformly as possible, which is a significant advance compared to an $r,\theta,\phi$ mesh. This eliminates the increased discretization errors stemming from singularities at the poles and the need for shorter time steps due to the concentration of mesh points at the poles.

In this work, we have carried out the 3D simulations of the magnetically channeled line driven winds for a template O star. For the first time, the simulations are done for a higher rotation rate ($0.5~{}\Omega_{crit}$) and larger magnetic field tilt angles ($0\arcdeg,45\arcdeg,75\arcdeg$); by utilizing the uniform meshing and state-of-the-art MHD algorithms of the Riemann Geomesh code. The simulations reach a quasi-steady state, as expected, and the mass-loss rate is observed to increase with rotation and decrease with increased magnetic field strength, due to the magnetic confinement of the wind. For the magnetized simulations, the wind exhibits the dynamics of magnetic channeling, as observed in previous 2D simulations. In addition, the results for the $45\arcdeg$ case demonstrates the interplay between the centrifugal force, which is maximum along the rotational equator, and the magnetic channeling, which is prominent along the magnetic equator.

Using these insights, we proceed to catalogue the mass-loss rate and angular momentum loss rate for different magnetic-field strengths, tilt angles and rotation rates. The spatial variations in total angular momentum flux ($dj_{\rm tot}/dt$) are examined for different magnetic tilt angles, with the help of images at the cross-section plane, as well as, at the outer surface of the simulation domain. Near the stellar surface, the maximal total angular momentum flux ($dj_{\rm tot}/dt$) is observed where the open magnetic field lines align with the plane of rotation. Farther from the star, the maximal angular momentum flux ($dj_{\rm tot}/dt$) is aligned with the confinement of magnetic field lines, emanating from either sides of the magnetic equator.

The gaseous, magnetic and total angular momentum loss rates ${\dot{J}}$ are computed as a function of radius, and plotted for different tilt angles and magnetic field strengths. With the corresponding mass-loss rates, the characteristic mass-loss time and spin-down time are tabulated for the simulations carried out in this work. The results show a longer mass-loss time with increasing magnetic field strength, which tracks the reduced mass-loss rate. We also find faster spin-down time with increasing magnetic field strength, showing the role of magnetic braking torque on the rotation of the star.

These simulations are the very first that have been reported with a code that can handle all possible tilt angles of the magnetic dipole with respect to the rotation axis. The use of a mapped mesh makes these simulations somewhat more expensive, but this comes with the tremendous advantage of full geometric generality. While we are not restricted to dipolar fields, we have focused on dipolar fields in this study. In order to develop a comprehensive understanding of the relations between, $\eta_{\star},~{}\zeta,~{}W$ and ${\dot{M}},{\dot{J}}$, a much larger set of data with more finer variations in $\eta_{\star},~{}\zeta,~{}W$ is necessary. Therefore, this is planned for the next set of study and it will be elaborated in a subsequent paper.

We acknowledge the extremely valuable inputs from Stanley Owocki. We also acknowledge the computer nodes, which are provided to us on i) Compute clusters at the Notre Dame Center for Research Computing, ii) Bridges-2 supercomputer at the Pittsburgh Supercomputing Center, and iii) Stampede-2 supercomputer at the Texas Advanced Computing Center. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562. Support for this work was provided by the National Aeronautics and Space Administration through Chandra Award Number TM1-22001 issued by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-03060. AuD acknowledges support from Pennsylvania State University Commonwealth Campuses Research Collaboration Development Program. Data was generated through this support from the Institute of Computational and Data Sciences. MG acknowledges summer support from a Provost’s Research Grant, and a work-release grant from the College of Science and Mathematics at West Chester University.

The data underlying this article are available in the article and in its online supplementary material.