Constraining the long-lived supramassive neutron stars by magnetar boosted kilonovae

Kilonovae are optical transients following the merger of neutron star binaries, which are powered by the r-process heating of merger ejecta. However, if a merger remnant is a long-lived supramassive neutron star supported by its uniform rotation, it will inject energy into the ejecta through spindown power. The energy injection can boost the peak luminosity of a kilonova by many orders of magnitudes, thus significantly increasing the detectable volume. Therefore, even if such events are only a small fraction of the kilonovae population, they could dominate the detection rates. However, after many years of optical sky surveys, no such event has been confirmed. In this work, we build a boosted kilonova model with rich physical details, including the description of the evolution and stability of a proto neutron star, and the energy absorption through X-ray photoionization. We simulate the observation prospects and find the only way to match the absence of detection is to limit the energy injection by the newly born magnetar to only a small fraction of the neutron star rotational energy, thus they should collapse soon after the merger. Our result indicates that most supramassive neutron stars resulting from binary neutron star mergers are short lived and they are likely to be rare in the Universe.


INTRODUCTION
Kilonovae (KNe, also called macronovae) are bright optical events that occur after the merger of a binary neutron star (BNS) systems (Li & Paczyński 1998;Metzger et al. 2010, see Rosswog 2015;Tanaka 2016;Fernández & Metzger 2016;Metzger 2019 for reviews), serving as the optical counterparts to gravitational wave (GW) sources.They arise from the thermal radiation emitted by the hot matter ejected during the BNS merger.The thermal energy of the ejected material originates from the radioactive decay of heavy elements produced through the r-process nucleosynthesis (Burbidge et al. 1957;Cameron 1957) which happens in a neutron-rich environment.To first order approximation, the evolution of a KN can be treated as an isotropic expanding hot ejecta.The ejecta is initially optically thick due to the bound-bound absorption (i.e., the line forest) by the r-process elements, but gradually gets transparent as it expands, resulting in a peak in the light curve.The spectrum of the emitted radiation, which can be approximated as thermalized emission, typically peaks in the optical or near-infrared wavelengths.Clear KN emission signatures were first observed as an electromagnetic counterpart of the notable event GW170817: the merger of a BNS system detected by the Laser Interferometer Gravitational-Wave Observatory (LIGO) (Abbott et al. 2017a,c).The observations mostly match with the theoretical modeling, and the recognition of Lanthanide el-★ E-mail: wang4145@purdue.edu† E-mail: dgiannios@purdue.eduements in the spectrum confirms the r-process heating as the energy source (Cowperthwaite et al. 2017).Together with the prompt GRB (Abbott et al. 2017d), its afterglow observation, and the host galaxy, GW170817 has been extensively applied in the research of physics and astrophysics, such as the neutron star matter equation of state (Abbott et al. 2018), GRB afterglow physics (e.g., Gill & Granot 2018;Lazzati et al. 2018;Margutti et al. 2018;Kathirgamaraju et al. 2019;Troja et al. 2019;Wu & MacFadyen 2019;Beniamini et al. 2020;Nathanail et al. 2020;Nakar & Piran 2021), cosmology (Abbott et al. 2017b;Hotokezaka et al. 2019;Wang & Giannios 2021) and fundamental physics (Wang et al. 2017).
While the brightness of a KN is inherently limited by the radioactive energy of the ejected material (approximately 10 46 erg, e.g.Metzger 2019), there is a possibility of augmenting their luminosity through a hypothesized energy source originating from a central remnant that remains active after the merger event (Yu et al. 2013;Metzger & Piro 2014;Kisaka et al. 2016).One such example is a millisecond magnetar.If a remnant of the merger persists due to rapid uniform rotation (rigid body rotation), its rotational energy could potentially reach levels up to a few 10 53 erg (Margalit & Metzger 2017;Radice et al. 2018) limited by the Keplerian rotation (also known as the mass-shedding limit).At this stage, the neutron star is referred to as a supramassive neutron star (SMNS) since its mass exceeds the maximum allowed mass of a static neutron star, known as the Tolman-Oppenheimer-Volkoff mass ( TOV ).It is believed that a hypothetical SMNS formed from a BNS merger is likely to also be a millisecond magnetar whose dipole magnetic field ranges from 10 14 to 10 16 G, where the upper limit is bounded by the stability of magnetized NS (e.g., Akgün et al. 2013), and the lower limit is caused by the amplification of magnetic fields during the differential rotation phase of the central remnant following the merger (e.g., Price & Rosswog 2006).A millisecond magnetar spins down and losses energy through magnetic dipole radiation.The majority of this released energy is transferred into the surrounding environment by the magnetar wind.If a fraction of this energy can be deposited into the ejecta as thermal energy, it has the potential to significantly enhance the luminosity of a KN -by more than two orders of magnitude (Metzger 2019), depending on the model.This enhanced luminosity enables detection at distances exceeding that for regular KNe by more than an order of magnitude, corresponding to a detectable volume more than three orders of magnitude greater than that of regular KNe.In this study, we refer to these exceptionally bright, and as of yet hypothetical transients, as magnetar-boosted KNe.Recently, works have argued that their luminosity can be reduced if the ejected material is Poynting flux dominated (Ai et al. 2022), or if the ejection is not isotropic (Wang et al. 2023).However, such a scenario is not considered in this work since the magnetic fields in the magnetar wind are mostly dissipated in our model (this will be explained in §2.3).In this paper, since we only care about a magnetar produced after a neutron star merger, we use the terms "magnetar" and "SMNS" interchangeably.Readers should not confuse it with the magnetars as remnants of single-star stellar evolution.
Despite the fact that the occurrence rate of magnetar-boosted KNe may constitute only a small fraction of the overall population of binary neutron star mergers, their detectability can still remain substantial due to the considerably larger detectable volume as compared with regular KNe.Numerous ground-based optical telescopes, such as the Zwicky Transient Facility (ZTF) (Bellm et al. 2019) and the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS) (Kaiser 2004), have been actively surveying the sky for rapidly evolving transients.Additionally, several upcoming optical telescopes, including the Vera C. Rubin Observatory (Ivezić et al. 2019), are ready to start their operations in the near future.However, over the past several years of sky surveys, no confirmed KNe have been reported (Andreoni et al. 2020(Andreoni et al. , 2021)).
The absence of detection provides a significant constraint on the characteristics and rates of magnetar-boosted KNe, specifically addressing the question of why they are so rare.One potential explanation lies in the formation rate of SMNS.It is possible that the occurrence of long-lasting SMNS is an exceptionally uncommon outcome of BNS mergers.The fate of a BNS merger remnant is determined by factors including the equation of state (EoS) of neutron star matter, the initial rotation speed during uniform rotation (if applicable), and the mass of the remnant.Depending on various conditions, four possible scenarios can arise, ranked here in order of decreasing remnant mass.Firstly, if the remnant is excessively massive, it will promptly collapse into a black hole without undergoing an intermediate stage.Secondly, if the remnant survives sudden collapse, its inner angular momentum will rapidly dissipate and redistribute through differential rotation.At this stage, the remnant is known as a hypermassive neutron star (HMNS).A HMNS may collapse into a black hole if the centrifugal force can't balance the gravity when it slows down.Thirdly, if the remnant remains stable against collapse after it enters uniformly rotating phase, it becomes a temporarily stable SMNS.Lastly, if the remnant's gravitational mass at rest remains below the Tolman-Oppenheimer-Volkoff mass ( TOV ), it becomes a indefinitely stable neutron star.The boundary between these scenarios relies on the aforementioned conditions, but the EoS and the statistical properties of the remnant's rotation and mass are still not well understood.Considering the lower limit of  TOV constrained by most massive pulsars (Antoniadis et al. 2013;Fonseca et al. 2021), assuming progenitors follow the mass distribution of Galactic neutron stars, and assuming that SMNS initially rotates at Keplerian speed, the recent work by Beniamini & Lu (2021) suggests that a non-negligible fraction of BNS remnants would result in long-lived SMNS.Consequently, the absence of detection should place stringent constraints on these assumptions.
Indeed, both observational and theoretical studies have indicated that the long-lived remnants are likely to be very rare.Late time radio observations of sGRBs have so far not shown evidence of a persisting radio source (Metzger & Bower 2014).Recently, Beniamini & Lu 2021 have found that the long-lived magnetar model are inconsistent with the signatures of X-ray plateaus found in sGRB afterglow, as well as the lack of bright sources in blind radio sources (for the latter point see also earlier predictions by Metzger et al. 2015).Margalit et al. 2022 performed numerical simulations of neutron star mergers and found that the core of the remnant will collapse into a black hole even if the remnant's total mass and angular momentum allows the formation of a temporarily stable SMNS, since the core slows down much faster than the "disk".Motivated by these studies, a similar constraint should be made by the aforementioned optical survey, provided that a boosted KNe model is well established.
To accurately predict the signatures of boosted KNe, one needs to carefully study the interaction between the magnetar wind and the ejecta.The energy injection efficiency should be calculated based on the interaction, rather than assuming a free efficiency parameter.A detailed calculation was carried out by Metzger & Piro (2014) (referred to as MP14).They considered the efficiency by incorporating a model involving a pulsar wind nebula (PWN) obstructed by an ejecta wall.In this model, the PWN is inflated by the magnetar wind, while the ejecta wall consists of the r-process elements ejected during the merger.The ejecta is photoionized and heated by the X-rays emitted from PWN.Within the PWN, ultra-relativistic pairs emit gamma-rays through synchrotron radiation and inverse Compton scattering.The gamma-rays subsequently annihilate with background photons and generate additional ultra-relativistic pairs, initiating the, so called, pair cascade.Due to the small size of the PWN constrained by the ejecta wall, the cascade becomes saturated, resulting in a fraction of ∼ 10% of spindown power turning to the rest mass of pairs in the PWN (Svensson 1987).Consequently, the PWN becomes highly opaque to Thomson scattering, and a significant amount of energy injection eventually turns to the kinetic energy of the ejecta through the pdV work (Metzger & Piro 2014).According to this model, the luminosity enhancement is considerably suppressed as compared with the energy input from the central engine.Nonetheless, they find a magnetar-boosted KN luminosity that is still more than two orders of magnitude brighter than a regular one.Correspondingly, the model predicts a detectable volume that is more than three orders of magnitude larger than that of regular KNe, which is in contrast with the lack of detections.It should be noted that in this model, the assumption is made that the magnetars are indefinitely stable.To further investigate the constraints of the rate of SMNS implied by observations, a more detailed investigation of this model, including the photoionization processes and a limited survival time for SMNS, may be necessary.
The model can also be improved by considering the Rayleigh-Taylor instability of the PWN-ejecta interacting surface, which arises due to the high acceleration and density difference.If this instability occurs, a significant portion of the matter in the PWN may escape from the ejecta, resulting in the formation of an ultra-relativistic blastwave.This blastwave propagates through the interstellar medium, accelerates electrons, amplifies microscopic magnetic fields, and generates synchrotron radiation, just as in the case of a GRB afterglow.However, unlike a GRB, the blastwave in this scenario is isotropic rather than confined to a narrow jet angle.Considering the substantial energy budget of the SMNS and the isotropic nature of the blastwave, such radiation might also be observed through sky surveys, and would be classified as a, so called, "orphan afterglow".Such events haven't been robustly identified, further constraining the formation rate of long lived rapidly rotating magnetized NS remnants of BNS mergers.
In this work, we present a refined model of magnetar-boosted KNe building upon the framework established by MP14, incorporating additional physical details.Specifically, we incorporate a limited survival time for the SMNS and a more detailed photoionization calculation of the ejecta.We also explore the potential occurrence of Rayleigh-Taylor instability and its afterglow-like radiation.Using this model, we perform an EoS-independent study to assess the observational potential of such remnants across the parameter space.We also perform a EoS-dependent simulation to study the detection rate, starting from a population of BNS mergers and incorporating the observations made by ground-based optical telescopes.Our results indicate that, in order to be consistent with observations, most SMNS can't be long-lived, suggesting that SMNS, as the merger remnant of BNS, are exceedingly rare in the Universe.
This paper is organized as follows.We describe the details of our model in section §2 and discuss its observational features and stability.In §3 we discuss the observational features and prospects of our model.In §4 we perform the EoS-independent and EoSdependent study, compare it with current observations, and make constraint on the merger remnants.In §5 we discuss the implications on related topics of our model.Finally, we summarize the points and conclude our study in section §6.

