The detection of strong Balmer breaks and absorption features in Little Red Dots (LRDs) suggests they host AGN embedded within dense gas envelopes, potentially powered by super-Eddington accretion. We present GLIMPSE-17775, a luminous ($L_{\rm bol}\sim10^{45}$ erg s$^{-1}$) LRD at $z=3.501$ behind Abell S1063 ($\mu\sim2$), observed with deep JWST/NIRCam and a $\sim$20 hr (80 hr de-lensed) NIRSpec/G395M spectrum. The data reveal 40+ emission and absorption features, including a rich forest of low-ionization FeII lines and numerous broad hydrogen recombination transitions. We use this depth to test the dense-gas interpretation through five independent diagnostics. Nearly all permitted lines show exponential wings with consistent FWHM, the signature of Thomson scattering requiring $n_e\gtrsim10^8$ cm$^{-3}$. Adopting this width yields $M_{\rm BH}\sim10^{6.7}M_\odot$, a factor of ten lower than Gaussian fits, and $\lambda_{\rm Edd}\sim1.8$. Additional diagnostics support the same picture: a pronounced Balmer break ($f_{\nu,4050}/f_{\nu,3670}=2.0\pm0.1$), enhanced HeI $\lambda7065$ and $\lambda10830$ with P-Cygni absorption, Bowen-fluorescent OI $\lambda8446$-$\lambda11290$ emission requiring Ly$\beta$ pumping, and 16 FeII lines matching fluorescence models. These features indicate a dense ($n\sim10^8$ cm$^{-3}$), partially ionized cocoon where scattering and fluorescence dominate line formation, providing strong evidence that at least some LRDs are powered by super-Eddington black-hole growth in the early Universe.

We present a sample of 341 "little red dots" (LRDs) spanning the redshift range $z\sim2-11$ using data from the CEERS, PRIMER, JADES, UNCOVER and NGDEEP surveys. Unlike past use of color indices to identify LRDs, we employ continuum slope fitting using shifting bandpasses to sample the same rest-frame emission blueward and redward of the Balmer break. This enables the detection of LRDs over a wider redshift range and with less contamination from galaxies with strong breaks that otherwise lack a rising red continuum. The redshift distribution of our sample increases at z&lt;8 and then undergoes a rapid decline at $z\sim4.5$, which may tie the emergence of these sources to the inside-out growth that galaxies experience during this epoch. We find that LRDs are $\sim1$ dex more numerous than X-ray and UV selected AGN at z~5-7. Within our sample, we have identified the first two X-ray detected LRDs. An X-ray spectral analysis confirms that these AGN are moderately obscured with $\log\,(N_{\rm H}/{\rm cm}^{2}$) of $23.3^{+0.4}_{-1.3}$ and $22.72^{+0.13}_{-0.16}$. Our analysis reveals that reddened AGN emission dominates their rest-optical light, while the rest-UV originates from their host galaxies. We also present NIRSpec observations from the RUBIES survey of 17 LRDs that show broad emission lines consistent with AGN activity. The confirmed AGN fraction of our sample is 71\% for sources with F444W<26.5. In addition, we find three LRDs with blue-shifted Balmer absorption features in their spectra, suggesting an outflow of high-density, low-ionization gas from near the central engine of these faint, red AGN.

The inverse cascade in MHD turbulence plays a crucial role in various astrophysical processes such as galaxy cluster formation, solar and stellar dynamo mechanisms, and the evolution of primordial magnetic fields in the early universe. A standard numerical approach involves injecting magnetic helicity at intermediate length scales to generate a secondary, time-dependent spectral peak that gradually propagates toward larger scales. Previous simulations have already suggested a resistive dependence of inverse transfer rates and demonstrated the significant influence of magnetic helicity flux density $\epsilon_\mathrm{H}$ on this process. On dimensional grounds, we have $E_\mathrm{M}(k,t)=C_\mathrm{H} \epsilon_\mathrm{H}^{2/3} k^{-1}$ where $C_\mathrm{H}$ represents a potentially universal dimensionless coefficient analogous to the Kolmogorov constant. We present a summary of the 25 distinct simulations conducted with the \textsc{Pencil Code}, systematically varying the forcing wavenumber $k_\mathrm{f}$, magnetic Prandtl number $Pm$, grid resolution $N^3$, and Lundquist number $Lu$. We obtained $C_\mathrm{H}$ and corresponding error bars by calculating the compensated spectrum and investigated its dependence with $Lu$ and $k_\mathrm{f}$. For the $C_\mathrm{H}$ - $Lu$ relationship, we observe strong correlations with power-law exponents of 1 and 2/3. In contrast, we find no significant correlation between $C_\mathrm{H}$ and $k_\mathrm{f}$.

