© ESO, 2015

1. Introduction

Magnetic cataclysmic variables are close binary systems containing a highly magnetized (B ~ 10−200 MG) white dwarf that accretes matter from a low-mass companion. For a subclass of those objects, called AM Her stars or polars, the magnetic field of the white dwarf is strong enough to lock the system into a synchronous rotation and to prevent the formation of an accretion disk (Cropper 1990; Warner 1995). Matter leaves the companion, follows a ballistic trajectory, and is then captured by the magnetic field lines of the white dwarf. Its falls down to the magnetic pole of the compact star at supersonic free-fall velocity (v_ff ~ 5−10 × 10⁸ cm/s), forming an accretion column structure.

The collision of the infalling flow with the white dwarf photosphere leads to forming a reverse shock, which propagates upstream and heats the plasma to temperature T_s ≃ 10−50 keV. These temperatures are due to the high accretion velocity (${\mathit{T}}_{\mathrm{s}}\mathrm{\propto }{\mathit{v}}_{\mathrm{ff}}^{\mathrm{2}}$), which is directly linked to the mass (and radius) of the white dwarf. Those extreme temperatures, as well as the strong magnetic field, imply that the post-shock region radiates through several radiative processes from the optical to X-ray domains. The radiation arising from these binary systems comes mainly from the accretion region. As a consequence, these systems can be used as a perfect probe for testing accretion models in high-energy regimes. The energy losses shape the density and temperature profiles in the post-shock region and induce a strong gradient of density, temperature, and velocity near the white dwarf surface. Since the radiation properties depend on the mass of the white dwarf, modelling the emitting region can be used to determine the physical properties of the compact object (see King & Lasota 1979; Cropper et al. 1998). Numerous studies have shown, in fact, that modelling of the post-shock region physics strongly affects any determination of the white dwarf properties (Ishida et al. 1991; Wu et al. 1995; Yuasa et al. 2010; Hayashi & Ishida 2014).

A stationary regime is frequently assumed to be established to help in interpreting observations (Fisher & Beuermann 2001). However, radiative instabilities that can potentially develop in the post-shock region can induce dynamical effects and lead to the unsteady behaviour of the post-shock accretion column (PSAC). Indeed, the efficiency of the bremsstrahlung process can trigger the cooling instability. The shock thus oscillates around the stationary position. First numerical studies of such accretion flows have highlighted the presence of the cooling instability in the column (Langer et al. 1981) supported by theoretical studies (Chevalier & Imamura 1982). A comprehensive review of the different stationary and instationary models was given later by Wu (2000).

The oscillation of the front shock leads to a variation in the luminosity due to the induced variation of the emission volume, as well as of the conditions of density and temperature inside the column (Imamura et al. 1991; Wu et al. 1992). The characteristic timescale of this phenomenon is the cooling timescale associated to the bremsstrahlung emission, i.e. t_brem ~ 1 s. This timescale value is consistent with the periods of the quasi-periodical oscillations (QPOs) observed in the optical light curves of a few polars (Middleditch 1982; Mason et al. 1983; Bonnet-Bidaud et al. 1985, 1996; Larsson 1987, 1989). These rapid oscillations are characterized by an amplitude of 1−5% rms and a period of 1−2.5 s. From the model of cooling instability, an X-ray counterpart of optical oscillations is theoretically expected. Thus QPOs have been extensively searched in X-ray light curves. New results are presented in a related paper (Bonnet-Bidaud et al. 2015, hereafter referred as Paper I).

In this context, various analytical and numerical studies have been performed to determine the sensitivity of the instability to physical parameters, such as the reverse shock Mach number (Strickland & Blondin 1995; Ramachandran & Smith 2006) or the influence of different possible physical processes taking part in the post-shock region, like two-temperature effects (Imamura et al. 1996; Saxton & Wu 1999; Saxton et al. 2005) or transverse magnetic field (Tòth & Draine 1993; Ramachandran & Smith 2005). A few studies have also described the effect of the coupling between several cooling processes (Wu et al. 1992; Saxton et al. 1998). In particular, recent astronomical observations (see Paper I) show that several polars are in regimes where either the bremsstrahlung dominates the losses or both bremsstrahlung and cyclotron have to be taken into account.

The two main conclusions of those previous works are the following: firstly, the cyclotron losses damp the oscillation. Thus, a critical magnetic field amplitude can be found for which the oscillation disappears after a few periods. To sustain an oscillation for a longer time, a noisy accretion rate has been invoked (Wolff et al. 1991; Wood et al. 1992; Wu et al. 1992). Secondly, the oscillations of the emission in the X-ray domain due to the bremsstrahlung process and in the optical domain due to the cyclotron one are expected to present the same frequencies and the same order of magnitude of the relative amplitudes (Wu et al. 1992). To complement those results and to better understand the physics of this specific regime, we present new numerical simulations to study the unsteady behaviour of the accretion column. Those simulations are compared to observational results (see Paper I).

In Sect. 2, the standard model is reviewed in order to build the parameter space and the relevant accretion regimes of the column as a function of the three leading parameters (specific accretion rate, white dwarf mass, and magnetic field). In Paper I, four parameters are considered because the specific accretion rate is not reached directly by observations but through the combination of the total accretion rate and the column cross-section. In Sect. 3, the general behaviour of a column in an unsteady state is studied with particular attention to the additional instabilities that can develop in the post-shock region. Then in Sect. 4, the dependency of the observational QPO parameters on the three leading parameters is developed further and discussed.

2. Models of radiation and magnetic processes in high-energy accretion regions

2.1. Brief review of the accretion column standard model

In the standard model of accretion column, the matter coming from the companion is trapped by the gravitational field of the compact object and falls at free-fall velocity, as given by v_ff = (2GM_WD/R_WD)^1/2, onto the magnetic poles of the white dwarf, with G, M_WD, and R_WD being the gravitational constant, the white dwarf mass, and radius, respectively. It should be noted that the matter is actually captured at the Alfven radius, R_A, so the free fall velocity should be corrected by a factor (1 − R_WD/R_A). However, for polars, the Alfven radius is much larger than the white dwarf radius so this effect has been neglected. Considering the Nauenberg mass-radius relation (Nauenberg 1972), the free-fall velocity only depends on the white dwarf mass. Typically, this velocity is v_ff ~ 8.75 × 10⁸ cm/s for a 1 M_⊙ white dwarf. The incoming flow is thus highly supersonic with a Mach number of ℳ ~ (10−100).

Owing to the magnetic collimation, matter spreads over a small fraction of the white dwarf surface, denoted f and defined by $\mathit{f}\mathrm{=}\mathit{S}\mathit{/}\mathrm{\left(}\mathrm{4}\mathit{\pi }{\mathit{R}}_{\mathrm{WD}}^{\mathrm{2}}\mathrm{\right)}$ where S is the area of the column cross-section. This quantity is not well constrained by observations. As an example, this quantity has been evaluated in the range f ~ 10^-5−3 × 10^-4 by Imamura (1984) for a sample of polars. The impact of the accreting flow onto the white dwarf surface leads to the formation of a reverse shock, which counter-propagates in the incoming plasma. The shock transforms the kinetic energy into thermal energy, and the PSAC is strongly heated to high-temperature T_s. Using the Rankine-Hugoniot conditions and assuming that the accreted medium is described by the perfect gas approximation, the post-shock temperature scales as ${\mathit{v}}_{\mathrm{ff}}^{\mathrm{2}}$ and is given by ${\mathit{T}}_{\mathrm{s}}\mathrm{=}\frac{\mathit{\mu }{\mathit{m}}_{\mathrm{H}}}{{\mathit{k}}_{\mathrm{B}}}\frac{\mathrm{2}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}}{\mathrm{(}\mathit{\gamma }\mathrm{+}\mathrm{1}{\mathrm{)}}^{\mathrm{2}}}{\mathit{v}}_{\mathrm{ff}}^{\mathrm{2}}\mathrm{=}\mathrm{1.1}\mathrm{\times }{\mathrm{10}}^{\mathrm{7}}\mathrm{K}{\left(\frac{{\mathit{v}}_{\mathrm{ff}}}{{\mathrm{10}}^{\mathrm{8}}\mathrm{cm}\mathrm{/}\mathrm{s}}\right)}^{\mathrm{2}}\mathit{,}$(1)with μ, m_H, k_B, and γ being the mean molecular weight, the hydrogen mass, the Boltzmann constant, and the adiabatic index, respectively. We use μ = 0.5 (fully ionized hydrogen) throughout the paper and γ = 5/3.

