  We study the eigenvalue problem for a one-dimensional Dirac operator with a
spatially varying ``mass'' term. It is well-known that when the mass function
has the form of a kink, or \emph{domain wall}, transitioning between strictly
positive and strictly negative asymptotic mass, $\pm\kappa_\infty$, at
$\pm\infty$, the Dirac operator has a simple eigenvalue of zero energy
(geometric multiplicity equal to one) within a gap in the continuous spectrum,
with corresponding \emph{zero mode}, an exponentially localized eigenfunction.
We prove that when the mass function has the form of \emph{two} domain walls
separated by a sufficiently large distance $2 \delta$, the Dirac operator has
two real simple eigenvalues of opposite sign and of order $e^{- 2
|\kappa_\infty| \delta}$. The associated eigenfunctions are, up to $L^2$ error
of order $e^{- 2 |\kappa_\infty| \delta}$, linear combinations of shifted
copies of the single domain wall zero mode. For the case of three domain walls,
there are two non-zero simple eigenvalues as above and a simple eigenvalue at
energy zero. Our methods are based on a Lyapunov-Schmidt reduction strategy and
we outline their natural extension to the case of $n$ domain walls for which
the minimal distance between domain walls is sufficiently large.
  The class of Dirac operators we consider controls the bifurcation of
topologically protected ``edge states'' from Dirac points (linear band
crossings) for classes of Schr\"odinger operators with domain-wall modulated
periodic potentials in one and two space dimensions. The present results may be
used to construct a rich class of defect modes in periodic structures modulated
by multiple domain walls.
