The Optically Thick Rotating Magnetic Wind from a Massive White Dwarf Merger Product. II. Axisymmetric Magnetohydrodynamic Simulations

We numerically construct a series of axisymmetric rotating magnetic wind solutions, aiming at exploring the observation properties of massive white dwarf (WD) merger remnants with a strong magnetic field, a fast spin, and an intense mass loss, as inferred for WD J005311. We investigate the magnetospheric structure and the resultant spin-down torque exerted to the merger remnant with respect to the surface magnetic flux Φ*, spin angular frequency Ω* and the mass-loss rate Ṁ . We confirm that the wind properties for σ≡Φ*2Ω*2/Ṁvesc3≳1 significantly deviate from those of the spherical Parker wind, where v esc is the escape velocity at stellar surface. For such a rotating magnetic wind sequence, we find (i) a quasiperiodic mass eruption triggered by magnetic reconnection along with the equatorial plane and (ii) a scaling relation for the spin-down torque T≈(1/2)×ṀΩ*R*2σ1/4 . We apply our results to discuss the spin-down evolution and wind anisotropy of massive WD merger remnants, the latter of which could be probed by a successive observation of WD J005311 using Chandra.


INTRODUCTION
Consequences of a merger of massive white dwarfs (WDs) are of great astrophysical importance.It may explode as a Type Ia supernova in particular when the binary constitutes of carbon-oxygen WDs with a total mass exceeding the Chandrasekhar limit (Webbink 1984;Iben & Tutukov 1984).Instead, if a super-Chandrasekhar oxygen-neon core is synthesized after the merger, it may collapse into a neutron star (NS) (Nomoto & Iben 1985;Saio & Nomoto 2004).Such a merger induced collapse has gotten attention as a scenario for the formation of peculiar type of neutron stars, e.g., sources of fast radio bursts (e.g., Kashiyama & Murase 2017;Kremer et al. 2021;Kirsten et al. 2022;Lu et al. 2022).
Recently, a candidate for a significantly younger merger product, WD J005311, was fortuitously discovered within an infrared nebula (Gvaramadze et al. 2019).The most remarkable characteristic of this WD is unveiled through optical spectroscopy, revealing an optically-thick wind emanating from it.This wind is enriched with carbon burning ashes and exhibits a remarkable velocity of v ∞ = 16, 000 ± 1, 000 km s −1 , accompanied by a mass loss rate of Ṁ = (3.5 ± 0.6) × 10 −6 M ⊙ yr −1 .While the direct measurement of the central WD's physical properties remains elusive, the presence of such a fast and intense wind strongly suggests that it is a rapidly rotating and strongly magnetized WD, potentially possessing a super-or near-Chandrasekhar mass (Gvaramadze et al. 2019;Kashiyama et al. 2019).
The mass and composition loaded on the WD J005311 wind is likely from the near-surface carbon burning.The launch of such a wind can be triggered by the Kelvin-Helmholtz contraction of the oxygen neon core of the merged WD, that can happen ∼ 1, 000-10, 000 yr after the merger (Schwab et al. 2016;Yao et al. 2023;Wu et al. 2023).The timing can be consistent with the post-merger age of the system estimated based on both the expansion velocity of the surrounding nebula and the ancient records on a historical Galactic SN, SN1181, which happened in the direction of WD J005311 ∼ 850 yr ago and is likely associated with the merger of the progenitor binary (Ritter et al. 2021;Lykou et al. 2022;Ko et al. 2023).
On the other hand, the expansion velocity of the wind observed in WD J005311 significantly surpasses the escape velocity of a WD with a typical mass.This suggests that the wind is either thermally driven, originating from a superor near-Chandrasekhar mass WD, or magnetically driven due to the rapid rotation and strong magnetic field of the WD.In the former case, the wind velocity will be (Parker 1965): while in the latter case, the maximum wind velocity along the equatorial plane is (Weber & Davis 1967;Michel 1969): The wind is so fast that it catches up and clashes into the surrounding supernova ejecta, forming a wind termination shock, which is observed as an inner X-ray nebula (Oskinova et al. 2020;Ko et al. 2023).The X-ray nebula is still in its infancy; given the observed angular size, it is only a few tens of years old (Ko et al. 2023).Subsequent observations may reveal the time variability and anisotropy of the wind, which is generally expected for a rotating magnetic wind but has not been explored in this context.These properties of the wind can also be linked to the mass-loss and spindown rates of the central WD, which are important in determining the fate of the central WD: whether it eventually collapses into a neutron star, and if so, how rapidly rotating and strongly magnetized the neutron star would be.
Here we model a system like WD J005311 by numerically constructing a 2D axisymmetric wind solution driven by rotating dipole, with implementing a wind launching region that mimics the near-surface carbon burning region.We investigate the wind structure together with its time evolution (i.e., how the mass, energy and angular momentum loss rate from the system evolves with time), and the scaling of the spin-down torque with respect to system parameters such as surface magnetic field, rotation frequency and mass loss rate.This paper is organized as follows.We introduce our setup in Sec. 2, including numerical details.In Sec. 3, we show our results on wind structure, time evolution and scaling of spin-down torque.Finally, we discuss several implications and applications on observational results in Sec. 4.