At such high temperatures, matter radiates strongly through optically thin bremsstrahlung, which dominates the hard X-ray emission (Frank et al. 2002). Moreover, for the highest magnetized polars, the energy can also be transported efficiently by the cyclotron process (Saxton et al. 1998), which emits mainly in the optical domain for the range of polar magnetic fields. As infalling material cools and slows down, it also gets denser with a sharp gradient near the surface. Thanks to the energy losses by radiation, the reverse shock position eventually stabilizes. Consequently, when using the steady state approximation and a cooling dominated by bremsstrahlung losses, the height of post-shock accretion region can be estimated and is given by Wu et al. (1994) as ${\mathit{x}}_{\mathrm{s}}\mathrm{=}\mathrm{3}\mathrm{\times }{\mathrm{10}}^{\mathrm{7}}\mathrm{cm}{\left[\frac{\mathit{ṁ}}{\mathrm{1}{\mathrm{g}\mathrm{/}\mathrm{cm}}^{\mathrm{2}}\mathrm{/}\mathrm{s}}\right]}^{-1}{\left[\frac{{\mathit{M}}_{\mathrm{WD}}}{\mathrm{0.5}{\mathit{M}}_{\mathrm{\odot }}}\right]}^{\mathrm{3}\mathit{/}\mathrm{2}}{\left[\frac{{\mathit{R}}_{\mathrm{WD}}}{{\mathrm{10}}^{\mathrm{9}}\mathrm{cm}}\right]}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}\mathit{,}$(2)with ṁ = ρ₀v_ff the specific accretion rate. A more complex but exact analytic expression can be found in Falize et al. (2009). With the Nauenberg relation, the relation (2) can be written as a function of ṁ and M_WD. Given the order of magnitude of the accretion shock height, the local white dwarf gravitational field g does not vary enough to strongly modify the post-shock structure (gx_s ≪ v_ff). This approximation is valid as long as the white dwarf mass is below 1 M_⊙ (Cropper et al. 1999).

On the whole, three key parameters determine the post-shock structure: the white dwarf mass, M_WD, determining the free-fall velocity and the post-shock temperature; the specific accretion rate, ṁ, determining the post-shock density; and the local magnetic field strength, B, which determines the relative importance of cyclotron over bremsstrahlung losses. By comparing the characteristic timescales of the different physical processes at the front shock, the possible post-shock regimes are defined (Lamb & Masters 1979; Langer et al. 1982).

Four main regimes are obtained and can be explicitly defined in the (M_WD, ṁ, B) diagram. Figure 1 displays these regimes for a 0.8 M_⊙ white dwarf in the (ṁ, B) plane. They are delimited by lines corresponding to an equal balance between two processes such as bremsstrahlung-cyclotron, electron-ion collisions, and ion-ion collisions. In the following, the characteristic timescales of these processes are recalled and used to derive the (ṁ, B) relations corresponding to the limits of the regimes.

Among all the radiative effects occurring in dense and hot magnetized plasmas, the dominant processes that shape the PSAC structure are the bremsstrahlung and cyclotron emissions (Lamb & Masters 1979; Wu 2000). Through the estimate of optical depth of the post-shock medium to bremsstrahlung photons (Frank et al. 2002), it can be deduced that radiation escapes the matter freely. Thus the loss of energy can be modelled by a cooling term that can be analytically expressed (Rybicki & Lightman 1979) as ${\mathrm{\Lambda }}_{\mathrm{brem}}\mathrm{=}\frac{\mathrm{32}\mathit{/}\mathrm{3}}{\mathrm{(}\mathrm{4}\mathit{\pi }{\mathit{\epsilon }}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{3}}}\sqrt{\frac{\mathrm{2}{\mathit{\pi }}^{\mathrm{3}}\mathit{k}}{{\mathit{m}}_{\mathrm{e}}^{\mathrm{3}}}}\frac{{\mathit{Z}}^{\mathrm{2}}{\mathit{e}}^{\mathrm{6}}}{\mathit{h}{\mathit{c}}^{\mathrm{3}}}{\mathit{n}}_{\mathrm{e}}{\mathit{n}}_{\mathrm{i}}{\mathit{g}}_{\mathrm{B}}\mathrm{\left(}\mathit{Z,T}\mathrm{\right)}\sqrt{\mathit{T}}\mathrm{\equiv }{\mathrm{\Lambda }}_{\mathrm{0}}{\mathit{\rho }}^{\mathrm{2}}{\mathrm{T}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}$(3)where g_B, m_e, e, Z, n_e, n_i, and ρ are the Gaunt factor, the electron mass, the electron charge, the mean charge number, the electronic and ionic density, and the gas density, respectively. Assuming that the plasma consists of fully ionized hydrogen, the parameter Λ₀ in Eq. (3)is estimated to Λ₀ ~ 5 × 10²⁰ erg cm³ g^-2 K^−1/2 s^-1. According to its definition in Rybicki & Lightman (1979), the Gaunt factor is taken as g_B ≃ 1 (cf. Saxton & Wu 1999).

The bremsstrahlung cooling timescale, noted t_brem is given by$\begin{array}{ccc}{\mathit{t}}_{\mathrm{brem}}& \mathrm{=}& \frac{\mathit{\rho e}}{{\mathrm{\Lambda }}_{\mathrm{brem}}}\\ {\mathit{t}}_{\mathrm{brem}}& \mathrm{\simeq }& \end{array}$(4)with e the internal energy per unit of mass. Parameters in relation (4)have been normalized to the typical values encountered in polars. Thus, for a white dwarf of 0.8 M_⊙ with ṁ ~ 1 g/cm²/s, the cooling timescale associated to the bremsstrahlung process is t_brem ~ 1 s.

When the local magnetic field close to the white dwarf surface is strong enough, the cyclotron process can also be an efficient radiative process for cooling the medium. Depending on the radiation frequency, the plasma is optically thick or optically thin to photons which strongly complexifies the modelling (Langer et al. 1982; Kalomeni et al. 2005). Theoretical simplifications and dimensional considerations allow the cyclotronic transport to be reduced to a cooling function (Langer et al. 1982; Saxton 1999), expressed as a function of the column cross-section area (S), the magnetic field, and the local thermodynamical conditions as $\begin{array}{ccc}{\mathrm{\Lambda }}_{\mathrm{cycl}}& \mathrm{=}& \mathrm{1.2}\mathrm{\times }{\mathrm{10}}^{\mathrm{8}}{\mathrm{erg}\mathrm{/}\mathrm{cm}}^{\mathrm{3}}\mathrm{/}\mathrm{s}{\left(\frac{\mathit{S}}{{\mathrm{10}}^{\mathrm{15}}{\mathrm{cm}}^{\mathrm{2}}}\right)}^{\mathrm{-}\mathrm{17}\mathit{/}\mathrm{40}}{\left(\frac{\mathit{B}}{\mathrm{10}\mathrm{MG}}\right)}^{\mathrm{57}\mathit{/}\mathrm{20}}\\ & & \end{array}$(5)For a given accretion column surface and a given white dwarf magnetic field, Eq. (5)is reduced to $\begin{array}{ccc}{\mathrm{\Lambda }}_{\mathrm{cycl}}& \mathrm{\equiv }& \end{array}$(6)where Λ_0,c is thus a constant.

