Impact of hyperon mixing on neutron star structure based on Skyrme-type equations of state: Systematic analysis of $\Lambda NN$ and $\Lambda\Lambda N$ three-body forces with Bayesisan inference

We study hyperonic density-dependent three-body effects in cold neutron-star matter using a Skyrme energy-density-functional framework. In beta-equilibrated $npe\mu\Lambda$ matter, the effective $\Lambda NN$ and $\Lambda\Lambda N$ terms are varied separately in the $(\beta,A_3)$ and $(\gamma,C_3)$ planes, and each tabulated equation of state is used in Tolman--Oppenheimer--Volkoff calculations. The calculated $P$--$\varepsilon$ branches are classified by monotonicity and extremum structure. The $\Lambda\Lambda N$ term does not affect the $\Lambda$-onset condition, but modifies the finite-$\Lambda$ post-onset EOS: increasing $C_3$ generally stiffens the post-onset branch and raises $M_{\max}$ in mechanically admissible regions, whereas increasing $\gamma$ reduces this enhancement at fixed $C_3$. In contrast, the $\Lambda NN$ term shifts the $\Lambda$-onset density and modifies the post-onset EOS simultaneously, producing organized branch-limited and Maxwell-candidate regions for some reference interactions. Representative two-extrema cases are examined with Maxwell constructions. We also perform an exploratory Bayesian analysis using neutron-star mass--radius information alone and apply XGBoost--SHAP surrogate diagnostics to summarize parameter sensitivities. Within the adopted likelihood and prior ranges, the posterior weight tends to favor sizable hyperonic three-body repulsion, and the SHAP analysis identifies $A_3$ and $C_3$ as important controls of $M_\text{max}$ and $R_{2.0}$. These results show that maximum-mass recovery in hyperonic neutron stars is not a single mechanism: $M_\text{max}$ maps must be interpreted together with onset behavior, branch admissibility, and extremum-count diagnostics. *shortened due to the arXiv's word limit.