HYPERCRITICAL ACCRETION ONTO A NEWBORN NEUTRON STAR AND MAGNETIC FIELD SUBMERGENCE

Most observed neutron stars show clear evidence for the presence of strong magnetic fields. In the case of magnetars (Woods & Thompson 2006) the estimated strength of the surface magnetic field is of the order of 10¹⁵ G, while for the grind of the mill radio pulsar 10¹² G is a typical value. Lower magnetic fields are, however, found in millisecond pulsars (Phinney & Kulkarni 1994) and in neutron stars in low-mass X-ray binaries (Psaltis 2006), but in both cases past or present, respectively, accretion is thought to be the cause of the magnetic field reduction. There is, however, a small group of neutron stars, found in young supernova remnants and dubbed CCOs ("Central Compact Object"; Pavlov et al. 2002; Ho 2013), which exhibit little or no evidence for the presence of a magnetic field.

The case of the supernova SN 1987A is also intriguing since no evidence for the presence of a compact object has yet been found, despite extensive searches (see, e.g., Manchester 2007). The presence of a young energetic pulsar is clearly excluded, but a slowly rotating neutron star, with period P > 100 ms, and with a not too strong dipolar magnetic field, B_d < 10¹² G, is still compatible with all current data (Manchester 2007; Ng et al. 2009; Larsson et al. 2011). Such a period and low magnetic field are within the range of values inferred, when possible, in some of the CCOs (Gotthelf & Halpern 2008; Gotthelf et al. 2013). The possibility that this supernova may have produced a black hole, as proposed, e.g., by Brown & Bethe (1994) on the basis of a very soft dense matter equation of state and, consequently, a low neutron star maximum mass around 1.6 M_☉, is now very slim in view of the existence of a 2 M_☉ pulsar (Demorest et al. 2010).

The origin of neutron star magnetic fields is still an unsolved problem (for recent reviews, see Spruit 2008, 2009). Two main mechanisms, a fossil field from the progenitor compressed during the core collapse and a proto-neutron star dynamo, are still competing and are likely both needed to explain the large variety of observed field strengths. In these scenarios, the magnetic field generation and/or adjustment process terminates within a minute after the neutron star's birth.

After this early field development the story is not necessarily over. The supernova shock is still pushing its way through the outer layers of the progenitor, and if it encounters a density discontinuity, a reverse shock can be generated. Depending on its strength and on how far out it was generated, this reverse shock can induce strong accretion onto the newborn neutron star on timescales of hours. Notice that this delayed accretion has to be distinguished from the initial fall-back that occurs seconds after the initial core collapse. A particularly favorable configuration for such late accretion is present in the cases for which the progenitor had a tenuous H/He envelope surrounding a dense He core, as in Type IIb supernovae or in SN 1987A (Smartt 2009).

Chevalier (1989) argued in favor of such late accretion in the case of SN 1987A and developed a simple analytical model for it (see also Brown & Weingartner 1994 for a similar model and Figure 1 of Bernal et al. 2010 for a depiction of these scenarios). Highly super-Eddington accretion results, in which the photons are trapped within the accretion flow and the energy liberated by the accretion is lost through neutrino emission close to the neutron star surface. Such a regime has been termed "hypercritical accretion" (Blondin 1986; Chevalier 1989; Houck & Chevalier 1991) and requires an accretion rate $\dot{m}$ higher than about $10^3 \, \dot{m}_\mathrm{Edd}$, where $\dot{m}_\mathrm{Edd} \sim 10^5$ g cm⁻² s⁻¹ is the Eddington rate. We note that this regime is relevant not only for SNe, but for gamma-ray bursts as well (Narayan & Quataert 2005; Nakar 2007; Lee & Ramirez-Ruiz 2007; Gehrels et al. 2009), allowing rapid mass accretion onto newborn black holes to produce the required power to account for isotropic equivalent luminosities that can reach above 10⁵² erg s⁻¹ in the prompt emission (Narayan et al. 2001; Kohri & Mineshige 2002; Lee et al. 2005), and possibly extended emission episodes as well (Rosswog 2007; Lee et al. 2009; Metzger et al. 2010). In the case of core collapse events, after the reverse shock hits the neutron star surface a third shock develops and starts moving outward against the infalling matter. Once this accretion shock stabilizes it will separate the infalling matter from an extended envelope in quasi-hydrostatic equilibrium.

Following the suggestion of Muslimov & Page (1995), Geppert et al. (1999) presented simple 1D ideal MHD simulations of the effect of this post-supernova hypercritical accretion on the newborn neutron star magnetic field. The result was a rapid submergence of the field into the neutron star. After accretion stopped, the field could diffuse back to the surface and result in a delayed switch-on of a pulsar (Michel 1994; Muslimov & Page 1995). Depending on the amount of accreted matter, the submergence could be so deep that the neutron star may appear and remain unmagnetized for more than a Hubble time (Geppert et al. 1999). This scenario was recently revisited by Ho (2011) and Viganò & Pons (2012) and applied to study the field evolution of the CCOs. Notice that, because of the violence and short duration of the accretion phase, an ideal MHD treatment of it is justified. In contradistinction, the back-diffusion of the field after accretion stops is due to the finite electrical conductivity of the neutron star crust matter. As a result, the back-diffusion becomes dependent on the thermal evolution of the star because of the temperature dependence of the electrical conductivity.

