Longhua Qin, Jiancheng Wang, Chuyuan Yang, Huaizhen Li, Ao Wang, Weiwei Na, Hushan Xu, Xiaogu Zhong, Zunli Yuan, Yubin Li, Guangbo Long

We construct a multi-epoch joint profile-likelihood framework to constrain axion-like particles (ALPs), employing a Gaussian copula ($ρ=0.03$) to account for weak correlations among distinct observational periods. Different from traditional independent likelihood addition, our scheme globally optimizes nuisance parameters under the no-ALP hypothesis prior to profiling ALP mass $m_a$ and coupling $g_{aγ}$, enhancing numerical stability and statistical consistency. Using five sets of simultaneous Fermi-LAT and LHAASO data of Mrk 421, we implement both two-zone hybrid and single-zone hadronic jet models. In the CAST-sensitive low-mass regime $m_a\lesssim1$ neV, the 95\% CL upper limit of $g_{aγ}$ is $8.33\times10^{-14}\,\mathrm{GeV}^{-1}$. Over $0.1$-$500$ neV, conservative bounds read $g_{aγ}<6.50\times10^{-12}\,\mathrm{GeV}^{-1}$ (two-zone) and $g_{aγ}<7.30\times10^{-12}\,\mathrm{GeV}^{-1}$ (single-zone). These results benefit from LHAASO's broadband very high energy coverage and long term monitoring of Mrk 421; the developed framework provides a robust and scalable pipeline suitable for future multi-messenger ALP searches with next-generation gamma-ray observatories.