It is easy to see that Eq. (14) gives the anticipated result for our model. Consider first the situation when ϵ > 6 t > 0 which according to Table I is a trivial insula- tor. In this case all 8 masses m$_{ℓ}$ are positive. Since four of the products v ℓ $_{x}$v ℓ $_{y}$v ℓ z are positive and four are nega- tive Eq. (14) gives θ = 0, as expected. Now suppose we tune ϵ so that its value drops below 6 t . This reverses the sign of a single Dirac mass at ℓ = (0 , 0 , 0). The corre- sponding θ ℓ also reverses sign and we obtain θ = π , as expected for a STI. For all other cases the values of θ are listed in Table I. These results confirm the one-to- one correspondence between the Z$_{2}$ invariant ν$_{0}$ and the axion parameter θ expressed by θ = πν$_{0}$ , as expected on very general grounds. 16 The system described by Hamiltonian (9) is inversion symmetric and we can thus employ the method devised in Ref. 9 to determine the topological class of its insulat- ing phases when the negative-energy bands are filled with electrons. This straightforward method requires comput- ing the eigenvalues of the parity operator at the 8 TRIM for the occupied bands. We find four distinct phases, de- pending on parameters ϵ and t , two of which are STI. The complete results are listed in Table I. According to general considerations 14,15 the STI phases should exhibit the axion term with θ = π . We evaluate θ for our model using the non-abelian Berry con- nection A αβ i = − i 〈 α k | ∂$_{i}$ | β k 〉 and the formula 14,15 We now consider a magnetic monopole in the interior of a STI. We model this situation by the Hamiltonian H defined by Eqs. (7) and (8) with a monopole posi- tioned at the center of a cubic unit cell. The magnetic field of the monopole couples to both the electron charge and the electron spin through the orbital and Zeeman couplings, respectively. The form of the orbital coupling is dictated by gauge invariance and is thus universal; in our lattice model it is implemented by the Peierls substi- tution, which attaches factors e iϑ$_{ij}$ to all hopping terms connecting sites i and j . Here ϑ$_{ij}$ = (2 π/ Φ$_{0}$) ∫ j i A · d l and A is the magnetic vector potential. The Zeeman coupling is of the form − gµ$_{B}$ B · S / ℏ where µ$_{B}$ = e ℏ / 2 m$_{e}$c is the Bohr magneton and S denotes the electron spin. For free electrons g is close to 2 but in solids the effective g can be substantially larger. The Zeeman coupling thus leads to an additional term in the Hamiltonian, where | β k 〉 is an eigenstate of H $_{k}$, the trace extends over occupied states and ∂$_{i}$ ≡ ∂/∂k$_{i}$ . The integral indicated in Eq. (13) is generally difficult to evaluate and numerical methods must be used to obtain θ for an arbitrary band structure. For a model in which the band touching is described by Dirac Hamiltonians, however, a simple ana- lytical evaluation of Eq. (13) is possible by noticing that in the limit of a small Dirac mass the entire contribution to the integral comes from the Dirac points. We show in the Appendix that each Dirac point contributes where B$_{j}$ is the magnetic field at site j of the lattice. This term is non-universal and its importance will depend on the ratio of gµ$_{B}$ | B | to the other relevant energy scales in the Hamiltonian set by λ , ϵ and t . We solve the Hamiltonian H = H$_{SO}$ + H$_{cd}$ + H$_{Z}$ in a cube containing L 3 sites by exact numerical diago- nalization. The monopole is placed inside the central unit cell (at the origin), so that the magnetic field of the monopole is B = (Φ$_{0}$ / 4 πr $^{2}$)ˆ r . We choose a gauge in which the system retains the four-fold rotational sym- metry around the z axis, 26 A = − Φ$_{0}$(1 + cos θ ) ∇ ϕ , with ( θ, ϕ ) the spherical angles. Exploiting this symmetry we are able to simulate system sizes up to L = 20, which requires diagonalizing a complex valued Hermitian ma- trix of size 1 $_{4}$(4 × 20 $^{3}$) = 8 , 000. In order to calculate the charge density at half filling we require knowledge of all the occupied eigenstates. to the total θ = ∑ ℓ θ ℓ mod 2 π . Although Eq. (14) has been derived for the case | m$_{ℓ}$ | ≪ | λ | we expect it to be more generally valid for all TRI Hamiltonians that can be deformed into the form of Eq. (9). This is because the value of θ in a TRI insulator is quantized and can only change when a band crossing closes the gap. As long as band crossings occur only at Γ$_{ℓ}$ points and are described by Dirac Hamiltonians, θ will be determined by Eq. (14), even when | m$_{ℓ}$ | is not small. We diagonalize the Hamiltonian, once with the mag- netic monopole and once without, obtaining charge den- sities ρ$_{1}$ and ρ$_{0}$ , respectively. The monopole-induced H$_{Z}$ = − gµ$_{B}$ 1 2 ∑ j B$_{j}$ · $^{(}$Ψ † $_{j}$σ Ψ$_{j}$ $^{)}$, (15) θ = 1 4 π ∫ BZ d $^{3}$kϵ $^{ijk}$Tr$^{[}$ A $_{i}$∂$_{j}$ A k + 2 i 3 A i A j A k ] , (13) θ ℓ = − π 2 sgn $^{(}$v ℓ $_{x}$v ℓ $_{y}$v ℓ $_{z}$m$_{ℓ}$ ) (14) E ℓ k = ± $^{√}$4 λ $^{2}$k2 + m 2 $_{ℓ}$. (12) v ℓ µ ≡ − 2 λ ( − 1) n$_{µ}$ are the Cartesian components of the Dirac velocities at Γ$_{ℓ}$ and m$_{ℓ}$ ≡ m$_{Γ}$$_{ℓ}$ are the correspond- ing Dirac masses. The Hamiltonians H ℓ k show spectra IV. CHARGE BOUND TO A MONOPOLE 3 TABLE I: Z$_{2}$ indices ( ν$_{0}$ ; ν$_{1}$ν$_{2}$ν$_{3}$ ) calculated according to Ref. 9 for our model insulator at half filling and the corresponding values of the axion parameter θ . ν$_{0}$ is the important ‘strong’ invariant while ν$_{i}$$_{=1}$$_{,}$$_{2}$$_{,}$$_{3}$ are the so called ‘weak’ invariants 7,8 which do not play a role in the present study but we list them here for completeness. It is assumed that λ, t > 0. WTI denotes a ‘weak’ topological insulator. Parameters Z$_{2}$ class Insulator type Axion θ | ϵ | > 6 t (0;000) trivial 0 − 6 t < ϵ < − 2 t (1;111) STI π − 2 t < ϵ < 2 t (0;111) WTI 0 2 t < ϵ < 6 t (1;000) STI π