The 1D MHD simulations of Geppert et al. (1999) were, however, carried out with some considerable simplifications. In a previous work (Bernal et al. 2010, Paper I hereafter) we presented results of 2D ideal MHD simulations of hypercritical accretion onto a magnetized neutron star surface, with simple magnetic field configurations. In that paper we simulated a hypercritical accretion flow in a rectangular domain, with the neutron star surface lying at its base, in which an initial uniform, horizontal or vertical, magnetic field is present, as well as a horizontal field with strength decreasing with height. For highly hypercritical accretion rates, $\dot{m} \gtrsim 10^{10} \, \dot{m}_\mathrm{Edd}$, complete submergence of the magnetic field was observed in a timescale of a few tens to a few hundreds of milliseconds.

In the present work we extend the scenarios considered in Paper I and consider an initial magnetic field configuration that consists of a magnetic loop protruding from the neutron star surface. Such a field configuration is intended to be more realistic than the ones we studied in Paper I: the feet of the field loop will be frozen into the neutron star surface while the rest of the loop will follow the motion of the matter since, with the field strengths we employ, matter and ram pressures still strongly dominate over the magnetic pressure. Further, the dimensions of the accretion column are significantly larger in the present study, allowing for greater freedom in the solution regarding the direction parallel to the stellar surface. We mostly consider the 2D case in a cylindrical box and one example of a 3D model. Moreover, we also extend our simulations to much lower accretion rates, down to $10^{6} \, \dot{m}_\mathrm{Edd}$, but still assume the accreting matter is unmagnetized.

In Section 2 we describe the numerical method and a simple analytical model. In Section 3 we describe our results, of which an interpretation is proposed in Section 4. Finally, we conclude in Section 5.

The present work is a direct extension of our previous work presented in Paper I. As described in Paper I, we use a customized version of the hydrodynamic code AMR FLASH2.5 (Fryxell et al. 2000), with the Split 8 wave solver, which solves the whole set of MHD equations. We work in the ideal MHD regime with only numerical resistivity and viscosity. The matter equation of state is an adaptation of FLASH's HELMHOLTZ package, which includes contributions from the nuclei, e⁻ − e⁺ pairs, and radiation, plus the Coulomb correction. Neutrino energy losses are dominated by the e⁻ − e⁺ annihilation process, but we also include the photo-neutrino, plasmon decay, e-ion bremsstrahlung, and synchrotron processes, which are implemented in a customized module (see Paper I). No nuclear reactions are taken into account.

Since the MHD equations can be solved by FLASH only in Cartesian coordinates, we consider wide columns in 2D and 3D, with a base of Δx = 2 × 10⁶ cm centered on x = 0 in 2D, or Δx × Δz = (2 × 10⁶) × (2 × 10⁶) cm² centered on (x, z) = (0, 0) in 3D, and a height Δy of several times 10⁶ cm, with y = 10⁶ being the neutron star surface. The mapping of the surrounding cone above the neutron star into a column is illustrated in Figure 1. The gravitational acceleration is taken as g_y = −GM/y², and we assume a neutron star mass of 1.4 M_☉.

Zoom In Zoom Out Reset image size

We consider as magnetic initial condition a magnetic field loop, in the shape of a hemi-torus. On the central hemi-circle of the loop the field has a strength B₀ = 10¹² G and about it is shaped as a Gaussian, i.e., with strength B(d) = B₀exp [ − (d/R_L)²)], d being the distance to the loop central hemi-circle and R_L = 1 km. The two feet of the loop are centered at x = −5 and x = +5 km, and z = 0 in the 3D model. The initial condition for matter will be described below.

As boundary conditions, we impose mass inflow along the top edge of the computational domain, and periodic conditions along the sides. At the bottom, on the neutron star surface, we use a custom boundary condition that enforces hydrostatic equilibrium (see Paper I). For the magnetic field, the two lateral sides are also treated as periodic boundaries, while at the bottom the field is frozen from the initial condition, i.e., the two feet of the loop are anchored into the neutron star and no field can be pushed into the star by the accretion. On the top boundary the magnetic field is set to zero, i.e., we assume the accreting matter to be non-magnetized.

We are interested in following the evolution of the magnetic field under the heavy hypercritical accretion from a reverse shock (see Figure 1 in Paper I) that reaches the neutron star surface hours after the core collapse. The various accretion rates, $\dot{m}$, we consider will be expressed in terms of a fiducial rate

$$\begin{equation} \dot{m}_0 = 1.75 \times 10^{15} \; \mathrm{g \, cm^{-2} \, s^{-1}} \end{equation} \tag{ 1 }$$

and always assumed to be constant during the entire simulations. This corresponds to a total accretion rate onto the neutron star of

$$\begin{equation} \dot{M}_0 \sim 350 \; M_\odot \, \mathrm{yr}^{-1} \end{equation} \tag{ 2 }$$