SETUP
We conduct a series of numerical simulations of a rotating magnetic wind from a massive WD merger product with a stable nuclear burning occurring at the near surface region.We first describe the general numerical setup including the governing equations, the Riemann solver, the mesh decomposition, and the boundary conditions in Sec.2.1.We then describe the source term that represents the injection of mass and internal energy at the near surface nuclear burning region in Sec.2.2.Finally, we elaborate on setups related to magnetic fields.

Magnetohydrodynamic (MHD) equations
We numerically integrate ideal MHD equations with central gravity; in the two dimensional spherical coordinate using Athena++1 (Stone et al. 2020).Here S ρ and S e are source terms that we use to mimic the matter and energy injection into the computational domain, which will be described in detail in Sec.2.2; the velocity vector v, magnetic field B, stress tensor T, total energy density ε tot , energy flux s are given as where ρ is the density, p is the pressure, I is the identity dyadic tensor, and ϕ = −GM * /r is the gravitational potential, where G is the gravitational constant, M * is the mass of the central WD.To close Eqs.(3)-( 6), we use the adiabatic equation of state with an index of γ = 4/3.The above ideal MHD equations are scale-free; we use a unit of G = M * = R * = 1 for the numerical calculations, where R * is the radius of the WD.When estimating quantities in a physical unit, we transform to the cgs unit with setting M * = 1 M ⊙ and radius R * = 0.009 R ⊙ .We note that this is consistent with the mass-radius relation of degenerate oxygen neon cores with an angular frequency of Ω * ≲ 0.5 s −1 (Kashiyama et al. 2019).
We use the HLLD approximate Riemann solver for the MHD equations (Miyoshi & Kusano 2005) with the secondorder piecewise linear reconstruction method (PLM).The time integration is carried out by the second-order Runge-Kutta method with Courant-Friedrich-Lewy number of 0.1.The computational domain is resolved with the mesh number of 128 for [0.9, 30] R * in the radial direction and 128 for [0, π] in the polar direction.We employ a non-uniform mesh in the radial direction, where the radial grid size is proportional to the radius.The fiducial value of the grid size ratio ∆r(i + 1)/∆r(i) is 1.02 so that the smallest cell size is 0.05, where i stands for the grid index.
At the outer boundary of the computational domain, we impose the zero-gradient boundary condition for the radial direction and connect the domain across the axes for the polar direction.On the other hand, we impose the zero gradient boundary condition for the inner boundary (r = r in = 0.9R * ), and set the velocity to be compatible with the rigid rotation of the central WD; In this paper, we consider the cases with Ω * = [0.05,0.07, 0.12, 0.16, 0.23, 0.35, 0.46] s −1 , which correpsonds to ∼ 5-50 % of the mass shedding limit.In terms of the inner ghost cell's density and pressure, we carefully prescribe their values to achieve a specific thermally-driven wind mass loss rate (see Sec. 2.2).Initially, we distribute a cold and homogeneous gas throughout the entire computational domain and inject the thermally-driven wind from the designated launching region.As the thermally-driven outflow reaches the outer boundary, we initiate an aligned dipole field at the inner boundary, facilitating the transformation of the wind into a rotating magnetic wind (see Sec. 2.3).