MODELING THE MAGNETAR-BOOSTED KILONOVAE
In this work, we consider a system consisting of two distinct regions: the inner PWN and the outer ejecta.A schematic representation of the system is illustrated in Figure 1.The PWN, predominantly composed of electron-positron pairs and X-rays, is inflated by the spindown power of the magnetar.It's surrounded and trapped by the ejecta wall that consists of r-process elements.Initially both PWN and ejecta are optically thick, and most of the internal energy converts to kinetic energy of the ejecta through pdV expansion rather than being radiated away.The x-rays diffuse from the PWN, photoionize and heat the ejecta, and are able to break out once the ejecta is fully ionized.The hot ejecta produces the observed thermal radiation, i.e., the KN.When the magnetar collapses, the PWN loses its energy supply and rapidly disappears due to pair annihilation, leaving an expanding thermal ejecta.We describe the details of the above process in the following parts.

The basic assumptions
To simplify the calculation, we build a toy model based on the following assumptions: (i) The ejecta has a uniform (but time evolving) density  ej .
(ii) To balance the pressure at the interface between PWN and ejecta, we assume a uniform and radiation dominated pressure throughout the PWN-ejecta system.
(iii) The expansion is homologous.In other words, the velocity  at the radius  follows  ∝ .

X-ray ionization front
Figure 1.Illustration of the structure of the post-merger system.The system is composed by an inner pulsar wind nebula (PWN) inflated by a magnetar, and an outer ejecta shell.The PWN is composed by X-ray photons and electronpositron pairs.The X-ray radiation ionize and heat the ejecta, leading to a boosted luminosity of KN.The evolution of the X-ray opacity in the ejecta can be characterized by a approximate ionization front (red-dashed line).The X-rays breakout from the ejecta once the ionization front reaches the ejecta surface.
The uniform pressure and the homologous condition result in following relations where  n and  ej are internal energy of nebula and ejecta, respectively.The radii of the PWN  n and the ejecta  ej are measured from the magnetar.Because most of energy is trapped in the PWN due to the high opacity, the above equation implies that  n ∼  ej , i.e., the ejecta shell must be thin.In fact, this is a natural result of pressure balance.The shell thickness is Δ sh =  ej −  n .Assuming a uniform density  ej in the ejecta, the kinetic energy of the ejecta is where  ej is the mass of the ejecta.The kinetic energy of the PWN is not considered because it's much lighter than ejecta.We also don't consider relativistic effects for the bulk motions in our model.In an extreme case, if the magnetar has a rotational energy  ini = 10 53 erg and the ejecta has a mass of  ej = 0.01 ⊙ , the ejecta is accelerated to mildly relativistic speed.However, we don't expect a significant modification to our results because: (i) In the case a SMNS forms, most of the ejecta mass comes from disk wind ejecta with a mass of approximately 0.1 ⊙ (e.g., Margalit & Metzger (2019)), which cannot be easily accelerated to relativistic speeds.(ii) The rotational For the readers' convenience, we list the symbols of geometric and thermodynamic variables in Table 1.

Magnetar
Assuming magnetic dipole radiation, the spindown power of the magnetar is (3) Note the spindown power index may not be precisely -2 when considering the variation of the moment of inertia during the early evolution of a fast spinning SMNS.However, this approximation provides a reasonable approximation of the spindown process.The initial spindown power is assumed to follow the magnetic dipole power which depends on the magnetar dipolar magnetic field  at the pole, the initial angular velocity Ω 0 and the neutron star radius  NS .In this work we parameterize the power by the initial rotational energy  ini instead of the angular velocity, because this determines the energy budget available to the KN system. ini is given by  ini = 1 2 Ω 2 0 with the moment of inertia estimated by  = 2 5  NS  2 NS .Assuming  NS ≈ 3 ⊙ and  NS ≈ 10 km, the initial spindown power becomes  sd,0 ≈ 10 50  ini 10 53 erg 2  10 15 G 2 erg/s. (5) Provided the initial rotational energy  ini , the spin-down time scale  sd can be calculated by  sd =  ini / sd,0 .Inserting the above equations, we have If the magnetar is indefinitely stable, or if the spin-down timescale is shorter than the KN peak time, its entire rotational energy is extracted and becomes available for enhancing the KN.However, in cases where the initial rotation speed at the onset of the uniform rotation stage is significantly smaller than the Keplerian speed, or if the magnetar is so massive that it collapses on a finite timescale  c , only a fraction of the energy will be accessible for the KN.We define this fraction as the energy extraction ratio .Its value is calculated as follows: The above discussion leads to a parametrization of spindown luminosity by the magnetic field , initial rotational energy  ini , and the energy extraction ratio .We list these symbols related with the magnetar in Table 2.

