Identifying latent variables and the causal structure involving them is essential across various scientific fields. While many existing works fall under the category of constraint-based methods (with e.g. conditional independence or rank deficiency tests), they face common empirical challenges such as testing-order dependency, error propagation, and the difficulty in choosing an appropriate significance level. These issues can potentially be mitigated by properly designed score-based methods, such as Greedy Equivalence Search (GES) (Chickering, 2002) in the specific case without latent variables. Yet, formulating score-based methods with latent variables is highly challenging. This work is, to the best of our knowledge, the first score-based method that is capable of identifying a causal structure containing causally-related latent variables with identifiability guarantees. Specifically, we show that a properly formulated BIC score can achieve score equivalence and consistency with latent variables in structure learning, though it was originally designed without considering latent variables. We further provide a rigorous characterization of the degrees of freedom for the marginal over the observed variables under multiple structural assumptions considered in the literature, and accordingly develop both exact and continuous score-based methods that can asymptotically identify the true Markov equivalence class. This offers a unified view of several existing constraint-based methods with different structural assumptions.