Wind launching region
We initialize our simulation with a cold, homogeneous, isotropic and non-magnetized atmosphere, and set up a "wind launching region"2 with a width of D near the WD surface, where the mass is injected to the computational domain to mimic the mass loading due to the carbon burning around the surface of massive WD merger product.To do that, we implement an isotropic relaxation function for both matter and energy source terms to update density and pressure in wind launching region: where ρ * and p * correspond to the density and pressure at the outer edge of the wind launching region.The actual value of the relaxation timescale τ is chosen to satisfy the condition, where v s, * ≡ γp * /ρ * is the adiabatic sound velocity and v A, * ≡ B 2 * /4πρ * is the Alfvén velocity with B * being the surface magnetic field strength at the equator (see Sec. 2.3).This condition is needed to stably inject mass to the computational domain by suppressing fluctuations associated with hydrodynamic and/or MHD waves in the wind launching region.The above source terms can self-consistently produce a thermal pressure-driven wind with ρ ∝ r −2 , v r ≈ v esc and a stable mass loss in the steady state, where v esc is the surface escape velocity.In this paper, we set the width of the wind launching region as D = 0.6, R * as our fiducial value and check the convergence of our results with respect to the value of D. We use a fixed value of τ , with which Eq. ( 15) is satisfied for the most strongly magnetized case.Then we set ρ * and p * so that the mass loss rate by the thermal pressure-driven wind becomes Ṁ = 10 −6 M ⊙ yr −1 .

Rotating magnetic wind
After the thermal pressure-driven wind settles down, we turn on a dipole magnetic field that is embedded on the rotating stellar surface, with magnetic moment µ ≡ B * R 3 * ẑ aligned with the rotation axis.We use the following vector potential to ensure that the divergence of magnetic field vanishes.We consider the cases with B * = [1.5 × 10 6 , 2.3 × 10 6 , 2.6 × 10 6 , 3.0 × 10 6 ] G, for which β * ≡ v s, * /v A, * = 10 −(2−3) ≪ 1 so that the magnetic pressure dominates in the near surface region.In order to numerically solve the MHD equations in such a low plasma beta gas, we implement the dual energy formalism (see Appendix A).The launched gas will then corotate with the rotating magnetic field, and the magnetic torque, which depends on the resultant magnetospheric structure and the polar angle, can also contribute to the wind acceleration in addition to the thermal pressure gradient.As a result, we expect a rotating magnetic wind to start blowing in an angle dependent manner, and relax to a quasi-steady state when it reaches to the outer boundary.We simulate the rotating magnetic wind for a few 10 × the spin period after turning on the magnetic field.
wind luminosity and spindown torque estimated at the outer boundary.As we show later, the strength of rotating magnetic winds can be characterized by a dimensionless parameter where Φ * ≡ 2π π/2 0 B r r 2 sin θdθ| r=R * is the half hemisphere magnetic flux and v esc = 2GM * /R * is the escape velocity at the WD surface 3 .With using σ, the Michel velocity (Eq. 1) can be described as v M,max ≈ σ 1/3 v esc .Our simulations cover the range of 1 ≲ σ ≲ 500.
Hereafter we take B1.5e6Ω0.23 with σ = 34.3 as the fiducial model, and first show the multi-dimensional structure of the rotating magnetic wind in Sec.3.1.We then investigate the time variability of the system primarily focusing on the impacts of quasi-periodic eruption along with the equatorial plane in Sec.3.2.Finally, we show how the time-averaged spin-down torque scales with system parameters in Sec.3.3.