PWN
The magnetar continuously injects ultra-relativistic pairs ( ≫ 10 4 ) into the system, leading to the formation of a PWN.The pairs quickly cool down through synchrotron radiation and inverse Compton scattering, producing high energy gamma-ray photons.The gamma-ray photons are able to annihilate with background photons which further produces ultra-relativistic pairs, resulting in a pair cascade.The extent of this cascade is characterised by the compactness parameter ℓ ≡  T  sd /( n  e  3 ).When ℓ ≫ 1, the pair cascade is saturated.As a result, around  ≈ 0.1 of spindown power turns to rest mass energy of the pairs.Consequently, the PWN is predominantly composed of low-energy electrons and non-thermal photons.The spectrum of the photons extends from the background photon energy to the pair annihilation threshold ∼ 1 MeV with a power-law index of -1 (Svensson 1987;Ghisellini 2013;see Vurm & Metzger 2021 for a more detailed calculation in the case of superluminous supernova).Since the PWN is trapped behind the ejecta wall, in this work, we approximate the background photon energy by the typical thermal photons energy of the ejecta, i.e., 3   ej , where   is the Boltzmann constant and  ej is the temperature of the ejecta.Note that because the energy density of the system is uniform, the PWN and the ejecta should share a common temperature where  is the Stefan-Boltzmann constant.
The pair density  ± (counting both electrons and positrons) is estimated by balancing the pair production and pair annihilation rate.The pair production rate can be estimated by saturated pair cascade  + =   sd /( e  2  n ).Balancing with the pair annihilation rate  − = 3 T  2 ± /16, the pair density is calculated by The timescale to reach this equilibrium is   ≃ 16 n /3 n (Metzger & Piro 2014), which is typically shorter than the evolution timescale.
We have tried to incorporate a dynamical pair density considering  ± / =  + −  − , but we find no practical difference as compared with the balanced value.
The radiation of the PWN can be estimated using the photon diffusion timescale The diffusion timescale,  d n , is a smooth interpolation between the optically thin and thick cases.Now the luminosity of the PWN can be estimated by Considering the shape of the spectrum, we can estimate the frequency dependent luminosity After the survival timescale  c , the magnetar collapses into a black hole.The sudden termination of the ultra-relativistic pair supply ceases the pair cascade.However, it may take some time before this information propagates to the nebula surface, whose speed can be approximated by the sound speed.Since the PWN material is a relativistic fluid, the sound speed is   = / √ 3, which can only be well-modeled by considering the relativistic effects.This is out of the scope of our work.To simplify the scenario, we simply assume that this information is instantaneously transmitted across the PWN, so the pair density directly starts dropping following the annihilation rate.At the same time, we also turn the spindown power  sd to 0. In this approximation the PWN quickly disappears after the collapse because it quickly becomes transparent.Although our treatment exaggerates the effects of the collapse, it is unlikely to strongly affect the peak luminosity, because  c is generally much earlier than the time at which the ejecta becomes transparent, at which stage the physical evolution of the PWN is hidden by the surrounding ejecta wall.From the observer's perspective, the central engine is active for a very short timescale, so the ejecta appears as if it undergoes an instantaneous energy injection, where the specific details of the injection process are no longer important.
The symbols related with the radiation of PWN (as well as ejecta to be discussed next) are listed in Table 3.

Ejecta
The composition of the ejecta is rather complicated and is subject to numerical study.The main challenge here is the modelling of photoionization of the ejecta, which requires the knowledge of the bound-free cross sections of r-process elements.However, due to the relatively short half-decay timescales of these elements and their isotopes, it is difficult to measure these values in ground-based laboratories.Currently, the available atomic data for the heaviest elements is iron-56.Moreover, in the situation where a long-lived magnetar is present, it will strongly irradiate the disk outflows with neutrinos, which tends to increase the electron fraction of the material to values ≳ 0.3 (Lippuner et al. 2017).As a result, the ejecta will be mostly composed by light r-process elements whose electron structure is similar to iron-56.Therefore, in this work, we follow MP14 and assume the ejecta is iron like.While this is a crude estimation, we will demonstrate in the following section that the process of thermalization does not play a dominant role in determining the luminosity of magnetar-boosted KNe.
Similar to PWN, the ejecta is initially optically thick.The radiation in the X-ray band suffers bound-free absorption, and the optical rays suffer bound-bound absorption.The heating efficiency of the ejecta and its resulting luminosity sensitively depend on the photoionization process, which we will discuss in detail below.Similar to previous sections, we summarize the frequently used symbols in Table 3.

Ionization
The X-rays radiated into the ejecta lead to the photoionization of the elements.As mentioned above, we approximate the ejecta by matter composed of iron-56 which is initially neutral.The ionization is assumed to be balanced throughout the whole time and evolves as a quasi-static process.The ionization balance in a photon bath is where   is the number density of the   ℎ excited state of iron,   is the radiative intensity of the X-rays,    is the photoionization cross section,   rec is the recombination rate (Woods et al. 1981) which depends on ejecta temperature  ej , and  e is the number density of free electrons.In our definition,  starts from 1 to 27 for iron, where 1 correspond to neutral iron and 27 corresponds to fully ionized iron.The approximations of cross section for different ions are picked from Verner & Yakovlev (1995) and Verner et al. (1996), which are analytical interpolations of atom data.In principle, The above approximations are effective from the ionization threshold energy and up to 0.5MeV, above which relativistic corrections need to be considered.Here we simply extrapolate them to 1MeV, assuming the relativistic effects don't strongly impact our overall results.The relation between cross section and photon energy plays a crucial role in calculating energy absorption and X-ray radiation, as it directly determines the optical depth of X-rays.Roughly speaking, the cross section peaks at the threshold energy and decreases approximately following a power-law relation.The threshold energy tends to be higher for high excited states compared to low excited states, while the cross section at the threshold energy tends to be lower at high excited states.This is due to the difficulty of ionizing inner shell electrons that have greater binding energies than outer shell electrons.
Because the ejecta is initially neutral, all free electrons come from ionized atoms, so we have The X-ray intensity can be estimated by the luminosity of the PWN at the interface We normalize the ion number density to   =   / Fe where  Fe =  ej /(56 p ) is the iron atom number density.The degree of ionization can be expressed by the fraction of free electrons to total electrons, i.e.,  e =  e /(26 Fe ), or 1 −  e which is a measure of the optical depth of the ejecta.Now   can be solved by the following equations Given the solutions, we can calculate the bound-free opacity In addition to bound-free absorption, the hard X-rays can also be absorbed by down scattering, the corresponding effective opacity can be estimated by where the scattering opacity is Note that X-ray photons can also be scattered by bound electrons.
The total X-ray opacity of the ejecta is then In MP14, the propagation of photoionization is approximated by an ionization front at which optical depth is equal to 1. X-rays are not allowed to escape from the ejecta until the ionization front reaches the ejecta surface.The ionization fronts are associated with the ion species which dominate the photoionization.In contrast to this approach, we consider a different way of modeling the ionization.Using the opacity estimated above, we are able to track the optical depth at each frequency.As an X-ray photon traverses through the ejecta, it can either pass through, be absorbed through ionization, or reflected back due to scattering.The fate of this X-ray photon depends on the optical depth of absorption and scattering at its frequency.Thus, we do not employ an assumption of ionization front, but rather a frequency dependent transmission rate.Our approach does not deviate significantly from the "front" approach, but it results in a smoother X-ray light curve during the breakout time.In addition, since the transmission rate is frequency-dependent rather than associated with specific ion species, we are able to model the evolution of the X-ray spectrum.The details of this calculation will be shown in §2.4.2.
Although our model does not require an ionization front, we can still define an effective penetration depth using a similar approach to MP14.This depth can be considered as an indicator of the transparency of X-ray at a specific frequency.When this depth reaches the surface it provides an estimate for the time of X-ray breakout.
Following MP14, we approximate the depth by assuming the effective optical depth is equal to 1.The effective optical depth  eff, is the absorption optical depth corrected by the path length factor due to scattering The absorption and scattering optical depth can be calculated by where Δ  is the frequency dependent penetration depth.Equating eq.21 to unity, this value can be analytically solved It's maximum value is limited to the ejecta thickness Δ sh .

Scattering
Because of the scattering effect, X-ray photons not only have the possibility of being absorbed or passing through the ejecta but also a chance of being reflected back to the PWN.As mentioned above, the overall effect can be described by frequency-dependent rates of reflection ( ref  ), absorption ( abs  ), and transmission ( tra  ), which satisfy the normalization  ref  +  abs  +  tra  = 1.The values of these rates directly depends on the frequency dependent absorption optical depth ( abs ) and scattering optical depth ( sca ).However, their relation is very complicated, as the joint process of absorption and scattering is highly non-linear.The only way to determine the relation is through a Monte Carlo simulation.Unlike the approach taken by MP14, who only simulated cases where the ejecta is optically thick due to both scattering and absorption (resulting in a dependence solely on  abs / sca ), we aim to cover all possible combinations of  abs and  sca .The simulation is described as follows.
We consider a slab with a width of unity and an infinite area.The normal direction of the slab is represented by the z-axis.The slab has an optical depth of  abs due to absorption and  sca due to scattering.Before injecting a photon, we generate a random variable,   , which represents the maximum path length the photon can travel before being absorbed.The probability of a photon being absorbed after traveling an accumulated path length  follows a distribution depending on the mean free path l: () = exp(−/ l)/ l.In our setup, the mean free path is 1/ abs .So   follows the probability distribution () =  abs exp(− abs ).After generating this variable, we inject photons from one side with random directions, but with a positive z-component of velocity.These photons then start a 3D random walk within the slab due to scattering.Since scattering can be regarded as an absorption and emission process, the length of each step follows the same probability distribution with mean free path 1/ sca : () =  sca exp(− sca ).We terminate the photon's walk if (i) the photon is absorbed, i.e., the cumulative path length exceeds   , (ii) the photon is reflected, i.e.,  < 0, or (iii) the photon passes through the slab, i.e.,  > 1.For every pair of ( abs ,  sca ), we inject one million photons and calculate the reflection, absorption, and transmission rates.These rates, as functions of  abs and  sca , are presented in a contour plot shown in Fig. 2. The results are then interpolated to obtain a smooth function.