As for the bremsstrahlung, the cyclotron cooling timescale is given by $\begin{array}{ccc}{\mathit{t}}_{\mathrm{cycl}}& \mathrm{=}& \frac{\mathit{\rho e}}{{\mathrm{\Lambda }}_{\mathrm{cycl}}}\\ {\mathit{t}}_{\mathrm{cycl}}& \mathrm{\simeq }& \mathrm{2}\mathrm{\times }{\mathrm{10}}^{\mathrm{2}}\mathrm{s}\hspace{0.17em}{\left(\frac{\mathit{S}}{{\mathrm{10}}^{\mathrm{15}}{\mathrm{cm}}^{\mathrm{2}}}\right)}^{\mathrm{17}\mathit{/}\mathrm{40}}{\left(\frac{\mathit{B}}{\mathrm{10}\mathrm{MG}}\right)}^{\mathrm{-}\mathrm{57}\mathit{/}\mathrm{20}}\\ & & \end{array}$(7)As expected, the cyclotronic effects are more efficient when the magnetic field is stronger and the plasma is hotter. The optically thick property of the medium can be retrieved from Eq. (7), since t_cycl increases with the surface and the density of the column. Indeed, the wider and denser the column, the longer it takes the radiation to escape the column.

To determine the prevailing mechanism between bremsstrahlung and cyclotron processes, the two cooling timescales computed right below the shock front are compared through a dimensionless ratio, noted ϵ_s (Wu et al. 1994), and given by $\begin{array}{ccc}{\mathit{ϵ}}_{\mathrm{s}}& \mathrm{\simeq }& \frac{{\mathit{t}}_{\mathrm{brem}}}{{\mathit{t}}_{\mathrm{cycl}}}\mathrm{=}\mathrm{2}\mathrm{\times }{\mathrm{10}}^{-4}{\left(\frac{\mathit{S}}{{\mathrm{10}}^{\mathrm{15}}{\mathrm{cm}}^{\mathrm{2}}}\right)}^{\mathrm{-}\mathrm{17}\mathit{/}\mathrm{40}}\\ & & \end{array}$(8)The bremsstrahlung is the dominating process when ϵ_s< 1 near the front shock, while the cyclotron is the dominant one when ϵ_s> 1. Since ϵ_s ∝ B^57/20ṁ^{− 37/20}, the first regime is obtained for high accretion rates and weak magnetic fields (Regime 1 in Fig. 1), and the second regime corresponds to strong magnetic fields and low accretion rates (Regime 2 in Fig. 1).

From Eq. (8)a critical magnetic field depending on ṁ, S, and v_ff can be obtained at the limit ϵ_s = 1, which separates the two regimes and is computed as $\begin{array}{ccc}{\mathit{B}}_{{\mathit{ϵ}}_{\mathrm{s}}\mathrm{=}\mathrm{1}}& \mathrm{\simeq }& \mathrm{200}\mathrm{MG}{\left(\frac{\mathit{ṁ}}{\mathrm{1}{\mathrm{g}\mathrm{/}\mathrm{cm}}^{\mathrm{2}}\mathrm{/}\mathrm{s}}\right)}^{\mathrm{37}\mathit{/}\mathrm{57}}{\left(\frac{{\mathit{v}}_{\mathrm{ff}}}{{\mathrm{10}}^{\mathrm{8}}\mathrm{cm}\mathrm{/}\mathrm{s}}\right)}^{\mathrm{-}\mathrm{39}\mathit{/}\mathrm{19}}\\ & & \mathrm{\times }{\left(\frac{\mathit{S}}{{\mathrm{10}}^{\mathrm{15}}{\mathrm{cm}}^{\mathrm{2}}}\right)}^{\mathrm{17}\mathit{/}\mathrm{114}}\mathrm{·}\end{array}$(9)With the free-fall velocity expression and the white dwarf mass-radius relation, B_{ϵs = 1} can be written as only a function of ṁ, S, and M_WD . It should be emphasized that the bremsstrahlung losses increase towards the surface, while the cyclotron losses decrease (see Sect. 4.2). In fact the radiative losses lead to the cooling of the plasma and its densification. This evolution implies an increase in the cyclotron cooling timescale (t_cycl ∝ ρ^17/20T^{− 3/2}) and a decrease in the bremsstrahlung scale (${\mathit{t}}_{\mathrm{brem}}\mathrm{\propto }{\mathit{\rho }}^{-1}\sqrt{\mathit{T}}$). Thus, even if the cyclotron dominates at the front shock, the bremsstrahlung will become the dominant process lower in the column, inside the high-energy region (see Sect. 4.2).

Fig. 1 Different regimes of accretion as a function of the accretion rate and the magnetic field strength, considering a white dwarf of M = 0.8 M⊙ and an accretion section of S = 1015 cm2.

The one-temperature approximation is valid as long as the Coulombian collisions between ions and electrons are sufficient to insure that the two populations have the same temperature despite some energy losses by radiation. When it is not the case, the electrons and ions are thermally decoupled and two-temperature effects appear. The characteristic timescale of electron-ion collisions is given by (Spitzer 1965) ${\mathit{t}}_{\mathrm{ei}}\mathrm{=}\mathrm{6.22}\mathrm{\times }{\mathrm{10}}^{-5}\mathrm{s}{\left(\frac{\ln {\mathrm{\Lambda }}_{\mathrm{C}}}{\mathrm{20}}\right)}^{-1}{\left(\frac{\mathit{ṁ}}{\mathrm{1}{\mathrm{g}\mathrm{/}\mathrm{cm}}^{\mathrm{2}}\mathrm{/}\mathrm{s}}\right)}^{-1}{\left(\frac{{\mathit{v}}_{\mathrm{ff}}}{{\mathrm{10}}^{\mathrm{8}}\mathrm{cm}\mathrm{/}\mathrm{s}}\right)}^{\mathrm{4}}\mathit{,}$(10)where lnΛ_C is the Coulomb logarithm. These effects are significant when t_ei is greater than the timescale of the most efficient cooling process. In the accretion column case, it happens when t_ei ~ t_cycl (Regime 3 in Fig. 1). A critical magnetic field can be estimated, as for the previous processes, and is given by $\begin{array}{ccc}\mathit{B}\mathrm{\simeq }{\mathit{B}}_{\mathrm{0}}{\left(\frac{\ln {\mathrm{\Lambda }}_{\mathit{c}}}{\mathrm{20}}\right)}^{\frac{\mathrm{20}}{\mathrm{57}}}{\left(\frac{\mathit{ṁ}}{\mathrm{1}{\mathrm{g}\mathrm{/}\mathrm{cm}}^{\mathrm{2}}\mathrm{/}\mathrm{s}}\right)}^{\frac{\mathrm{37}}{\mathrm{57}}}{\left(\frac{{\mathit{v}}_{\mathrm{ff}}}{{\mathrm{10}}^{\mathrm{8}}\mathrm{cm}\mathrm{/}\mathrm{s}}\right)}^{\frac{-157}{\mathrm{57}}}{\left(\frac{\mathit{S}}{{\mathrm{10}}^{\mathrm{15}}{\mathrm{cm}}^{\mathrm{2}}}\right)}^{\frac{\mathrm{17}}{\mathrm{114}}}& & \end{array}$(11)with B₀ = 1.9 × 10³ MG.

