Interactions of large numbers of spiking neurons give rise to complex neural dynamics with fluctuations occurring at multiple scales. Understanding the dynamical mechanisms underlying such complex neural dynamics is a long-standing topic of interest in neuroscience, statistical physics and nonlinear dynamics. Conventionally, fluctuating neural dynamics are formulated as balanced, uncorrelated excitatory and inhibitory inputs with Gaussian properties. However, heterogeneous, non-Gaussian properties have been widely observed in both neural connections and neural dynamics. Here, based on balanced neural networks with heterogeneous, non-Gaussian features, our analysis reveals that in the limit of large network size, synaptic inputs possess power-law fluctuations, leading to a remarkable relation of complex neural dynamics to the fractional diffusion formalisms of non-equilibrium physical systems. By uniquely accounting for the leapovers caused by the fluctuations of spiking activity, we further develop a fractional Fokker-Planck equation with absorbing boundary conditions. This body of formalisms represents a novel fractional diffusion theory of heterogeneous neural networks and results in an exact description of the network activity states. This theory is further implemented in a biologically plausible, balanced neural network and identifies a novel type of network state with rich, nonlinear response properties, providing a unified account of a variety of experimental findings on neural dynamics at the individual neuron and the network levels, including fluctuations of membrane potentials and population firing rates. We illustrate that this novel state endows neural networks with a fundamental computational advantage; that is, the neural response is maximised as a function of structural connectivity. Our theory and its network implementations provide a framework for investigating complex neural dynamics emerging from large networks of spiking neurons and their functional roles in neural processing.