Comparisons of MHD Propeller Model with Observations of Cataclysmic Variable AE Aqr

We have developed a numerical MHD model of the propeller candidate star AE Aqr using axisymmetric magneto-hydrodynamic simulations. We suggest that AE Aqr is an intermediate polar-type star, where the magnetic field is relatively weak and some kind of an accretion disc may form around the white dwarf. The star is in the propeller regime, and many of its observational properties are determined by the disc-magnetosphere interaction. Comparisons of the characteristics of the observed versus modelled AE Aqr star show that the model can explain observational properties of AE Aqr if the magnetic field of the star is B\approx (2.8-4.2)x10^5 G. In a representative model, the magnetic field is B\approx 3.3x10^5 G and the time-averaged accretion rate in the disc is 5.5x10^{16} g/s. The star is in the strong propeller regime, where 85 per cent of the disc mass is ejected into conically-shaped winds. The numerical model explains the rapid spin-down of AE Aqr through the outflow of angular momentum from the surface of the star into the wind. The model can explain the low accretion rate onto the star if the radiative efficiency of accretion eta\approx 0.07. The energy budget in the outflows, 2x10^{33} erg/s, is sufficient for explaining the observed flaring radiation in different wavebands. The time scale of ejections into the wind matches the short time scale variability in the light-curves of AE Aqr. Most of the matter is ejected into the wind by the magnetic force, and matter flows into the wind at relatively low velocities, 100-700 km/s. This can explain the predominantly low outflow velocities observed in the tomograms of AE Aqr.