Finally for extremely low accretion rates, as observed at some epochs in most polars (Cropper 1990), the flow becomes collisionless (t_ii>t_cycl, where t_ii is the ion-ion collision timescale) and no shock is formed at the impact (regime 4 in Fig. 1). This regime has been identified observationally among different polars referred to as LARPs (low accretion rate polars, see Schwope et al. 2002). By using the expression of t_ii given in Kylafis & Lamb (1982), this specific regime turns out to be crucial when the magnetic field is above a critical value that has the same dependency on the accretion rate, free-fall velocity, and accretion surface as in Eq. (11)but with a different value of the constant B₀ = 8.5 × 10³ MG.

Considering more typical conditions commonly found in polars (see Paper I), we focus our numerical studies on the first two regimes (Regimes 1 and 2), which are the only ones considered in the following sections.

2.2. Radiation hydrodynamic numerical code HADES

The post-shock region is structured by the cooling processes, thus the position of the reverse shock can be evaluated using Eq. (2). For typical values of density and velocity of the incoming flow, x_s ~ 3 × 10⁷ cm considering bremsstrahlung cooling alone¹. As a consequence, the curvature of the magnetic field line close to the pole can be neglected, and the accretion column can be assimilated to a cylinder. Moreover, considering the spatial extension of the column, the magnetic field gradient can be neglected along the post-shock region, and the terms linked to the magnetic field in the ideal MHD description vanished so an hydrodynamical description is sufficient to model the base of the column (Cropper 1990; Warner 1995; Wu 2000).

The code Hades is dedicated to the resolution of radiative hydrodynamics equations (Michaut et al. 2011). Considering our system, we treat radiative losses as a source term in the energy conservation equation $\begin{array}{ccc}& & \frac{\mathit{\partial \rho }}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathit{.}\mathrm{\left(}\mathit{\rho }{v}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}\\ & & \frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{\left(}\mathit{\rho }{v}\mathrm{\right)}\mathrm{+}\mathrm{\nabla }\mathit{.}\left(\mathit{\rho }{v}\mathrm{\otimes }{v}\mathrm{+}\mathit{P}\mathit{.}\mathit{I̅}\right)\mathrm{=}{0}\mathit{,}\\ & & \frac{\mathit{\partial }}{\mathit{\partial t}}\left(\mathit{\rho e}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{\mathit{v}}^{\mathrm{2}}\right)\mathrm{+}\mathrm{\nabla }\mathit{.}\left(\mathit{\rho e}{v}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{v}{\mathit{v}}^{\mathrm{2}}\mathrm{+}\mathit{P}{v}\right)\mathrm{=}\mathrm{-}\mathrm{\Lambda }\mathrm{\left(}\mathit{\rho ,T}\mathrm{\right)}\mathit{,}\end{array}$where t, ρ, P, v, e, Λ(ρ,T), and are the time, density, pressure, velocity, internal energy per unit of mass, cooling function of the system, and the identity matrix, respectively. Equations (12)−(14)come from the conservation of the mass, impulsion, and energy of the system, respectively.

The perfect gas approximation is assumed to close the system of equations. The code is designed to solve high-Mach number flows (ℳ > 10) thanks to a MUSCL-Hancock scheme (Van Leer 1979) using a HLLE solver to compute the Riemann problem (Einfeldt 1988). The cooling function Λ is expressed as a sum of the power laws of density and of pressure, thus modelling the various cooling processes at stake in the shocked medium (Busschaert et al. 2013a). In this numerical model, the initial conditions are the mass density (ρ₀), the free-fall velocity (v_ff), and the temperature of the infalling flow, respectively. In the presented simulations, the infalling flow velocity evolves from 3 × 10⁸ cm/s (M_WD = 0.4M_⊙) to 5.5 × 10⁸ cm/s (M_WD = 0.8M_⊙), and the infalling flow temperature has been set at T = 10⁷ K to ensure a Mach number in the flow higher than ten to obtain a strong shock but also to ensure numerical robustness since the code cannot sustain Mach numbers that are too high.

Moreover, a specific boundary condition has been developed to mimic the accretion onto the white dwarf at the bottom of the column. Indeed, since the cooling processes stop when the plasma is sufficiently cooled, a layer of cold and dense material accumulates at the white dwarf surface. On the astrophysical scale, this matter is assimilated by the star, so to avoid this unphysical accumulation, the mass flow is conserved along the column. Therefore at the bottom of the column, the boundary condition satisfies ρ_max × v_min = ρ₀ × v_ff. The parameter v_min is determined so that the velocity contrast does not exceed three orders of magnitude, insuring the calculation’s robustness. We verified that the minimum temperature reached resembles temperature at which the cooling is no longer efficient. More details can be found in Busschaert et al. (2013a) concerning the implementation and validation of this boundary condition.

3. Non-stationary phenomena in the post-shock region

The structure and dynamics of the PSAC can be strongly modified by unsteady phenomena. In particular the cooling instability (Chevalier & Imamura 1982; Mignone 2005) can develop, and the post-shock matter can radiate so strongly that it again becomes supersonic close to the white dwarf surface, which leads to the formation of a secondary shock (Falle 1981). These unsteady phenomena imply that no stationary regime can be established. These time variabilities are used to explain observations of QPOs in the optical luminosity curves of some polars. The effects of these two mechanisms onto astronomical observations are discussed. Our aim is to quantify the impact of these effects to evaluate the relevancy of the stationary model. These conditions are critical since they are used to determine the white dwarf mass (Wu et al. 1995).

3.1. The cooling instability

Thanks to the bremsstrahlung cooling, the post-shock region can developed the cooling instability (Chevalier & Imamura 1982). Indeed, according to the Field criterium (Field 1965), a radiating fluid characterized by a cooling function in the form Λ ∝ ρ²T^β is instable if β< 1. The linear analysis has been the subject of various studies where the stationary solution is perturbated. The non-linear regime can be explored in details with numerical hydrodynamical simulations (Imamura et al. 1984; Mignone 2005).

Figure 2 shows the dynamics of the accretion shock simulated using the Hades code with a cooling function Λ ∝ ρ²T^1/2, relevant for a pure bremsstrahlung cooling (see Eq. (4)). The position of the front shock with time is plotted, showing a quasi-periodic oscillation with a period of the order of t_brem (see Eq. (5)), in accordance with previous works (Mignone 2005). This shock height variation induces a modification of the post-shock volume of the hot and highly radiating material, as well as a modification of the post-shock temperature. As a consequence, the integrated luminosity is also modified.

Fig. 2 Evolution of the front shock position measured from the white dwarf surface (top blue line) and the X-ray (0.5−10 keV) luminosity (bottom red line), for a PSAC with an accretion rate of ṁ = 10 g/cm2/s onto a 0.8 M⊙ white dwarf and a pure bremsstrahlung cooling.

The X-ray luminosity and spectrum are computed using the local bremsstrahlung emissivity combined with the density and temperature profiles extracted from the simulations. To compare with observations, the integrated X-ray luminosity is evaluated by integrating the radiative losses both over the volume of the column and over photon energies between 0.5 and 10 keV that correspond to the XMM-Newton detector range used in Paper I. The resulting integrated X-ray light curve is also shown in Fig. 2 and clearly displays oscillations at the same main frequency than those observed for the shock height but with a more complex structure.

Fig. 3 Density (top) and temperature (bottom) profiles of the column at successive times across the oscillation (shown in the inset) for an accretion rate of ṁ = 10 g/cm2/s, a white dwarf of 0.8 M⊙, and a pure bremsstrahlung cooling. Starting after t = 6.35 s, a secondary shock clearly forms at the basis of the accretion column and propagates upstream while the main shock falls towards the white dwarf surface.

3.2. Formation of secondary shock: Falle criterion