$N$-body simulations are computationally expensive, so machine-learning (ML)-based emulation techniques have emerged as a way to increase their speed. Although fast, surrogate models have limited trustworthiness due to potentially substantial emulation errors that current approaches cannot correct for. To alleviate this problem, we introduce COmoving Computer Acceleration (COCA), a hybrid framework interfacing ML with an $N$-body simulator. The correct physical equations of motion are solved in an emulated frame of reference, so that any emulation error is corrected by design. This approach corresponds to solving for the perturbation of particle trajectories around the machine-learnt solution, which is computationally cheaper than obtaining the full solution, yet is guaranteed to converge to the truth as one increases the number of force evaluations. Although applicable to any ML algorithm and $N$-body simulator, this approach is assessed in the particular case of particle-mesh cosmological simulations in a frame of reference predicted by a convolutional neural network, where the time dependence is encoded as an additional input parameter to the network. COCA efficiently reduces emulation errors in particle trajectories, requiring far fewer force evaluations than running the corresponding simulation without ML. We obtain accurate final density and velocity fields for a reduced computational budget. We demonstrate that this method shows robustness when applied to examples outside the range of the training data. When compared to the direct emulation of the Lagrangian displacement field using the same training resources, COCA's ability to correct emulation errors results in more accurate predictions. COCA makes $N$-body simulations cheaper by skipping unnecessary force evaluations, while still solving the correct equations of motion and correcting for emulation errors made by ML.

Very recent work has identified a new satellite galaxy, Ursa Major III/UNIONS I, which is the faintest such system ever observed. Dynamical considerations indicate that if the system is in equilibrium, it is likely to be highly dark matter dominated. This, in combination with its proximity, predicts that it may be the preeminent dwarf spheroidal galaxy target for dark matter indirect detection searches. We utilize 15 years of Fermi-LAT data to search for $\gamma$-ray emission from Ursa Major III. Finding no excess, we set strong constraints on dark matter annihilation. Intriguingly, if the high J-factor of Ursa Major III is confirmed, standard thermal dark matter annihilation to $b\bar{b}$ final states would be ruled out for dark matter masses up to 4 TeV. The discovery of Ursa Major III, combined with recent tentative measurements of other high J-factor systems, suggests the exciting possibility that near-future data could produce transformative constraints on thermal dark matter.

In this work, we emphasize that it is necessary to take into account one-loop corrections of $2.0\%$ to the neutrino matter potential in the precision measurements of neutrino oscillation parameters and in the experimental searches for new physics beyond the Standard Model. With the numerical simulation of the DUNE experiment, we study how radiative corrections to the matter potential affect neutrino oscillation probabilities, and thus, the event rates in the presence of neutrino non-standard interactions (NSIs). We find that neglecting one-loop corrections may lead to wrong conclusions for the discovery of NSIs. The implications for the determination of neutrino mass ordering and constraints on the NSI parameters in future long-baseline accelerator neutrino experiments are explored in a quantitative way.

The recent discovery of a large number of massive black holes within the first two billion years after the Big Bang, as well as their peculiar properties, have been largely unexpected based on the extrapolation of the properties of luminous quasars. These findings have prompted the development of several theoretical models for the early formation and growth of black holes, which are, however, difficult to differentiate. We report the metallicity measurement around a gravitationally lensed massive black hole at redshift 7.04 (classified as a Little Red Dot), hosted in a galaxy with very low dynamical mass. The weakness of the [OIII]5007 emission line relative to the narrow H$\beta$ emission indicates extremely low metallicity, about $4\times 10^{-2}$ solar, and even more metal poor in the surrounding few 100 pc. We argue that such properties cannot be uncommon among accreting black holes around this early cosmic epoch. Explaining such a low chemical enrichment in a system that has developed a massive black hole is challenging for most theories. Models assuming heavy black hole seeds (such as Direct Collapse Black Holes) or super-Eddington accretion scenarios struggle to explain the observations, although they can potentially reproduce the observed properties in some cases. Models invoking "primordial black holes" (i.e. putative black holes formed shortly after the Big Bang) may potentially explain the low chemical enrichment associated with this black hole, although this class of models also requires further developments for proper testing.