INTRODUCTION
AE Aqr is a nova-like cataclysmic variable (CV) in which a white dwarf star rotates rapidly, with a period of P = 33.08 s (Patterson 1979).It consists of a magnetic white dwarf and a late-type companion star with a spectral type of K3-K5.The binary rotates with a relatively long period of 9.88 hours.It is widely believed that the companion star's atmosphere fills its Roche lobe and matter flows out of the companion's Roche lobe to the white dwarf (Casares et al. 1996).
AE Aqr could be classified as a typical intermediate polar (DQ Her) type CV, except that it has a number of unusual properties: (1) It shows flaring oscillations in optical, ultraviolet and X-ray bands, which are all correlated with each other (e.g., Patterson 1979 and references therein; see, also Fig. 1 from Mauche et al. 2012).(2) Flaring radiation in the radio band shows a non-thermal spectrum typical for electrons radiating in a magnetized plasma (Bastian et al. 1988).Authors compared the radio flares in AE Aqr with those of the microquasar Cen X-3.(3) Pulsations with a period of 33.08 s were detected in optical (Patterson 1979), UV (e.g., Eracleous et al. 1994), soft X-ray (Patterson et al. 1980), and hard X-ray (Kitaguchi et al. 2014) wavebands.(4) The star is spinning down rapidly, at a rate of Ṗ = 5.64 × 10 −14 s s −1 , which corresponds to a high spin-down luminosity of Ėsd ≈ 6 × 10 33 erg/s (De Jager et al. 1994).( 5) The estimated accretion luminosity, Ėacc 10 32 erg/s, is much lower than the spindown luminosity.(5) The Hα Doppler tomograms of AE Aqr are different from the tomograms of other intermediate polars, which typically indicate the presence of an accretion disc (e.g., Marsh et al. 1990).No signatures of a disk were observed in the AE Aqr tomograms (e.g., Wynn et al. 1997;Welsh et al. 1998).Instead, the tomograms indicate that matter flows away from the star.Different models were proposed to explain the observational properties of AE Aqr.In one type of model, it is suggested that the magnetic field of AE Aqr is very large, Bs ≈ 5 × 10 7 G, and the white dwarf spins down due to the magneto-dipole radiation, as in the case of pulsars (e.g., Ikhsanov 1998Ikhsanov , 2006)).In the second type of model, it is suggested that the magnetic field of the star is typical for that of intermediate polars, B ∼ 10 5 − 10 6 G (Warner 1995), and the star spins down due to the interaction of the rapidly-rotating magnetosphere with the surrounding matter in the propeller regime (e.g., Eracleous & Horne 1996).
The apparent lack of an accretion disk in the Doppler tomograms led to the suggestion that an accretion disc does not form; instead, matter flows from the secondary star in a stream of blobs, which interact with the magnetosphere of the WD directly without forming a disc (e.g., King 1993;Wynn & King 1995;Wynn et al. 1997).In their model, blobs of matter interact with the WD ballistically, and a diffusive term has been added to the equations to describe the interaction of the blobs with the magnetosphere of the star.The flaring radiation in different wavebands is explained by the collisions of blobs during their exit path from the white dwarf (Welsh et al. 1998).Meintjes & de Jager (2000) suggested that a stream of matter flowing from the secondary star interacts with the external layers of the magnetosphere, loses its angular momentum and forms a ring or a small disc around the white dwarf.The matter distribution in this ring is strongly inhomogeneous.Recent 3D MHD simulations of matter flow in AE Aqr have shown that a stream of matter flowing from the secondary star collides with itself after the 1st turn around the white dwarf, forming a ring around it that subsequently forms a turbulent disc (Isakova et al. 2016).Next, the matter of the disk is ejected by the magnetosphere of the WD due to the propeller mechanism.In their model, the magnetic field of the WD is large, 50MG, which leads to a rapid ejection of the disk matter by the large, rapidly-rotating magnetosphere.
It is interesting to note that some members of the IP class (SW Sex type) that do show an accretion disk in their quiescent states also show broad, single-peaked flaring phase-shifted emission lines similar to those observed in AE Aqr (Horne 1999).The authors suggested that the propeller mechanism may operate in these stars as well, and possibly in all nova-like IPs during their high-rate accretion states.The authors also concluded that the propeller mechanism is probably connected with the interaction of the magnetosphere with the inner disk, while the flaring emission lines originate in the outflows that form streams above the disc (Horne 1999).In another set of observations, Echevarría et al. (2008) performed high-dispersion absorption-line spectroscopy of AE Aqr and concluded that a weak accretion disk is present, which appears as either a set of clumps or as a more ordered feature.
Overall, it is reasonable to suggest that, in the cases of IP-scale magnetic fields, some kind of disc may form around the white dwarf in AE Aqr because the expected size of the magnetosphere (a few WD radii) is much smaller than the size of the Roche lobe, R0 ≈ 2.2R ≈ 1.5 × 10 11 cm (e.g., Patterson 1979).The disc may have unusual properties compared with the discs around non-propelling stars: it may be small and/or temporary due to the propeller action; it may also have different thermal properties from the standard discs around non-propelling stars.
The propeller regime is expected in AE Aqr, and different observational properties may be connected with this regime.However, the disc-magnetosphere interaction in the propeller regime is a complex phenomenon, and has not been sufficiently studied.It has been investigated in a few theoretical works on spherical accretion (e.g., Illarionov & Sunyaev 1975) and disc accretion (Lovelace et al. 1999).However, only restricted numerical models of the propeller regime have been developed so far.Wang & Robertson (1985) studied the propeller regime in two-dimensional numerical simulations using polar coordinates, and thus studied the processes in the equatorial plane of the disc.They observed that the matter of the inner disc interacts with the magnetosphere of the star due to the magnetic interchange instability (see also Arons & Lea 1976).No outflows were observed due to the twodimensional, polar geometry of their coordinate system.
The propeller regime has also been studied in axisymmetric simulations (Ustyugova et al. 2006;Romanova et al. 2005Romanova et al. , 2009;;Lii et al. 2014).These simulations have shown that the disc-magnetosphere interaction is a strongly non-stationary process, where the inner disc oscillates, most of the matter is ejected into the outflows from the disc-magnetosphere boundary, while a smaller amount of matter accretes onto the star.The WD spins down due to the outflow of angular momentum into the matter-and magnetically-dominated winds.This model can potentially explain the different observational properties of CV AE Aqr.In particular, the model provides: (1) The rate of stellar spin-down; (2) The values of matter fluxes In Sec. 2, we describe our model of AE Aqr.In Sec. 3, we compare our model with the observations.In Sec. 4, we summarize our results.In Appendix A, we provide the details of our numerical model.
We suggest that a small disc forms around the WD.The luminosity associated with accretion onto the stellar surface is low, and is estimated to be Ės 10 32 erg/s (e.g., Mauche 2006).It is connected with the matter flux onto the star as where Ṁs15 = Ṁs/10 15 g/s and η is the efficiency of conversion of the gravitational energy to radiation.The value of η is not well known, particularly in propelling stars, where accretion is expected to be strongly non-stationary.The corresponding accretion rate onto the star is: where Ės32 = Ės/10 32 erg/s is the normalized energy flux.The accretion disc is stopped by the magnetosphere of the star at the radius of rm, where the matter stress in the disk is equal to the magnetic stress in the magnetosphere.This condition can be approximately described by the formula for the Alfvén radius, obtained for non-rotating stars (e.g., Pringle & Rees 1972;Ghosh & Lamb 1978): where µs is the magnetic moment of the WD and Ṁd is the accretion rate in the disc (which, in the propeller regime, is larger than the accretion rate onto the star, Ms); km ∼ 1 is a dimensionless coefficient.B5 = Bs/10 5 G is the normalized magnetic field of the WD, and Ṁd17 = Ṁd /10 17 g/s is the normalized accretion rate in the disc.
If the magnetic field of the WD is typical for that expected in the IPs, then the magnetospheric radius is rm ≈ 3.60Rs in the case of Bs = 10 5 G, and rm ≈ 13.43Rs in the case of Bs = 10 6 G.One can see that, in both cases, the magnetospheric radius is larger than the corotation radius, and the star is in the propeller regime.However, the strengths of the propellers are different.
The strength of a propeller is often measured by the fastness parameter, ωs (e.g.Ghosh 2007), which is the ratio of the angular velocity of the star, Ωs, to the angular velocity at the inner disk, Ω d : ωs = Ωs/Ω d .In the case of a Keplerian disc, Ω d = GMs/r3 m , and the fastness parameter can be re-written in the following form:

Axisymmetric MHD Simulations of AE Aqr
As a base, we use the axisymmetric model of the propeller regime developed in our group (e.g., Romanova et al. 2018).
We start with the initial conditions that are expected in discaccreting propelling stars: the simulation region consists of (a) a star with mass Ms, radius Rs and magnetic field Bs; (b) an accretion disc that is cold and dense, and has an aspect ratio of H/r ≈ 0.15 (H is the initial half-thickness of the disc); (c) a low-density and high-temperature corona.
The disk is turbulent.The turbulence is driven by the magneto-rotational instability (MRI, e.g., Balbus & Hawley 1991;Stone et al. 1996;Armitage 1998;Hawley 2000), which is initiated by a weak poloidal magnetic field placed inside the disk and is subsequently supported by the inward flow of matter with a weak magnetic field from the outer boundary1 .In our model, the accretion rate corresponds to an effective α−parameter of αmag ≈ 0.02 − 0.04.A diffusivity term has been added to the code, with the coefficient of diffusivity constructed in analogy with α−viscosity (Shakura & Sunyaev 1973): η d = α d csH, where cs is the sound speed and α d is a dimensionless parameter.
In our models, we suggest that the 3D instabilities should provide high diffusivity, and we use α d = 1 at radii r < 7Rs, where the disk typically interacts with the magnetosphere, and α d = 0 in the rest of the disc (to avoid suppressing the MRI-driven turbulence).The star rotates with an angular velocity of Ωs such that the magnetosphere rotates more rapidly than the inner disc, and the star is in the propeller regime.
We use a Godunov-type code to solve the MHD equations in cylindrical coordinates.The models are calculated in dimensionless variables.The conversion procedure for dimensionalization and other details of the model are provided in Appendix A3.
To model AE Aqr, we take the corotation radius rcor = 2Rs (see Eq. 1) or rcor = 2 in dimensionless units (see Appendix A3), and calculate all models for this corotation radius.Since the magnetic field of AE Aqr is not known, we model propellers of different magnetospheric radii, rm.To do so, we vary the parameter μ, which we call the dimensionless magnetic moment of the star2 .We have developed numerical models for μ = 30, 60, 100 and obtained different magnetospheric radii (and therefore different fastness values).
We observed that, in all three models, the inner disc strongly oscillates and most of the inner disc matter is redirected into the outflows.The top panels of Fig. 2 show the matter flux density and sample field lines in the model μ60.One can see that the magnetic field lines inflate and open and some of the matter is ejected into conically-shaped outflows.
Matter is ejected at approximately 45 − 50 degrees relative to the rotational axis.Ejections above and below the equatorial plane alternate. 3 The disc-magnetosphere interaction occurs in a cyclic manner, where: (1) The inner disc penetrates diffusively through the external layers of the magnetospehere, (2) The field lines inflate and open after a sufficient amount of matter has accumulated at the inner disc, making the outflows possible; (3) The magnetosphere expands, pushing the inner disc out to larger radii; (4) The inner disc moves inward, diffusively penetrates through the magnetosphere, and the cycle repeats (see analysis in Romanova et al. 2018, also Lii et al. 2014) 4 .
We find the dimensionless magnetospheric radius rm from the simulations using the balance of magnetic and matter stresses, B2 /8π = P + ρṽ 2 , where B, P and ρ are the magnetic field, pressure and density of matter in the equatorial plane, respectively (Romanova et al. 2018).This formula provides the instantaneous value of rm.The middle row of Fig. 2 shows that the inner disc strongly oscillates and the magnetospheric radius, rm, varies with time, which is why  we calculate the time-averaged magnetospheric radius, The dashed line in the middle row of Fig. 2 shows the timeaveraged value rm .We also calculated the time-averaged fastness parameter, Tab. 1 shows the time-averaged magnetospheric radius, rm , the ratio rm /rcor and the fastness parameter ωs for our three models.We obtained that the time-averaged size of the magnetosphere increases with μ, and is rm ≈ 5, 5.9, 6.9 for models μ30, μ60, and μ100, respectively.The fastness parameter also increases systematically: ωs = 3.9, 5.1, 6.5 for the same models, respectively.All three values of ωs are within the interval of fastnesses, 2.3 ωs 16.6, expected for the magnetic fields of IPs5 .
Simulations, performed in dimensionless form, provide us with different values that can be used for the subsequent development of the dimensional model of AE Aqr and for comparisons with the observations.The bottom row of Fig. 2 shows the matter fluxes onto the star, Ṁ s, and into the wind, Ṁ w .One can see that the fluxes strongly oscillate.Matter accretes onto the star in rare bursts, because most of time accretion is blocked by the centrifugal barrier of the rapidlyrotating magnetosphere.Ejections into the wind (blue line) are more frequent than accretion events (red line).To characterize each model, we introduce the time-averaged matter fluxes, which are calculated using a formula similar to Eq. 6.The dashed lines show the time-averaged values of the matter fluxes to the star, Ṁ s , and to the wind, Ṁ w .Simulations show that most of the inner disc matter is ejected into the winds.
The dimensionless models also provide the values of the angular momentum fluxes from the star, Lsd (which are used to calculate the spin-down of the star), and the values of the energy fluxes from the disc to the wind carried by matter, Ėm , and by the magnetic field, Ėw .Tab. 1 shows the corresponding time-averaged values for our models.
Next, we use the data obtained in our dimensionless models to develop the dimensional models of AE Aqr.To obtain a dimensional value A, we take the dimensionless value Ã, obtained in the simulations, and multiply it by the reference value A0 : A = ÃA0.We derive the reference values A0 using our standard procedure described in Appendix A3.