and is the accretion rate originally estimated by Chevalier (1989) for SN 1987A. We consider two different scenarios regarding the initial conditions for the accreting matter.

In the first case, corresponding to initial free fall, we start the simulation just before the reverse shock reaches the magnetic loop, which is initially immersed in a low density, 10² g cm⁻², medium, and above it the falling-back matter has a free-fall velocity and a density obtained by mass conservation:

$$\begin{equation} v_\mathrm{ff} = \sqrt{2GM/y} \;\;\;\; \mathrm{and} \;\;\;\; \rho _\mathrm{ff} = \dot{m}/v_\mathrm{ff}. \end{equation} \tag{ 3 }$$

As we will see, after a few tens of milliseconds an accretion envelope develops that is in a quasi-hydrostatic equilibrium, separated from the continuous accretion inflow by an accretion shock.

The structure of this quasi-hydrostatic equilibrium envelope (QHEE hereafter) can easily be obtained analytically (adjusting the treatment of Chevalier 1989 to our box geometry). The pressure, density, and velocity profiles are simply given by

$$\begin{equation} P = P_\mathrm{sh} (y_\mathrm{sh}/y)^4,\quad \rho = \rho _\mathrm{sh} (y_\mathrm{sh}/y)^3,\quad v = v_\mathrm{sh} (y/y_\mathrm{sh})^3, \end{equation} \tag{ 4 }$$

where the first two values come from imposing hydrostatic equilibrium for a polytropic equation of state, P∝ρ^γ with γ = 4/3, and the velocity is fixed by mass conservation. Once the shock position y_sh is known (see below for its determination), P_sh and ρ_sh are determined by the strong shock condition and v_sh by mass conservation as

$$\begin{equation} P_\mathrm{sh} = \frac{49}{8} \rho _\mathrm{sh} v_\mathrm{sh}^2, \quad \rho _\mathrm{sh} = 7 \rho _\mathrm{ff}(y_\mathrm{sh}), \quad v_\mathrm{sh} = \frac{1}{7} v_\mathrm{ff}(y_\mathrm{sh}). \end{equation} \tag{ 5 }$$

For given M and R, and a fixed $\dot{m}$, the vertical location of the accretion shock, y_sh, is controlled by energy balance

$$\begin{equation} \frac{GM\dot{m}}{R} = \int _R^\infty \epsilon _\nu (y) dy \end{equation} \tag{ 6 }$$

between the accretion power and the integrated neutrino losses, per unit neutron star surface area. Neutrino losses are dominated by e⁻ − e⁺ pair annihilation, for which a simple rate is (Dicus 1972)

$$\begin{equation} \epsilon _\nu = 1.83 \times 10^{-34} P^{2.25} \; \mathrm{erg \, cm^{-3} \, s^{-1}} \;. \end{equation} \tag{ 7 }$$

Due to the resulting strong y dependence of _ν, y⁻⁹, the value of the upper limit in the integral of Equation (6) is not important, and this fixes the height of the accretion shock as

$$\begin{equation} y_\mathrm{sh} \simeq 7.4 \times 10^6 \, (\dot{m}_0/\dot{m})^{10/63} \; \mathrm{cm} \,. \end{equation} \tag{ 8 }$$

A stationary envelope, in quasi-hydrostatic equilibrium, will expand, or shrink, so that the physical conditions at its base allow neutrinos to carry away all the energy injected by the accretion. Once emitted, neutrinos will act as an energy sink provided the material is optically thin to them. Under the present conditions of density and temperature at the base of the flow, the fluid consists mainly of free neutrons, protons, and electrons. The main sources of neutrino opacity are then coherent scattering off neutrons and protons and pair annihilation. For example, the corresponding cross section for coherent scattering (Tubbs & Schramm 1975; Shapiro & Teukolsky 1983) is σ_N = (1/4)σ₀[E_ν/(m_ec²)]², where σ₀ = 1.76 × 10⁻⁴⁴ cm². As these are thermal neutrinos, their energy is E_ν ∼ k_BT, and with temperatures T ≲ 10¹¹ K ∼10 MeV as we will find, we have σ_N ≲ 7 × 10⁻⁴² cm². The maximum densities reached at the bottom of the envelope will be below 10¹¹ g cm⁻³, and in such conditions, the neutrino mean free path is l_ν = (n_Nσ_N)⁻¹ ≳ 2.5 × 10⁶ cm, which is safely larger than the depth of the dense envelope, of the order of a few km. Above this dense region the envelope density decreases rapidly, Equation (4), and the whole envelope is practically transparent to neutrinos. We will, hence, ignore neutrino absorption and heating.

To compare the magnetic pressure P_mag = B²/8π ∼ 10²³(B/10¹² G)² to the matter pressure in the QHEE, from the above we deduce

$$\begin{equation} P = 3.26 \times 10^{28} (\dot{m}/\dot{m}_0)^{4/9} (10 \,\mathrm{km}/y)^4 \, \mathrm{dyne\,cm^{-2}} \,. \end{equation} \tag{ 9 }$$

For all accretion rates we will consider we will always be in the regime P_mag ≪ P.