Anisotropic wind structure
Fig. 1 shows a snapshot of our fiducial model (B1.5e6Ω0.23)after the wind structure reaches a quasi-steady state.As explained in Sec. 2, mass and internal energy are continuously injected into the wind launching region, as indicated by the lightly shaded area around the WD surface.An aligned rotating magnetic dipole is situated within the WD, and the resulting magnetic field lines are represented by the solid lines.Since the plasma beta (top-right panel of Fig. 1) at the WD surface is significantly smaller than unity, the injected gases co-rotate with the magnetic field up to approximately the Alfvén radius, r A , shown with the dotted line; we determine r A from the condition ρ As shown in the bottom-right panel of Fig. 1, the poloidal component of the magnetic field dominates inside the Alfvén radius, maintaining the dipolar structure.In this region, the gases acquire azimuthal velocities due to the magnetic centrifugal force.On the other hand, the gases are also accelerated by the thermal pressure gradient at the outer edge of the wind launching region, causing them to expand radially.As the magnetic field strength decreases more rapidly with radius than the inertia of the expanding gases, the magnetic field structure undergoes modification, and the toroidal component dominates outside the Alfvén radius.
In the quasi-steady state, magnetic fields are fully open in directions away from the equatorial plane (θ ≲ 80 • and θ ≳ 100 • ), where the wind is primarily accelerated by the pressure gradient at the outer edge of the wind launching region and becomes supersonic at r ∼ 2 R * .In Fig. 1, the sonic radius r s is depicted with the dashed line; where we determine r s based on the condition |v(r s )| = c s (r s ).Note that the terminal velocity is comparable to the escape velocity (as shown in the bottom-right panel of Fig. 1), and the azimuthal velocities are at most a few percent of the radial velocities.Therefore, the properties of the wind in these directions are broadly consistent with the nonmagnetized spherical Parker wind, even though the plasma beta at small radii is significantly less than unity.a Surface magnetic field strength (Eq.16); b Spin angular frequency of the white dwarf; c Mass loss rate (Eq.17); d Spin-down luminosity (Eq.18); e Spin-down torque (Eq.19); f Magnetization parameter (Eq.20); † Time averaged values.
In the equatorial direction (80 • ≲ θ ≲ 100 • ), magnetic fields are closed at small radii, forming a corotating magnetosphere.Beyond the last closed loop, the magnetic field lines are open with a predominant toroidal component, having opposite polarities with respect to the equatorial plane.The transition of the magnetic field configuration is mediated by reconnection occurring at around the tip of the last closed loop, or the Y point.As can be observed from the top-right panel of Fig. 1, the plasma beta in this transition region is higher than those along the open magnetic fields, implying that the gas is trapped mainly by magnetic tensions.In this high plasma-beta region sandwiched by low plasma beta regions, gases are pinched and radially accelerated in the reconnection region, eventually become supersonic at around r ∼ 5 R * .
We confirm that the terminal velocity of the fastest portion becomes comparable to the Michel velocity, v M,max ∼ σ 1/3 v esc .The azimuthal velocity at the Y point is 30% of v esc , which roughly corresponds to the corotation velocity at that location.After becoming ballistic, it gradually decreases as ∝ 1/r due to angular momentum conservation.
The latitudinal angle dependence of the wind at the outer boundary is illustrated in Fig. 2, where the gray shaded regions represent the dynamical range of radial velocity (top-left), luminosity (top-right), torque (bottom-left), and mass loss rate (bottom-right).As described in the previous paragraphs, radial velocities in the near-equatorial direction can reach and transiently exceed the Michel velocity (represented by the horizontal dotted line) during eruptions caused by reconnection at the Y-point.Consequently, the wind luminosity, dominated by the radial kinetic term, also exhibits a sharp peak in the equatorial direction.On the other hand, the torque is primarily exerted by corotation with the magnetic fields and reaches its peak slightly off the equatorial plane (θ ∼ 85 • ), corresponding to the edge of the concave shape of the Alfven radius.