COMPARISONS OF MODEL WITH OBSERVATIONS OF AE AQR
Below, we compare the different values calculated in our model with the values observed in AE Aqr.

Magnetic field of AE Aqr derived from comparisons of modelled and observed spin-down rates
In this section, we compare the spin-down rates obtained in our numerical model with the observed spin-down rate of P sd = 5.64 × 10 −14 s/s, and derive the possible values of the magnetic field.
The angular momentum of the star is Ls = IsΩs, where the value of the moment of inertia of the white dwarf is Is = kMsR 2 s = 1.57× 10 50 k0.2M0.8R 2 7000 gcm 2 .The time-averaged spin-down rate can be estimated from the relation Lsd = Is Ωsd , and can be re-written in the following form: Table 2 shows the spin-down rates obtained in our models.
To derive the magnetic field of AE Aqr, we equate the spin-down rates Ṗsd obtained in the simulations (see Eq. 8) with the observed spin-down rate.We obtain the magnetic field in the following form: Tab. 2 shows the values of the magnetic field, B Psd , which are Bs ∼ (2.8, 3.3, 4.2) × 10 5 G, for μ30, μ60, μ100, respectively.One can see that the strength of the magnetic field slightly increases with μ.It does depend on the moment of inertia coefficient, k0.2, which can vary by a factor of ∼ 2 (depending on the model of the white dwarf).Overall, these values are in the range of the magnetic field values typical for IPs.
We can also calculate the spin-down luminosity of the Values of the magnetic field, B Psd , derived from the comparisons of the spin-down rates, Ṗsd .Ėsd is the calculated value of spin-down luminosity.
modelled AE Aqr.The time-averaged spin-down luminosity of the star is: Substituting in the values of μ60 and Lsd , we obtain spindown luminosities that are approximately the same in all three models (see Tab. 2).These values are comparable to those derived from the observations, Ėsd ≈ 6 × 10 33 erg/s.These comparisons show that the propeller model offers a good explanation for the spin-down properties of the star.However, we cannot yet select one specific model over another, because all of them can explain the observations, although at slightly different values of the stellar magnetic field.

Matter fluxes to the star and to the wind
We can now calculate the dimensional values of the matter fluxes onto the star and into the wind.We take the reference value, Ṁ0, from Eq. A4 and take the dimensionless values from Tab. 1 to obtain the time-averaged matter fluxes to the star (subscript 's') and to the wind (subscript 'w'): Substituting in the values of the magnetic field from Tab. 2, we obtain the values of the matter fluxes (see Tab. 3).The table also shows the time-averaged total matter flux through the disc: Ṁd = Ṁs + Ṁw .One can see that the matter flux in the disc decreases when parameter μ increases.The matter flux in the disc obtained in our models, 4.45× 10 16 g/s Ṁd 1.0 × 10 17 g/s, is comparable to the mass transfer rates expected in nova-like CVs with orbital periods of ∼ 10 hours: ∼ 1.0 × 10 17 g/s (e.g., Dhillon 1996).Therefore, the model is in general agreement with the suggestion that AE Aqr is a nova-like IP.
Tab. 3 also shows the efficiency of the propeller: One can see that propeller efficiency is f eff = 0.82, 0.85, 0.72 for μ = 30, 60, 100, respectively.That is, 72 − 85 per cent of the disc matter is ejected into the wind, while 15% − 28% accretes onto the star.We use Eq. 2 and our values of Ṁs to estimate the luminosity associated with accretion onto the star, Ės .Tab. 3 shows the values of luminosity under the assumption that all matter falling onto the white dwarf is converted to radiation, that is, the efficiency of conversion η = 1.In this case, the luminosity Ės is more than 10 times larger than the accretion luminosity deduced from the observations: Ėacc 10 32 erg/s.
To explain this inconsistency, we can suggest that only a small part of the gravitational energy of the falling matter is converted to radiation.Equating the values of Ės obtained from the simulations with the observed value of 10 32 erg/s, we obtain the values of radiation efficiency: η = 0.033, 0.076, 0.029 for models μ30, μ60 and μ100, respectively.We can also suggest that, in the propeller regime, the processes of accretion and conversion to radiation are different from the processes of quasi-stationary accretion in other IPs: accretion is non-stationary, and the accreting gas may have higher temperatures (due, e.g., to the heating processes at the disc-magnetosphere boundary and/or in the wind).Recent observations of AE Aqr with the NuSTAR and SWIFT telescopes in the hard X-ray band have shown pulsing variability in AE Aqr and a thermal-like spectrum, which has been interpreted as originating in the accretion columns (Kitaguchi et al. 2014;Rodrigues et al. 2017).In their interpretive model, the columns should be tall, taller than the radius of the white dwarf.In this case, the radiation efficiency may be lower than in the cases of quasi-steady accretion observed in nonpropelling IPs.
Another possible explanation for the unusually high accretion rate onto the star may be connected with the axisymmetry of our models and our suggestions with respect to the diffusive layer at the disc-magnetosphere boundary.In our axisymmetric simulations, we suggest that the diffusivity is determined by some 3D instabilities (which we cannot model) and is high, such that the α−coefficient of diffusivity in the code is large: α d = 1.We suggest a high diffusivity in the entire region of r < 7Rs, and place α d = 0 in the rest of the simulation region.It is possible that our diffusivity is higher than that expected in nature, and that our accretion rate onto the star is too high.This issue should be investigated in future work. 6.
On the other hand, we can suggest that our models imply a higher magnetic field of the star, and therefore the corresponding values of fastness, ωs , should be higher.Indeed, our earlier simulations do show that the efficiency of a propeller increases with fastness.Fig. 3 shows the distribution of efficiencies for models with μ = 30, 60, 100 and different values of rcor, resulting in a wide range of fastness values: 1.3 < ωs < 10.2 (from Romanova et al. 2018).One can see that efficiency f eff increases with fastness and approaches 1 at ωs ≈ 14 for models μ = 30, 607 .One could suggest that we should simply take the model with a higher fastness value and obtain a much lower matter flux onto the star.However, the analysis of velocities shows that, in the models with high values of fastness, the velocities in the outflows are too high and may not agree with the distribution of velocities in the tomograms (see analysis of velocities in Sec.3.4).
In summary, we think that the origin of the low luminosity associated with accretion onto the star may either be connected with low radiation efficiency, η, or with a lower rate of matter penetration through the magnetosphere.Future 3D simulations with a proper treatment of diffusivity and heating processes associated with the 3D instabilities may help answer this question.