We revisit the possibility that Dark Matter is composed of stable scalar glueballs of a confining dark ${\rm SU}(3)$ gauge theory coupled only to gravity. The relic abundance of dark glueballs is studied for the first time in a thermal effective theory accounting for strong-coupling dynamics. An important ingredient of our analysis is the use of an effective potential for glueballs that is fitted by lattice simulations. We predict the relic abundance to be in the range $0.12\zeta_{T}^{-3}\Lambda/(137.9 {\rm eV}) \lesssim \Omega h^{2}\lesssim 0.12\zeta_{T}^{-3}\Lambda/(82.7 {\rm eV})$, with$\Lambda$ being the confinement scale, $\zeta_{T}$ the visible-to-dark sector temperature ratio and the uncertainty is coming from the fit to lattice data. This prediction is an order of magnitude smaller than the existing glueball abundance results in the literature. Our framework can be easily generalised to different gauge groups and modified cosmological histories paving the way towards consistent exploration of strongly-coupled dark sectors and their cosmological implications.

We present spatially resolved analysis of stellar populations in the brightest $z>6$ galaxy known to date (AB mag 23), the strongly lensed MACS0308$-$zD1 (dubbed the ``Cosmic Spear'') at $z_{\rm spec}=6.2$. New JWST NIRCam imaging and high-resolution NIRSpec IFU spectroscopy span the rest-frame ultraviolet to optical. The NIRCam imaging reveals bright star-forming clumps and a tail consisting of three distinct, extremely compact star clusters that are multiply-imaged by gravitational lensing. The star clusters have effective radii of $R_{\rm{eff}} \sim 5$ pc, stellar masses of $M_{*} \sim 10^{6}-10^{7}\,M_{\odot}$, and high stellar mass surface densities of $\Sigma_{*} > 10^{4}\,M_{\odot}~\rm{pc}^{-2}$. While their stellar populations are very young ($\sim 5-9$ Myr), their dynamical ages exceed unity, consistent with the clusters being gravitationally bound systems. Placing the star clusters in the size vs.~stellar mass density plane, we find they occupy a region similar to other high-redshift star clusters within galaxies observed recently with JWST, being significantly more massive and denser than local star clusters. Spatially resolved analysis of the brightest clump reveals a compact, intensely star-forming core. The ionizing photon production efficiency ($\xi_{\rm{ion}}$) is slightly suppressed in this central region, potentially indicating a locally elevated Lyman continuum escape fraction facilitated by feedback-driven channels.

The Zwicky Transient Facility (ZTF) was expected to detect more than one strong gravitationally-lensed supernova (glSN) per year, but only one event was identified in the first four years of the survey. This work investigates selection biases in the search strategy that could explain the discrepancy and revise discovery predictions. We present simulations of realistic lightcurves for lensed thermonuclear (glSNIa) and core-collapse supernova (glCCSN) explosions over a span of 5.33 years of the survey, utilizing the actual observation logs of ZTF. We find that the magnitude limit in spectroscopic screening significantly biases the selection towards highly magnified glSNe, for which the detection rates are consistent with the identification of a single object by ZTF. To reach the higher predicted rate of detections requires an optimization of the identification criteria for fainter objects. We find that around 1.36 (3.08) Type Ia SNe (CCSNe) are identifiable with the magnification method per year in ZTF, but when applying the magnitude cut of m < 19 mag, the detection rates decrease to 0.17 (0.32) per year. We compare our simulations with the previously found lensed Type Ia SNe, iPTF16geu and SN Zwicky, and conclude that considering the bias towards highly magnified events, the findings are within expectations in terms of detection rates and lensing properties of the systems. In addition, we provide a set of selection cuts based on simple observables to distinguish glSNe from regular, unlensed, supernovae to select potential candidates for spectroscopic and high-spatial resolution follow-up campaigns. We find optimal cuts in observed colours $g-r$, $g-i$, and $r-i$ as well as in the colour SALT2 fit parameter. The developed pipeline and the simulated lightcurves employed in this analysis can be found in the $LENSIT$ github repository.

