We show the existence of all the 36 eightfold way mesons and determine their masses and dispersion curves exactly, from dynamical first principles such as directly from the quark-fluon dynamics. We also give a proof of confinement below the two-meson energy threshold. For this purpose, we consider an imaginary time functional integral representation of a 3+1 dimensional lattice QCD model with Wilson action, SU(3)f global and SU(3)c local symmetries. We work in the strong coupling regime, such that the hopping parameter κ>0 is small and much larger than the plaquette coupling β>1/g02⩾0 (β⪡κ⪡1). In the quantum mechanical physical Hilbert space H, a Feynman-Kac type representation for the two-meson correlation and its spectral representation are used to establish an exact rigorous connection between the complex momentum singularities of the two-meson truncated correlation and the energy-momentum spectrum of the model. The total spin operator J and its z-component Jz are defined by using π∕2 rotations about the spatial coordinate axes, and agree with the infinitesimal generators of the continuum for improper zero-momentum meson states. The mesons admit a labelling in terms of the quantum numbers of total isospin I, the third component I3 of total isospin, the z-component Jz of total spin and quadratic Casimir C2 for SU(3)f. With this labelling, the mesons can be organized into two sets of states, distinguished by the total spin J. These two sets are identified with the SU(3)f nonet of pseudo-scalar mesons (J=0) and the three nonets of vector mesons (J=1,Jz=±1,0). Within each nonet a further decomposition can be made using C2 to obtain the singlet state (C2=0) and the eight members of the octet (C2=3). By casting the problem of determination of the meson masses and dispersion curves into the framework of the the anaytic implicit function theorem, all the masses m(κ,β) are found exactly and are given by convergent expansions in the parameters κ and β. The masses are all of the form m(κ,β=0)≡m(κ)=−2lnκ−3κ2/2+κ4r(κ) with r(0)≠0 and r(κ) real analytic; for β>0,m(κ,β)+2lnκ is jointly analytic in κ and β. The masses of the vector mesons are independent of Jz and are all equal within each octet. All isospin singlet masses are also equal for the vector mesons. For each nonet and β=0, up to and including O(κ4), the masses of the octet and the singlet are found to be equal. But there is a pseudoscalar-vector meson mass splitting given by 2κ4+O(κ6) and the splitting persists for β>0. For β=0, the dispersion curves are all of the form w(p⃗)=−2lnκ−3κ2∕2+(14)κ2∑j=132(1−cospj)+κ4r(κ,p⃗), with ∣r(κ,p⃗)∣⩽const. For the pseudoscalar mesons, r(κ,p⃗) is jointly analytic in κ and pj, for ∣κ∣ and ∣Impj∣ small. We use some machinery from constructive field theory, such as the decoupling of hyperplane method, in order to reveal the gauge-invariant eightfold way meson states and a correlation subtraction method to extend our spectral results to all He, the subspace of H generated by vectors with an even number of Grassmann variables, up to near the two-meson energy threshold of ≈−4lnκ. Combining this result with a previously similar result for the baryon sector of the eightfold way, we show that the only spectrum in all H≡He⊕Ho (Ho being the odd subspace) below ≈−4lnκ is given by the eightfold way mesons and baryons. Hence, we prove confinement up to near this energy threshold.