Energy and angular momentum fluxes in the winds
The wind that flows from the inner disc carries energy and angular momentum out of the star and the disc.The timeaveraged energy fluxes carried by matter (m) and by the magnetic field (f) through the surface S(r = 10, z = ±10) are: Ėm,f erg/s .(13) Table 4 shows these energy fluxes in our models.One can see that, in all models, the amounts of energy carried by matter and by the magnetic field are comparable, although the energy carried by matter dominates slightly in the weaker propeller, μ30, while the energy carried by the field is somewhat larger in the stronger propellers, μ60 and μ1008 .Both the magnetic and centrifugal forces are important in driving the outflows.The magnetic force is particularly important in driving the low-velocity matter into the outflows, and is a dominant driving force in the weaker propellers.In the stronger propellers, the role of the centrifugal force becomes larger in driving the outflows, and the outflows typically have higher velocities.
Energy is ejected into the wind in bursts, and its variation resembles that of the matter flux into the wind (see Fig. 4).These non-stationary ejections may lead to a formation of shock waves and flares of radiation in these shock waves.We suggest that the main flaring radiation observed in AE Aqr may be connected with this mechanism.The total energy flux of Ėtot ≈ 2 × 10 33 erg/s obtained in our models is more than sufficient for explaining the flaring radiation in optical, UV and X-ray bands.The ejected plasma is magnetized, and therefore a part of the energy can be released in the radio band due to the synchrotron mechanism.
We should note that, in Tab. 4, we estimated the timeaveraged energy fluxes.Fig. 4 shows that the energy flux in the individual ejections can be 8-10 times larger than its timeaveraged value, and can be as large as Ėtot ≈ 2 × 10 34 erg/s.
Part of the energy flux carried by the magnetic field flows Table 3.The time-averaged matter fluxes onto the star, Ṁs , into the wind, Ṁw , and in the disc, Ṁd .Ės is the energy flux associated with accretion onto the star and calculated using Eq.2; η is the efficiency of radiation during accretion onto the star. in the axial region in the form of a Poynting flux jet, where a small amount of matter is accelerated rapidly by the magnetic force.The Poynting flux outflows are non-stationary, and particles may be additionally accelerated at the magnetic shocks (e.g., Romanova & Lovelace 1997), potentially up to very high energies.However, it is not clear which part of the energy will propagate in the form of an invisible, magnetic jet and which part will be converted to accelerated particles and radiation9 .