In the second set of scenarios, we start with this QHEE and the magnetic loop immersed within it. In the cases of highly hypercritical accretion rates the accretion shock lies within the integration domain and we complement the QHEE initial condition with free-fall matter above it. For this purpose we use the ρ, P, and v profiles given by the above analytical solution but adjust them within the first few kilometers above the neutron star surface to take into account the piling up of matter onto the surface, as shown in Figure 2. We adjust both ρ and P as rising above the analytical solution and v dropping to zero within these 2 km. (These adjustments mock up the numerical results of the QHEE obtained from a free-fall initial condition that will be shown in Figure 3.)

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Besides being used as an initial condition in the QHEE scenarios, this analytical envelope model will be useful in the free-fall scenarios for comparison with the numerical results when the system has reached a stationary state.

In Paper I we performed detailed comparisons of the hydro solver and the MHD solver with zero magnetic field and obtained excellent agreement, as well as agreement with the above described analytical envelope models when a stationary state had been reached. More details can be found in Paper I.

We will here present results from two series of simulations, with different accretion rates. In the first series, we consider highly hypercritical accretion rates, $\dot{m} = (1, 10, 100) \cdot \dot{m}_0$, while in the second series weakly hypercritical rates, $\dot{m} = (10^{-4}, 10^{-3}, 10^{-2}, 10^{-1}) \cdot \dot{m}_0$, are studied. We will describe the scenario of matter initially in free-fall, with the resulting formation of an accretion shock and the development of a QHEE, and the simplified scenario where we start with a previously formed QHEE and follow its evolution. The time step in the FLASH code is adaptive and depends on local conditions. Typically, the time resolution of the simulations is dt ≃ 10⁻⁷ s.

With the high accretion rates of the first series we are able to follow the complete evolution of the accretion shock and the development of the QHEE in the free-fall scenario and compare the results with the simplified purely QHEE scenario. In contradistinction, in the presence of the lower accretion rates of the second series, the accretion shock would move outward to very large distances and leave the computational domain and we are, hence, forced to restrict ourselves to the initial QHEE scenario. Our distinction between highly and weakly hypercritical rates is thus purely numerical and does not imply different physical conditions and/or results. At the highest accretion rate, $100 \, \dot{m}_0$, since the accretion shock expands to small heights we will present results of a 3D MHD simulation while for the other cases only 2D simulations have been performed.

If we ignore convection effects, the timescale required for the quasi-stationary solution to set in is a few sound crossing times, t_cross ≃ r_shock/c_s. For a shock radius of ≃ 50 km and c_s ≃ c/10, this is t_cross ≃ 1–2 ms. The simulations presented here run for hundreds of ms, so this is established quite rapidly. The timescale for convection is of course much longer and will depend on the equilibrium between infall and cooling at the base of the envelope.

3.1. Highly Hypercritical Accretion Rates

We describe here results for highly hypercritical accretion rates, $\dot{m} = (1, 10, 100) \cdot \dot{m}_0$, for which we are able to follow the evolution of the accretion shock. In Figure 3 we plot the radial profiles of density, pressure, velocity, and temperature for 2D simulations in the free-fall scenario at the end of the simulation when the accretion shock and the P and ρ profiles have become stationary. In the density and pressure profiles the piling up of matter close to the neutron star surface is notorious (an effect not accounted for in the analytical approach). Notice that, although the system has reached a quasi-stationary state in all cases, significant noise remains in the velocity profile: because of the periodic boundary condition on the vertical sides, matter can freely flow in the horizontal direction, as well as bounce off the neutron star surface, which prevents a full stationary state from being reached. Nevertheless, the mean velocity profile is close to the analytical one, indicating that the system has largely relaxed despite these fluctuations. In these free-fall simulations, for our fiducial accretion rate $\dot{m}_{0}$ the system reaches the quasi-stationary state in about 600 ms, whereas for higher accretion rates this time is substantially reduced: 300 ms for $10\, \dot{m}_{0}$ and 100 ms for $100\, \dot{m}_{0}$. The $\dot{m}_0$ model of Figure 3 should be compared with Figure 2, which exhibits the analytical profiles and shows a very good agreement between the two. The only visible difference between the numerical and analytical results is the exact position of the shock, Equation (8), in that the numerical results do not exactly follow the $\dot{m}^{-10/63}$ scaling: since y_sh is determined by the energy balance of Equation (6) with the neutrino losses, the larger P just above the neutron star surface in the numerical results implies higher energy losses and, hence, a smaller y_sh.