Gravitational waves (GWs) can provide crucial information about the central engines of core-collapse supernovae (CCSNe). In order to unveil the nature of GW emission in CCSNe, we apply perturbative analyses with the same underlying equations as simulations to diagnose oscillations of the proto-neutron star (PNS) during $\sim$1 s postbounce. In the pseudo-Newtonian case, we find that radial profiles of GW emission match well between the perturbative analysis with $l=2$ and simulations inside the PNS at \emph{any} frequency and time. This confirms that the GW emission of CCSNe arises from the global PNS oscillations in the perturbative regime. Based on this, we solve for the discrete eigenmodes with a free PNS surface and tentatively identify a set of $g$ modes and the $f$ mode contributing to the peak GW emission. We also offer a possible explanation for the power gap in the GW spectrum found in simulations that lies at the frequency with vanishing cumulative emission of the PNS. Our results enhance the predictive power of perturbative analyses in the GW signals of CCSNe.

We discuss the low-frequency gravitational wave signals from three state-of-the-art three-dimensional core-collapse supernova models produced with the \textsc{Chimera} supernova code. We provide a detailed derivation of the gravitational wave signal sourced from the anisotropic emission of neutrinos and provide the total (fluid sourced and neutrino sourced) gravitational waves signal generated in our models. We discuss the templatablity of this low-frequency signal, which is useful for future work involving matched filtering for signal detection and parameter estimation.

We study the impact of rotation on the multimessenger signals of core-collapse supernovae (CCSNe) with the occurrence of a first-order hadron-quark phase transition (HQPT). We simulate CCSNe with the \texttt{FLASH} code starting from a 20~$M_\odot$ progenitor with different rotation rates, and using the RDF equation of state from \textit{Bastian} 2021 that prescribes the HQPT. Rotation is found to delay the onset of the HQPT and the resulting dynamical collapse of the protocompact star (PCS) due to the centrifugal support. All models with the HQPT experience a second bounce shock which leads to a successful explosion. The oblate PCS as deformed by rotation gives rise to strong gravitational-wave (GW) emission around the second bounce with a peak amplitude larger by a factor of $\sim10$ than that around the first bounce. The breakout of the second bounce shock at the neutrinosphere produces a $\bar{\nu}_e$-rich neutrino burst with a luminosity of serveral 10$^{53}$~erg~s$^{-1}$. In rapidly rotating models the PCS pulsation following the second bounce generates oscillations in the neutrino signal after the burst. In the fastest rotating model with the HQPT, a clear correlation is found between the oscillations in the GW and neutrino signals immediately after the second bounce. In addition, the HQPT-induced collapse leads to a jump in the ratio of rotational kinetic energy to gravitational energy ($\beta$) of the PCS, for which persistent GW emission may arise due to secular nonaxisymmetric instabilities.

The elusive nature of Dark Matter (DM) remains a mystery far from being solved. A vast effort is dedicated to search for signatures of feeble DM interactions with Standard Model particles. In this work, we explore the signatures of axion DM boosted by interactions with Supernova neutrinos: Neutrino-Boosted Axion DM ($\nu$BADM). We focus on $\nu$BADM converting into photons in the Galactic magnetic field, generating a peculiar gamma-ray flux. This signal falls in the poorly explored MeV energy range, that will be probed by next generation gamma-ray missions. Once more, astrophysical searches might act as a probe of fundamental physics, unveiling the nature and properties of DM.

Early studies of the AMS-02 antiproton ratio identified a possible excess over the expected astrophysical background that could be fit by the annihilation of a weakly interacting massive particle (WIMP). However, recent efforts have shown that uncertainties in cosmic-ray propagation, the antiproton production cross-section, and correlated systematic uncertainties in the AMS-02 data, may combine to decrease or eliminate the significance of this feature. We produce an advanced analysis using the DRAGON2 code which, for the first time, simultaneously fits the antiproton ratio along with multiple secondary cosmic-ray flux measurements to constrain astrophysical and nuclear uncertainties. Compared to previous work, our analysis benefits from a combination of: (1) recently released AMS-02 antiproton data, (2) updated nuclear fragmentation cross-section fits, (3) a rigorous Bayesian parameter space scan that constrains cosmic-ray propagation parameters. We find no statistically significant preference for a dark matter signal and set strong constraints on WIMP annihilation to $b\bar{b}$, ruling out annihilation at the thermal cross-section for dark matter masses below $\sim200$~GeV. We do find a positive residual that is consistent with previous work, and can be explained by a $\sim70$~GeV WIMP annihilating below the thermal cross-section. However, our default analysis finds this excess to have a local significance of only 2.8$\sigma$, which is decreased to 1.8$\sigma$ when the look-elsewhere effect is taken into account.