Velocities of matter in the wind and comparisons with tomograms
Fig. 5 shows the tomograms of AE Aqr obtained in the Hα line (from Welsh et al. 1998).The tomograms show that there is a significant amount of matter flowing at velocities of 100-500 km/s (darker regions in the tomograms) and a smaller amount of matter flowing at higher velocities (lighter regions in the tomograms).The tomograms show that only a small amount of matter flows above 1,000 km/s.Below, we analyze whether the distribution of velocities in our models can explain the observed distribution of velocities in the tomograms.
In our propeller model, the matter flowing into the wind has different poloidal velocities.A significant part of this matter is ejected at relatively low velocities, which are much lower than the local escape velocity, vesc = 2GMs/r, where r is the distance from the star.This is because the outflows are mainly magnetically-driven and only occur during the inflation events, when the stellar magnetic field lines threading the inner disc inflate and open.The centrifugal force also contributes to driving the outflows, and becomes the dominant force in very strong propellers.
We calculated the matter fluxes to the wind through the surface S(r = 30, z = ±30) using different criteria for minimum velocity.First, we calculated the total matter flux that satisfies the condition of the poloidal velocity vp being larger than some minimum velocity: vp > vmin, where vmin = k × (0.1vesc), k = 1, 2, 3....20.To calculate the dimensional values of vmin, we took into account the fact that, at rp = 30, vesc = 2/rpv0 ≈ 0.258v0 and v0 = 3, 900km/s, so that in dimensional units the condition for minimum velocity is vmin = k × 100.7km/s,where k = 1, 2, 3, ....20.Fig. 6 (left panel) shows that, in all our models, a significant amount of matter flows at relatively low velocities, vmin < 500 km/s, while a smaller amount of matter flows at high velocities, vmin > 1, 000 km/s.One can see that, in models μ30 and μ60, the amount of matter flowing at vmin 1, 000 km/s is very small, while in the model μ100 there is more matter flowing at velocities vmin 1, 000 km/s.We conclude that the first two models explain the tomograms somewhat better than model μ100.
To refine our analysis, we calculated the matter fluxes in different velocity intervals.Namely, we divided all the velocities into equal intervals of width ∆v k = 0.1vesc = 100.7 km/s, such that the velocities corresponding to the middles of these intervals are: The right panel of Fig. 6 shows the matter flux distributions for the intervals of velocities ∆v k in our models.One can see that, in all models, a significant amount of matter is ejected at low velocities, in the interval of 101-202 km/s (see leftmost point on the plot at v k ≈ 150).At larger values of v k , the matter fluxes decrease.One can see that a smaller, yet still significant amount of matter flows with velocities centered at v k = 250, 350, 450, 550, 650, 750 km/s.In all the models, the matter fluxes are much lower at higher velocities, v k > 800 km/s, and continue to decrease at higher velocities.One can see that, in the model μ100, the amount of matter flowing at high velocities, v k 1, 000 km/s, is larger than in the other two models.The tomograms do not show velocities greater than ∼ 1, 000 km/s, and therefore we believe that models μ30 and μ60 are more favorable than the μ100 model.
In our earlier work on propellers, we calculated the timeaveraged maximum velocities vmax in the outflows from propellers of different strengths, from ωs = 1.2 to ωs = 10.2 (Romanova et al. 2018).We observed that the time-averaged maximum velocities only depend on the fastness parameter ωs, and there is no dependence on the size of the magnetosphere (parameter μ).Fig. 7 shows this dependence.One can see that the time-averaged maximum velocities corresponding to our models of AE Aqr are: vmax ≈ 800 km/s for ωs = 3.9, vmax ≈ 1, 000 km/s for ωs = 5.1, and vmax ≈ 1, 500 km/s for ωs = 6.510 .Looking at Fig. 7, we can see that the models with lower values of ωs (models μ30 and μ60) are preferable, due to their lower maximum velocities, over model μ100, where the maximum velocities are too high.Moreover, the models with even higher values of ωs will provide even higher velocities and are even less desirable.
We can also estimate the number of ejected shells associated with the strong outbursts in the field of view of the tomograms at any given moment in time.The size of our simulation region in the radial direction is r ≈ 34Rs ≈ 2.3 × 10 10 cm, which is 6.5 times smaller than the size of the Roche lobe of the white dwarf, RR ≈ 1.5 × 10 11 cm.The binary separation is A = 3.02R ≈ 2.1 × 10 11 cm.Matter, ejected at, e.g., 500 km/s will propagate to a distance of R outb = 2 × A = 4.2 × 10 11 cm during time t outb ≈ 8400s ≈ 140min.Fig. 2 shows that the outbursts occur every ∆t ≈ 10 minutes.Therefore, the number of outbursts observed simultaneously in this region is N outb = t outb /∆t ≈ 14.These outbursts can explain the inhomogeneous structure observed in the tomograms.
The tomograms also show asymmetry, which may reflect the asymmetric ejections of matter in the propeller regime.In our model of AE Aqr, most of the matter is ejected at approximately Θw = 45 − 50 • about the rotational axis of the disc (see top panels of Fig. 2).This asymmetry may provide the asymmetric picture in the tomograms (e.g., left panel of Fig. 5).