Thermal Radiation
Similar to §2.3, the thermal luminosity of ejecta can be estimated by where  ej =  ej  ej Δ sh is the bound-bound optical depth.The opacity  ej is hard to estimate due to the poor knowledge of the line forest   (Metzger 2019).As previously mentioned, the ejecta is likely composed of iron-like elements with an electron fraction  e ≳ 0.3 (Lippuner et al. 2017).Such a composition likely results in a relatively low opacity (Tanaka et al. 2020).Here we follow the common simplification in literature and set it to a constant value  ej = 1cm −2 g −1 (though we consider a broader parameter space in §4.2).
We assume that all photons absorbed are turned to thermal energy, and assume that the ejecta is in thermal equilibrium.The radiation then follows the black body formula.The black body temperature of the ejecta is the effective temperature at the surface The flux of the KN is where   is Planck's formula and  L is the luminosity distance.We have considered cosmological effects here, since later we will show that the most optimal magnetarboosted KNe can be observed up to a few Gpc.In this study we assume the following cosmology parameters:  0 = 68 km s −1 Mpc −1 , Ω  = 0.286 and Ω Λ = 0.714.

Evolution Equations
Now that we have all the necessary ingredients, we are ready to derive the evolution equations.The system losses its internal energy due to pdV work and radiation.Applying Eq. 1 and Eq. 2, the evolution equations can be summarized as follows (32) The system can be solved given the parameter set ,  ini , ,  ej .In our calculations, the ejecta is initially taken to be a sufficiently small spheroid with an initial velocity of 0.1c.The results are independent of the initial ejecta radius.
In this study we omit the r-process heating of the ejecta.This is because the radioactive power (e.g., Korobkin et al. 2012) is many orders of magnitudes smaller than the dipole power of the magnetar, and has no practical effect on our results.

Rayleigh-Taylor Instability
While solving the evolution equations, we also test the Rayleigh-Taylor (RT) instability of the system.This is very likely to occur in this scenario, because the heavy ejecta is accelerated by the light PWN matter at early times.If the ejecta breaks apart before the light curve peaks, the PWN matter will leak away forming an ultra-relativistic blastwave, and the rest of the energy will be insufficient to boost the KN.Fully capturing the dynamics of the RT instability would require hydrodynamical simulations, which is beyond the scope of our work.
Here we simply consider the linear growth rates and provide a rough test.
The growth timescale of RT instability is roughly estimated by where A is the Atwood number where  is the acceleration of the ejecta and  is the wave number of the instability.Here we approximate it by  ≈ 2/Δ sh .
We calculate this value throughout the evolution of the system.The system is considered unstable if the growth timescale is shorter than the dynamic timescale, i.e.,  RT <  dyn .The dynamic timescale is defines to be  dyn =  ej / ej , where  ej  =  ej / is the velocity of the ejecta.The evolution of the factor  RT / dyn will be shown together with other variables in §3, where we will see that the system is generally unstable to RT instability.

LIGHT CURVE SIGNATURES
In this section, we investigate the behavior of the PWN-ejecta system over the parameter space and the predicted optical and X-ray light curves.1 f e L abs / L sd sh / R ej t RT / t dyn / sh @ 0.1 keV / sh @ 1 keV / sh @ 10 keV / sh @ 100 keV In the upper panel we show the bolometric luminosity of the optical radiation by the black solid line, and the X-ray luminosity in given bands by colored dashed lines.Together with the luminosity, we also show the spin down power in a green solid line for reference.In the lower panel, the solid colored lines show the evolution of some typical parameters, including the reverse of ionization degree 1 −  e (an indicator of X-ray opacity), heating efficiency  abs / sd , the growth rate of RT instability  RT / dyn , and the ejecta shell thickness Δ  / ej as indicated in the figure.The dashed and dotted colored lines show the evolution of the approximate ionization front of X-rays at typical frequencies.

Typical temporal behavior of a indefinitely stable kilonova
First, in Figure 3, we show the temporal evolution of some critical parameters in an indefinitely stable magnetar boosted KN.The parameters are described in the figure caption.The upper panel shows the luminosity of thermal (optical) and X-ray radiation.The KN (optical) peaks after a few days with a luminosity ∼ 10 45 erg/s.The X-rays show a (frequency-dependent) sharp breakout as expected.
The peak luminosity of the X-rays is rather complicated, since it sensitively depends on the photoionization process.Our results are in general agreement with MP14.In MP14, the X-ray trend after the peak follows the same power index as the magnetic dipole formula.This is because, at that stage, the ejecta is fully ionized, allowing X-rays from the PWN to freely escape.However, in our model, with the exception of the relatively low mass of the ejecta, we find that the X-ray power index is initially slightly steeper than the dipole radiation.This is caused by the recombination of ions that increases the optical depth of ejecta, as demonstrated in the lower panel of the figure.However, at an even later stage, when 1−  e becomes constant, we still anticipate the same asymptotic trend following the magnetar spindown power.
In the lower panel of the figure, we present the evolution of several key variables of the system.The 1 −  e term represents the degree of ionization, which serves as an indicator of the X-ray opacity.It can be seen that the ejecta is highly ionized around the peak of the light curve, indicating that the radiation of the PWN is sufficient to fully ionize the ejecta, allowing X-rays to pass through.This result is compatible with previous studies using photoionization code CLOUDY (e.g., Margalit et al. 2018).As the X-ray intensity decreases after the peak, the ionization degree decreases and the opacity increases due to recombination.The heating efficiency, characterized by  abs / sd , varies significantly throughout the evolution, ranging from 0.01 to 0.5.Unlike some studies that assume a constant value, we find that the heating efficiency is dynamic in our model.The evolution of the variable Δ sh / ej confirms that the ejecta is compressed into a thin shell due to the high pressure of the PWN.The test of the Rayleigh-Taylor (RT) instability is shown by the ratio  RT / dyn .This ratio is generally less than 1 in the early stages, indicating that the system is prone to RT instability.It is important to note that our modeling only accounts for the linear stages of the RT instability, and may not fully capture the entire process.Numerical hydrodynamics is needed for a more accurate study, but is out of the scope of our work.Furthermore, even if the KN is disrupted by the instability, the existence of a SMNS can be revealed by the non-thermal signatures of the resulting blastwave.We will explain this case in §3.4.

Typical features of optical radiation
The most important features of an optical light curve are its peak luminosity (or flux) and peak time.The peak luminosity is difficult to analytically estimate in our model because it depends on the details of ionization, but the peak time is relatively easier.Before showing our numerical result, we first provide an analytical estimation of peak time  peak , which can be useful for direct comparisons with observations.
The peak time is roughly the time when the diffusion time scale of the ejecta  d ej reduces to the dynamical time scale  dyn .In the optically thick case, the diffusion time scale is We can roughly estimate the radius by  ej ∼  ej .Matching the diffusion time scale and the dynamical time scale, we have The value of  ej depends on whether the ejecta's kinetic energy around the peak time is dominated by the initial kinetic energy or the injected energy.Based on the two scenarios, we have the following derivations: i) Initial energy dominated.In this case the velocity maintains its initial value  ej =  ej,0 , so we can simply use eq.37 with the proper scaling  peak = 4.87 ii) Injected energy dominated.In this case most of the injected energy (i.e.,  ini ) is transformed to the kinetic energy of the ejecta, because the system remains optically thick before the peak time.We can calculate the velocity by  ej = √︃ 2 ini / ej  2 .Insert it into eq.37, we have To determine which case is relevant, we can first try case 1, get  1 peak and perform a consistency check, i.e., we can calculate the total injected energy up to this time and then calculate the corresponding velocity.If it's smaller than   ,0 , the result is valid, otherwise we can move to the second case.
From the above estimation we find the peak time is mostly dominated by ejecta mass and opacity, while other parameters have moderate impacts.This is in agreement with our anticipation, since these are the dominating parameters that determine when the ejecta becomes transparent.
The numerical optical light curves for different combinations of parameters are shown in Fig. 4, where we have compared the impact of different parameters by varying them against the previous "typical" case.We find the above peak time estimation agrees well with our numerical result within a deviation smaller than a factor of 2. Additionally, the fact that the above two cases provides similar estimation indicates that the peak time of the magnetar-boosted KN is relatively universal, typically occurring within a few days after the merger.
Besides the peak time, we can also roughly examine the impact of different parameters on the peak luminosity, as shown in fig. 4. We first examine the light curves produced by a indefinitely stable magnetar, where only the parameters ,  ini , and  ej vary.We find that a large magnetic field can significantly suppress the luminosity due to increased pair production in the PWN, leading to higher opacity.The ejecta mass, on the other hand, doesn't strongly affect the results.This finding is mostly consistent with MP14.
Surprisingly, the peak luminosity appears to be insensitive to the initial rotational energy of the magnetar, which is essentially the available energy budget of the system.We have explored various parameter ranges and consistently found this result in our model.One explanation is that a larger rotational energy also corresponds to a higher spindown power, which leads to increased pair production similar to the case with a large magnetic field.Consequently, the increased energy budget is counterbalanced by the corresponding pair suppression.This can be understood by revisiting equation 5 and equation 10, where the  ini boosts the spindown power which increases the optical depth due to pair suppression.Another possible reason is that the spindown power is no longer dependent on  ini when  ≫  sd .This can be seen by taking this limit for eq. 3, 5, and 6.This finding suggests that for indefinitely stable magnetars as merger remnants, the luminosity of the boosted KN is primarily determined by the magnetic field, regardless of the initial rotation of the magnetar.
In contrast to other parameters,  has a much bigger impact on the light curve.Except for the early small bump caused by the collapse of magnetar, the rest of the light curve evolves like that of an adiabaticaly expanding shell, and fades away when it becomes optically thin.The peak luminosity drops by 3 orders of magnitude even if  drops by just 1 order of magnitude, corresponding to  ini amount of total energy input.The effect of  is not same as trivially reducing  ini accordingly, because the termination of central engine completely changes the acceleration process of the ejecta.