The careful examination of the temperature and density profiles of the postshock region and their evolution with time clearly show the formation of a secondary shock (see Fig. 3). This secondary shock, which has been somewhat overlooked by previous works, is associated with the development of the catastrophic instability, or Falle instability, that can take place for a post-shock region cooled mainly by the bremsstrahlung process. This instability was first evoked by Falle (1981) in the context of the filamentation of supernovae remnants. Secondary shocks in the post-shock region can arise since right below the main shock, the fluid velocity is subsonic, but owing to the cooling, the flow can again become supersonic close to the white dwarf surface. This kind of phenomenon has already been noted in several studies concerning accretion shocks (Blondin & Cioffi 1989; Strickland & Blondin 1995; Wolff et al. 1989; Mignone 2005), but its impact onto the main shock dynamics or onto the observables has not been discussed yet.

To determine the favourable conditions for developing the catastrophic cooling instability, a local characteristic cooling length, L_cool, is defined as the product of the sound velocity, c_s, and the cooling timescale, t_cool as ${\mathit{L}}_{\mathrm{cool}}\mathrm{=}{\mathit{c}}_{\mathrm{s}}\mathrm{\times }{\mathit{t}}_{\mathrm{cool}}\mathrm{\propto }{\mathit{T}}^{\mathrm{3}\mathit{/}\mathrm{2}\mathrm{-}\mathit{\beta }}{\mathit{\rho }}^{\mathrm{1}\mathrm{-}\mathit{\alpha }}\mathit{,}$(15)where a general cooling function Λ = Λ₀ρ^αT^β has been considered. Then, an acoustic wave propagates on a length of the order of L_cool before its dynamics is strongly impacted by the cooling processes. It can be shown that to develop the catastrophic cooling instability, the local cooling length has to be smaller than the local distance from the surface of the white dwarf (Falle 1981) $\frac{{\mathit{L}}_{\mathrm{cool}}}{\mathit{x}}\mathit{<}\mathrm{1}\mathit{,}$(16)and considering a cooling function characterized by a parameter α = 2, as for the bremsstrahlung process, then a second criterion is for the second parameter to satisfy (Falle 1981) $\mathit{\beta }\mathit{<}\mathrm{3}\mathit{/}\mathrm{2.}$(17)The evolutions of density and temperature in the post-shock region are presented at four representative times in Fig. 3. Considering a system cooled by the bremsstrahlung process, the condition (17)is satisfied since in this case β = 1/2. The first condition is more complicated to check out. We calculated the ratio L_cool/x over the column, which decreases from the shock towards the surface. When the secondary shock appears, the ratio tends towards 1 by superior value close to the surface. To ensure that the structure observed is truly a shock, we computed the Mach number associated with this structure over an oscillation (see Fig. 4) and checked that the value is around 1.2−2. Thus, assuming a system cooled only by bremsstrahlung, a secondary shock appears in the post-shock region (see Fig. 3, t = 6.4 s). The secondary shock is formed near the surface of the white dwarf and progresses upstream up to the main reverse shock (see Fig. 3, t = 6.43 s) until the two shocks collide.

Fig. 4 Evolution with time of the position of the main shock (blue filled points), the secondary shock (red filled points) and the rarefaction wave (grey dashed line) in the same condition as in Fig. 3. The Mach number value for the two shocks are also shown by the open symbols. Different critical times through the oscillation are labelled by letters and discussed in the text.

The positions of the main and secondary shocks have been tracked and are presented in Fig. 4, together with the Mach number value for the two shocks. The interaction of the two shocks across the oscillation can be fully described. At maximum extension, the main shock progressively slows down and then falls back towards the white dwarf surface due to the cooling instability (Fig. 4, a’). Meanwhile, the secondary shock develops from the white dwarf surface (Fig. 4, a) and propagates upwards until it collides with the main shock (Fig. 4, b). After both shocks have collided, the main shock propagates upstream, and a rarefaction wave (Zel’dovitch & Raizer 2002) propagates through the shocked medium towards the white dwarf surface, as expected (see Fig. 3, t = 6.48 s and Fig. 4, c’). In the meantime, a new secondary shock has formed close to the white dwarf surface (Fig. 4, c). On the temperature profile (see Fig. 3), both the rarefaction wave and the secondary shock are clearly visible at t = 6.48 s near x ~ 5−6×10⁶ cm and x ~ 2 × 10⁶ cm, respectively. As the rarefaction wave encounters the newly formed secondary shock (see Fig. 4, d), the shock accelerates further and then collides again with the main shock. Then this cyclic behaviour starts again. The secondary shock always appears close to the white dwarf surface when the extension of the PSAC is close to its maximum. After each collision between the two shocks, only the main shock remains and is pushed away from the white dwarf surface.

Secondary shocks are also present in previous studies, which considered at least four different boundary conditions (see Fig. 5 in Mignone 2005), though such shocks were not fully discussed in this work. To ensure that this physical phenomena truly exists, we verified that its presence and this cyclic behaviour are independent of boundary conditions at the basis of the accretion column and initial conditions.

To study the development of the catastrophic cooling instabilities, several cooling functions have been tested with the parameters α = 2 and various β values. For values of β higher than 0.7, this instability does not develop. Indeed, even if the second criterion (see Eq. (17)) is fulfilled, the cooling process acts too quickly after the initial collision to allow the main shock going far enough from the white dwarf surface to fulfil the first criterion (see Eq. (16)). Also, it should be noted that when both cyclotron and bremsstrahlung are effective, considering the stabilizing effect of cyclotron losses, the secondary shock does not appear for a high magnetic field for the same reason.

From the simulations, the detailed X-ray spectrum emitted by the PSAC over an oscillation can be computed and is shown in Fig. 5. The spectrum evolves strongly through the different phases of the oscillation. When only the primary shock is present (t = 6.35 s), the emitted spectrum is softer than when the secondary shock has appeared (t = 6.4−6.43 s). The relative difference between the mean spectrum calculated over a couple of oscillations and the spectrum emitted by the theoretical stationary column is presented at the bottom of Fig. 5. The difference between the two spectra is quite negligible at low energy, but at high energy, the mean spectrum is significantly lower than the stationary one. Thus, if we attempt to evaluate the white dwarf mass through the mean spectrum, we expect to find a lower mass.

Fig. 5 Top: evolution of the emitted bremsstrahlung spectrum through the oscillation for the same parameters and the times indicated in the inset of Fig. 3. The mean spectrum calculated over a couple of oscillations is represented by the thick grey line. Bottom: relative difference between the mean spectrum and the spectrum emitted by the theoretical stationary column.

Considering both cooling instabilities, the column structure, when cooled mainly by the bremsstrahlung process, is quite different from the steady state solution. This strongly affects both light curves and emitted spectra.

4. Influence of system parameters on luminosity variations

In Paper I, observational X-ray light curves of a set of polars were analysed to search for oscillations as predicted by the cooling instability model. To directly compare these observational constraints with the theoretical expectations from the shock model, we used the Hades simulations to compute the amplitudes and frequencies of the expected oscillations more precisely for different sets of parameters that are relevant to the polars, considering both bremsstrahlung and cyclotron cooling.

4.1. Bremsstrahlung dominated columns

The temporal characteristics of the front shock position and the computed PSAC X-ray flux emerging from the simulations were analysed using fast Fourier transforms (FFT) to better understand the properties of the quasi-periodic phenomena. Representative power spectra of the shock front position and X-ray luminosity are presented in Fig. 6 where we consider a white dwarf of 0.8 M_⊙ with an accretion rate of ṁ = 10 g/cm²/s and only taking the bremsstrahlung losses into account. The FFT covers simulation data obtained at 0.01 s resolution over 2048 time steps at long times between t = 19.76 s and t = 30 s.

