Recently, the LIGO (Aasi et al. 2015) and Virgo (Acernese et al. 2015) scientific collaborations (LVC) reported the detection of one of the most enigmatic gravitational wave (GW) mergers till date (Abbott et al. 2020 ). This event, named GW190814, has been associated with a compact object binary of mass-ratio, q = 0.112 +0.008 −0.009 , and primary and secondary masses m 1 = 23.2 +1.1 −1.0 M and m 2 = 2.59 +0.08 −0.09 , respectively. Since, an electromagnetic (EM) counterpart has not been found for this particular event and the tidal deformability has not been measurable from the GW signal, the secondary component might well be the lightest BH ever found. However, EM emissions are expected for only a fraction of neutron star binaries, and tidal deformabilities are known to be small for massive NSs, hence the secondary in this case cannot be ruled out as a neutron star. In the latter scenario, it would become the heaviest neutron star (NS) observed in a binary system, given its well-constrained mass. Either hypothesis deserves a deep study owing to its far-reaching implications on the formation channels of such objects and the nature of the densest form of matter in the universe.
Discoveries of massive pulsars in past decades have severely constrained the equation of state (EoS) of supranuclear matter inside their cores (Demorest et al. 2010; Antoniadis et al. 2013; Fonseca et al. 2016; Arzoumanian et al. 2018; Cromartie et al. 2019) . These observations provided a very strong lower bound of ∼ 2M on the maximum mass of nonrotating NSs that all the competing EoS models from nuclear physics must satisfy. Furthermore, GW170817 (Abbott et al. 2017 ) has prompted several studies predicting an upper bound of ∼ 2.2 − 2.3M on M max of nonrotating NSs, based on the mass ejecta, kilonova signal and absence of a prompt collapse (Shibata et al. 2017; Margalit & Metzger 2017; Ruiz et al. 2018; Rezzolla et al. 2018; Shibata et al. 2019) . While the simultaneous mass-radius measurements of PSR J0030+0451 by NICER collaboration (Riley et al. 2019; Miller et al. 2019 ) indicate a tilt towards slightly stiffer EoS (Raaijmakers et al. 2020; Landry et al. 2020) , the distribution of m 2 would require even higher M max . Possible formation channels of GW190814-type binaries have also been studied in some recent works (Zevin et al. 2020; Safarzadeh & Loeb 2020; Kinugawa et al. 2020) . While there is a general consensus that the fallback of a significant amount of bound su-pernova ejecta on the secondary compact remnant leads to its formation in the lower mass-gap region, the nature of its state at the time of the merger being a BH or a NS remains unclear. Nevertheless, GW190814 has motivated experts to reevaluate the knowledge of dense matter and stellar structure to determine the possible scenarios in which one can construct such configurations of NSs while satisfying relevant constraints (Most et al. 2020; Zhang & Li 2020; Fattoyev et al. 2020; Tsokaros et al. 2020; Tews et al. 2020; Lim et al. 2020; Dexheimer et al. 2020; Sedrakian et al. 2020; Godzieba et al. 2020; Huang et al. 2020 ). Most of these works suggest rapid uniform rotation with or without exotic matter, such as hyperons or quark matter, exploiting the caveat that the spin of m 2 is unconstrained. Other possibilities such as m 2 being a primordial black hole (Vattis et al. 2020; Jedamzik 2020; Clesse & Garcia-Bellido 2020) , an anisotropic object (Roupas 2020) [see also (Biswas & Bose 2019 ) for a detailed study on anisotropic object] or a NS in scalar-tensor gravity (Rosca-Mead et al. 2020) have also been considered.
In this Letter, we investigate the possibility of the GW190814 secondary being a NS within a generic parameterized nuclear EoS, and study its releated properties under assumptions of it being both slowly and rapidly rotating. We also constrain its spin using a universal relation developed in (Breu & Rezzolla 2016 ).
The mass of the secondary object in GW190814 measured by the LVC falls into the so called "mass gap" region (Bailyn et al. 1998; Özel et al. 2010 ) and, therefore, demands a careful inspection of its properties before it can be ruled out as a BH or NS. A non-informative measurement of the tidal deformability or the spin of the secondary, or the absence of an electromagnetic counterpart associated with this event, have made it difficult to make a robust statement about the nature of this object. We begin by examining if the GW mass measurement along with a nuclear-physics informed EoS alone can rule it out as a NS. In Fig. 1 the posterior distribution of secondary mass m 2 is plotted, in blue, by using publicly available LVC posterior samples 1 . In orange, the posterior distribution of M max is overlaid from the empiri-1 LVK collaboration, posterior samples for GW190814 are available here cal EoS analysis by Biswas et al. (2020) using combined GW170817 2 , GW190425 3 and NICER 4 data. (Biswas et al. 2020) , is shown in orange. The distribution shown in green is obtained with the same EoS samples as for the orange one, but considering uniform NS rotation at 716 Hz. These two distributions are compared with the probability distribution of the secondary's mass m2 (in blue) deduced from the GW190814 posterior samples in Ref. Abbott et al. (2020) .
Given these two distributions -both for nonrotating stars -we calculate the probability of m 2 being greater than M max , i.e., P (m 2 > M max ) = P (m 2 − M max ). This probability can be easily obtained by calculating the convolution of the m 2 and −M max probability distributions, which yields P (m 2 > M max ) = 0.91. Therefore, the mass measurement implies that the probability that the secondary object in GW190814 is a NS is 8.9%. However, this type of analysis is highly sensitive to the choice of EoS parameterization as well as on the implementation of the maximum-mass constraint obtained from the heaviest pulsar observations. Biswas et al.
(2020) used a hard 2M cut-off for maximum mass, whereas, the LVC analysis (Abbott et al. 2020) which is based on the spectral EoS parameterization (Lindblom 2010), obtained ∼ 3% probability for the secondary to be a NS using GW170817-informed EoS samples from Abbott et al. (2018) . The addition of NICER data might increase this probability. added NICER data in their analysis of GW observations based on a nonparametric EoS and also examined the impact of different assumptions about the 2 LVK collaboration, https://dcc.ligo.org/LIGO-P1800115/public 3 LVK collaboration, https://dcc.ligo.org/LIGO-P2000026/public 4 Released mass-radius samples released by Miller et al. (Miller et al. 2019) , https://zenodo.org/record/3473466#.XrOt1nWlxBc compact object mass distribution. The P (m 2 > M max ) probabilities technically depend on the mass prior assumed for the secondary, but showed that, regardless of assumed population model, there is a less than ∼ 6% probability for the GW190814 secondary to be a NS. In the discovery paper, LVC also reported an EoS-independent result using the pulsar mass distribution, following Farr & Chatziioannou (2020) , which suggests that there is less than ∼ 29% probability that the secondary is a NS. Despite the differences inherent to these studies, they all suggest that there is a small but finite probability of the secondary object in GW190814 to be a neutron star. It is also important to note that they all assumed the NS to be either nonrotating or slowly rotating (χ < 0.05).
Another possibility is that the secondary object is a rapidly rotating NS (Most et al. 2020; Tsokaros et al. 2020) . It is known that uniform rotation can increase the maximum mass of a NS by ∼ 20% (Friedman & Ipser 1987; Cook et al. 1992 Cook et al. , 1994 . Therefore, rapid rotation may improve the chances that the GW190814 data are consistent with a neutron star.
From pulsar observations, we know that NSs with spin frequencies as high as ν obs max = 716 Hz exist in nature (Hessels et al. 2006) . Using this value for the spin frequency and the EoS samples of Biswas et al. (2020) we can deduce the maximum improvement in probability that the GW190814 secondary is a NS. We used this information in the RNS code (Stergioulas & Friedman 1995) and obtained a corresponding distribution of maximum mass denoted as M 716Hz max . The superscript "716 Hz" emphasizes that all configurations here are computed at that fixed spin frequency. In Fig. 1 , the distribution of M 716Hz max is shown in green. From the overlap of this distribution with P (m 2 ), we find there is ∼ 20% probability that m 2 is a rapidly rotating NS. However, this probability is obtained on the basis of two sets of choices: (a) EoS samples used here were obtained using combined GW-NICER data and a hard upper cut-off of 2M for nonrotating NS. (b) All the NSs are assumed to be rotating at fixed spin frequency of 716 Hz.
Alternatively, if the GW190814 secondary were indeed a NS, then the LVC mass measurement sets a lower limit on the maximum NS mass for any spin at least up to ν obs max . We next relax this constraint by considering all theoretically allowed values of the spin frequency, which for some masses and EoSs may exceed the maximum observed value. In the next two sections, we investigate the properties of NSs -for various rotational frequencies -using a Bayesian approach based on an empirical EoS parameterization.
For slowly rotating NS, a Bayesian methodology was already developed in Biswas et al. (2020) by combining multiple observations based on an empirical EoS parameterization. In this Letter, we relax the assumption of the 2M constraint used in Biswas et al. (2020) and instead use the m 2 distribution of GW190814 as the maximum-mass threshold. We use Gaussian kerneldensity to approximate the posterior distribution of m 2 . The resulting posteriors of radius (R 1.4 ) and tidal deformability (Λ 1.4 ) obtained from this analysis are plotted in Fig. 2 . We find that R 1.4 = 12.7 +0.6 −0.7 km and Λ 1.4 = 616 +195
−177 , at 90% CI, which are in good agreement with previous studies (Abbott et al. 2020; Tews et al. 2020) .
Interestingly, the posterior distributions of R 1.4 and Λ 1.4 obtained in this study are quite similar with those obtained previously in Biswas et al. (2020) . This is most likely because our empirical EoS parameterization already supports stiff EoS, i.e., it is in some tension with previous astrophysical observations. The ddition of GW190814 makes the EoS stiffer, especially in the high density region since now a very small subspace of the EoS family can support a ∼ 2.6M NS. In the right panel of Fig. 2 the posterior of the pressure inside the neutron star is plotted as a function of energy density. This plot clearly shows that the addition of GW190814 places a very tight constraint on the high-density part of the EoS.
In this Letter, for the first time, we develop a Bayesian formalism to constrain the EoS of NS that allows for rapid rotation. At first, given an EoS we create a sequence of mass-radius-tidal deformability up to M TOV max for nonrotating stars. Then, we use a universal relation found by Breu & Rezzolla (2016) which relates the maximum mass of a uniformly rotating star (M rot rmax ) with the maximum mass of a nonrotating star (M TOV max ) for the same EoS,
(1) where a 1 = 0.132 and a 2 = 0.071. χ is the dimensionless spin magnitude of a uniformly rotating star and χ kep is the maximum allowed dimensionless spin magnitude at the mass-shedding limit. Given a χ/χ kep value, we calculate M rot rmax using this universal relation. Its use makes our computation much faster but can cause up to ∼ 2% deviation from the exact result, as noted by Breu & Rezzolla (2016) . We assume that the error is constant throughout the parameter space; we took it to be distributed uniformly in [−2%, 2%] and marginalized over it to get an unbiased estimate of the properties of the object.
We combine data from two other binary neutron stars, namely GW170817 and GW190425, as well as NICER data assuming nonrotating NS following Biswas et al. (2020) . We assume that the Cromartie et al. pulsar (Cromartie et al. 2019 ) is the heaviest possible nonrotating NSs and use a Gaussian likelihood of 2.14 M median and 0.1 M of 1σ credible interval to place a maximum-mass threshold on M TOV max . Then, the m 2 distribution of GW190814 is used for the maximummass threshold of a uniformly rotating star, i.e., M rot rmax . We use a nested sampler algorithm implemented in Pymultinest (Buchner et al. 2014 ) to simultaneously sample the EoS parameters and χ/χ kep . These posterior samples are then used in the RNS code (Stergioulas & Friedman 1995) to calculate several properties of the secondary object associated with GW190814.
In the upper left and middle panel of Fig. 3 posterior distributions of equatorial radius (R e ) and ellipticiy (e) are plotted, respectively. Within the 90% CI we find R e = 15.7 +1.0 −1.7 km and e = 0.66 +0.03 −0.08 . Such high values of equatorial radius and ellipticity imply a considerable deviation from a spherically symmetric static configuration. From the distribution of χ shown in the upper left panel of Fig. 3 we find its value to be χ = 0.64 +0.05 −0.10 . Most et al. (2020) have also obtained a similar bound on χ with simpler arguments. In this Letter, we provide a distribution for χ employing a Bayesian framework as well as place a more robust bound on this parameter. In the lower left panel of Fig. 3 , the posterior distribution of rotational frequency is plotted in Hz. We find its value to be f = 1143 +194 −155 Hz. As noted above, till date PSR J1748-2446a (Hessels et al. 2006 ) is known as the fastest rotating pulsar, with a rotational frequency of 716 Hz. Therefore, if the secondary of GW190814 is indeed a rapidly rotating NS, it would definitely be the fastest rotating NS observed so far. In the lower-middle and right panels, the posterior distributions of the moment of inertia and quadrupole moment of the secondary are shown, respectively.
EoS constraints derived from the observation of nonrotating NSs also provide an upper bound on the maximum spin of a NS. The maximum spin frequency is given empirically as f lim 1 2π (0.468 + 0.378χ s ) GMmax Rmaxc 2 , with M max and R max being the maximum mass and its corresponding radius of a nonrotating NS, respectively. We use M max − R max posterior samples that were deduced in Biswas et al. (2020) by using combined GWs and NICER data to calculate f lim . In the left panel of Fig. 4 , its distribution is shown by the grey shaded region. We overlay that distribution with distributions of frequencies of the secondary object of GW190814 and those of a few hypothetical rotating NSs of various masses -all Gaussian distributed, but with medians of 2.4 M , 2.8 M and 3.0 M , respectively, and each having a measurement uncertainty of 0.1M . We also assume the primary component of GW190425 to be a rapidly rotating NS, since by using a high-spin prior LVC determined its mass to be 1.61M − 2.52M . In our calculations, for GW190425 we used the publicly available high-spin posterior of m 1 obtained by using the PhenomPNRT waveform. We find observations like m 1 of GW190425 and simulations like N (2.4M , 0.1M ) correspond to posteriors of rotational frequency that are comparatively lower than limiting values of rotational frequencies. However, as the mass increases, the posterior of frequency eventually almost coincides with f lim . Therefore, if the secondary of GW190814 was a rapidly rotating NS, it would have to be rotating rather close to the limiting frequency.
Any rotating star is generically unstable through the Chandrasekhar-Friedman-Schutz (CFS) mechanism (Chandrasekhar 1970; Friedman & Schutz 1978) . This instability occurs when a certain retrograde mode in the rotating frame becomes prograde in the inertial frame. For example, the f -modes of a rotating NS can always be made unstable for a sufficiently large mode number m (not to be confused with component masses m 1,2 ) even for low spin frequencies, but, the instability timescale increases rapidly with the increase of m. Numerical calculations have shown (Stergioulas & Friedman 1998; Morsink et al. 1999) , that for maximum mass stars m = 2 mode changes from retrograde to prograde at T /|W | ∼ 0.06, where T is the rotional energy and W the gravitational potential energy of the NS. We computed this ratio for all the cases considered in this section and plot the distributions in the right panel of Fig. 4 . From this analysis we find that the secondary of GW190814 should be f −mode unstable as for most of the allowed EoSs T /|W | is significantly larger than 0.06. The CFS instability is even more effective for r−modes (Lindblom et al. 1998; Andersson et al. 1999 ) as they are generically unstable for all values of spin frequency. However, an instability can develop, only if its growth timescale is shorter than the timescale of the strongest damping mechanism affecting it. A multitude of damping mechanisms, such as shear viscosity, bulk viscosity, viscous boundary layer, crustal resonances and superfluid mutual friction (each having each own temperature dependence) have been investigated (see (Kokkotas & Schwenzer 2016; Paschalidis & Stergioulas 2017; Andersson 2019) and references therein). The spin-distribution of millisecond pulsars in accreting systems (Papitto et al. 2014) can be explained, if the r-mode instability is effectively damped up to spin frequencies of ∼ 700 Hz, (Ho et al. 2011 ) and operating at higher spin rates. This would not allow for the secondary in GW190814 to be a rapidly rotating NS at the limiting spin frequency.
On the other hand, if the secondary of GW190814 was a rapidly rotating NS at the limiting frequency, then the f -mode and r-mode instabilities must be effectively damped both during the spin-up phase in a low-mass-X-ray binary, where it acquires rapid rotation, as well as during its subsequent lifetime up to the moment of merger. This might be possible, if both the f -mode and the r-mode instabilities are damped by a particularly strong mutual friction of superfluid vortices below the superfluid transition temperature of ∼ 10 9 K (see (Lindblom & Mendell 2000; Gaertig et al. 2011 ) and in particular the case of an intermediate drag parameter R ∼ 1 in (Haskell et al. 2009) ). If this is the case, then the limiting frequency observed in the spin-distribution of millisecond pulsars must be explained by other mechanisms, see (Gittins & Andersson 2019 The distribution of f lim is plotted considering three simulated rapidly rotating NS whose mass measurements are Gaussian distributed with median 2.4 M , 2.8 M and 3.0 M , respectively and each having a measurement uncertainty of .1M . The same has been overlaid using the secondary component of the GW190814 and the primary of the GW190425 events. In the right panel, the corresponding ratio of rotational to gravitational potential energy T /|W | is shown.
would have strong implications on the physics of superfluidity in neutron star matter (in particular constraining the drag parameter R of mutual friction) and on the astrophysics of accreting systems. So far, we have analyzed the impact on NS EoS properties arising from the hypothesis that the secondary object in GW190814 is a NS. On the other hand, if that secondary object is a BH, then again novel information about the NS EoS can be obtained, since it will set an upper bound on the NS maximum mass, but only if one were to assume that the NS and BH mass distributions do not overlap. In our analysis, we take this value to be 2.5M , which is the lowest possible value of the secondary object within 90% CI. Then, using Bayesian inference for nonrotating stars, we combine GW and NICER data and also use the maximum mass constraint from massive pulsar observations (in a manner similar to what was done in Sec. 4 for nonrotating stars) to place further constraints on the NS EoS. In Fig. 5 , the distribution of the maximum mass for nonrotating NSs is shown in blue using the EoS samples obtained from this analysis. Within the 90% CI we find M max = 2.12 +0.16 −0.15 M , which is the most conservative bound on NS maximum mass obtained so far in this work.
Assuming a high-spin prior, the mass of the primary component of GW190425 is constrained between 1.61 − 2.52M . In Fig. 5 , its distribution is over-plotted in orange. From the overlap with the newly obtained M max distribution and the m 1 distribution of GW190425, we find that there is ∼ 49% probability that the primary of GW190425 is a BH.
Based on the maximum mass samples obtained from Biswas et al. (2020) , we find that there is ∼ 9% probability that the secondary object associated with GW190814 is a nonrotating NS. However, such an estimation depends on the choice of EoS parameterization and the maximum mass threshold. Nevertheless, the possibility of the secondary being a nonrotating NS is not inconsistent with the data. Based on our empirical EoS parameterization, we find that the addition of GW190814 as an upper limit of maximum mass for nonrotating stars provides a very stringent constraint on the EoS specially in the high density region. We also discussed the alternative that the secondary is a rapidly rotating NS. Based on the assumption that J0740 + 6620 (Cromartie et al. 2019 ) provides the maximum mass threshold for nonrotating NSs, we find that in order to satisfy the secondary mass estimate of GW190814, its spin magnitude has to be close to the limiting spin frequency for uniform rotation. In fact, it would be the fastest rotating NS ever observed. However, this could be the case, only if gravitational-wave instabilities are effectively damped for rapidly rotating stars, which opens the possibility of constraining physical mechanisms, such as mutual friction in a superfluid interior.