Time variability
The acceleration of the rotating magnetic wind in the equatorial direction occurs in a time-variable manner, associated with magnetic reconnection at the Y-point.Consequently, the overall flux of mass, energy, and angular momentum from the central WD can also vary with time.Fig. 3 displays the time evolution of the mass loss rate ( Ṁ ), luminosity (L), and torque (T ) in our fiducial model.All quantities are normalized by their time-averaged values.
The lower panel of Fig. 3, which provides a long-term perspective, reveals a recurrent eruptive behavior.A recurrent cycle consists of pre-eruption, reconnection, post-eruption phases: In the pre-eruption phase, gases injected into the near-equatorial plane become trapped within the closed field lines.Due to the centrifugal force, the gases accumulate at the tip of the last closed loop, resulting in a continuous decrease in plasma beta in that region.When the centrifugal force acting on the accumulated gases exceeds the tension of the closed magnetic fields, the tip of the closed zone starts expanding radially and is subsequently ejected as a plasmoid through reconnection.Such a plasmoid can be observed at r ∼ 4 -5 R * in Fig. 1.Afterwards, the cycle returns to the pre-eruption phase and restores the gases within the closed magnetic field lines.This type of recurrent eruptions has been known as slingshot prominence in the context of magnetically-active rapidly-rotating stars (e.g., Ferreira 2000;Townsend & Owocki 2005;Jardine & Collier Cameron 2019).
The upper panel of Fig. 3 focuses on the fluxes during a rotation period (t = 15-16 [2π/Ω * ]), as also depicted in Fig. 2. In comparison to the pre-and post-eruption phases, represented by the purple and red lines, respectively, the observed fluxes of mass, energy, and angular momentum consistently increase as the erupted plasmoids reach the outer boundary, indicated by the green lines.Notably, the magnetic torque significantly contributes to the overall torque increase during the eruption phase and plays a dominant role in the central WD's spin-down.
We note that the specifics of reconnection dynamics, such as the frequency of recurrent eruptions and the resulting time evolution of mass, energy, and angular momentum fluxes, may be influenced by our numerical parameters, including spatial resolution (which governs numerical resistivity) and the width of the wind launching region.However, we have verified that the time-averaged values of wind velocities, mass loss rate, luminosity, and torque have all reached convergence concerning the spatial resolution in our simulations and the width of the wind launching region (see Appendix B).

Scaling relation of the spin-down torque
Here we consider how the spin-down torque of rotating magnetic wind depends on the system parameters based on our numerical results.As we show in the previous sub-sections, the time-averaged torque is essentially determined by the magnetic torque exerted on the gases at the tip of the last closed field lines, or the Y-point, where the field configuration is still roughly compatible with the rotating dipole.In this case, the (electro)magnetic torque at the Y point can be estimated as where r Y represents the Y-point radius.Eq. ( 21) is based on the analogy with the force-free limit (e.g., Contopoulos & Spitkovsky 2006).In the force-free limit, the last closed field line corresponds to the light cylinder, r Y ≈ r lc = c/Ω * , and the spin-down torque of a rotating dipole is roughly given as Based on our numerical results, the Y-point radius is determined by the balance between the centrifugal force and the magnetic tension force excerted on the gases, which can be described as  where ρ Y is the density, B Y is the magnetic field streangth, and κ Y is the curvature of the magnetic field at the Y point.In our case, the mass injection from the wind launching region is designed to be spherical, allowing us to describe the density at the Y point as For the latter equation, we take into account that the gases at the Y point is quasi-hydrostatic in the radial direction and v r,Y ≈ v esc is satisfied in all cases examined in Table 3.On the other hand, given again that the magnetic field configuration at around the Y point is still roughly compatible with the rotating dipole, the strength and curvature of the magnetic field can be estimated as respectively.By substituting Eqs.(23-25) into Eq.( 22), the Y-point radius can be obtained as Using Eq. ( 19) and the dimensionless parameter σ, the torque of a rotating magnetic wind can be expressed as * Ω 2 * / Ṁ v 3 esc .Each color of the points corresponds to an angular frequency ranging from Ω * = [0.05,0.07, 0.12, 0.16, 0.23, 0.35, 0.46] s −1 .The dashed line represents a fitting formula, T = 0.5 × σ 1/4 .Fig. 4 shows the relation between the dimensionless parameters σ and T .While the above derivation of Eq. ( 27) is crudely approximate, the derived scaling relation is broadly consistent with our simulation results.The data points, regardless of their angular frequencies, can be effectively fitted by a single relation: T = 0.5 × σ 1/4 .With restoring the physical dimensions, we obtain a fitting formula for the time-averaged spin-down torque of the rotating magnetic wind as ) which can be applicable at least to the cases with 1 ≲ σ ≲ 10 3 .

