Digital Twins (DTs) are computational models that reflect real-world dynamical systems with high fidelity, playing a crucial role in simulations, understanding system dynamics, and supporting decision-making across diverse domains. However, existing DT approaches often struggle to generalize to unseen conditions in data-scarce settings, a crucial requirement for such models. To further understand these challenges, our work begins by outlining the essential desiderata for an effective DT. Subsequently, we introduce hybrid digital twins (HDTwins), a novel approach that symbolically represents a system as a composition of components, specifying the dynamic functions using both mechanistic and neural components. This hybrid model simultaneously embeds (partial) domain knowledge symbolically and capitalizes on the expressiveness of neural networks, while its modular design facilitates enhanced evolvability. However, automatically learning HDTwins from data poses significant challenges due to the complex search space and the flexible integration of domain priors. To address this, we propose an evolutionary multi-agent algorithm (HDTwinGen), augmented by Large Language Models, to efficiently and autonomously optimize HDTwins. Our empirical results reveal that HDTwinGen produces generalizable, sample-efficient, and evolvable models, significantly advancing DTs' efficacy in real-world applications.