We present the results of low-resolution spectroscopic and densely sampled multiband simultaneous optical imaging ($ugi$ and $vrz$ bands) follow-up of supernova (SN) 2024aecx. The photometric data is supplemented with $Swift$/UVOT and ATLAS survey observations. The SN was discovered in the spiral galaxy NGC 3521 (distance $\sim$11 Mpc) within a day after the explosion. The early spectra of SN 2024aecx show a weak signature of hydrogen lines, which disappeared in $\sim$30 days after the explosion. Light curves in all bands show a distinct feature of two peaks, and the first peak is likely due to the shock cooling emission. The early phase light curve evolution of SN 2024aecx has similarity with the typical Type IIb events, but the decay rate in different bands (e.g., $\rm \Delta m_{15}$ = 1.60 $\pm$ 0.05 mag, $g$-band) is significantly faster in the post-peak phase. It attained the secondary maximum in $\sim$19 days ($g$-band) with a peak absolute magnitude of M$_{g}$= -17.94 $\pm$ 0.10 mag. The color evolution of SN 2024aecx is displaying a red-blue-red trend between days $\sim$8 to 40. The analytical model fitting to the light curves reveals an envelope mass and progenitor radii in the range of $\sim$0.03 - 0.24 $M_\odot$ and $\sim$169 - 200 $R_\odot$, respectively. Modeling of the pseudo-bolometric light curve suggests that synthesized $^{56}$Ni in the explosion was $\sim$0.15 M$_{\odot}$ with ejecta mass and kinetic energy of $\sim$0.7 M$_{\odot}$ and $\sim$0.16 x 10$^{51}$ erg, respectively. The observational properties and modeling indicate that the SN 2024aecx progenitor belongs to the extended progenitor category.

The stellar mass-star formation rate ($\mathrm{M_*}$-$\mathrm{SFR}$) plane is a fundamental diagnostic for distinguishing galaxy populations. However, the evolutionary pathways of galaxies within this plane across cosmic time remain poorly understood. This study aims to observationally characterize galaxy migration in the $\mathrm{M_*}$-$\mathrm{SFR}$ plane using reconstructed star formation histories (SFHs) of galaxies at z &lt; 4. Our goal is to provide insights into the physical processes governing star formation and quenching. We analyze a sample of 299,131 galaxies at z &lt; 4 from the COSMOS-Web NIRCam survey (m_{\mathrm{F444W}} &lt; 27, 0.54 deg$^2$). Using non-parametric SFH modeling with CIGALE, we derive physical properties and reconstruct SFHs. To trace galaxy evolution, we define migration vectors, quantifying their direction ($\Phi_{\mathrm{dt}}$ [deg]) and velocity norm ($r_{\mathrm{dt}}$ [dex/Gyr]) on the $\mathrm{M_*}$-$\mathrm{SFR}$ plane. The reliability of these vectors is assessed using the Horizon-AGN simulation. We find that main-sequence galaxies exhibit low-amplitude migration with scattered directions, suggesting oscillations within the main sequence. Their progenitors predominantly lie on the main sequence 1 Gyr earlier. Starburst galaxies show rapid mass assembly ($50\%$ within 350 Myr) and originate from the main sequence, while passive galaxies display uniformly declining SFHs. Massive passive galaxies emerge as early as 3.5 &lt; z &lt; 4, increasing in number density over time. Only &lt;20\% of passive galaxies were starbursts 1 Gyr prior, indicating diverse quenching pathways. By reconstructing SFHs to z &lt; 4, we present a coherent picture of galaxy migration in the $\mathrm{M_*}$-$\mathrm{SFR}$ plane, linking evolutionary phases to their star formation signatures.