We show in Figure 4, for a $100\, \dot{m}_{0}$ accretion rate, color maps of density with magnetic field intensity contours¹ superimposed, for the 2D free-fall scenario. The initial panel at t = 0 shows the reverse shock just before it reaches the magnetic field loop. When the reverse shock hits the neutron star surface and bounces, a violent, convective layer appears, and the accretion shock develops and moves rapidly upward against the infalling matter (t = 1 ms panel). The flow smoothens gradually as the accretion shock stabilizes. Instabilities of the Rayleigh–Taylor type are present in this regime, but they disappear when the system reaches equilibrium (see the panel at t = 100 ms). The magnetic loop is immediately torn apart by the reverse shock and the initial strong convection that is dragging the field with it within the developing envelope (panel at t = 1 ms). Although it has a turbulent dynamics, the magnetic field, having its feet frozen into the neutron star surface, always remains confined below the accretion shock. After 50 ms, the magnetic field begins to be submerged and is trapped into the material that is piling up onto the neutron star surface. The matter outside the accretion shock is still falling with constant accretion rate, but it is the fluid inside the envelope, with its much lower downward velocity, that is nevertheless responsible for the magnetic field submergence. After 100 ms, the submergence is completed, with a maximum magnetic field strength ≃ 5 × 10¹² G. At t = 0 the magnetic loop was immersed in a medium of density 10² g cm⁻³, the infalling matter just behind the reverse shock had a density ∼10⁷ g cm⁻³, but after 100 ms the magnetic field has been submerged into matter at density ≳ 5 × 10⁹ g cm⁻³.

Zoom In Zoom Out Reset image size

For comparison, in Figure 5 we show the results of a similar simulation but within the QHEE scenario. At t = 0 the magnetic field loop is immersed within the QHEE. The initial, strongly convective phase seen in the free-fall scenario is absent. Some weak convection is nevertheless present in the upper part of the envelope, whose height initially oscillates (t = 1 ms panel in Figure 5). Due to the smoother evolution, compared to the free-fall scenario, even at this early time the magnetic loop has already been compressed by the infalling envelope but is not strongly disrupted as in the free-fall scenario. After 100 ms, the QHEE has recovered its initial shape and height, with an equilibrium shock radius r_sh ≃ 4 × 10⁶ cm, and the magnetic field is seen to be submerged. Although the accretion velocity within the envelope is much lower than above the shock, it is enough to completely submerge the magnetic loop.

Zoom In Zoom Out Reset image size

In the free-fall scenario the initial smashing of the magnetic loop by the reverse shock is immediately counteracted by the bouncing of the matter off the neutron star surface that drags the field with it, upward, within the forming envelope. The field submergence in the QHEE scenario at early times shows that it is the slow settling of the envelope material, under continuous accretion, that is actually responsible for the field submergence. In the free-fall scenario the field submergence follows the same route as in the QHEE scenario once the envelope stabilizes and its matter slowly settles onto the neutron star surface. The final geometry of the magnetic field is seen to be very similar under both free-fall and QHEE initial conditions.

For the lower accretion rates, 10 and $1 \, \dot{m}_{0}$, the evolution is very similar but the accretion shock reaches much larger heights. We choose to illustrate only the $100\, \dot{m}_{0}$, in Figures 4 and 5 since the lower rate cases involve much higher columns. We show in Figure 6 the radial profiles of the magnetic field for our three highly hypercritical accretion rates, in the free-fall scenario, corresponding to the times at which the quasi-hydrostatic equilibrium has been reached. In the three cases the submergence of the magnetic field is clearly seen. The field is confined in layers where there is still a significant fluid motion (see Figure 3), which explains why it is not totally smashed onto the neutron star surface but is rather hovering above the surface (our surface boundary condition does not allow penetration of the field into the neutron star).

Zoom In Zoom Out Reset image size

It is thus clear from the above results that for strongly hypercritical mass accretion rates, both approaches, the free fall and QHEE initial conditions, lead to the same conclusion, namely, that the field is promptly submerged by the sheer force of the settling fluid in the envelope. Thus we believe this to be a robust result reflecting the actual fate likely to befall an initial magnetic field anchored on the neutron star once a quasi-stationary settling envelope has been established.

3.1.1. 3D Simulation

In all the high accretion rates described above we observe the submergence of the magnetic field. These were, however, 2D geometries and, since the magnetic field behavior is an intrinsically 3D phenomenon, one may wonder how much the addition of the third dimension degree of freedom would alter the results.² To tackle this issue, we performed a 3D MHD simulation for the free-fall scenario. From the available computing power we were limited to simulate only a short column, hence we consider only our highest accretion rate, $100 \, \dot{m}_{0}$.

In Figure 7 we show four snapshots of the time evolution of the magnetic field iso-surfaces until the quasi-stationary state is reached. Notice that the fluid reaches the stationary state in 60 ms, i.e., in a shorter time than in the corresponding 2D simulation. Overall, the evolution is very similar to the corresponding 2D simulations. After 1 ms the magnetic field loop has been smashed by the accretion shock and the field is beginning to be spread by the violent fluid motion. The accretion shock rapidly progresses upward and reaches its final height, r_sh ≃ 4 × 10⁶ cm as in the 2D case, in about 20 ms. Notice, at that time, that the magnetic field is also being dragged upward by the fluid motion but is not able to reach the location of the shock. In no moment will the field be able to rise above the forming envelope. However, when the QHEE has built up and the initial strongly convective regime has died away the magnetic field has been forced downward by the slow settling of matter in the envelope and is submerged inside the denser layers onto the neutron star surface. Notice that the magnetic field has been spread all over the whole bottom area of the simulation box.

Zoom In Zoom Out Reset image size