SUMMARY AND DISCUSSION
We have conducted a series of axisymmetric MHD simulations for rapidly rotating and strongly magnetized WDs, taking into account a near-surface carbon burning process as observationally inferred for WD J005311.We systematically investigated the wind anisotropy, time variability, and the spin-down evolution with respect to the dimensionless parameter σ (Eq.20).We have confirmed that a co-rotating magnetosphere forms beyond the wind launching region and inside the Alfvén radius for σ ≳ 1, which leads to an anisotropic wind structure.In the near-equatorial directions there happens recurrent eruptions of plasmoids that are triggered by reconnections near the Y point.These plasmoids are accelerated to a radial velocity compatible with the Michel velocity, while the wind properties remain broadly consistent with the Parker wind away from the equatorial plane.We found a scaling relation for the spin-down torque (Eq.27) that can be consistently explained by the criteria for reconnections to happen around the Y point, based on the numerical facts we obtained.Our results complement previous studies on solar-like stars with relatively slow rotation (e.g., Ud-Doula et al. 2009;Matt et al. 2012;Raives et al. 2023), and can be applied to not only massive WD merger remnants, but also various stellar objects with 1 ≲ σ ≲ 10 3 .
We now discuss implications of our numerical results on the properties of WD J005311.To reproduce the observed maximum wind velocity of v ∞ = 16, 000 ± 1, 000 km s −1 by the rotating magnetic wind near the equatorial plane, and considering the carbon burning to be occurred near the WD surface, the WD paramters are constrained to be M * ∼ 1.1-1.3M ⊙ , B * ∼ (2-5) × 10 7 G, and Ω ∼ 0.2-0.5 s −1 (Kashiyama et al. 2019).Consequently, the dimensionless parameter ranges from σ ∼ 2-3.
• Given the Michel velocity to be v M,max ≈ σ 1/3 v esc , the contrast in radial velocity between the equatorial and polar directions is ∼ σ 1/3 ∼ 1.2-1.4.Such an anisotropic velocity profile could manifest in the optical spectrum.
To identify this signature, a multi-dimensional line transfer calculation based on our optically-thick rotating magnetic wind solution is necessary.
• Assuming the mass injection from the carbon-burning region to be spherical, the time-averaged mass loss is also presumed to be spherical.In other words, the quantity ρv r remains relatively constant concerning the latitudinal angle.Consequently, the difference in wind ram pressure, which is proportional to ρv 2 r , between the equatorial and polar directions is estimated to be ∼ σ 1/3 , roughly within the range of 1.2-1.4.This can result in a nonspherical expansion of the wind termination shock.The wind nebula of WD J005311 has recently been shown to have an extended structure by Chandra (Ko et al. 2023).Continued observations might identify any asymmetry or non-spherical characteristics.
• The reconnection around the Y-point occurs in a time-dependent manner, which makes the wind acceleration and the resultant non-thermal radiation also time-variable.However, given that the Y-point is well within the photosphere in the case of WD J005311 (with r ph ∼ 0.15R ⊙ ) and the light crossing time at the wind termination shock is much longer than the expected reoccurrence time of reconnection, any time variability induced by reconnection may become smeared out and is difficult to detect.This can potentially explain the absence of apparent variabilities in WD J005311.
• Using the scaling relation for the spin-down torque (Eq.28), we can estimate the spin-down timescale of WD J005311 as ) Hence, even if the currently observed wind of WD J005311 is a rotating magnetic one and continues to blow for a Kelvin-Helmholtz timescale of the central WD, which is ∼ 1, 000-10, 000, yr, the spin-down will be negligible.When the carbon burning in the near-surface region ceases, the mass loss rate will significantly decrease, which increases the dimensionless parameter σ.The rotating magnetic wind will then become relativistic and eventually enter the force-free regime without significantly spinning down the WD.In this case, the remnant WD may serve as a non-thermal radiation source, or or the so-called WD pulsar (e.g., Kashiyama et al. 2011).
Finally, we address some caveats in our numerical simulations.We have implemented a simple prescription for the near-surface carbon burning region as source terms (Eqs.13 and 14), referred to as the wind launching region.However, the actual near-surface carbon burning region should be convective, and can be affected by the strong magnetic field.The structure of the convective region, the resulting wind launch, and its chemical composition would also be influenced by the radiative transfer.For accurate multi-wavelength spectrum calculations, it is desirable to conduct a comprehensive radiative MHD simulation that covers from the carbon burning layer to the photosphere radius.Also, we only investigate the aligned rotating dipole magnetic fields in this paper, while a more complicated field configuration such as oblique or off-centered dipole may be realized for the remnant WD system.Finally, the deformation of the central WD due to its rapid rotation and anisotropic carbon burning can alter the observed properties as well.We save the investigations into the above topics for our future work.
YZ is supported by the International Graduate Program for Excellence in Earth-Space Science (IGPEES) at the University of Tokyo.

