We study optimal synchronization in networks of heterogeneous phase oscillators. Our main result is the derivation of a synchrony alignment function that encodes the interplay between network structure and oscillators' frequencies and can be readily optimized. We highlight its utility in two general problems: constrained frequency allocation and network design. In general, we find that synchronization is promoted by strong alignments between frequencies and the dominant Laplacian eigenvectors, as well as a matching between the heterogeneity of frequencies and network structure. PACS numbers: 05.45.Xt, 89.75.Hc A central goal of complexity theory is to understand the emergence of collective behavior in large ensembles of interacting dynamical systems. Synchronization of networkcoupled oscillators has served as a paradigm for understanding emergence [1][2][3][4], where examples arise in nature (e.g., flashing of fireflies [5] and cardiac pacemaker cells [6]), engineering (e.g., power grid [7] and bridge oscillations [8]), and at their intersection (e.g., synthetic cell engineering [9]). We consider the dynamics of N network-coupled phase oscillators θ i for i = 1, . . . , N , whose evolution is governed bẏHere ω i is the natural frequency of oscillator i, K > 0 is the coupling strength, [A ij ] is a symmetric network adjacency matrix, and H is a 2π-periodic coupling function [10]. We treat H(θ) with full generality so long as H ′ (0) > 0. The choices H(θ) = sin(θ) and H(θ) = sin(θ − α) with the phase-lag parameter α ∈ (−π/2, π/2) yield the classical Kuramoto [10] and Sakaguchi-Kuramoto models [11].Considerable research has shown that the underlying structure of a network plays a crucial role in determining synchronization [12][13][14][15][16][17][18][19][20][21], yet the precise relationship between the dynamical and structural properties of a network and its synchronization remains not fully understood. One unanswered question is, given an objective measure of synchronization, how can synchronization be optimized? One application lies in synchronizing the power grid [22], where sources and loads can be modeled as oscillators with different frequencies. To this end, we ask: what structural and/or dynamical properties should be present to optimize synchronization?We measure the degree of synchronization of an ensemble of oscillators using the Kuramoto order parameterHere re iψ denotes the phases' centroid on the complex unit circle, with the magnitude r ranging from 0 (incoherence) to 1 (perfect synchronization) [10]. In general, the question of optimization (maximizing r) is challenging due to the fact that the macroscopic dynamics depend on both the natural frequencies and the network structure. To quantify the interplay between node dynamics and network structure, we derive directly from Eqs. (1) and (2) a synchrony alignment function which is an objective measure of synchronization and can be used to systematically optimize a network's synchronization. We highlight this result by addressing two classes of op...