In Figure 8 we show the radial profiles of the density, pressure, velocity, and magnetic field in the last stage. The density and pressure profiles are essentially identical to the 2D ones shown in Figure 3 with the clear piling up of matter onto the neutron star surface. The velocity profile, however, is much less noisy in 3D than in 2D. These three, ρ, P, and v, profiles are also in close agreement with the analytical solution. The magnetic field strength profile clearly exhibits the submergence within the high-density region. We thus conclude that at this level of mass accretion, no significant differences in morphology are apparent by considering the 3D nature of the problem, and our previous statements concerning the submergence of the field from the 2D calculations remain valid.

Zoom In Zoom Out Reset image size

As a further energy check, the neutrino luminosity integrated over the whole domain, at t = 60 ms when the quasi-stationary state has been reached, is L_ν ≃ 2.3 × 10⁵⁰ erg s⁻¹, but the emission is strongly concentrated within the first two km above the stellar surface due to the strong temperature dependence of the emissivity. The gravitational energy liberated by the accretion is ${\sim } 0.3 \times \dot{m} c^2 \simeq 2.1 \times 10^{50}$ erg s⁻¹ and is hence emitted essentially in its entirety in neutrinos.

Once the quasi-stationary state has been achieved, at t = 60 ms, the adiabatic indices connected with the sound velocity and with the system energy are γ_c = 1.35 and γ_e = 1.34. Then the adiabatic and radiative gradients are ∇_ad = 1–1/γ_c ≃ 0.26 and ∇_rad = (dln T/dln P) ≃ 0.25. The system is, hence, convectively stable, but only marginally.

3.2. Weakly Hypercritical Accretion Rates

We can now investigate what happens with the magnetic loop anchored onto the neutron star surface for lower accretion rates. To answer this, we carry out 2D simulations for the cases $\dot{m} = (10^{-4}, 10^{-3}, 10^{-2},10^{-1}) \cdot \dot{m}_0$. For these lower accretion rates it is computationally impossible to follow the expansion of the accretion shock because it reaches to extremely large heights, leaving the computational domain, and, moreover, the relaxation time of the envelope becomes exceedingly long. Attempting to follow the evolution of a free-fall scenario by allowing the accretion shock to leave the integration domain would be physically, and numerically, inconsistent with our setting: it is not possible to impose an upper boundary condition that incorporates in any reasonable manner the strongly turbulent motion of the initial phase of envelope formation. Due to this, we are forced to restrict ourselves to the QHEE scenarios. The results of our study of the highly hypercritical accretion rates in Section 3.1 showed that the field submergence is occurring in a very similar manner in both types of scenarios and that, in the long term, it is mainly due to the settling of matter onto the neutron star surface. This gives us confidence that limiting ourselves to the QHEE scenario for the weakly hypercritical rates will give us representative results.

We simulated an accretion column with height Δy = 4 × 10⁶ cm, using the upper boundary condition with matter injection at the same $\dot{m}$ and the corresponding ρ and v as given by the analytical solution for a QHEE. The initial magnetic field configuration was that of a magnetic loop, as described in Section 2, immersed within a QHEE envelope. Each simulation exhibited a short initial relaxation phase, but the system rapidly reached stationary density and pressure profiles. In Figure 9 we show color maps of density, with magnetic field intensity contours superimposed, for these four accretion rates at times where a stationary state was well established in the density profile. At the highest accretion rate, $\dot{m} = 0.1 \, \dot{m}_0$, in 100 ms the magnetic loops are clearly submerged within the high-density bottom of the envelope. This result is similar to what we obtained at the highly hypercritical rates in Section 3.1.

Zoom In Zoom Out Reset image size

At lower rates, $\dot{m} \ll 0.1 \, \dot{m}_0$, however, the situation becomes more complex. As seen in Figure 9, the magnetic field is mostly confined in the high-density region but is able to extend to slightly larger heights than before. The evolution of our lowest accretion rate scenario, $\dot{m} = 10^{-4} \, \dot{m}_0$, is depicted in greater detail in Figure 10, a simulation we could pursue for 400 ms (after that time the turbulence reached the upper boundary of our integration domain and conflicted with the inflow condition, so the simulation was stopped). The left panels show that, after an initial strong compression, the field is slowly rising: at 100 ms it is essentially confined within the high-density region but, as matter piles up onto the neutron star surface and this high-density region expands upward, the magnetic field follows it. The right column in Figure 10 illustrates the pattern of convective cells with very high velocities in the upper low density regions that contrast with the much slower motion within the much denser matter accumulating onto the neutron star surface. The size of the convective cells is expected to be of the order of the pressure scale height, H_P = −dr/dlog P, which is ∼10 km in the envelope but only ∼1 km in the forming crust. These two distinct patterns and their length scales are clearly seen in Figure 11, where the absolute value of the vertical velocity is plotted as a function of radius, and are conserved during the evolution in spite of a more than one order of magnitude decrease in velocity. It is this clear difference in the two convection patterns, large, low density cells with high velocities on top of small, high-density cells with low velocities, that keeps the crushed magnetic field confined to the high-density region of the forming new crust.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have continued our study of hypercritical accretion onto newborn neutron stars, in particular with the aim of studying the possible submergence of the magnetic field for various accretion rates. This is particularly relevant in the context of making the neutron star eventually invisible as a pulsar following the supernova explosion. Extending our exploration of parameter space to wide columns spanning a significant fraction of the stellar surface, we have paid special attention to the formation of the quasi-steady atmosphere that forms after the reverse shock has reached the hard surface of the star, and how it affects an initial loop of field anchored in it. We find that whether an initial condition in free fall or an analytical atmosphere in quasi-hydrostatic equilibrium is initially used, the field is rapidly submerged into the forming high-density new crust for strongly hypercritical rates, above our fiducial value of $\dot{M}=\dot{M}_{0}=350$ M_☉ yr⁻¹. The comparison between results for different initial conditions is important because it is extremely challenging to follow the evolution of the system, in particular the shock front, for lower accretion rates. This is due to the scaling of the equilibrium position of the shock, given in Equation (8), that shows the shock expands to very large heights at small accretion rates. We have further carried out a calculation in three dimensions, to investigate whether the flow exhibits qualitative changes, in particular regarding the submergence of the field, as compared to the 2D simulations. Finding none, we conclude that 2D studies are able to correctly capture the main features of the field's behavior in this respect.

