  At low energies, the strong interaction is governed by the Goldstone bosons
associated with the spontaneous chiral symmetry breaking, which can be
systematically described by chiral perturbation theory. In this paper, we apply
this theory to study the $\theta$-vacuum energy density and hence the QCD axion
potential up to next-to-leading order with $N$ non-degenerate quark masses. By
setting $N=3$, we then derive the axion mass, self-coupling, topological
susceptibility and the normalized fourth cumulant both analytically and
numerically, taking the strong isospin breaking effects into account. In
addition, the model-independent part of the axion-photon coupling, which is
important for axion search experiments, is also extracted from the chiral
Lagrangian supplemented with the anomalous terms up to $\mathcal{O}(p^6)$.