Time intervals between flares
Next, we compare the variability time scales obtained in our numerical models with the time scales observed in the lightcurves of AE Aqr.Fig. 1 shows that there are two major types of flares: some flares occur more rarely, with a time interval of a 1-3 hours between each flare (see individual large flares marked by vertical dashed lines); other flares occur more frequently, on time intervals of 10-30 minutes (see multiple flares on the left-hand side of Fig. 1, and also the subflares that occur on other time intervals).We can compare the observed time scales with the time scales obtained in our models.
Fig. 2 shows the variability obtained in model μ60.The bottom panel shows the matter fluxes onto the star, Ṁs, and into the wind, Ṁw.The blue line in the top panel shows that the strongest ejections into the wind occur every 10-20 minutes.This time scale is similar to the time scales of the frequent flares in the observed light-curve.Ejections into the wind are slightly more frequent (∼ 5 − 10 min) in the model with the smaller magnetosphere, μ30, and are less frequent (∼ 10 − 50 min) in the model with the larger magnetosphere, μ100.Overall, the variability time scales match the short time scale oscillations in the observed light-curve.
In our simulations, ejections into the wind occur in brief bursts, because the outflows are only possible if the field lines connecting the star with the inner disc inflate and open.Inflation becomes possible when a sufficient amount of matter accumulates at the inner disc and diffuses through the field lines of the outer magnetosphere.Most of the time, matter is accumulated at the inner disc, while the outflow events are relatively brief.This is why the variability curve looks spiky (see bottom panel of Fig. 2).However, the small scale oscillations observed in AE Aqr are not spiky.We suggest that a proper inclusion of radiation in our model may change the shape of the variability curve, possibly making the oscillations less spiky.
On the other hand, the variability curve associated with the variation of the inner disc radius, rm (see middle row of Fig. 2), is not spiky.It strongly resembles the observed lightcurve for short time scale flares.In some theoretical models, it is suggested that the flaring radiation is connected with the processes at the disc-magnetosphere boundary and in the inner disc, such as the heating in the turbulent layer of the boundary (Papitto & Torres 2015), or the acceleration of particles in the magnetized plasma of the inner disc (Meintjes & de Jager 2000).If the radiation really does originate at the inner disc, then the oscillations of the inner disc (observed in our model) will lead to the variations in the observed lightcurves.
The origin of the less frequent flares in the light-curves of AE Aqr may be connected with the non-stationary accretion expected in the propeller regime.On the other hand, the disc may be inhomogeneous and non-stationary due to the nonstationary accretion of matter from the secondary star, which may have the form of inhomogeneous streams, or blobs (e.g., Wynn & King 1995).More detailed analysis of these processes is beyond the scope of this paper.

SUMMARY
In this work, we have developed a propeller model of AE Aqr using axisymmetric simulations.We suggested that some type of an accretion disc forms around the white dwarf and interacts with the magnetosphere of the star in the propeller regime.We compared the results of our three models, μ30, μ60 and μ100, with the observations.In these models, the time-averaged magnetospheric radii are rm ≈ 5, 6 and 7, and the fastness parameter values are ωs = 3.9, 5.1 and 6.5, respectively.Our conclusions are the following: 1.All three models can explain the rapid spin-down of AE Aqr, although at slightly-different values of the magnetic field of the white dwarf: Bs ≈ (2.9, 3.3, 4.2) × 10 5 G for μ = 30, 60, 100, respectively.
2. In all models, the disc-magnetosphere interaction is a strongly non-stationary process, where the inner disc oscillates.The total time-averaged accretion rate in the disc is Ṁd = (1.0,0.55, 0.44) × 10 17 g/s in the above three models, respectively.

3.
The models predict that 80, 85 and 72 per cent of the inner disk matter is ejected into the conically-shaped winds, for μ = 30, 60 and 100, respectively.

4.
The main flaring variability can be explained through the processes in the non-stationary outflows (possibly by the radiation in shocks).The time-averaged total energy budget in the outflows, 2 × 10 33 erg/s, is sufficient for explaining the flares observed in different wavebands.

5.
The models predict that 20, 15 or 28 per cent of the disc matter accretes onto a star (in the above models, respectively), which is relatively high.The models can explain the low accretion luminosity of AE Aqr (Eacc 10 32 erg/s) through the process of matter accretion onto the stellar surface, if the efficiency of matter conversion to radiation is η ≈ 0.033, 0.076, 0.029, respectively.The accretion rate onto the star is sensitive to the value of diffusivity at the discmagnetosphere boundary.We accepted a high rate of diffusivity in our model.If the diffusivity is lower, then the accretion rate onto the star will be lower than that obtained in our model.

6.
Most of the outflowing matter has relatively low velocities, v = 100 − 700 km/s, because the outflows are predominantly driven by the magnetic force.This property of the outflows may explain the relatively low velocities of matter, observed in the tomograms of AE Aqr.

7.
Both accretion and ejections are non-stationary and occur in brief episodes.The ejections of matter into the outflows and the oscillations of the inner disc occur on a time scale of ∆t ≈ 10 − 20 minutes.This variability matches the short time scale variability observed in the light-curves of AE Aqr.The longer time scale variability may be connected with the non-stationary accretion from the disc.
Overall, the developed models of AE Aqr do not have any major contradictions with the observational data.One of the inconsistencies is the presence of an accretion disk in our model and no observational evidence of an accretion disk in the Doppler tomograms (e.g., Wynn et al. 1997).We should note that, in our models of the strong propeller regime, the inner disk strongly oscillates.This may lead to a variable, non-ordered disk, which can be difficult to detect using the Doppler tomography technique (which suggested a steady flow of matter in the equatorial plane of the binary, e.g., Echevarría 2012).Our propeller models work even in the cases where the disk is variable, relatively small, or when it forms as a temporary feature.Additionally, in our model, most of the disc matter is ejected into the outflows, which can distort information about the accretion disk in the Doppler tomograms.
for the B φ component.In addition, we adjust the matter velocity vectors to be parallel to the magnetic field vectors.This models the frozen-in condition on the star.Top and bottom boundaries: all variables have free boundary conditions along the top and bottom boundaries.In addition, we implement the outflow boundary conditions on velocity to prohibit matter from flowing back into the simulation region once it leaves.Outer side boundary: the side boundary is divided into a "disk region" (|z| < z disk ) and a "coronal region" (|z| > z disk ), with where Rout is the external simulation radius.Matter along the disk boundary (|z| < z disk ) is allowed to flow inward at a low radial velocity, vr = −δ 3 2 p ρvK (Rout) , δ = 0.02, and a poloidal magnetic field corresponding to the calculated magnetic field at r = Rout.The remaining variables have free boundary conditions.The coronal boundary (|z| > z disk ) has the same boundary conditions as the top and bottom boundaries.