As the accretion rate is lowered, reaching down to 10⁻⁴ $\dot{M}_0 \sim 0.035\, M_\odot$ yr⁻¹, the field is progressively less affected by the infall, and the initial loop reaches upward of the neutron star surface after the initial transient, to a height of a few kilometers. Intriguingly, the pattern of convective cells helps explain why the field is not dragged to greater heights. Close to the surface of the star, a high-density region where matter accumulates is present, and this is also where most of the magnetic field remains confined. As described in Section 3.2, there is a sharp difference in the convection patterns in the high-density region, with cell sizes ∼1 km, compared to the overlaying envelope, with cell sizes ∼10 km. Thus, the envelope plasma circulations are unable to reach into the region of high field and drag it to larger altitudes. Notice that even at the lowest accretion rate we considered, the matter pressure is much larger than the magnetic pressure and the field is continuously dragged by the motion of the plasma.

Our distinction between highly and weakly hypercritical accretion rates is purely due to our computational limitations. However, we do find a physical distinction in the field submergence that is simply due to the evolution of the convection in the envelope and the forming crust. In the highly hypercritical cases of Section 3.1 and the highest rate, $0.1 \dot{m}_0$, of Section 3.2 our simulation lasted long enough for convection in the envelope to essentially disappear, while at the lower rates significant convection was always present. In the former cases the only remaining fluid motion is the slow settling of the accreting matter onto the neutron star and, since the magnetic field is frozen into the plasma, the field is completely submerged into the high-density region of the forming crust. This result is still clearly seen in the $0.1 \dot{m}_0$ case of Figure 9. In contrast, in the latter cases continued convection maintains mixing of the magnetic field, but the different patterns of convection in the forming crust compared to the envelope trap most of the magnetic field in the high-density regions, and the field present in the envelope is strongly reduced compared to its initial value. The evolution displayed in the left panels of Figure 10 clearly shows the slow submergence of the field in the high-density region as matter piles up onto the neutron star and the continuous presence of a residual field in the bottom of the envelope in the interface region between the two regions of different convection patterns. Notice that the apparent rise of the magnetic field in Figure 10 is actually just the rise of the height of the building new crust while the residual field never extends into the envelope much beyond the interface region. If this accretion phase last for several hours, most of the magnetic field should eventually be pushed deep into the crust.

The critical factor separating a total versus a partial submergence of the initial magnetic field hence seems to be the timescale for convection in the accretion envelope to disappear. Once convection stops, the magnetic field is unavoidably submerged into the neutron star curst. As described by Fryer et al. (1996), at low accretion rate the development of a stationary non-convective envelope may take more time than the duration of the hypercritical accretion phase. Supernovae in which a strong reverse shock was produced and in which the accreting envelope reached a convectively stable profile are likely to produce a neutron star with a vanishingly small surface magnetic field. On the other hand, supernovae that produced a weak reverse shock whose accreting envelope remained convective can be expected to produce weakly magnetized neutron stars. In the case where no reverse shock developed, the neutron star will simply maintain the magnetic field it had at its birth.

We have extended our previous study (Paper I) of the effect of post-supernova hypercritical accretion on the magnetic field of the newborn neutron star. Our study is still restricted to cases where the accreting matter pressure dominates over the magnetic pressure, so that the magnetic field is only a "victim" of the accretion. We found that the critical factor determining the fate of the magnetic field is whether the developing quasi-hydrostatic equilibrium accreting envelope will reach a state of convective stability or not. Once convection stops the magnetic field is unavoidably pushed into the forming new crust and a non-magnetized neutron star will result. If the accreting envelope remains convective, most of the initial magnetic field is still pushed into the forming crust, but the fluid motions, due to the different patterns of convection in the crust and the surrounding envelope, maintain a residual field protruding out of the star. In this latter case, the resulting external field certainly has a very complicated structure, and anisotropy in the accretion, as, e.g., from rotation, will also increase the complexity of the final field. The dipolar component of the final field may be several orders of magnitude smaller than the initial one, but localized spots with a much stronger field seem unavoidable.