C. CHANGE OF THE MASS LOSS RATE IN MHD REGIME
As we claimed in Sec. 2, the mass loss rate is controlled to be the same for the pressure driven wind, based on the prescription of our wind launching region.However, for the rotating magnetic wind solutions we obtained, we found that time-averaged mass loss rate in MHD regime is slightly altered by increasing B * and Ω * , as shown in Fig. 6.For all the cases we explored (σ ∼ 10 0−3 ), ṀMHD / ṀHD varies ∼ 2 times at most, which is generally caused by the recurring eruption events (see the peak in the top right panel of Fig. 2).Though minor in the regime we explored, it may get significant when beta further decreases, and the magnetic effects become increasingly important.
a summary of our simulation.A model BxΩy corresponds to the case with B * = x and Ω * = y in the cgs unit.When the rotating magnetic wind becomes quasi-steady, it can be characterized by the mass loss rate Ṁ = 2π π 0 ρv r r 2 sin θdθ,

Figure 1 .
Figure 1.Snapshot of the rotating magnetic wind of B1.5e6Ω0.23 showing the gas density (left top), the plasma beta (right top), the radial velocity normalized by the escape velocity from the WD, and the ratio of the strength of toroidal field over poloidal field (right bottom).The solid, dashed, and dotted lines indicate the poloidal magnetic field lines, the position where the radial velocity of the wind exceeds the adiabatic sound velocity and the Alfvén velocity, respectively.The pale shaded region around the WD surface corresponds to the wind launching region.

Figure 2 .
Figure 2. Latitudinal angle dependence of the rotating magnetic wind of B1.5e6Ω0.23 showing the radial velocity vr (top left panel), the wind luminosity L (top right panel), the torque T (bottom left panel), and the mass loss rate Ṁ (bottom right panel) at the outer boundary.A time sequence during a rotation period (t = 15-16 [2π/Ω * ]) is represented with colors across the gray shaded region, indicating the entire dynamic range during the simulation.The thick purple, green, and red lines highlight the timings of pre-eruption, eruption, and post-eruption, respectively.These timings are marked with vertical dashed lines in the upper panel of Fig. 3.The orange dotted line in the top left panel indicates the Michel velocity (Eq.2).

Figure 3 .
Figure 3.Time evolution of mass loss rate Ṁ , luminosity L, and torque T of the rotating magnetic wind of B1.5e6Ω0.23 estimated at the outer boundary of the computational domain.The quantities are normalized by the time-averaged values.The upper panel displays a close-up view of a rotation period (t = 15-16 [2π/Ω * ]), where the vertical dashed lines indicate the timings of pre-eruption, eruption, and post-eruption highlighted in Fig. 2.

BBFigure 5 .
Figure 5.Time evolution of spin-down torque before and after (a) increasing resolution (b) decreasing the size of wind launching region in our fiducial case (B1.5e6Ω0.23).

Figure 6 .
Figure 6.Parameter dependence of the ratio between time-averaged mass loss rate after turning on the magnetic fields ṀMHD and the value for pressure-driven wind ṀHD on surface magnetic field B * and WD spin frequency Ω * .

Table 1 .
Summary of our simulations for rotating magnetic winds from white dwarfs