A2 Grid and code description
Grid description The axisymmetric grid is in cylindrical (r, z) coordinates with mesh compression towards the equatorial plane and the z-axis, so that there is a larger number of cells in the disk plane and near the star.In the models presented here, we use a non-uniform grid with 190 × 306 grid cells, corresponding to a grid that is 36 by 66 stellar radii in size.At r = 20, the number of grid cells that cover the disk in the vertical direction is about 60.

Code description
We use a Godunov-type numerical method with a five-wave Riemann solver similar to the HLLD solver developed by Miyoshi & Kusano (2005).The MHD variables are calculated in four states bounded by five MHD discontinuities: the contact discontinuity, two Alfvén waves and two fast magnetosonic waves.Unlike Miyoshi & Kusano (2005), our method solves the equation for entropy instead of the full energy equation.This approximation is valid in the cases (such as ours) where strong shocks are not present.We ensure that the magnetic field is divergence-free by introducing the φcomponent of the magnetic field potential, which is calculated using the constrained transport scheme proposed by Gardiner & Stone (2005).The magnetic field is split into the stellar dipole and the calculated components, B = B dip +B (Tanaka 1994).No viscosity terms have been included in the MHD equations, and hence we only investigate accretion driven by the MRI turbulence.Our code has been extensively tested and has been previously utilized to study different MHD problems (see Koldoba et al. 2016 for tests and some astrophysical examples).

A3 Reference values
We take the reference mass and radius to be M0 = Ms = 0.8M and R0 = Rs = 7000 km, respectively.The reference velocity is the Keplerian orbital velocity at r = R0,

Figure 1 .
Figure 1.The light-curves of cataclysmic variable AE Aqr in different spectral bands.Chandra X-ray (upper histogram); GALEX NUV (blue points); Skinakas & Mount Laguna (black points); CBA & AAVSO (red points) B-band optical; and VLA radio (lower histogram).The times of some of the more prominent flares are marked with dashed blue vertical lines.From Mauche et al. (2012).

Figure 2 .
Figure 2. Top row: A two-dimensional picture of matter flow in the model μ60, shown at several moments in time.The colour background shows the matter flux density and the lines are sample field lines.Middle row: Variation of the inner disc radius, rm.Bottom row: Variation of matter fluxes at the stellar surface, Ṁs, and into the wind, Ṁw.Dashed vertical lines show the moments in time corresponding to the top panels.Dashed lines show the time-averaged values.
values of AE Aqr, the values of the fastness parameter are ωs ≈ 2.3 and ωs ≈ 16.6 for Bs = 10 5 G and Bs = 10 6 G, respectively.These estimates are helpful in restricting the values of fastness in the modelling AE Aqr.
three models of the strong propeller regime, where the corotation radius rcor = 2Rs.Parameter μ determines the final values of the time-averaged magnetospheric radii, rm , and the fastness parameter, ωs.Here, Ṁ s is the time-averaged matter flux onto the star; Ṁ w is the matter flux into the wind, calculated through the surface S(r = 10, z = ±10) at condition v min > 0.1vesc; f eff is the propeller efficiency; Lsd is the angular momentum flux from the stellar surface; Ėm and Ėf are the energy fluxes carried into the outflows through the surface S(r = 10, z = ±10) by matter and by the magnetic field, respectively.

Figure 4 .
Figure 4. Temporal variation of energy fluxes in model μ60.Ėm (blue line) and Ėf (green line) are the energy fluxes carried by matter and the magnetic field, respectively.Dashed lines show the time-averaged values.

Figure 5 .
Figure 5. Doppler tomograms of AE Aqr, computed with the filtered backprojection method in Hα spectral line during 3 nights.vx and vy are inertial velocities plotted with respect to a corotating coordinate system, where the x−axis is along the line of centers of the binary, and the y−axis is in the direction of instantaneous orbital motion.Also drawn in the velocity space are the centre of mass of the binary (+), the position of the white dwarf (dot at Vy = −102 km/s), the Roche lobe of the secondary star, and the gas stream.The darker regions in the image denote stronger emission.FromWelsh et al. (1998).

Figure 6 .
Figure 6.Left panel: Fraction of matter flowing into the wind through surface S(r = 30, z = ±30) at the condition that the minimum velocity in the outflow is v min in models µ30, µ60 and µ100.Right panel: Matter fluxes in the intervals of velocities ∆v k as a function of velocity.

Figure 7 .
Figure 7.The time-averaged maximum velocity , vmax, as a function of fastness.Vertical lines show values of fastness in our models of AE Aqr.Horizontal lines show the maximum time-averaged velocities, corresponding to these models.Figure is based on the plot from Romanova et al. (2018).

Table 1 .
Dimensionless values obtained in

Table 4 .
Time-averaged energy fluxes to the wind carried by matter, Ėm , and by the magnetic field, Ėf .Ėtot is the total energy flux into the wind.