X-ray Spectrum evolution
Another interesting behavior of this typical case is the frequency dependent evolution of the X-ray opacity, as presented in the lower panel of Figure 3.It is important to note that the breakout time of X-rays is not a monotonically increasing function with photon energy up to 1 MeV.Instead, we find that photons with energies around 10 keV are the last to pass through the ejecta (see dotted pink line in the lower panel of Fig. 3).This behavior is attributed to the nature of the ionization cross section as a function of photon energy.As we have explained in §2.4.1, the cross section extends from the ionization threshold energy to 1 MeV, which approximately follows a decaying power-law formula.Ions at low excited states generally have lower threshold energies, so the low energy photons can only ionize ions of low excited states.After the ejecta is highly ionized, only rare ions remain in low excited states.Consequently, the ejecta becomes nearly transparent to low-energy photons since they are unable to reach the threshold energy of the existing ions.On the other hand, high-energy X-rays can also easily pass through the ejecta due to their small cross section.It is the photons in the intermediate energy that experience the most absorption as they are capable of reaching the threshold energy while still having a relatively larger cross section.Therefore, our findings indicate a potential evolution of the X-ray spectrum.This behavior is not specific to iron-like elements but rather a common characteristic of heavy elements, as they generally follow a similar rule for cross sections.While our results may offer interesting observational predictions in the x-ray band, this is not the primary focus of our work, and we leave it for a future study.

The limit of forming an ultra-relativistic blast wave
As mentioned above, we find the system is generally unstable to the Rayleigh-Taylor instability.Though our estimation based on the linear growth rate may not be precise enough, it's still worthwhile to study the extreme limit in which the PWN matter completely leaks out from the ejecta.
In this case, we assume the leaked energy forms an ultra-relativistic blastwave propagating into the surrounding environment.The blastwave accelerates electrons which produce synchrotron radiation, just like the case of a GRB.However, unlike regular GRBs, once the PWN leaks from the ejecta shell, there is no mechanism to confine the material into a narrow cone.As a result, the dynamics of the blastwave will be quasi-isotropic rather than jet-like.The light curve should be very similar to a GRB afterglow, except that it can be observed from all directions and exhibits no jet-break.We consider a simple analytic model of the blastwave (see, e.g.Kumar & Zhang 2015).To simplify the calculation, we only consider the slow cooling case where the synchrotron injection frequency,  m , is less than than the cooling frequency  c and where the observed frequency falls within the range between the self-absorption frequency  a and the cooling frequency  c .This frequency range is generally enough to encompass the optical waveband.The observed flux is where  is the electron energy distribution power-law index.The synchrotron peak frequency  m is calculated by where  e and  B are the fractions of internal energy of the blastwave converted to non-thermal electrons and magnetic fields, respectively.The peak flux is iso 10 52 erg where  0 is the ambient number density.Typically, there should be hydrodynamic breaks other than the spectral break, such as transition from coasting to deceleration, the lateral expansion (if the blast-wave is sufficiently anisotropic), and the transition from relativistic motion to Newtonian velocities.However, the first one is much earlier than the time of interest (tested for initial Lorentz factor of the blastwave Γ > 100), while the later two happen when the flux has significantly dropped off, thus not impacting our result.For simplicity, we don't include them here.
The microphysical parameters of this afterglow-like light curve should be similar to those of short GRBs.The ambient density should be relatively low with  0 = 10 −3 -10 −2 cm −3 .We consider a fixed  e = 0.1 since it is observationally constrained to be narrowly distributed between different GRBs (Beniamini & van der Horst 2017).As for  B , which is observationally less well-determined and may vary more from burst to burst, we consider a range between 10 −3 to 10 −2 .Moreover, to make a fair comparison, we assume the energy of the blast wave matches with our "fiducial" case, i.e., we fix the isotropic energy to 3 × 10 52 erg.The waveband is the same as for the KN.
To compare the light curve in this scenario with the boosted-KNe model, we present them together in fig 5.The flux of the KN is calculated using eq.27.The unit of flux is converted to AB magnitude for convenience.The KN parameters remain the same as those in Figure 4, and the observable features of the flux match with the luminosity, so we don't repeat their description here.As we can see, the afterglow is in general comparable to or even brighter than the KN, though it peaks at a much earlier time.Moreover, since the blastwave is quasi-isotropic, the prospect of detection is not limited by the jet opening angle, thus significantly increasing the potential number of observable sources.If this scenario were correct, the detection rate of orphan afterglow would be expected to be very high, which is not evident in current observations (e.g.Ho et al. 2022).Therefore, the disruption of the KN ejecta because of the Rayleigh-Taylor instability is not a plausible explanation for the absence of detections of the transients following neutron star mergers.A quantitative simulation to explore the observed rate of the blast wave signatures is not pursued in this work.

CONSTRAINT BY ABSENCE OF OBSERVATION
In this section, we aim to compare our model with the current observations.As described in the previous section, the bright nature of magnetar-boosted KNe implies that they should be observable over a large cosmic volume.This suggests that the detection rate of these events could be comparable to, or even higher than, that of regular KNe, despite having a lower intrinsic event rate.The lack of confirmed detections puts strong constraints on either the parameters of the model or the characteristics of the remnants from neutron star mergers.
In this work, we perform Monte Carlo simulations to estimate the observational rates.However, certain model parameters, such as the initial rotational energy  ini and energy extraction efficiency , are not free to setup but depend on factors like the progenitor mass of the binary system and EoS.Furthermore, there is still a significant uncertainty in the EoS, which would introduce model dependency to the EoS based simulation.To address these considerations, we employ two approaches: a model-independent study and a model-dependent study.In the model-independent study, we generically evaluate the maximum detectable distance based on various parameter combinations, which serves as an indicator of the event rate.On the other hand, the model-dependent study begins with a population of binary systems, evolves them to remnants, and generates model parameters based on an assumed EoS.We then simulate observations using a sky survey strategy and collect the observed events to determine the detection rate.
The model-independent study provides an overview of our model's predictions regardless of EoS, while the model-dependent study offers a quantitative result.The details of the two approaches are described below.