Several discrete frequencies emerge from theses analyses, and the same frequencies are present in both signals but with different relative amplitudes. The different amplitudes are explained by considering that not only does the variation in the emission volume modify the luminosity but also the physical conditions of temperature and density throughout the post-shock region. For instance, the post-shock temperature depends on the pre-shock velocity in the shock frame. Thus, when the shock is going upstream (resp. downstream), the infalling flow goes faster (resp. slower) than the free-fall velocity, and the post-shock temperature is higher (resp. lower). This phenomenon is clearly visible in the temperature profiles presented in Fig. 3. Amongst the most significant frequencies, the lowest is found at 2.24 Hz (period of 0.45 s) (mode f). This frequency is very visible in the first instants (see Fig. 2) and remains present later. It is correlated with the setting up of the PSAC structure: when the flow impacts the surface, the shock goes away from the surface and progressively slows down. At t = 0.24 s, the shock falls back and reaches the minimum position at t = 0.42 s before bouncing against the secondary shock. The second strongest mode is present at a frequency of 6.7 Hz (period of 0.15 s) (mode a). This typical timescale corresponds to the cooling timescale that is calculated using Eq. (4). Modes of higher frequencies are also present and can make a significant contribution to the oscillations. A mode (b) has been emphasized in Fig. 6, since it has a strong contribution for higher magnetic field strength. For the regime with pure bremsstrahlung cooling, the FFT exhibits the same overall behaviour for the different white dwarf masses tested (M_WD = 0.4,0.6,0.8 M_⊙). However, the relative contribution of each mode can be slightly different, and it can be demonstrated how the frequency depends on the accretion rate and the white dwarf mass.

Fig. 6 Fast Fourier transform of the front shock position (left) and the X-ray luminosity (right) respectively normalized by their mean values for the same parameters as in Fig. 3. The frequencies are identical for both signals but with different amplitudes. The main modes discussed in the text are referred as (f), (a), and (b).

Considering the bremsstrahlung cooling as dominant (see Regime 1, Fig. 1), the mechanism that induces the oscillations is the cooling instability with the oscillation frequency expected to evolve as the inverse of the cooling time (see Eq. (4)) $\mathit{\nu }\mathrm{\propto }{\mathit{t}}_{\mathrm{brem}}^{-1}\mathrm{\propto }\mathit{ṁ}\mathit{/}{\mathit{v}}_{\mathrm{ff}}^{\mathrm{2}}\mathit{.}$(18)Introducing the expression of the free-fall velocity coupled to the white dwarf mass-radius relation (Nauenberg 1972), the frequency depends on the specific accretion rate and the white dwarf mass as $\mathit{\nu }\mathrm{=}{\mathit{A}}_{\mathrm{0}}\mathit{ṁ}{\left[{\left(\frac{\mathrm{1.44}{\mathit{M}}_{\mathrm{\odot }}}{{\mathit{M}}_{\mathrm{WD}}}\right)}^{\mathrm{8}\mathit{/}\mathrm{3}}\mathrm{-}{\left(\frac{\mathrm{1.44}{\mathit{M}}_{\mathrm{\odot }}}{{\mathit{M}}_{\mathrm{WD}}}\right)}^{\mathrm{4}\mathit{/}\mathrm{3}}\right]}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}$(19)where A₀ is a constant. In Fig. 7, we represent the dominant oscillation frequency of the X-ray luminosity (mode a), extracted from simulations for three different white dwarf masses and accretion rates. The dependence from Eq. (19)is plotted as solid lines and gives a good estimate of the oscillation frequency. The constant A₀, which provides the best fit, is A₀ = 0.9,0.9, and 0.7 cm²/g for ṁ = 0.1,1, and 10 g/cm²/s, respectively, while the expected theoretical value is A_0,th = 0.5 cm²/g. These values are slightly different since the theoretical value is evaluated just behind the shock and does not take the strong variations along the column into account.

These results imply that, when QPOs are observed, their frequencies should allow constraining the value of the accretion rate and white dwarf mass (see Paper I).

Fig. 7 Variation in the main oscillation frequency (mode a) of the X-ray light curve as a function of the white dwarf mass, for different accretion rates in the case of pure bremsstrahlung cooling. The curves represent the best fits following the theoretical tendency given by the relation (19).

4.2. PSAC efficiently cooled by the cyclotron process

When the magnetic field is strong enough, the cyclotron cooling becomes efficient and a two-process cooling has to be considered. The losses for the bremsstrahlung and cyclotron processes are presented in Fig. 8 as a function of the position in the column for a typical configuration. The cyclotron losses are higher than the bremsstrahlung ones at the shock front in this case since ϵ_s = 1.2, which corresponds to a magnetic field of 30.3 MG, assuming a column section of 10¹⁵ cm². Because the bremsstrahlung mechanism is more sensitive to variation in density than temperature (see Eq. (3)), the losses progressively increase from the shock to the surface of the white dwarf along with the density. In contrast, since cyclotron mechanism is more sensitive to variations in temperature than in density (see Eq. (6)), it is more efficient at the shock where matter is hotter than closer to the white dwarf surface. Thus, ϵ_s, defined as the ratio of cooling times (see Eq. (8)), is always maximum at the shock and then decreases along the column towards the white dwarf surface.

Fig. 8 Variation in the bremsstrahlung (solid blue line) and cyclotron (dotted red line) losses along the column assuming a white dwarf of 0.8 M⊙, an accretion rate of ṁ = 10 g/cm2/s and ϵs = 1.2 (corresponding to a magnetic field of 30.3 MG) at t = 100 s. The abscisses have been normalized to the front shock position (xs ~ 7.7 × 106 cm).

The analysis in Paper I gives upper limits for the amplitude of QPOs in X-rays. Up to now, oscillations in light curves have only been observed in the optical range for a few selected sources. In the context of the shock instability model for polars, it should be emphasized that optical variations associated to cyclotron are expected together with the X-ray ones linked to bremsstrahlung. From the simulations, we extracted both components separately and compare them here in terms of amplitude and frequency.

In Fig. 9, the luminosity in the optical range due to the cyclotron losses is compared to the X-ray emission due to the bremsstrahlung process, considering a white dwarf with M_WD = 0.8 M_⊙ and an accretion rate of 10 g/cm²/s, for a typical magnetic field of 18.7 MG (ϵ_s = 0.3). The behaviour in the two energy ranges is very similar. We present in Fig. 10 the corresponding Fourier transforms of the cyclotron and bremsstrahlung emissions, where the power shown is the oscillation amplitude normalized by the mean luminosity. They present the exact same frequencies of oscillations but with an oscillation amplitude in the dominant mode (a), which is lower for the cyclotron signal (~18%) than the bremsstrahlung one (~28%). This behaviour has been systematically observed in all simulations. Thus, we expect to observe stronger oscillations in the X-ray light curve than in the optical domain. As already noted by Wu et al. (1992), we confirm that the two signals are slightly out of synchronization.

Fig. 9 Variation in the (0.5−10 keV) bremsstrahlung (upper blue line) and cyclotron (lower red line) luminosity considering a white dwarf of 0.8 M⊙, an accretion rate of ṁ = 10 g/cm2/s, and a magnetic field of 18.7 MG (ϵs = 0.3).

Fig. 10 Fourier transform of the (0.5−10 keV) bremsstrahlung luminosity and cyclotron total luminosity with the same parameters as in Fig. 9.

Considering the perturbative analysis, the system is expected to stabilize as cyclotron becomes a strong contribution in the radiative losses (Saxton et al. 1998; Saxton 1999; Saxton & Wu 1999). In the results presented in this work, for a white dwarf of 0.8 M_⊙, the oscillations have been considered as completely damped when the front shock oscillation is only a few meshes in length. It corresponds to a luminosity variation of less than ~0.1%. We observed that for the accretion rate of ṁ ~ 1 g/cm²/s, the oscillations are effectively damped for ϵ_s ~ 1.5; and for ṁ ~ 10 g/cm²/s, the oscillations are almost damped for the same value of ϵ_s.