Non-convective envelopes are rapidly obtained only at extremely high accretion rates, larger than ∼10 M_☉ yr⁻¹ in our study. It is not clear how many newborn neutron stars undergo such events, and how many of them would survive without collapsing into a black hole (Fryer et al. 1996). Our conclusion is that there may be only very few young neutron stars that have been totally de-magnetized by hypercritical accretion and they are likely very massive ones.

At lower accretion rates, field submergence is only partial and the fraction of newborn neutron stars having passed through such events is likely significant. Fryer et al. (1996) were very careful in the treatment of neutrino heating and found that when $\dot{M} > 0.1\, M_\odot$ yr⁻¹ a neutrino-driven explosion may result, expelling the envelope, while at lower rates steady accretion is achieved with, eventually, a non-convective envelope. As the accretion rate decreases one may reach a point where the residual magnetic field may begin to play a dynamical role. The final outcome of such evolution is impossible to treat numerically and can only be guessed at now, but a weakly magnetized neutron star, with a complicated field geometry, seems very likely.

The CCOs are the primary candidates for neutron stars having undergone such a field submergence, and they constitute half of the known population of young neutron stars: Popov & Turolla (2012) list a dozen pulsars with age <10 kyr, three of them being CCOs, i.e., a total of nine non-CCOs pulsars, while Gotthelf et al. (2013) present eight confirmed CCOs, all having ages below 10 kyr, plus three candidates, to which one may add the yet undiscovered neutron star in the remnant of SN 1987A. The three CCOs with measurements of spin periods P, and also derivatives $\dot{P}$ from which surface dipolar field strengths B_s can be inferred, are

• 1.  
  PSR J1852+0040 in Kesteven 79, with a period P = 105 ms (Gotthelf et al. 2005) and an inferred B_s = 3.1 × 10¹⁰ G (Halpern & Gotthelf 2010);
• 2.  
  PSR J1210−5226 in PKS 1209-51/52, with P = 424 ms (Zavlin et al. 2000) and B_s = 9.8 × 10¹⁰ G (Halpern & Gotthelf 2011) as confirmed by Gotthelf et al. (2013); and
• 3.  
  PSR J0821−4300 in Puppis A with P = 112 ms (Gotthelf & Halpern 2009) and B_s = 2.9 × 10¹⁰ G Gotthelf et al. (2013).

The unavoidable prediction of a submerged magnetic field is the subsequent growth of the surface field by back-diffusion of the internal field (Muslimov & Page 1995). Hence, burial of the field in a newborn neutron star may manifest itself as a delay in the star turning on as a classical radio pulsar. Possible evolution scenarios were studied in 1D models in Geppert et al. (1999) and Ho (2011) and recently by Viganò & Pons (2012) in a 2D scenario. A growing magnetic field also naturally manifests itself through a braking index inferior to 3 (Muslimov & Page 1996, 1999) and all reliably measured braking indices actually are <3 (Espinoza et al. 2011). Along this line, Güneydaş & Ekşi (2013) contemplate the possibility that a neutron star that receives a large kick at birth may accrete less matter than a slow-moving one. If this is the case, they argue there may be a correlation between the field growth timescale, estimated from the measured value of the braking index, and the observed transverse velocity of the pulsar and find some marginal support for it.

In a different approach, Popov & Turolla (2013) consider high-mass X-ray binaries (HMXBs) that are binary systems in which a neutron star accretes from a massive companion and which have ages inferior to a few millions years, the life expectancy of the massive companion. However, these authors find no robust candidates for low-field neutron stars in HMXBs, in sharp contradistinction to the young (<10⁴ yr) neutron star population. In the case where the initial distribution of magnetic field of neutron stars in HMXBs is similar to the one in the isolated neutron star populations, Popov & Turolla (2013) conclude that the magnetic field of the CCOs must grow on a timescale of 10⁴–10⁵ yr.

In the present work we were still dealing with magnetic field strengths, and accretion rates, in which the magnetic forces are negligible, even though the FLASH code does take them into account. More interesting effects can be expected in the case of a magnetar-sized field (Thompson & Murray 2001) and will be considered in a future paper.

Finally, among the many simplifications of our model, the assumption of non-magnetized accreted matter is clearly disputable. If the neutron star magnetic field is of fossil origin, the accreting flow is likely bringing in some significant magnetic field that may just replace the submerged one. One may conceive that accretion will hence complicate the geometry of the surface field without necessarily reducing its strength, but possibly reducing the strength of the dipolar component. However, turbulence in the accretion envelope can act as a diamagnetic medium (Vaĭnshteĭn & Zel'dovich 1972), which may also result in a strong suppression of the surface magnetic field. Such an evolution certainly deserves further study.

We are grateful to DGTIC-UNAM for allowing us to use its KanBalam Cluster, where all the simulations were performed. The software used in this work was in part developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago. This work was supported in part by CONACyT grants CB-2009-1 #132400 and CB-2008-1 #101958. Finally, we thank the referee for her/his very useful comments.