Model independent study
In this model-independent study we examine the observation potential of a parameter set indicated by its maximum detection distance.The maximum detection distance is roughly set by requiring the peak flux reaches a telescope's detection threshold.We assume the threshold is  AB =20.5, matching with the performance of ZTF (Smith et al. 2014;Dekany et al. 2020).Similar to the previous section, the filter is set to r-band.We notice that a realistic confirmation of a detection requires more than one convincing data point, thus the peak flux should be slightly higher than detection threshold.We leave this effect to the model dependent study.
The optical peak flux of the magnetar-boosted KN is primarily dependent on the values of magnetic field  and energy extraction efficiency , as explained in the previous section.In order to simplify our study, we fix the value of  ej to 0.1 ⊙ motivated by numerical simulations (Margalit & Metzger 2019).To illustrate the impact of both  and , we consider four cases: B=10 16 G, B=10 15 G, B=10 14 G, and B=3 × 10 12 G.The first three cases account for the possible varying domain of magnetic field, while the last one serves as a selfconsistency check, which will be described below.For each case, the maximum detection distance should be a function of , which is solved by requiring the peak flux is equal to the detection threshold, i.e.  AB = 20.5.To account for the minor influence of different  ini , each case is represented by a shaded region, encompassing  ini values ranging from 10 52 erg to 10 53 erg.
Our results are shown in Figure 6.As we can see, the peak luminosity of the KN increases with decreasing magnetic field strength, leading to the increase of maximum detectable distance.This is because reducing the magnetic field leads to a smaller spin down power, which produces less pairs in the PWN that suppress the energy injection.On the other hand, it is expected that the luminosity will eventually decrease as we keep decreasing the magnetic field, since the spin-down timescale could be significantly longer than the peak time, resulting in less energy injection while the ejecta is optically thick.To provide a self-consistency check, we verify this phenomenon by the B=3 × 10 12 G case.As shown in the blue shadow region, the maximum detectable distance indeed decreases as compared with the B=10 14 G case.
The implications of this result will be discussed in §4.3.
10 2 10 1 10 0 10 1 10 2 10 3 10 4 D l, max [Mpc] B=3 × 10 12 G B=10 14 G B=10 15 G B=10 16 G Figure 6.The maximum detectable distances of magnetar-boosted KN as a function of  and with magnetar field strengths of  = 3 × 10 12 G,  = 10 14 G,  = 10 15 G, and  = 10 16 G.The shaded region is bracketed by the limits: 10 52 erg <  ini < 10 53 erg.The lowest  case serves as a self-consistency check, where the peak luminosity should eventually reduce for such low magnetic fields, because the spin down timescale is much longer than the diffusion timescale.In this case the peak luminosity is also independent of  since the magnetar never collapses in the time of interest.

Model dependent study
In this model dependent study, we simulate the whole process from a given population of BNS systems to the detection rate of a boosted-KNe.Our method is explained in the following steps.
In the first step, we consider a population of BNS systems and let them merge and generate remnants.We assume that the population follows the Galactic neutron star mass distribution, which is approximated by a normal distribution with a mean value of 1.33 ⊙ and a standard deviation of 0.09 ⊙ (Antoniadis et al. 2016;Özel & Freire 2016).The properties of the remnants (which will be determined in the next step) depend on two parameters that are determined during this step: the baryonic mass and the initial rotational energy, i.e.,  ini .The baryonic mass is obtained by summing the masses of the progenitor system and subtracting the ejecta mass.As previously mentioned, we anticipate the ejecta mass to be around 0.1 ⊙ if the merger remnant is a SMNS (Margalit & Metzger 2019).To account for potential uncertainties, we consider a uniform distribution ranging from 0.01  ⊙ to 0.1  ⊙ , with an emphasis on the higher end of the range.The initial rotational energy is influenced by energy losses of the system during the differential rotation phase, which is not yet fully understood.We treat it as a free parameter scaled by the maximum allowed rotation, i.e., the Kepler rotation.Specifically, this free parameter is denoted as  ini / kep .
In the second step, we calculate the fate of the merger remnant.A remnant can either collapse under its own gravity or survive for a period of time, depending on whether it reaches the threshold mass of uniform rotation.The survival of a remnant is highly sensitive to the equation of state (EoS), which can be roughly characterized by the Tolman-Oppenheimer-Volkoff (TOV) mass, denoted as  TOV .In our study, we consider two EoS: UU (Wiringa et al. 1988) and SLY (Douchin & Haensel 2001), which correspond to  TOV of 2.2 ⊙ and 2.05 ⊙ , respectively.These values mostly cover the lower limit set by the most massive pulsars (Antoniadis et al. 2013;Fonseca et al. 2021) and the upper limit constrained by the GW170817 (e.g., Mar-galit &Metzger 2017 andMa et al. 2018).To determine the evolution and status of a remnant based on its initial conditions, we employ the RNS code (Stergioulas & Friedman 1995).Specifically, we calculate the threshold mass required for the formation of a SMNS, which depends on  ini / kep , and compare it with the remnant mass calculated in the previous step to determine whether they can survive.We also calculate the critical rotational energy  crit just before the collapse of the SMNS.This critical energy represents the end state of the SMNS.By comparing the initial energy  ini with  crit , we can derive the energy extraction efficiency in our KN model  = ( ini −  crit )/ ini .If the magnetar is indefinitely stable, we set this parameter to 1.Note that due to truncation errors in the RNS code, the obtained  TOV values slightly deviate from the theoretical values.To avoid inaccuracies that could be caused by artificially scaling the results, we adopt the values provided by the code, which are approximately 2.17 ⊙ for SLy and 2.07 ⊙ for UU.This implementation doesn't change our conclusion.The precise value of  TOV for UU implies a higher rate of SMNS formation than what we've considered, imposing a more stringent constraint on this EoS.On the other hand, the exact value of  TOV for SLy suggests a shorter survival timescale for SMNS, which aligns with our conclusion.
In the third step, we generate a large set of event parameters.The ejecta mass  ej , initial rotational energy  ini , and energy extraction efficiency  have been determined in the previous steps.To account for a range of magnetic fields, in this work, we assume a lower limit of 10 14 G and upper limit of 10 16 G with uniform distribution in logarithmic space.This assumption is based on the consideration of magnetic field amplification during the previous differential rotation phase of the magnetar, and the maximum magnetic field allowed for a stable magnetar configuration.The distribution of distances at which these events occur depends on the evolution of the BNS merger rate () with redshift.For simplicity, we assume that () is proportional to the sGRB rate and adopt the analytical approximation derived by Wanderman & Piran (2015) (see equation 9 therein).The rate is then scaled to match the local BNS merger rate.Given other uncertainties in our model, we believe this assumption is sufficient.Using these parameter distributions, we generate a set of millions of parameter combinations.The total number of events is denoted by  tot .
Finally, we calculate the light curves based on our model.To address the uncertainty in ejecta opacity, we assume a uniform distribution of  in logarithmic space, ranging from 1 cm −2 g −1 to 10 cm −2 g −1 , with an emphasis on the lower values.We then proceed to calculate the detection rate by selecting the light curves that can be detected.In order to compare with the ZTF observation, we designed a similar sky survey strategy.For each generated light curve, we select a series of time points in the r-band.The time interval between neighboring points is fixed at 3 days, in order to mimic the approximate cadence of the ZTF survey strategy (e.g., Andreoni et al. 2020).To introduce some variability, the time series is randomly shifted.A data point is considered "observed" if its flux exceeds the threshold value.Similar to the model-independent study, we assume the threshold flux to be  AB = 20.5.We confirm the detection of an event if at least three data points are observed.We count all detected events and get the total number  det .The expected yearly detection rate is then calculated by where  is comoving distance corresponding to redshift  and Ω fov = 47 • is the field of view of ZTF (Andreoni et al. 2020).In this simulation, the maximum distance  max is set to be sufficiently large SLy.The black dashed line is the approximate observational upper limit set by 1/ ZTF .In the lower panel we show the fraction of BNS mergers leading to long-lived magnetar formation as a reference.To clarify the population of these events, we separate them to two classes.In both panels, the dashes lines are the events produced by indefinitely stable magnetars, and the the dotted lines are the temporarily stable magnetars.
(corresponding to z∼1.5) in order to cover the most optimal case (see fig. 6).The BNS merger rate is scaled to match the local rate (0) = 300 Gpc −3  −1 , motivated by the recent constraint (Mandel & Broekgaarden 2022).Different local BNS rates can easily be accommodated for by scaling the results given here.
The results are shown in fig. 7. To demonstrate the population of detected events, the total detection rates (thick solid lines) are divided into two populations: events produced by indefinitely stable magnetars (thin dashed lines) and events produced by temporarily stable magnetars (thin dotted lines).The rates for each EoS are plotted in the upper panel 1 .In the lower panel, we show the formation rate of the two populations as a reference.Furthermore, considering the lack of confirmed detections, we can derive an approximate upper limit of 1/ ZTF yr −1 , where  ZTF represents the effective operational time of the ZTF so far.Based on the reported effective operational time by the ZTF team in 2020 (approximately 2 years), we adopt the value of 4 years for this study.It is important to note that this upper limit is only an approximate estimation, and we will provide a more strict discussion below.