Abridged - Stars with ZAMS masses between 140 and $260 M_\odot$ are thought to explode as pair-instability supernovae (PISNe). During their thermonuclear runaway, PISNe can produce up to several tens of solar masses of radioactive nickel, resulting in luminous transients similar to some superluminous supernovae (SLSNe). Yet, no unambiguous PISN has been discovered so far. SN2018ibb is a H-poor SLSN at $z=0.166$ that evolves extremely slowly compared to the hundreds of known SLSNe. Between mid 2018 and early 2022, we monitored its photometric and spectroscopic evolution from the UV to the NIR with 2-10m class telescopes. SN2018ibb radiated $>3\times10^{51} \rm erg$ during its evolution, and its bolometric light curve reached $>2\times10^{44} \rm erg\,s^{-1}$at peak. The long-lasting rise of$>93$ rest-frame days implies a long diffusion time, which requires a very high total ejected mass. The PISN mechanism naturally provides both the energy source ($^{56}$Ni) and the long diffusion time. Theoretical models of PISNe make clear predictions for their photometric and spectroscopic properties. SN2018ibb complies with most tests on the light curves, nebular spectra and host galaxy, potentially all tests with the interpretation we propose. Both the light curve and the spectra require 25-44 $M_\odot$ of freshly nucleosynthesised $^{56}$Ni, pointing to the explosion of a metal-poor star with a He-core mass of 120-130 $M_\odot$ at the time of death. This interpretation is also supported by the tentative detection of [Co II]$\lambda$1.025$\mu$m, which has never been observed in any other PISN candidate or SLSN before. Powering by a central engine, such as a magnetar or a black hole, can be excluded with high confidence. This makes SN2018ibb by far the best candidate for being a PISN, to date.

Energetic nuclear outbursts have been discovered in luminous and ultraluminous infrared galaxies (U/LIRGs) at unexpectedly high rates. To investigate this population of transients, we performed a search in mid-IR data from the Wide-field Infrared Survey Explorer (WISE) satellite and its NEOWISE survey to detect and characterise luminous and smoothly evolving transients in a sample of 215 U/LIRGs. We report three new transients, all with $\Delta L > 10^{43}$ erg s$^{-1}$, in addition to two previously known cases. Their host galaxies are all part of major galaxy mergers, and through radiative transfer model fitting we find that all have a significant contribution from an active galactic nucleus (AGN). We characterised the transients through measurements of their luminosities and resulting energetics, all of which are between 10$^{50.9}$ erg and 10$^{52.2}$ erg. The IR emission of the five transients was found to be consistent with re-radiation by the hot dust of emission at shorter wavelengths, presumably originating from an accretion event, onto the supermassive black hole. The corresponding transient rate of (1.6-4.6)$\times$10$^{-3}$ / yr / galaxy is over an order of magnitude higher than the rate of large amplitude flares shown by AGN in the optical. We suggest that the observed transients are part of a dust-obscured population of tidal disruption events (TDEs) that have remained out of the reach of optical surveys due to the obscuring dust. In one case, this is supported by our radio observations. We also discuss other plausible explanations. The observed rate of events is significantly higher than optical TDE rates, which can be expected in U/LIRG hosts undergoing a major galaxy merger with increased stellar densities in the nuclear regions. Continued searches for such transients and their multi-wavelength follow-up is required to constrain their rate and nature.

Gravitational waves (GWs) from binary neutron stars (BNSs) are expected to be accompanied by electromagnetic (EM) emissions, which help to identify the host galaxy. Since GW events directly measure their luminosity distances, joint GW-EM observations from BNSs help to study cosmology, particularly the Hubble constant, unaffected by cosmic distance ladder systematics. However, detecting EM counterparts from BNS mergers is not always possible. Additionally, the tidal deformations of BNS components offer insights into the neutron star (NS) equation of state (EoS). In such cases, the tidal parameters of NSs, combined with the knowledge of the NS EoS, can break the degeneracy between mass parameters and redshift, allowing for the inference of the Hubble constant. Several efforts have aimed to infer the Hubble constant using the tidal parameters of BNSs, without EM counterparts, termed dark BNSs. Moreover, some studies have focused on the joint estimation of population and NS EoS for unbiased NS EoS estimation. However, none of the works consistently combined the uncertainties of population, cosmology, and NS EoS within a Bayesian framework. In this study, we propose a novel Bayesian analysis to jointly constrain the NS EoS, population, and cosmological parameters using a population of dark BNSs detected through GW observations. This method can well constrain the Hubble constant with as few as $5$ BNS observations using current-generation detectors. This level of precision is unattainable without incorporating the NS EoS, especially when observing BNS mergers without EM counterpart information.