As the magnetic field strength increases, the low-order modes are expected to be rapidly damped. This type of behaviour is found in our numerical simulations. The typical evolution of the X-ray light curve FFT is presented in Fig. 11 for various values of ϵ_s that correspond to a magnetic field strength ranging from 0 to 28 MG, considering an M_WD = 0.8 M_⊙ white dwarf and an accretion rate of ṁ = 10 g/cm²/s. For weak magnetic fields, the dominant mode is around (ν ≃ 6−7 Hz), but for higher fields, this mode decreases strongly and disappears, while higher modes (ν ≃ 20−25 Hz) appear.

Fig. 11 Changes in the FFT of the X-ray light curve for a growing magnetic field characterized by ϵs, for a white dwarf of 0.8 M⊙ and an accretion rate of ṁ = 10 g/cm2/s. The range in ϵs from 0 to 1 corresponds to a magnetic field from 0 to 28 MG. The different frequency modes are identified as (a), (b), and (c) and listed in Table 2.

In Fig. 11, it is clear that the distribution of the various modes and their identification are not obvious. However, a few main modes can be tracked when the magnetic field grows. The frequencies and relative amplitudes of the three main modes are reported in Tables 1 and 2 and plotted in Fig. 12. We see that mode (a) decreases monotonically as the magnetic field increases. Its amplitude becomes negligible for a magnetic field around 23 MG. Mode (b) increases initially and then drops suddenly around the same critical magnetic field as mode (a). Mode (c) has the same behaviour as mode (b), but its amplitude keeps increasing when the first two modes decrease strongly. Its amplitude ultimately decreases for higher magnetic field amplitude. The different behaviours exhibited by these first three modes is interesting since it reveals the great richness of the non-linear dynamics of the PSAC. Moreover, it should be emphasized that progressively, as the magnetic field grows, the power is distributed in higher order modes. It should be noted that this non-linear behaviour is coherent with the linear stability studies published by Saxton et al. (1998). In Tables 1 and 2, the total relative amplitude is also given as a quadratic sum of the amplitudes in the different modes when the power is split into different frequencies.

Fig. 12 Variation in the relative amplitude of the X-ray oscillations for the three main modes observed in the FFT for a white dwarf of 0.8 M⊙ and an accretion rate of ṁ = 10 g/cm2/s. The total relative amplitude is given by the full black line.

For a given mode, a frequency variation with the magnetic field strength is observed in the simulations that can be explained analytically. Under the assumption that oscillations are induced by the cooling instability, the expected frequency is directly connected to the cooling time by the relation $\mathit{\nu }\mathrm{\propto }\mathrm{1}\mathit{/}{\mathit{t}}_{\mathrm{cool}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{t}}_{\mathrm{brem}}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{t}}_{\mathrm{cyc}}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{t}}_{\mathrm{brem}}}\left[\mathrm{1}\mathrm{+}\mathit{ϵ}\mathrm{\left(}\mathit{B}\mathrm{\right)}\right]\mathit{,}$(20)where ϵ(B) is the ratio between the bremsstrahlung and cyclotron cooling time in each spatial position. If we suppose that the typical timescales calculated behind the shock are representative of the value along the accretion column, then ϵ(B) ≃ ϵ_s, and the relation can be written as $\mathit{\nu }\mathrm{=}{\mathit{\nu }}_{\mathrm{brem}\mathit{,}\mathrm{s}}\left[\mathrm{1}\mathrm{+}{\mathit{C}}_{\mathrm{0}}{\left(\frac{\mathit{B}}{\mathrm{10}\mathrm{MG}}\right)}^{\mathrm{57}\mathit{/}\mathrm{20}}\right]\mathit{,}$(21)where ν_brem,s is the oscillation frequency at the front shock in the asymptotic bremsstrahlung dominated regime, and C₀ is a constant that depends on the accretion rate and the free fall velocity as$\begin{array}{ccc}{\mathit{C}}_{\mathrm{0}}& \mathrm{\simeq }& \mathrm{2}\mathrm{\times }{\mathrm{10}}^{-4}{\left[\frac{\mathit{S}}{{\mathrm{10}}^{\mathrm{15}}{\mathrm{cm}}^{\mathrm{2}}}\right]}^{\mathrm{-}\mathrm{17}\mathit{/}\mathrm{40}}{\left[\frac{\mathit{ṁ}}{\mathrm{1}{\mathrm{g}\mathrm{/}\mathrm{cm}}^{\mathrm{2}}\mathrm{/}\mathrm{s}}\right]}^{\mathrm{-}\mathrm{37}\mathit{/}\mathrm{20}}\\ & & \end{array}$(22)The specific behaviour of the oscillation mode frequencies predicted here is different to the results presented by Wu et al. (1992) in the context of a stochastic accretion. In Wu et al. (1992), the suggested simplified evolution law ν ∝ (1 + C₀B²) is only based on phenomenological assumptions and not justified by physical considerations.

The evolution of the frequency for a range of magnetic fields from 0 to 22 MG is presented in Fig. 13 for the two main modes (a) and (b), assuming a 0.8M_⊙ white dwarf and an accretion rate of 10 g/cm²/s. It should be noted that for strongest magnetic fields, the amplitude of those modes decreases drastically so no oscillation frequency can be found. The curves in Fig. 13 represent the best fits using the analytical trend presented in Eq. (21). The constant determined from these fits is C₀ = 0.022 and is the same for the two modes. In this specific configuration, the expected analytical value of the constant as derived from Eq. (22)is C_0,th = 0.061. The difference with the fitted value from the numerical simulations comes from the structure of the post-shock region. Indeed, to theoretically infer the characteristic timescale, only the values just behind the shock have been taken into account (see also 4.1).

Fig. 13 Evolution of oscillation frequency with magnetic field for the first two modes (a and b) for a white dwarf of 0.8 M⊙ and an accretion rate of ṁ = 10 g/cm2/s. The lines are the best fits according to the predicted dependency of Eq. (21).

5. Conclusions

In this work, a numerical study of the dynamics of PSAC has been presented. Radiation hydrodynamic simulations have shown that the coupling between magnetic fields, radiation processes, and extreme hydrodynamics leads to a complex dynamics and a wide variety of dynamical behaviours. To validate the numerical simulations, they have been compared to theoretical stability criteria, analytical results, and dimensional tendencies. Good agreement has been obtained, in general. To link simulation results to observational data presented in Paper I, a post-treatment of numerical results were performed to extract and study the X-ray and optical luminosities and their dynamics, the X-ray spectral luminosity, and the frequency spectrum. Given the observational data (see Paper I), we have focused our analysis on the two most important regimes of PSAC: the bremsstrahlung-dominated regime (weak magnetic field and/or high accretion rate) and the coupled bremsstrahlung-cyclotron regime (strong magnetic fields and/or lower accretion rate).

In the bremsstrahlung-dominated regime, one of the most important results is the development of the catastrophic instability, which leads to the formation of secondary shocks. According to radiation hydrodynamic simulations, those additional shocks modify the dynamics of the primary shock, as well as density and temperature profiles. It induces modifications on the total and spectral emissivities and plays an important role in the oscillations of the shock height. The complex coupling between the secondary shock dynamics and the cooling instability development leads to QPOs in the light curves. Specifically, the secondary shock implies a harder X-ray spectrum than in the stationary regime, which can affect the mass determination from the classical X-ray approach. The presence and role of the secondary shocks in PSAC have not been studied much as yet. In principle, the presence of the secondary shock can also be revealed through variations induced in the X-ray spectrum across the oscillations. However, the sensitivity of the present X-ray telescopes is not yet high enough to detect such variations. Laboratory experiments, using powerful lasers already in use or near completion such as NIF or LMJ, are a promising way to directly investigate these shocks.