Constraints on SMNS and the neutron star EoS
Our results place a strict constraint on the fate of SMNS.
In Figure 6, our results reveal that if the magnetar is indefinitely stable, the optimal detectable distance can reach several to 10 Gpc.Even with strong magnetic field suppression (B=10 16 G), the maximum distance can still reach approximately 1 Gpc.Considering a neutron star merger rate of approximately 300 Gpc −3 yr −1 , this implies that more than one thousand merger events can enter the detectable volume since the start of the optical survey.The absence of detections therefore suggests that the fraction of BNS mergers leading to SMNS formation must be smaller than 10 −3 , which is in contrast to the estimation based on the current constraints of EoS assuming  = 1 (i.e., assuming Keplerian rotation at the birth time of SMNS).
A plausible explanation is that either most of merger remnants collapse during the differential rotation phase or, at the onset of uniform rotation, their energy is very close to the critical energy required to support them against collapse.The latter scenario corresponds to cases with  ≪ 1 as shown in Figure 6.
This result can be more qualitatively seen in the model dependent study.In fig 7 .we can see that if the SMNS starts from Kepler rotation (i.e.,  ini / kep = 1), both UU and SLy predicts detection number above the observational limit.For the EoS of UU, the expected detection rate is above the limit even with  ini / kep = 10 −2 .Therefore, our result tend to prefer SLy over UU.In other words, the  TOV should be close to the lower limit of the current observational constraints.However, even with SLy, the expected detection rate is able to match with observation only when  ini / kep ≲ 1.Note if  ini / kep is too small, it will be impractical for the formation of SMNS, since their survival timescale would be negligible.If these events are the majority of the population, the events will be dominated by indefinitely stable neutron stars, instead of temporarily stable SMNS.
To provide a rigorous statistical analysis, we can perform a hypothesis test by assuming that the number of detections follows a Poisson distribution.The mean value for each  ini / kep case is the expected detection number over a period of  ZTF = 4 years.The probability of observing 0 detection in a Poisson distribution is simply given by  − , where  represents the mean value.We can then convert this probability into equivalent  levels in a normal distribution, which serves as an indicator of the "rejection level".The result is shown in fig.8.
In this result, we find that the scenario where the BNS merger remnants following the UU EoS forms a long-lasting SMNS is rejected at a significance level of > 3 (see the blue line in fig.8) It is important to note that UU EoS itself is not ruled out, but the SMNS with high initial rotational energy is ruled out.
On the other hand, for the SLy EoS, the long-lasting SMNS scenario is rejected at a significance level of ≳ 1 if the SMNS is initially rotating at the Keplerian speed.From these results, we can conclude that regardless of the EoS, when the merger remnant transitions into the uniformly rotating phase, it is unlikely to have a Keplerian rotation speed.This suggests that there must be some energy and angular momentum loss during the differential rotation phase.This result is consistent with some recent studies and numerical simulations, where it is argued that the a merger remnant could collapse into a black hole before it becomes a SMNS, even if its mass allows to stabilize itself at Keplerian speed.This could happen if the remnant losses significant fraction of its angular momentum before it enters such phase (Beniamini & Lu 2021).Alternatively, if the remnant rearranges its angular velocity profile during the differential rotation, such that the core slows down faster and initiates the collapse before reaching a uniform rotating configuration (Margalit et al. 2022).
Our results are based on the assumption that the masses of merging neutron star binaries are similar to those observed in our Galaxy.Our constraint on EoS and initial rotation speed of SMNS can be relieved if the BNS mass distribution in the Universe is heavier than the Galactic distribution.However, this assumption leads to even more strict condition for SMNS formation.Therefore, the conclusion that SMNS are rare and short-lived objects is unchanged.

Implications on potential boosted kilonova candidates
Besides the non-detection of boosted KN from sky surveys, there are some candidates (regular KNe) found in sGRB afterglow, such as GRB 130603B (Berger et al. 2013;Tanvir et al. 2013), GRB 060614 (Yang et al. 2015;Jin et al. 2015), GRB 050709 (Jin et al. 2016) and GRB 080503 (Perley et al. 2009).There is no clear evidence suggesting that the optical excess of these events require an additional energy source in the form of long-lived magnetars.In some of the events (e.g., GRB 130603B and AT2017gfo) there are strong limita- tions on the presence of long-lived magnetars, since their associated kilonovae show no signs of boosting.
To better demonstrate how the candidate events are in tension with the boosted KN model, we present a contour plot in Figure 9, showing the peak luminosity of boosted KNe as predicted by our model.This plot shows the variation of the peak luminosity with respect to the two dominant parameters,  and .Additionally, we also calculate the corresponding survival timescale  c and present it on the same figure.We find that a long-lived magnetar that survives longer than the spindown timescale ( ≳ 0.5) and lasts for ≳ 1000 seconds generally boosts the KN to a luminosity  peak > 10 43 erg/s.Such a luminosity exceeds any confirmed or candidate KN known to date.In other words, to explain the lower luminosity of these events compared to our model predictions, one must assume a short survival timescale of the central magnetars, or rule out the magnetar explanation.
For instance, consider the only confirmed KN, AT2017gfo with a bolometric luminosity of  bol ∼ 10 42 erg/s.If we assume a magnetar origin, in Figure 9 we can find that even in the most extreme cases ( ∼ 10 16 G), it still requires  < 0.5, corresponding to a magnetar surviving less than the spindown timescale.In addition, the survival timescale in the allowed parameter region is less than an hour.This result further indicates that late time X-ray activity of GW170817 is unlikely to originate from a magnetar remnant.A similar conclusion can also be drawn regarding the other aforementioned candidate events in GRB afterglows.Furthermore, the blast wave kinetic energy in those events (e.g.,  k ≲ 10 50 erg for AT2017gfo, Balasubramanian et al. 2022 and  k ≲ 10 51 erg even for the beaming-corrected energy of the GRB afterglow blast wave in the same event, see Margutti & Chornock 2021 and references therein) is much lower than would be expected in case the nebula material had leaked out of the KN ejecta (10 52 erg -10 53 erg).Therefore, even if magnetars have ever existed in these events, they are unlikely to have been long-lived or responsible for late time activity.
The presence of a long-lived magnetar as the merger remnant of GW170817 is also disfavored by various other studies (e.g., Margalit & Metzger 2017;Granot et al. 2017;Murase et al. 2018).Recently, a luminous candidate was found in GRB 200522A (Fong et al. 2021), which may serve as a candidate of boosted KN (see however O'Connor et al. 2021 for a different interpretation).The luminosity (≳ 10 42 erg/s) of this event, although brighter than the KN of GW170817, is still much fainter than the typical luminosity of a stable magnetar predicted in our model.The small energy injection can also be explained if there are additional processes to slow down the magnetar except for dipole radiation.Possible candidates are gravitational waves and neutrino cooling.Previous works have argued that gravitational wave losses may dominate the spindown process if the ellipticity of the SMNS is sufficiently large (Fan et al. 2013a,b;Li et al. 2018;Ai et al. 2018).Since these energy dissipation processes do not inject additional energy into the ejecta, modeling these processes is similar to using a shorter remnant survival timescale, i.e., a smaller .This is because the small survival timescale compared with evolution timescale will hide the details of energy injection so that the result is insensitive to the spin down power index.However, in order to achieve the desired level ( ≪ 1) as constrained by observations, these energy extraction processes must contribute approximately 10 times more energy than the dipole power.It remains uncertain whether this can be accomplished through gravitational waves or neutrinos.In particular, if such amount of loss is due to gravitational wave radiation, it will lead to unstable magnetic field configuration in the magnetar.Even if it happens, such powerful energy lose also implies a short survival timescale of the SMNS.In fact, recent studies (e.g.Sarin et al. 2022) that incorporate gravitational wave emission have reached conclusions consistent with our work.

Implications on Orphan GRB afterglow
The lack of observational evidence for powerful KNe does not necessarily imply the absence of a powerful energy injection by the SMNS.It is plausible that the PWN energy is not transferred into the ejecta.Instead, the nebular-ejecta interface may turn out to be violently unstable to the RT instability which produces holes through the ejecta where PWN matter can escape.As we have mentioned before, the escaped energy will form a blastwave which produces powerful emission similar to the GRB afterglow, which implies an overwhelming number of orphan afterglow in contrast to observations.However, as our simulation indicates, there is still a chance that some merger remnants form indefinitely stable magnetars with low initial rotation energy.If prone to the RT instability, they may serve as a potential population of orphan afterglows.

Implications on progenitors of Fast Radio Bursts
Fast rotating magnetars as BNS merger remnants are also considered as possible sources of fast radio bursts (FRB).There are two important constraints on this scenario.One is that the dispersion measure (DM) from the source can't exceed the total DM in FRBs, which is typically a few hundred pc cm −3 .In our model, the PWN is rich of pairs and the ejecta is highly ionized, implicating a very large DM, which may not be compatible with observations.The other is the free-free absorption of radio waves.
We first calculate the DM following where  e is the electron number density along the path and  0 is the site where FRB is generated.Depending on models, the FRB is either produced near the magnetar or in the magnetar wind, so the limit of  0 is 0 <  0 <  n .However, in our model, the PWN is optically thick before the light curve peak.This means that any FRB produced near the magnetar can't escape from the PWN before the diffusion timescale (around peak time) due to Thompson scattering.Therefore, we assume  0 =  n , implying that FRBs are produced at the outer layer of PWN, and the source of DM is purely from the ejecta.In the ejecta, the number density of free electrons is  e = 26  e  ej /(56 p ) which are produced by the photoionization.Note that if an FRB is produced near the magnetar after the light curve peak, it should likely be able to escape.However, since the pair density calculated by eq. 9 is overestimated in this stage (though it doesn't affect the luminosity calculation), it can't be used to calculate the DM of the PWN.Our estimate of the DM (i.e., considering the ejecta only) at this stage should be regarded as a lower limit.
The free-free absorption optical depth of the radio waves places another constraint.For a similar reason as above, we also only calculate the optical depth in the ejecta One may refer to Table 3 for the definition of the symbols.Also note that  = 1 corresponds to neutral atoms in our definition, meaning that ions at  th ionization state have charges of  − 1.In this work we assume the Gaunt factor is ḡff = 1.The frequency is set to =1 GHz.The radio waves are unable to escape if the ejecta is optically thick to free-free absorption.
Our results are shown in fig.10.The parameters of the model are taken to be the same as in fig. 3. We can clearly see that even if we only consider the DM from ejecta, it reduces to values consistent with observations after at least 10 to 100 days.Moreover, the situation becomes even more constraining when considering the free-free absorption optical depth, since the ejecta gets transparent only after ∼ a year (note that at this time the DM from the ejecta is ≲ 4pc cm −3 and as such is no longer constraining).These results indicate that FRBs can't escape from the ejecta for at least a year after the merger, thus challenging the merger model.
We caution that if the ejecta is disrupted due to the RT instability, the above calculation of DM no longer holds, since the pair cascade won't be triggered.In this case, we should also see a non-thermal multi-waveband radiation from the blastwave, and FRB radiation (if produced) may also escape more easily.We also caution that these arguments do not apply to the possibility of an FRB forming prior to the merger of BNS, e.g., a FRB produced by the interaction between the magnetospheres of the two neutron stars (Sridhar et al. 2021;Most & Philippov 2022).

