Causal structure learning (CSL) plays a pivotal role in causality and is often formulated as an optimization problem within score-and-search methods. Under the assumption of an infinite dataset and a predefined distribution, several well-established and consistent score functions have been shown to be both optimal and reliable for identifying ground-truth causal graphs. However, in practice, these idealized assumptions are often infeasible, which can result in CSL algorithms learning suboptimal structures. In this paper, we introduce L-SFE, a framework designed to automatically discover effective score functions by exploring the "score function space". L-SFE addresses this task from a bi-level optimization perspective. First, it leverages a Large Language Model (LLM) to interpret the characteristics of score functions and generate the corresponding code implementations. Next, L-SFE employs evolutionary algorithms along with carefully designed operators, to search for solutions with higher fitness. Additionally, we take the BIC as example and prove the consistency of the generated score functions. Experimental evaluations, conducted on discrete, continuous, and real datasets, demonstrate the high stability, generality and effectiveness of L-SFE.