In the coupled bremsstrahlung-cyclotron regime, numerical simulations show the stabilizing effect of the cyclotron process. The predicted amplitude of oscillations in the X-ray light curve considering a pure bremsstrahlung regime is ~30%. As the magnetic field grows, the amplitude of oscillation decreases and a critical magnetic field depending on the accretion parameters can be found where no oscillation persists. Thanks to numerical simulations, a more detailed study of the distribution of power in the different modes of oscillations has been explored. In particular, the low-order modes are the first to be damped as the magnetic field grows, and the power is then distributed in higher order modes. These two effects have been theoretically expected in the linear regime and persist with non-linear effects. An important correlation between the magnetic field and the oscillation frequency has been established. A comparison between the numerical results and the theoretical relations based on dimensional arguments leads to very good agreement.

Numerical simulations exhibit strong oscillations in the X-ray luminosity with the same frequencies as in the optical luminosity curve. They show the presence of QPOs with a wider range of frequencies in the X-ray domain (0.6−50 Hz) than commonly assumed. This wider range, which has never been fully explored observationally until now, has been investigated in Paper I with negative results. There is therefore no compelling evidence yet of significant X-ray oscillations, in contrast to what is expected from shock instabilities. These numerical results raise important questions on the origin of physical processes that lead to QPO.

The inconsistency between X-ray and optical observational data, on the one hand, and between numerical expectations and observations, on the other, highlights that the physics of QPO is still not completely understood. The situation is somewhat similar among classical T Tauri stars for which shock-heated accreting material is expected to give rise to QPOs (Koldoba et al. 2008; Sacco et al. 2008) that are not yet observed (Drake et al. 2009; Günther et al. 2010). In this context, 2D axisymmetric MHD simulations have shown that significant distortions of the accretion structure can occur as a result of the gas pressure. At the bottom of the column, accumulation of matter may lead to outflows parallel to the star surface that can increase the cooling, and the oscillations are expected to be reduced or suppressed (Orlando 2010, 2013). However, the analogy between polars and T Tauri stars is somewhat limited. Indeed, for the typical magnetic field and velocity and density regimes, the plasma parameter (noted β = P_th/P_mag, where P_th and P_mag are respectively the thermal and magnetic pressure) is significantly different. It is in the range β_{T Tauri} ~ 10^-2−10² for T Tauri (Matsakos et al. 2013), but only β_polar ~ 10^-4−10^-3 for polars. Therefore the magnetic collimation is much more efficient for polars, and 2D effects are expected to be less important. In T Tauri, low β values have been shown to lead to a fragmentation of the column into independent oscillating tubes or fibrils (Orlando et al. 2010; Matsakos et al. 2013). If out of phase, the oscillations of the different fibrils can cancel each other out. Such an effect can also occur in polar PSAC and would also induce a frequency widening, as well as a significant decrease in the overall QPO amplitudes. These effects will be the subject of an upcoming study where numerical simulations will take multi-dimensional effects into account, as well as non-homogeneous accretion and a more detailed treatment of the boundary conditions at the white dwarf photosphere.

Additionally, it should be emphasized that the PSAC presents interesting similarity properties that allow us to reproduce a scaled model in laboratory on measured scales (Falize et al. 2011b; Busschaert et al. 2013b). This is the main objective of a laboratory astrophysics project: the POLAR project (Falize et al. 2011a). The possibility of producing relevant experiments has been demonstrated on the LULI2000 facility (Falize et al. 2012). The experimental way is a complementary approach to increase our understanding of the radiation hydrodynamics of accretion shock. Thus the QPO physics will be explored using dedicated designed experiments with powerful lasers.

Acknowledgments

We wish to thank an anonymous referee for helpful comments. C.B. and C.M. acknowledge financial support from the PNPS of CNRS/INSU, France, as well as from the Scientific Council of the Observatoire de Paris (France). C.B., E.F., M.M., and C.M. acknowledge financial support from the PNHE, France. The authors also thank H.C. Nguyen and M. Mancini for the development of HADES.

References

All Tables

All Figures

	Fig. 1 Different regimes of accretion as a function of the accretion rate and the magnetic field strength, considering a white dwarf of M = 0.8 M⊙ and an accretion section of S = 1015 cm2.
In the text

	Fig. 2 Evolution of the front shock position measured from the white dwarf surface (top blue line) and the X-ray (0.5−10 keV) luminosity (bottom red line), for a PSAC with an accretion rate of ṁ = 10 g/cm2/s onto a 0.8 M⊙ white dwarf and a pure bremsstrahlung cooling.
In the text

Fig. 3 Density (top) and temperature (bottom) profiles of the column at successive times across the oscillation (shown in the inset) for an accretion rate of ṁ = 10 g/cm2/s, a white dwarf of 0.8 M⊙, and a pure bremsstrahlung cooling. Starting after t = 6.35 s, a secondary shock clearly forms at the basis of the accretion column and propagates upstream while the main shock falls towards the white dwarf surface.

In the text

	Fig. 4 Evolution with time of the position of the main shock (blue filled points), the secondary shock (red filled points) and the rarefaction wave (grey dashed line) in the same condition as in Fig. 3. The Mach number value for the two shocks are also shown by the open symbols. Different critical times through the oscillation are labelled by letters and discussed in the text.
In the text

	Fig. 5 Top: evolution of the emitted bremsstrahlung spectrum through the oscillation for the same parameters and the times indicated in the inset of Fig. 3. The mean spectrum calculated over a couple of oscillations is represented by the thick grey line. Bottom: relative difference between the mean spectrum and the spectrum emitted by the theoretical stationary column.
In the text

	Fig. 6 Fast Fourier transform of the front shock position (left) and the X-ray luminosity (right) respectively normalized by their mean values for the same parameters as in Fig. 3. The frequencies are identical for both signals but with different amplitudes. The main modes discussed in the text are referred as (f), (a), and (b).
In the text

	Fig. 7 Variation in the main oscillation frequency (mode a) of the X-ray light curve as a function of the white dwarf mass, for different accretion rates in the case of pure bremsstrahlung cooling. The curves represent the best fits following the theoretical tendency given by the relation (19).
In the text

	Fig. 8 Variation in the bremsstrahlung (solid blue line) and cyclotron (dotted red line) losses along the column assuming a white dwarf of 0.8 M⊙, an accretion rate of ṁ = 10 g/cm2/s and ϵs = 1.2 (corresponding to a magnetic field of 30.3 MG) at t = 100 s. The abscisses have been normalized to the front shock position (xs ~ 7.7 × 106 cm).
In the text

	Fig. 9 Variation in the (0.5−10 keV) bremsstrahlung (upper blue line) and cyclotron (lower red line) luminosity considering a white dwarf of 0.8 M⊙, an accretion rate of ṁ = 10 g/cm2/s, and a magnetic field of 18.7 MG (ϵs = 0.3).
In the text

	Fig. 10 Fourier transform of the (0.5−10 keV) bremsstrahlung luminosity and cyclotron total luminosity with the same parameters as in Fig. 9.
In the text

	Fig. 11 Changes in the FFT of the X-ray light curve for a growing magnetic field characterized by ϵs, for a white dwarf of 0.8 M⊙ and an accretion rate of ṁ = 10 g/cm2/s. The range in ϵs from 0 to 1 corresponds to a magnetic field from 0 to 28 MG. The different frequency modes are identified as (a), (b), and (c) and listed in Table 2.
In the text

	Fig. 12 Variation in the relative amplitude of the X-ray oscillations for the three main modes observed in the FFT for a white dwarf of 0.8 M⊙ and an accretion rate of ṁ = 10 g/cm2/s. The total relative amplitude is given by the full black line.
In the text

	Fig. 13 Evolution of oscillation frequency with magnetic field for the first two modes (a and b) for a white dwarf of 0.8 M⊙ and an accretion rate of ṁ = 10 g/cm2/s. The lines are the best fits according to the predicted dependency of Eq. (21).
In the text