Implications on future Multi-messenger observation
Magnetar-boosted KNe also serve as electromagnetic counterparts of gravitational waves.Our result (i.e., fig 6) has shown that their bright nature allows for a maximum detection distance of a few Gpc in the most optimal cases.This distance is much larger than the horizon distance of BNS mergers for LIGO's GW detectors in the current and future planed observing run.Our study indicates that such events are very rare, so we generally do not expect a detection as a counterpart of the gravitational waves.That being said, if such events indeed happen, considering their immense brightness, they 3.Note the typical DM observed in FRBs, which may be attributed to its local environment, is less than a few hundred.Since the estimated DMs are much greater than observed values for FRBs, and the ejecta is optically thick to radio waves, our result indicates that a FRB is not likely to form in the first year after the merger.The grey shaded region indicates the time period ruled out for FRB production owing to the free-free optical depth.
are very likely to be detected.Such a detection would be extremely useful for constraining SMNS formation.The detection of a booted KN will place a very strict constraint on the magnetar model.The fate of merger remnants can be sensitively constrained by the peak luminosity, which mostly depends on the survival timescale of the SMNS.The type of central engine (i.e., a stable magnetar, a SMNS, or a black hole) can be constrained by the slope of the light curve, which asymptotically follows the power of energy injection.
As mentioned above, a long-lived SMNS is also a source of gravitational waves if it has some ellipticity.If its spin is fast and its survival timescale is long enough, the collected cycles of waves may be enough for detection before it collapses.The detection prospect sensitively depends on the ellipticity, which further depends on the type of instability during the rotation and the EoS (see Lasky 2015 for a review).Our result indicates that SMNS can't be long-lived, thus the gravitational waves from this stage are unlikely detectable.The detection of gravitational waves in this stage will imply a completely different interpretation of the lack of boosted-KN.

SUMMARY AND CONCLUSION
In this study we built a magnetar-boosted KN model with a detailed description of relevant physical processes such as the collapse of the magnetar central engine and find its effects on the peak brightness of the light curve.To compare with observations, we conduct both a model dependent and model independent studies to estimate the detection rate.We find that the only way to match the non-detection result is to require an energy injection that is significantly smaller than the rotational energy of a maximally rotating remnant, which means that the SMNS (if formed) is likely to be short-lived in the vast majority of the cases.
In principle, the system under consideration, a light PWN confined by the heavier KN ejecta is prone to the Rayleigh-Taylor instability.Indeed the linear growth rate of a mode with wavelength comparable to the ejecta thickness is typically faster than the expansion timescale of the system.This leads to the possibility that the PWN breaks out from the KN ejecta driving a relativistic blast wave into the ambient gas.We find the blast wave produces an afterglow-like emission, which is quasi-isotropic and of comparable or brighter than the boosted KN signal.Therefore, the uncommon detection of orphan afterglows implies that this is not a common occurrence.
Our result has several implications.The statistically short-lived SMNS either means the initial rotation of a newly born SMNS is much slower than the Kepler rotation, or they collapse into a black hole before most of energy are released.One possible reason is that during the differential rotation phase, the angular momentum transfer from the core to the outer shell is fast, such that the SMNS collapse into a black hole before it establishes uniform rotation.Other possibilities are additional energy extraction mechanisms, such as the gravitational wave radiation or neutrino cooling, or alternatively, a heavier neutron star mass distribution.These all lead to early collapse.Due to their short survival timescale, we do not expect to see the magnetar-boosted KNe in LIGO's upcoming observation run.

Figure 2 .
Figure 2. The Monte Carlo simulation result of absorption, reflection, and transmission rate of a slab caused by both absorption and scattering processes.The red lines show the contour plot of  tra as a function of optical depth of two processes, and the blue lines are the  ref . abs can be calculated by 1 −  tra −  ref .

Figure 3 .
Figure3.The solution of a typical parameter set:  = 10 15 G,  ini = 3 × 10 52 erg,  = 1,  ej = 0.1  ⊙ .In the upper panel we show the bolometric luminosity of the optical radiation by the black solid line, and the X-ray luminosity in given bands by colored dashed lines.Together with the luminosity, we also show the spin down power in a green solid line for reference.In the lower panel, the solid colored lines show the evolution of some typical parameters, including the reverse of ionization degree 1 −  e (an indicator of X-ray opacity), heating efficiency  abs / sd , the growth rate of RT instability  RT / dyn , and the ejecta shell thickness Δ  / ej as indicated in the figure.The dashed and dotted colored lines show the evolution of the approximate ionization front of X-rays at typical frequencies.

Figure 4 .
Figure 4.The comparison of optical light curves given different parameter combinations.The black line is the same as in previous figure, which serves as a "fiducial" case.In the colored lines we change model parameters individually as indicated in the label corner.

Figure 5 .
Figure 5.The flux of magnetar boosted KN.The parameters are the same as the previous figure.The redshift is set to 0.1, corresponding to a luminosity distance of ∼ 474 Mpc.The wave band is chosen to the  band (∼ 634 nm).In addition, we also plot the light curve of a blastwave in a grey shaded region assuming it is afterglow-like.The parameters of the shown blastwave are:  iso = 3 × 10 52 erg,  = 10 −3 − 10 −2 cm −3 ,  e = 0.1 and  B = 10 −3 − 10 −2 .The waveband is the same as for the KN.

Figure 7 .
Figure7.The simulated detection rate in the model dependent study.In the upper panel we show the detection rate (solid lines) as a function of the initial condition  ini / kep of SMNS assuming the equations of state of UU and SLy.The black dashed line is the approximate observational upper limit set by 1/ ZTF .In the lower panel we show the fraction of BNS mergers leading to long-lived magnetar formation as a reference.To clarify the population of these events, we separate them to two classes.In both panels, the dashes lines are the events produced by indefinitely stable magnetars, and the the dotted lines are the temporarily stable magnetars.

Figure 8 .
Figure8.The rejection level of the given the initial rotation  ini / kep assuming equation of state of UU and SLy.The rejection probability is calculated assuming the number of detections follows a Poisson distribution.The probability is then converted to corresponding significance levels.Assuming EoS of UU, the scenario of SMNS is ruled out at a 3- level.Assuming EoS of SLy, the initial rotation of SMNS at Kepler speed is ruled out at a 1- level.

EFigure 9 .
Figure9.The peak luminosity and survival time of magnetar boosted KN with respect to two dominating parameters: the magnetic field  and the energy injection efficiency .The solid lines represents the KN peak luminosity for two values of the initial energy:  ini = 10 52 erg and  ini = 10 53 erg.The dashed lines represents the survival time  c assuming  ini = 10 53 erg. c will be 10 times larger for  ini = 10 52 .The ejecta mass and opacity are taken to be the typical values for a boosted KN:  ej = 0.1  ⊙ and  = 1.

Figure 10 .
Figure10.The dispersion measure (shown in the solid black line) and freefree optical depth (calculated at 1 GHz, shown in the dashed black line) of the ejecta shell along the evolution.The red dotted line is the peak time of the light curve for reference.The model parameters are taken to be same as fig.3.Note the typical DM observed in FRBs, which may be attributed to its local environment, is less than a few hundred.Since the estimated DMs are much greater than observed values for FRBs, and the ejecta is optically thick to radio waves, our result indicates that a FRB is not likely to form in the first year after the merger.The grey shaded region indicates the time period ruled out for FRB production owing to the free-free optical depth.

Table 3 .
where the optical depth  n Radiation process variables.
caused by Thomson scattering can be estimated by