Repulsion of fallback matter due to central energy source in supernova

The flow of fallback matter being shocked and repelled back by an energy deposition from a central object is discussed by using newly found self-similar solutions. We show that there exists a maximum mass accretion rate if the adiabatic index of the flow is less than or equal to 4/3. Otherwise we can find a solution with an arbitrarily large accretion rate by appropriately shrinking the energy deposition region. Applying the self-similar solution to supernova fallback, we discuss how the fate of newborn pulsars or magnetars depends on the fallback accretion and their spin-down power. Combining the condition for the fallback accretion to bury the surface magnetic field into the crust, we argue that supernova fallback with a rate of dM_{fb}/dt ~ 10^{-(4-6)} M_sun /s could be the main origin of the diversity of Galactic young neutron stars, i.e., rotation-powered pulsars, magnetars, and central compact objects.


Introduction
Compact objects such as neutron stars and black holes are formed in collapsing stars some of which eventually end up with supernova explosions. In general, a fraction of supernova ejecta falls back toward a newborn compact object (Colgate 1971;Zel'dovich et al. 1972;Michel 1988). The accretion rate and its temporal evolution can be determined by e.g., the strength of the supernova shock and the progenitor structure (e.g., ). On the other hand, newborn compact objects can continuously deposit energy into the fallback matter via e.g., neutrino emission, pulsar activity and/or accretion-disk wind (e.g., Piro & Ott 2011). Competition between the fallback accretion onto and the energy deposition from compact objects may result in a variety of outcomes and lead to a diversity of compact objects.
In particular, the dynamics of the fallback accretion may be related to the diversity of relatively young (∼ 1-10 kyr) neutron stars; there are three distinct populations, ordinary pulsars, magnetars, and central compact objects (CCOs) 1 . Based on the spectral and timing observations, they are considered to have different energy sources, the rotation, magnetic-field, and thermal energies, respectively. The most distinctive difference between them is the magnetic-field strength; young pulsars typically have dipole fields of B d < ∼ 10 13 G 2 . Magnetars are considered to be powered by decays of stronger fields of B * > a few × 10 13 G (Olausen & Kaspi 2014) while CCOs have considerably weaker dipole fields of B d ≪ 10 12 G (Pavlov et al. 2001;Park et al. 2006, and so on). Although the diversity of the magnetic field strength can be attributed to conditions before and during the neutron star formation, i.e., the magnetic field strength of the iron core of a progenitor (the fossil scenario; Ferrario & Wickramasinghe 2006) or the magnetic field amplification by a dynamo process (Duncan & Thompson 1992;Mösta et al. 2015), it may be also possible to attribute the diversity to conditions after the formation, i.e., various rates of the fallback accretion onto the neutron star.
If the accretion rate is sufficiently high, the fallback matter can compress and bury the magnetic field into the neutron star crust (e.g., Bernal, Lee, & Page 2010;Torres-Forné et al. 2016). Such a case may correspond to a CCO formation (the hidden magnetic field scenario; Muslimov & Page 1995;Young & Chanmugam 1995). For example, Torres-Forné et al. (2016) obtained the critical accretion rate for burying the magnetic field by considering the competition between the magnetic pressure of the surface field and the ram pressure of the fallback matter. In these studies, however, the energy deposition from the neutron star has been neglected. If the energy deposition rate is sufficiently large, the fallback matter can be repelled before reaching the surface. The critical conditions can be obtained by considering the competition between the fallback accretion onto and the energy deposition from the neutron star. This may define additional bifurcations of the neutron star population.
1 They sometimes show links to another population. For example, some pulsars underwent soft-gamma repeater like bursts (van der Horst et al. 2010), and a CCO was found to exhibit magnetar-like activity (D'Aì et al. 2016;Rea et al. 2016 To this end, we here investigate the dynamics of fallback matter being pushed back by an energy deposition from the central object with newly constructed spherically symmetric self-similar solutions. We consider fallback matter marginally bound by the gravitational field of the central object with a mass M c and a power-law energy deposition rate from the central object; Here L l and l are constants and t is the time measured from the onset of the energy deposition.
Although our solutions are described by a single dimensionless variable where r is the radial coordinate and G denotes the gravitational constant, a variety of density and velocity structures can be realized depending on the strength of the deposition and so on, which, to our knowledge, have not been seen in previous studies.
The structure of the paper is as follows. The next section describes our model. Section 3 presents our self-similar solutions with various values of parameters. In Section 4, we discuss some applications of our solutions to neutron star formation. Section 5 summarizes the results and discusses relations to previous works.

Model
We consider the flow of fallback matter affected by the energy deposition from a central object.
Following Masuyama et al. (2016), we divide the flow into three distinct regions (see Fig. 1): the innermost region (Region I) adjacent to the central object where the energy is deposited uniformly, Region II where the fallback matter is shocked, and the outermost region in which cold matter freely falls due to the gravity of the central object (Region III). Regions I and II are separated by the contact discontinuity specified by the dimensionless variable introduced above as ξ = ξ c . The fluid in Region II has a constant adiabatic index γ and turns into outflow at a certain point depending on ξ c and γ.
The shock front divides the fallback matter into Regions II and III at ξ = ξ s . Such a flow structure across Regions I and II has been reproduced by numerical simulations (e.g., Masuyama et al. 2016).

Energy deposition region (Region I)
The central object is assumed to deposit energy in a uniform spherical region with a radius R c = ξ c (GM c t 2 ) 1/3 . The internal energy E c and pressure P c in this region evolve with time t according to the first law of thermodynamics as Here L l and l are constant. The subscript l of L l indicates that the dimension of this quantity depends on the value of l. This region becomes very hot and the equation of state can be approximated by that for ultra-relativistic gas, that is, E c = 4πP c R 3 c .

Accreted matter (Regions II and III)
Accreted matter can be divided into the shocked fallback matter (Region II) and the surrounding accreted matter (Region III). The motion of the surrounding matter is not affected by the energy deposition but controlled solely by the gravity from the central object. We describe the flow in Region II by using a self-similar solution. We note that the problems treated in this paper involve two constant quantities with physical dimensions GM c and L l .
Here we assume that the spherically symmetric flow depends on time t and r only through t and ξ(r, t). Thus the density ρ, the velocity v, and the pressure p take the forms of as functions of ξ and t. Here we have introduced dimensionless functions D(ξ), V (ξ), and Pr(ξ).
The governing equations are described as ∂ρ(ξ, t) ∂t + ∂r 2 ρ(ξ, t)v(ξ, t) r 2 ∂r = 0, After some manipulations we obtain ordinary differential equations for dimensionless functions D(ξ), V (ξ), and Pr(ξ) as where ′ denotes the derivative with respect to ξ. Eqs. (10-12) can be rewritten as From the second term of the right hand side of equation (13), we can find that the derivative of the density diverges at ξ c where is satisfied while the other derivatives of the pressure and the velocity do not diverge. Thus this point ξ = ξ c defines the contact surface between Regions I and II. To obtain a solution, we numerically integrate Eqs.(13-15) from ξ = ξ s to ξ = ξ c .
The solution for the flow in Region III (ξ ≥ ξ s ) is obtained by ignoring pressure in equations (10) and (11) (or (13) and (14)) as Here a constant D fb denotes D(ξ s ) in the un-shocked flow (Region III). It follows from substitutions of these expressions into equations (4) and (5) that the distributions of the density and velocity in the fallback matter evolve as and that the density approaches a power-law function proportional to r (3l−1)/2 independent of t in the limit of r → ∞. Thus the flow approaches a stationary state for a large r (ξ). The velocity expressed by equation (20) implies that the fallback matter is initially at rest at large distances from the center. Thus the total energy of the flow in Region III is equal to zero while the energy of the flow in Region II becomes positive due to the energy supply from the central source. This indicates that all the solutions presented here describe the flow eventually repelled by the energy supply. Because the density distribution depends on the parameter l that specifies the nature of the central energy source, this density distribution is required to lead to a self-similar solution in which the shock front expanding with the radius proportional to t 2/3 for a given energy source. If the central source is more energetic, which corresponds to larger L l , the fall back matter can be denser, while if the gravity of the central object is stronger, which corresponds to a larger M c , then the fallback matter needs to be more sparse to have an expanding shock front.
We ignore physical processes like neutrino heating and cooling, radiative diffusion, photodisintegration of nuclei, and electron-positron pair production and annihilation, some of which might play crucial roles in the fate of the fallback matter. Instead, we have presented solutions with varying the adiabatic index and the exponent l ofQ in equation (1), which mimic some of the effects of these processes.

Contact surface
We require a condition that the flow has a continuous pressure distribution through the contact surface between Regions I and II where equation (16) holds. From the evolution of the pressure derived from equation (3), this requirement leads to This boundary condition is equivalent to the condition that the energy deposited by the central source is equal to the total energy of the flow inside the shock front. Note here that the accreted matter entering into the shock front carries no energy (the sum of the energy densities of the matter is v 2 (r, t)/2 − GM c /r = 0 from equation 20). We found that the self-similar flow around the contact surface behaves as where D c is a constant. From equation (22), the density distribution near the contact surface drastically changes if γ > 3/2 or not. If γ < 3/2, the density decreases to zero toward the contact surface. If γ > 3/2, solutions with l = (−3 + γ)/[3(γ − 1)] give constant densities at the contact surface between Regions I and II. Otherwise the density at the contact surface diverges to infinity (

Shock front
We assume that a strong shock propagates in the accreted matter. The location of the shock front is specified by ξ = ξ s , where ξ s is a constant. The Rankine-Hugoniot jump conditions at the shock front can be written as Here we have used the fact that the dimensionless velocity of the accreted matter in Region III at the shock front is − 2/ξ s (see eq. 17) and that the dimensionless density at the shock front is denoted by D fb . The value of D fb is determined to satisfy the boundary condition (21) at the contact surface for each given ξ s .
Though it might be possible to extend this shock condition to that with an arbitrary strength as was done in the context of failed supernovae by Coughlin et al. (2018), we take the strong shock limit to simplify the procedure to obtain solutions.

Results
We can find solutions of equations (10)-(12) satisfying the boundary conditions at the shock front and the contact surface for l > −1 and γ > 1. No solution exists in the other range of these parameters; the total deposited energy becomes infinite for l ≤ −1 and the density becomes negative or infinite at the shock front for γ ≤ 1 (see eq. 24). We obtain a solution for any positive value of ξ c . Note that ξ s monotonically increases with ξ c and the solutions in the limit of ξ c → 0 give the minimum ξ s (denoted by ξ s, 0 ). It follows that there are solutions with the maximum possible accretion rates at the shock front when γ ≤ 4/3 for each l. Since the central source is supposed to deposit energy at the rate ofQ, it is convenient to normalize the accretion rateṀ at the shock front withQ as where r s denotes the radius of the shock front andṀ = 4πr 2 s ρ(ξ s , t) 2GM c /r s is the mass accretion rate at the shock front. Another normalization is possible by taking account of the dimension as We show these dimensionless accretion rates, ξ s , and D fb as functions of ξ c for solutions with l = 0 and γ = 6/5, 4/3, 7/5, and 5/3 in Figure 2. Solutions with γ < 4/3 have maxima of the dimensionless accretion rates given by equation (27) or (28) at a finite ξ c . For γ = 4/3, solutions in the limit of ξ c → 0 give the maximum accretion rates. Solutions with γ > 4/3 can sustain any accretion rates. In other words, the accretion rate monotonically increases to infinity as ξ c → 0 in these solutions. At the same time, the value of D fb becomes infinite in this limit. These characteristic values are listed in Table 1 for some γ and l. Table 1 shows results of solutions with various l. The energy deposition rate usually decreases with time if a single physical process is involved and thus described with negative l. Such a deposition rate might mimic the declining phase of a sudden energy release due to a glitch in the crust followed by the transport of energy toward the surface (Eichler & Cheng 1989;van Riper et al. 1991). the contact surface has a sign opposite to that of the pressure gradient (see the left panel of Fig. 3).
This instability tends to bring dense matter toward the contact surface. Thus the central activity may fail to repel the fallback matter even if γ > 4/3.
Solutions with the maximum accretion rates are shown in Figure 4. These solutions have negative velocities in some part of Region II because of their small ξ s 's. From equation (25), the velocity immediately behind the shock front becomes negative when ξ s < [3(γ − 1)] 2/3 /2. If γ = 4/3 (γ = 6/5), this criterion becomes ξ s < 1/2 (ξ s < (5/3) 2/3 /2 ∼ 0.356). On the other hand, the solutions with γ = 5/3 presented in Figure 3 happen to have large ξ s 's that do not satisfy this criterion (ξ s < 1/2 1/3 for γ = 5/3), thus the shocked matter flows outward while solutions with γ > 4/3 and negative l have ξ s, 0 greater than the upper limit of the criterion as seen in Table 1. Shocked matter in a solution with a positive l tends to flow inward due to the weak central engine in the early phase. Because the contact surface always moves outward in all the solutions (see eq. (16)), the pressure increased by the central energy deposition repels the fallback matter somewhere in Region II even for these solutions with great accretion rates. In solutions with small accretion rates, the energy deposition can repel the fallback at the shock front.
The kinetic energy E k (t) of the repelled ejecta evolves with time t following the equation Since the value of the dimensionless integral is of the order of unity here, the factor in front of the integral gives a rough amount of the kinetic energy.
After a successful supernova shock propagates through the progenitor star, a bulk of the stellar material is ejected while a minor fraction falls back. The accretion rate and its temporal evolution can be determined by the strength of the supernova shock and the inner structure of the progenitor. For example,  estimated the fallback rate based on one-dimensional numerical simulations of neutrino-driven explosions. The fallback typically starts when the neutrino luminosity significantly decreases, i.e., t fb > ∼ 10 s after the core bounce, and the accretion rate subsequently decreases as ∝ (t/t fb ) −5/3 for t > ∼ t fb . The total mass M fb of the fallback matter ranges from ∼ 10 −4 M ⊙ to ∼ 10 −2 M ⊙ . More matter falls back in a more massive star. Correspondingly, the peak mass accretion rate is estimated to beṀ fb < ∼ 10 −(3-5) M ⊙ s −1 . In general, the maximum fallback accretion rate that can be repelled by the central energy source is described with M c , r s , (4πD fb √ ξ s ), andQ from equation (27). Since we consider a fallback accretion onto a neutron star, we take M c = M * ∼ 1.4 M ⊙ and set the radius of the neutron star as R * ∼ 12 km. When the fallback accretion starts, the position of the shock will be r s ≈ ξ s (GM * t fb 2 ) 1/3 , or r s ∼ 4.2 × 10 9 cm ξ s t fb 20 s In the critical situation corresponding to the maximum mass accretion rate, the radius of the inner contact surface should be set as the neutron star radius, i.e., ξ c,crit ≈ R * /(GM * t fb 2 ) 1/3 , or ξ c,crit ∼ 2.8 × 10 −4 t fb 20 s The temperature in the shocked region is typically high and the radiation pressure dominates so that the adiabatic index will be slightly larger than 4/3. Combining this fact with equation (31) Hereafter we will assume l = 0 which is typically valid for t ∼ t fb > ∼ 10s since the timescale t sd of the spin-down is estimated to be much longer than that of the fallback from the following formula: where I denotes the moment of inertia of the neutron star, Ω(= 2π/P ) the angular frequency (P the spin period), B * the magnetic field at the magnetic pole on the surface, and c denotes the speed of light. Finally, the maximum energy deposition rate from a rotating neutron star can be described aṡ where µ = B * R * 3 is the magnetic dipole moment and R lc = c/Ω is the light cylinder radius. We note that the above energy deposition rate is different from the classical dipole spin-down rate, The additional factor (R lc /R * ) 2 represents the enhancement of spin-down luminosity (Parfrey, Spitkovsky, & Beloborodov 2016); the magnetic fields are maximally open, like a split monopole, due to the accretion. From equations (30-34), we obtain the critical accretion rate aṡ IfṀ fb < ∼Ṁ crit,repul , a bulk of the fallback matter will be directly repelled by the spin-down power. Otherwise, it is accreted on the neutron-star surface.
IfṀ fb > ∼Ṁ crit,repul , the fallback matter reaches the neutron star surface. If the accretion rate is so large, then the surface magnetic field can be buried and the spin-down power is significantly reduced (e.g., Bernal, Lee, & Page 2010;Torres-Forné et al. 2016). A necessary condition for burying the magnetic field is set by the balance between the magnetic pressure and the ram pressure at the surface; In the case ofṀ crit,repul < ∼Ṁ fb < ∼Ṁ crit,bury , the situation will be more complicated. At first, the fallback matter can be accreted on the neutron star; the magnetic fields and spin-down energy are confined in a near surface region. As the fallback rate decreases with time, large-scale fields emerge and the spin-down power pushes back the fallback matter. In Figure 5, we show how the consequences of a fallback accretion depend on B * and P forṀ fb = 10 −6 M ⊙ s −1 (left), 10 −5 M ⊙ s −1 (center), and 10 −4 M ⊙ s −1 (right).
Let us now discuss possible connection between the diversities of neutron-star formation with fallback accretion and the observed young neutron stars. From Figure 5, the conditionṀ fb <Ṁ crit,repul 3 There is a typo in their eq. (25). is always satisfied for fast-spinning strongly-magnetized neutron stars with B * > ∼ 10 13 G and P < ∼ a few ms. Such cases have been proposed as a plausible central engine of extragalactic transients like gamma-ray bursts, superluminous supernovae, and fast radio bursts (see e.g., Metzger et al. 2015;Kashiyama et al. 2016;Kashiyama & Murase 2017;Margalit et al. 2018 and references therein) 4 . A similar range of B * and P has been also considered in the context of Galactic magnetar formation; the magnetic field amplification can be attributed to the proto-neutron-star convection coupled with a differential rotation (Duncan & Thompson 1992;Thompson & Duncan 1993) or the magneto-rotational instability (e.g., Mösta et al. 2015). Note, however, that so far there is no observational support to the dynamo scenario (e.g., Vink & Kuiper 2006).
Neutron stars with relatively weak magnetic fields (B * < ∼ 10 13 G (Ṁ fb /10 −5 M ⊙ s −1 ) 2/3 ) and slow spins (P > ∼ a few 10 ms) satisfy the conditionṀ fb >Ṁ crit,bury . Both the magnetic field and spin-down power are strongly suppressed by the accreted matter. As long as the situation continues, the thermal energy stored in the neutron star will be the main source of the emission. Such neutron stars will be observed as CCOs as proposed by Torres-Forné et al. (2016).
The boundary between the former and latter cases is set by equation (35), i.e., (B * /10 13 G) × (P/10 ms) −1 ∼ (Ṁ fb /10 −4 M ⊙ s −1 ) 1/2 , which might correspond to the boundary between rotationpowered pulsars and CCOs. For example, the initial magnetic field and spin period of the Crab pulsar have been estimated to be B * ∼ 10 13 G and P = a few 10 ms, respectively, and consistent with the above point ifṀ fb < ∼ 10 −5 M ⊙ s −1 . Such a relatively small fallback accretion rate might be also consistent with a relatively low-mass progenitor inferred for the Crab pulsar.
As above, we show that the observed diversity of Galactic neutron stars could be connected to (B * , P,Ṁ fb ) at their birth. We speculate that ordinary pulsars, CCOs, and magnetars originate from the blue, red, and yellow shaded regions in Figure 5, respectively. The fact that the three classes have a comparable population (e.g., Keane & Kramer 2008) can be naturally explained since the typical parameters at their birth, B * ∼ 10 13 G and P ∼ a few 10 ms, roughly coincide with the intersection of the boundaries of three regions. We should note, however, that the self-similar solutions in Sec.
3 and calculations by Torres-Forné et al. (2016) are one dimensional. Multi-dimensional effects in fallback accretion have to be taken into account consistently to make the above criterions more quantitatively accurate. For instance, as mentioned in section 3, solutions with l = 0 might be subject to the Rayleigh-Taylor instability. Nevertheless, multi-dimensional effects are not expected to change the above criterions because the energy of the shocked flow in Region II is always positive. Of course, one needs to perform multi-dimensional calculations to know detailed outcomes of this instability. We should also note that, in order to connect their states at birth to observational properties at t age > ∼ a few 100 yrs, the long-term evolution of spins and magnetic fields become important. We will investigate these topics elsewhere.

Summary and Discussion
We have presented series of self-similar solutions for fallback matter being shocked and repelled by the energy deposition from a central object. The behavior of the solutions changes depending on the adiabatic index γ in the shocked region. For γ > 4/3, we can find a self-similar solution for an arbitrarily high fallback accretion rate by taking the radius of the contact surface correspondingly small. On the other hand, for γ ≤ 4/3, there are upper-bounds to the accretion rate at the shock front The critical values of 4πD fb √ ξ s , ranging from ∼ 1−100, are shown in Table 1. Note that 4πD fb √ ξ s = 1 corresponds to the cases where the accretion luminosity GṀ M c /r s is equal to the energy deposition rate from the central objectQ. ForṀ >Ṁ crit , the fallback matter will plunge into the central object.
We have applied the self-similar solution to neutron star formation with fallback accretion when the neutron star deposits energy due to its spin-down power. The solution suggests that the critical mass accretion rate above which the fallback matter can be accreted on the surface is given as a function of the initial spin rate and surface magnetic filed. It is shown that this condition together with the criterion that the accreted matter can bury the magnetic field and suppress the spin-down power may determine the fate of the newborn neutron star to become a magnetar, a pulsar, or a CCO.

Relation to Previous Works
We here discuss relations between the self-similar solution constructed in this paper and previous studies on the dynamics of fallback accretion.
The shock propagates outward with the radius proportional to t 2/3 in our models. A similar behavior of the shock front was obtained from numerical computations for fallback presented in Zel'dovich et al. (1972). Thus our model can capture some features of realistic models that treat the energy source originating from gravitational energy of accreted matter. Zel'dovich et al. (1972) discussed the accretion of up to 10 −5 M ⊙ assuming a static power-law density distribution proportional to r −3/2 (r denotes the radial coordinate) as the initial condition, which results in a constant accretion rate. They computed hydrodynamics of the fallback matter including radiative diffusion of photons and neutrino emission. Colgate (1971) emphasized the role of neutrino emission in driving intense fallback motion. These investigations focused on the early phase affected by the unknown explosion mechanism of a core collapse supernova. Chevalier (1989) discussed the fallback phenomenon after a reverse shock wave reaches the neutron star. He argued that the accretion rate of the bound debris decreases with time t following a power law as t −5/3 provided that the mass fraction of the debris is uniformly distributed in binding energy. This power law evolution of the accretion rate had been already pointed out from a dimensional analysis (Michel 1988). Though Chevalier (1989) started his computations from the uniform motion of uniform matter to obtain this accretion rate, he pointed out that this situation can be realized if the density ρ d (r) of the debris initially has a static power-law distribution ρ d (r) ∝ r −4 (Sakashita & Yokosawa 1974). He also discussed the effects of the pulsar activity on the fallback and concluded that the effects are negligible if only the magnetic pressure is considered.
On the contrary, newborn spinning neutron stars should deposit energy in the surrounding matter as a result of the activities originating from rotating strong magnetic fields. This interaction has been discussed with self-similar solutions in 1D spherical symmetry (Chevalier 2005;Suzuki & Maeda 2017) and numerical simulations in 2D axi-symmetry (Chen et al. 2016;Suzuki & Maeda 2017;Blondin & Chevalier 2017) or 3D (Blondin & Chevalier 2017). All of these studies have dealt with interactions between pulsar winds and expanding supernova ejecta without fallback matter. In reality, the innermost part of once ejected matter falls back and the energy deposited by the pulsar wind should affect the motion of the fallback matter, which is important for determining the nature of young neutron stars as discussed in section 4. Our solutions cover some of these aspects though we simplified the situation by introducing the power law energy deposition rate and assuming a specific motion of the fallback matter above the shock front, which is different from the assumptions in the above previous work. Thus our solutions have a different temporal evolution of the mass accretion rate. Table 1. Critical dimensionless quantities in solutions with some l and γ. For γ = 5/3 and 7/5, values of ξ s, 0 of the solutions in the limit of ξc → 0 are listed. For γ = 4/3, in addition to ξ s, 0 , the corresponding eigenvalue of D fb (eq. (24)) as well as the corresponding maximum values of the two accretion rates (eqs. (27) and (28)) are listed. For γ = 6/5, values of ξs (ξ s, 1 and ξ s, 2 ) that give the maximum values of the two accretion rates and the corresponding values of D fb .