Model Realization and Numerical Studies of a Three-Dimensional Bosonic Topological Insulator and Symmetry-Enriched Topological Phases

We study a topological phase of interacting bosons in (3+1) dimensions which is protected by charge conservation and time-reversal symmetry. We present an explicit lattice model which realizes this phase and which can be studied in sign-free Monte Carlo simulations. The idea behind our model is to bind bosons to topological defects called hedgehogs. We determine the phase diagram of the model and identify a phase where such bound states are proliferated. In this phase we observe a Witten effect in the bulk whereby an external monopole binds half of the elementary boson charge, which confirms that it is a bosonic topological insulator. We also study the boundary between the topological insulator and a trivial insulator. We find a surface phase diagram which includes exotic superfluids, a topologically ordered phase, and a phase with a Hall effect quantized to one-half of the value possible in a purely two-dimensional system. We also present models that realize symmetry-enriched topologically-ordered phases by binding multiple hedgehogs to each boson; these phases show charge fractionalization and intrinsic topological order as well as a fractional Witten effect.


I. INTRODUCTION
The study of topological phases of matter has been a major component of condensed matter research in recent decades. Among the many phases studied, the topological insulator (TI) is one of the most prominent. 1,2 The TI is a three-dimensional phase of free fermions. Though it is insulating in the bulk, its topological behavior can be deduced from its unusual surface properties, in particular the odd number of Dirac cones it has on its surface. The topological insulator is an example of a symmetry-protected topological phase (SPT). Like all SPT's, it has short-ranged entanglement, which implies that it has only conventional excitations in the bulk and a unique ground state on any closed manifold. This is in contrast to intrinsically topologically ordered states like the fractional quantum Hall states. The relevant symmetries for the topological insulator are charge conservation and time-reversal, and if either of these symmetries is violated, the phase loses its topological properties.
One obvious extension of research into topological insulators is to consider the effects of interactions on their properties. This is however a difficult task. Many of the methods used to study TI's involve the properties of their band structure, and these methods obviously do not apply if interactions are strong. As an introduction to this difficult problem, one can try to study an analog of the topological insulator, constructed of interacting bosons instead of fermions. In bosonic systems we know that the non-interacting case would be a condensate, so we can be sure that the topological behavior is due to the interactions. In addition, certain theoretical techniques, like the Monte Carlo studies employed in this paper, work mainly for bosonic systems.
The study of topological phases of interacting bosons is relatively recent, but much progress has been made. [3][4][5][6][7][8][9][10][11][12][13] Chen, Liu, Gu, and Wen 3 have used group cohomology theory to determine which symmetries and dimensions can lead to nontrivial topological phases. However, this approach tells us little about the properties of these phases, which must be determined through other methods. One well-studied case is that of SPT phases with U (1) symmetry in two dimensions. Lu and Vishwanath 4 proposed a phenomenological Chern-Simons field theory to describe both the bulk and edges of these states and showed that such states have a Hall effect quantized to an even integer (in units of e 2 /h) and possess gapless counter-propagating edge modes. Therefore such states are called "bosonic integer quantum Hall phases". Senthil and Levin 5 proposed how to realize such phases in a quantum Hall-like setting by starting with two species of bosons in a magnetic field and using mutual flux attachment. One drawback of the flux attachment technique is that it is difficult to relate it precisely to a microscopic model. To address this, in our work 7 we provided an alternative exactly solvable model that replaces flux attachment by a (precisely formulated) dynamical binding of bosons to vortices. We studied this model in Monte Carlo and confirmed that it realizes the integer Quantum Hall phase of bosons, including observation of the gapless edge modes. We also showed that by binding multiple vortices to bosons we can realize Symmetry-Enriched Topological phases (SET) 14,15 with long-ranged entanglement.
In this paper we focus on understanding interacting bosonic topological phases which are analogs of the electronic topological insulator in three dimensions. Motivated by the formal cohomology results, 3 Vishwanath and Senthil in an inspirational paper 6 found effective field theories which can describe both the bulk and the surface of a three-dimensional "bosonic topological insulator" with charge conservation and time-reversal symmetry. They found exotic behavior on the surface which can be used to assert the topological behavior in the bulk.
In particular, Vishwanath and Senthil found three kinds of exotic phases on the surface of the bosonic TI. These surface phases cannot exist in a purely two-dimensional system, and their existence on the surface of a three-dimensional system shows that the system is topological. The first kind of phase is a superfluid, which spontaneously breaks charge conservation symmetry. There are actually several different types of superfluids which can exist on the surface. The gapped vortex excitations in these superfluids have properties which cannot exist in a purely two-dimensional system. Phase transitions between the different superfluids are predicted to be deconfined critical points. Another kind of surface phase appears when we break time-reversal symmetry on the surface. This phase has a Hall conductivity quantized to one-half of the elementary value possible in a purely two-dimensional bosonic system. Since the Hall conductivity in the bosonic integer quantum Hall effect is quantized to even integers, this surface phase is expected to have an odd integer Hall conductivity. Finally, it is possible to have a surface phase which breaks no symmetries but has a symmetry-enriched intrinsic topological order of a kind impossible in a purely two-dimensional system with these symmetries.
In both the two-and three-dimensional cases, the topological behavior can be thought of as coming from the binding of bosons to point topological defects, and the condensation of such objects. In our construction in two dimensions the topological defects were vortices, and the binding is essentially an exact dynamical realization of the flux attachment. In the three-dimensional systems that we will study here, the topological defects are hedgehogs which we introduce by adding an additional SO(3) "spin" symmetry, and we construct a bosonic TI by binding hedgehogs of this symmetry to the bosons. 6 In an independent work, Metlitski and Fisher also produced a construction which explicitly binds hedgehog-like "monopole" objects to bosons without enlarging the continuous symmetry and showed that it gives the three-dimensional bosonic TI, 16 while the present setting is a bit simpler to analyze and is amenable to Monte Carlo studies.
In this work, we construct explicit models which realize the interacting bosonic analog of the topological insulator. The models have both charge conservation U (1) symmetry and a time reversal Z T 2 symmetry to be discussed in the main text. We present two different models which realize this physics. In the first model, described in Section II, the spins are represented by SO(3) degrees of freedom with Heisenberg interactions. We introduce a term in our action which energetically binds hedgehogs to bosons, and we show that this term can lead to a phase (which we call the "binding phase") where these bound states are proliferated.
In order to show that the binding phase is indeed the bosonic topological insulator, we will attempt to find the phases predicted by Vishwanath and Senthil on its surface. Initially we find only superfluids on the surface. In our Monte Carlo simulations, we do not have simple access to properties of the gapped vortex excitations in these surface phases. Therefore we cannot determine whether the superfluids we observe are the exotic superfluids predicted by Vishwanath and Senthil, or more conventional superfluids. We find that the surface superfluids in our model are connected by a direct transition. If this transition were second-order, then it could be the predicted deconfined critical point. However, we cannot access large enough system sizes to determine the order of the observed phase transition.
We can try to find other exotic surface phases by explicitly breaking the Z T 2 symmetry of our model on the surface. If the bulk of the system is in a topological phase, this will lead to a Hall conductivity quantized to an odd integer. In order to measure this Hall conductivity in our system, we need to introduce external gauge fields. In the model in Section II, this is complicated by our definition of the hedgehog number, which is a discontinuous function of the Heisenberg spins. In this case, we do not know how to properly couple to the external gauge fields, which prevents us from measuring the surface Hall conductivity.
To remedy this, in Section III we introduce a second model which binds hedgehogs to bosons. In this model, we represent the spins with an easy-plane CP 1 model. We can identify hedgehogs with the monopoles of the internal compact gauge field of this CP 1 model. This formulation will allow us to make several measurements which demonstrate that our phase is a bosonic topological insulator.
One measurement that we can make is called the Witten effect. [16][17][18][19] This is the binding of one-half of a boson charge to external monopoles in the spin sector. We introduce such external monopoles into the bulk of our system and find that each binds precisely half of a boson charge.
We then determine the surface phases of our model. We again find superfluids connected by a direct transition. We can also break the Z T 2 symmetry on the surface. In the CP 1 model, we can measure the surface Hall conductivity and we find it to be quantized to odd integers as predicted. We reiterate that this Hall conductivity cannot be observed in a purely two-dimensional model, and so we must be measuring the Hall conductivity on the surface of a topological phase. Finally, we also find a surface phase which breaks none of the symmetries of the model. We suspect that this phase has symmetry-enriched intrinsic topological order of the kind predicted by Vishwanath and Senthil, 6 but we do not know how to test this using Monte Carlo.
In our studies of the (2+1)-dimensional "bosonic quantum Hall" phases, 7 we found that we can obtain symmetryenriched phases with intrinsic topological order by binding multiple topological defects (vortices) to bosons. In Section IV we consider such a binding in (3+1) dimensions. We find that binding of multiple hedgehogs to bosons leads to a bulk phase with intrinsic topological order, thus opening the study of bosonic SET phases in (3+1) dimensions. To help understand the properties of this SET phase, in the Appendix we consider a more abstract model where monopoles of a compact quantum electrodynamics (CQED) are bound to bosons; we analyze the resulting SPT-like and SET-like phases in such a CQED×boson theory, and also discuss how things change upon including matter fields coupled to the CQED, which is the situation in the CP 1 representation of the spins.

II. REALIZING THE TOPOLOGICAL INSULATOR BY BINDING BOSONS TO HEDGEHOGS OF SO(3) SPINS
We first study the binding between bosons and topological defects by using hedgehogs of a Heisenberg model as our topological defects. We demonstrate the existence of a phase which can be loosely viewed as a condensate of bound states of bosons and hedgehogs, and explore its properties. This model provides an intuitive introduction to the physics of such binding. Additional properties of the resulting SPT phase will be considered in Section III.

A. Model and its Bulk Phase Diagram
We study the following action, in (3+1)D Euclidean spacetime: S spin is an action which controls fluctuations in the SO (3) spins. The second term provides the binding interaction between bosons and hedgehogs. The bosons are represented by integer-valued conserved currents, J µ (r), defined on the links of a four-dimensional cubic lattice, where r is a site label on the lattice and µ ∈ (x, y, z, τ ) is a direction. These currents represent the world-lines of the bosons in the (3+1)dimensional space-time. The Q µ (r) variables represent the hedgehogs, which are also integer-valued conserved currents and will be defined shortly. When the real number λ is large, this term will bind bosons and hedgehogs together. We work with periodic boundary conditions and require no net boson charge and no net boson spatial currents, so that the J µ spacetime currents have zero total winding number; such conditions are automatically satisfied for the hedgehog currents Q µ defined below. Imposing such conditions on the boson currents is just a convenient choice, which does not affect the bulk physics, but allows a precise change of variables involving both J and Q currents [see Eq. (10) below]. In this section we represent the spin degrees of freedom by SO(3) unit vectors. For S spin , we take the following action, which describes a Heisenberg model: Here n(R) = (n a , n b , n c )(R) are three-component unit vectors which represent the spins. They reside on a different lattice from the J µ (r) variables above. This lattice has its sites labelled by R and located at the centers of the (hyper)cubes of the r-lattice in Eq. (1), i.e., the R and r lattices are dual to each other. From these n(R), we can define the hedgehog currents Q µ (r) using the prescription in the literature 20-23 generalized to four dimensions. We summarize this prescription here. We first define variables α µ (R), which reside on links connecting the spins n i ≡ n(R) and n j ≡ n(R +μ): Here N 0 is a reference vector which we can choose arbitrarily. We then define placket "fluxes" ω µν ∈ (−π, π] as follows: Here ∇ µ α ν (R) ≡ α ν (R +μ) − α ν (R). One can show that changing the reference vector N 0 corresponds to a gauge transformation of the α µ (R) variables so that ω µν (R) are independent of the reference vector. Finally, we can define the hedgehog current: with implied summation over repeated indices. Consider for example The right-hand-side is defined on a cube [R, R+x, R+ŷ, R+ z], and can also be associated with a point R + (x +ŷ +ẑ)/2 in the center of this cube. The ω variables are fluxes passing through the plaquettes of this cube, and the net flux out of the cube is guaranteed to be integer multiple of 2π. We then define hedgehog number to be the net outgoing flux divided by 2π. When all four dimensions are considered, the center of this cube can be equivalently associated with a link in the τ direction from a dual lattice site r = R + (x +ŷ +ẑ −τ )/2, and the hedgehog number becomes the τ -component, Q τ (r), of the hedgehog four-current.
We have studied the above action in Monte Carlo simulations. In addition to simple independent updates of n and J µ , when λ is large one needs to try updates which change both J µ and Q µ in such a way that the second term in Eq. (1) is unchanged. To do this we choose a spin and an amount to update it, check to see if any Q µ variables will change, and include matching changes in the J µ variables as part of the attempted Monte Carlo move.
In our simulations, we monitor the "internal energy per site," ǫ = S/Vol, where Vol ≡ L 4 is the volume of the system, which we take to have linear size L in all four directions. From this, we can determine the specific heat per site: We can locate phase transitions in our model by looking for singularities in the specific heat. We also monitor the magnetization per spin: When the spins are disordered the magnetization is proportional to 1/ √ Vol, while in the ordered phase the magnetization remains non-zero in the thermodynamic limit. Therefore we can use measurements of the magnetization at different sizes to determine if the spins are ordered.
To study the behavior of the boson currents, we monitor current-current correlators, defined as: where k is a wave vector, µ is a fixed direction, and In the space-time isotropic system, ρ J (k) is independent of the direction µ, and when we show numerical data we average over all directions to improve statistics. In an ensemble which would allow non-zero total winding number, ρ J (0) would be the familiar superfluid stiffness. In our model J(k = 0) = 0, so this measurment is not available in our simulations. Instead, we evaluate the correlators at the smallest non-zero k. For example, if µ is in the x direction, we can take k min = (0, 2π L , 0, 0), (0, 0, 2π L , 0), and (0, 0, 0, 2π L ) and average over these; we exclude ( 2π L , 0, 0, 0) because in our ensemble the net winding of the J x current is zero, so the J x (k) evaluated at this wavevector is also zero. We also monitor current-current correlators of the hedgehog currents, ρ Q (k), which are defined in the same way as for the boson currents.
In the phase where the J µ are gapped, only small loops contribute to the current-current correlators and ρ J (k min ) ∼ k 2 min ∼ 1/L 2 , while when the J µ proliferate ρ J is independent of the system size. Therefore we can use finite-size scaling of this quantity to determine the locations of phase transitions. For the hedgehog currents, ρ Q (k min ) ∼ k 2 min in all phases, so we cannot use finite-size scaling of this quantity to find phase transitions; this is similar to the properties of vortex currents in a (2+1)-dimensional XY model and originates from effective long-range interactions of these topological defects.
The bulk phase diagram of this model is shown in the inset in Fig. 1 and can be understood essentially analytically. Indeed, in the case where the entire system is in the same phase (i.e., the boson-hedgehog binding is applied everywhere), we can change to new variables in the partition sum which satisfy the same conditions as the original boson currents. Expressed in these variables, the first and second terms in Eq. (1) decouple, and we can study them independently. At small β, the n spins are disordered, which also implies that the hedgehog currents are proliferated. As β is increased, the spins order. This means that there is a large energy cost for hedgehog currents to exist. We say that the hedgehog currents are gapped, and only small loops of them are present. We can determine the location of the spin-ordering phase transition by finding singularities in the heat capacity, or by performing finite-size scaling on the magnetization as described above. The value found by our numerics agrees with the literature on the 4D classical Heisenberg model. 24 At small λ, theJ µ variables are proliferated. In our original variables, this means that the boson currents are effectively independent of the hedgehog currents and are condensed. At large λ, theJ µ variables are gapped. The boson currents are bound to the hedgehog currents. The result of this can be seen in Fig. 1, which shows the current-current correlators for the boson and hedgehog currents. For small λ, there is no relation between the bosons (red solid line) and the hedgehogs (blue dashed line). As λ is increased, the current-current correlators become essentially identical as the bosons and hedgehogs are bound together. We can determine the location of the phase transition in λ by studying singularities in the heat capacity, or by peforming finite-size scaling on ρJ as described above.  1) and (2) with bosons and Heisenberg spins. The phase diagram is mathematically equivalent to a system of decoupled currents and spins, hence the straight line boundaries. At λ = 0, the system has a paramagnetferromagnet transition as β is increased, while the bosons are superfluid throughout. As λ is increased, the boson currents bind to the hedgehog currents. The loop pictures in the phases show a "snapshot" of the phase. Red solid loops mean that boson currents are proliferated in the phase, while blue dashed loops indicate proliferated hedgehog currents. The phase of interest is the "binding" phase in the upper left corner where bosons are bound to hedgehogs. The main figure shows the current-current correlations of the bosons and hedgehogs as λ is increased while β = 0, for a system of linear dimension L = 6. We see that the correlators become essentially equal as the system enters the upper left phase, indicating that bosons have bound to hedgehogs.
We can now summarize the phase diagram shown in the inset in Fig. 1. At small λ and β both the boson and hedgehog currents are proliferated and are not bound, while at large λ and β all currents are gapped. At small λ and large β the boson currents are condensed but the hedgehog currents are gapped. Finally, when λ is large and β is small the system is in the "binding" phase with "proliferated" bound states, which we will argue is a (3+1)D topological phase protected by the appropriate symmetries.
It is also helpful to think about the "easy-plane" regime for the spin variables, in which the spins, n = (n a , n b , n c ), are roughly in the ab-plane, with only small c components; we will denote the corresponding global symmetry of spin rotations in the ab-plane as U (1) spin . In the easy-plane case we can define "vortices" of the XY spins (n a , n b ) (i.e., phase windings of the complex order parameter ∼ n a + in b ). The vortices are defined on the plackets of the cubic lattice. Therefore in the (3+1)D space-time they are represented as twodimensional "world-sheets." When dealing with vortices, one can gain intuition by thinking in terms of only the three spatial dimensions of our (3+1) dimensional space-time. In this picture the bosons and hedgehogs are represented by point particles, while the vortices are represented by lines. The next two paragraphs can be most easily understood by thinking in this picture.
Though ordinary XY spins are not defined at the core of a vortex, our spins have a c-component which can point either up or down at the vortex core. We can define two species of vortices, which we call the ↑ and ↓ species, depending on whether n c is positive or negative at the core. This description is useful since an ↑ vortex ending and continuing as a ↓ vortex is a hedgehog. Therefore our system can be thought of as a system of vortex lines having two different species. These vortex lines can change species, and the locations where this happens are hedgehogs. 25 This is illustrated in Fig. 2.
Hedgehog number can be either positive or negative, and this is determined by the properties of the vortex line the hedgehog is attached to. If, when looking along the line from the ↑ core to the ↓ core the vorticity is clockwise (counterclockwise), we define a positive (negative) hedgehog. This is pictured in Fig. 2.

Importance of Discrete Symmetry
The action in Eq. (1) has a U (1) symmetry which comes from the conservation of the J µ currents; this is boson charge conservation symmetry. It also has an SO(3) symmetry from the spins. Both of these symmetries participate in the protection of the topological phase, in the sense that if they are broken it is possible to continuously connect the topological phase to a trivial phase. In addition, the action has a Z 2 symmetry, which is obtained by reflecting the n spins in a plane in the spin space. To see how this affects the hedgehog current, we can examine Eq. (3), taking the reference vector N 0 to be in the plane of reflection. We see that reflecting the spins changes the sign of the imaginary part of e iαµ , and therefore the hedgehog current changes sign under such a reflection. For our entire action to be invariant, we therefore need to combine such reflections with an operation which changes the sign of the boson currents. For concreteness we will consider the Z 2 symmetry corresponding to reflections of the n variables in the ab plane of the spin space. In this case the Z 2 symmetry can be summarized as: Note that it is also possible to reflect the n spins around a different plane, but this is not a distinct symmetry since it is simply the product of the above Z 2 symmetry and an element of SO(3). By analogy with the electronic topological insulator, we would like the Z 2 symmetry described above to be a "timereversal" symmetry, i.e. it should be anti-unitary. The symmetry in Eq. (11) can be either a unitary or anti-unitary symmetry. Note that Eq. (1) is a real, (3+1)-dimensional action which is assumed to arise from the Trotter decomposition of the imaginary time propagator (i.e., Euclidean path integral) of a three-dimensional quantum Hamiltonian. The symmetry operations in Eq. (11) can be derived from the action of a symmetry operation on the quantum Hamiltonian. Therefore asking whether the symmetry in Eq. (11) is anti-unitary is the same as asking whether the symmetry of the quantum Hamiltonian which generates Eq. (11) is anti-unitary. This is a difficult question for us to answer as we do not know the quantum Hamiltonian which has Eq. (1) as its Euclidean path integral. Nevertheless, we take the perspective where we can check whether the original Hamiltonian has time reversal by complex-conjugating the action combined with the appropriate variable transformations, and in this way we can view the above symmetry of the action also as time reversal.
In the easy plane case, our system has U (1) boson ×U (1) spin in addition to this discrete symmetry. We can think of the U (1) spin as also coming from a boson, and n a + in b gives the phase degree of freedom of this boson. We can imagine allowing tunnelling between the two U (1) symmetries. If we do this, then the discrete symmetry should act the same way on the two species, and we can see that the above symmetry changes the number of the bosons but not their phase. In a system consisting of one species of bosons, the symmetry which acts in this way is anti-unitary. From now on the symmetry described above will also be treated as anti-unitary and denoted by Z T 2 . The action in Eq. (1) is invariant under several Z 2 symmetries, but for our purposes we will consider only the Z T 2 symmetry described above as important, as it protects the topological behavior. To see this we can break this symmetry and argue that the topological phase is destroyed. We break the Z T 2 symmetry by introducing a Zeeman field into our action: Here h is the strength of the Zeeman field, which points in the c-direction. Note that the Zeeman field does not break U (1) spin or U (1) boson . In our picture of two species of vortices (Fig. 2), the Zeeman field forbids one of the species. Since the vortex lines cannot change species, hedgehogs are forbidden and the binding phase is destroyed. This can be made more precise if we replace the binding term in Eq. (1) with the following term: where η is a real number. If we choose the parameters β and λ so that the system is initially in the binding phase, the introduction of η allows us to tune the system between the binding phase (η = 1) and the trivial insulator (η = 0). Without a Zeeman field, the system undergoes a phase transition as η is changed between 1 and 0. When a Zeeman field is applied, the hedgehogs are effectively forbidden, and so we can tune parameters h and η without going through a phase transition. Indeed, when η = 1, the change of variables in Eq. (10) leads to decoupledJ currents and spins, so there is no phase transition when making h arbitrarily large to align all spins. We can then tune η to zero and finally reduce h back to zero, all without undergoing a phase transition. Thus, in the presence of the Zeeman field, the phase with η = 1 is not distinct from the trivial insulator.

Binding of Multiple Bosons to a Hedgehog
The above methods also allow us to answer the question of what happens to the system if η in Eq. (13) is an integer larger than 1. In a U (1) × U (1) system in two dimensions, we studied the binding of multiple bosons to vortices (realized by taking integer η) and found that each number of bound bosons led to a different symmetry-protected topological phase. 7 There were therefore as many SPTs as there are integers, in agreement with the cohomology classification. 3,4 In the present three-dimensional case the classification of Chen et al. 3 for U (1) and Z T 2 symmetry predicts the existence of only a finite number of such symmetry protected topological phases, implying that not every value of η would lead to a distinct phase. Indeed, we find that all systems with η an even integer are topologically trivial, while when η is odd we have the same topological phase as η = 1.
We can justify this claim by showing that η = 2 can be continuously connected to η = 0 which is a trivial insulator. This argument can then be extended to show that any two systems where η differs by 2 are in the same phase. Our argument is inspired by Ref. 26, except that here we are working with a microscopic model rather than a topological field theory. We start by considering two copies of our action, each with its own bosons and Heisenberg spins and with η = 1: where the superscripts indicate which copy a variable is from. We now couple the two copies by adding the following terms: When A is large and positive, the first term above locks spins of the different copies together, n (1) ≈ n (2) . The hedgehog variables therefore take on the same values, and either can be viewed as the hedgehog number of the whole spin system in this case: Q ≈ Q (1) ≈ Q (2) . On the other hand, when A is large and negative the spins of different types are locked in opposite directions, n (1) ≈ − n (2) , and Q (1) ≈ −Q (2) . The Φ variables are 2π-periodic phases that can be thought of as conjugates to the J µ variables. More precisely, in our path integral we sum over only the configurations of J µ in which the currents are divergenceless. We can instead sum over all configurations of J µ , and include the following term in our path integral: which dynamically enforces the constraint that the boson currents be conserved. We introduce Φ variables for each copy of the boson currents, and the B term is tunnelling between the two copies. When B is large (of either sign) only the sum of J (1) and J (2) is conserved and can be identified as the current of the whole boson system (i.e., combining both copies). When A and B are large and positive, we can expand the terms on the second line of Eq. (14), and take J = J (1) +J (2) , Q = Q (1) = Q (2) . In this case we obtain coupling between J and Q which is effectively the same as in Eq. (13) with η = 2. On the other hand, when A is large and negative, but B is still large and positive, J is unchanged but its coupling to the hedgehogs vanishes as contributions from Q (1) and Q (2) cancel; this gives Eq. (13) with η = 0. We can continuously deform A = 0, B = 0 to A = ∞, B = ∞ without undergoing a phase transition. In addition, we can deform A = 0, B = 0 to A = −∞, B = ∞ without undergoing a phase transition. This implies that we can tune from η = 2 to η = 0 without undergoing a phase transition, and so both of these cases are in the trivial insulating phase.

B. Phase Diagram on the Boundary Between the Binding Phase and a Trivial Insulator
By analogy to the fermionic topological insulator, we expect that one way to investigate the topological nature of our phase is to study the physics of its surface. In particular, a number of interesting phases have been predicted on the surface of the bosonic TI, 6 and if we can identify these phases on the surface of our binding phase, it will be a powerful argument that the binding phase is a bosonic TI.
We introduce a surface between the binding phase and a trivial insulator by allowing η in Eq. (13) to vary spatially. We vary η in the z direction, so that: This leads to the binding phase in the region z L < z < z R , while the trivial phase occupies the rest of the space. Note that there are two surfaces of the binding phase: one at z L and one at z R . The above geometry in our Monte Carlo setup with periodic boundary conditions corresponds to a four-dimensional torus which is divided into two parts along the z-direction.
On the surface, we can measure all of the quantities which we measured in the bulk, but we now only sum over the sites near the surface (and when averaging over directions we only use thex,ŷ andτ directions).
The binding between hedgehogs and bosons in the bulk of the topological phase leads to exotic physics on the surface. In both the binding and trivial phases, hedgehog currents are proliferated, while boson currents are bound to the hedgehog currents in the binding phase and are absent in the trivial phase, as shown in Fig. 3. Consider what happens when a hedgehog loop tries to cross the boundary between the phases. Since the boson currents must form closed loops, the above conditions cannot both be satisfied without something interesting happening on the surface. For example, we could have unbound boson currents on the surface, or hedgehogs could be effectively forbidden from crossing the surface.
In order to study the surface physics and search for the exotic surface behavior predicted in Ref. 6, we first determine the surface phase diagram. To do this, we fix the values of β and λ in the bulk and tune them only on one of the surfaces, by setting: Here we focused on the surface at z R , and µ denotes link orientation for either β or λ terms. We show the surface phase diagram in the inset of Fig. 4. All data was taken with β bulk = 0, λ bulk = 5.2, parameters which put the bulk deep into the binding phase. The surface phase diagram contains three distinct phases. When λ surf is small the bosons are in a superfluid phase, breaking their U (1) symmetry. This is the scenario pictured in Fig. 3, where hedgehogs can cross the surface, and these crossings are connected by boson currents. If β surf is also small the SO(3) symmetry is unbroken as the spins are disordered. As β surf increases at small λ surf , the SO(3) symmetry breaks, and so both symmetries are broken. Finally, at large λ surf bosons see a large energy cost, and so the U (1) symmetry is unbroken. This forbids hedgehogs from crossing FIG. 3. A snapshot of the system, when it is spatially divided into a region which is a trivial insulator and a region which is in the binding phase. In the trivial phase, boson currents are gapped and hedgehog currents are proliferated. In the binding phase both currents are gapped individually, but their bound states are proliferated. Large hedgehog current loops can exist in either region; however, when such a loop tries to cross the boundary the conservation of boson current leads to a conflict between the two sides, and this can lead to exotic surface physics.
the surface, leading to a system of SO(3) spins with hedgehogs forbidden, which on the three-dimensional cubic latice is known to have magnetic order thus breaking the SO(3) symmetry. 21,23 The locations of the phases and phase transitions in Fig. 4 were determined by studying singularities in the specific heat, as well as by studying the surface magnetization and currentcurrent correlators. As an example of such data, in the main plot of Fig. 4 we show the magnetization, multiplied by the square-root of the volume of the surface (L 3/2 ). This quantity should be constant when the spins are disordered and should grow with system size when they are ordered. The data was taken with β surf = 0 and increasing λ surf . We can see that the spins order at λ surf ≈ 4.
In order to use the current-current correlators to detect the breaking of the U (1) boson symmetry on the surface, we must think through such measurements carefully. In a D-dimensional system of boson worldlines (space-time currents), the argument that ρ J (k min ) ∼ k 2 min ∼ 1/L 2 in the gapped phase relies on the conservation of the currents J. However, this is no longer true if we consider only the surface, as currents can enter and exit from the rest of the system. Therefore we instead measure the current-current correlators in the entire system (though only in thex,ŷ, andτ directions, which are parallel to the surface). When currents are gapped everywhere, this quantity will still be proportional to 1/L 2 . If we measured the current-current correlator only on the surface in the phase where the surface is a superfluid, we would expect the result to be independent of system size. Since we are measuring the correlator over the whole system, the result should The inset shows the phase diagram of the surface of the binding phase, without a Zeeman field. It was obtained by tuning the bulk parameters deep into the binding phase and then varying β and λ only on the surface. We find that in this model our surface always spontaneously breaks a symmetry. At small β surf and λ surf the boson U (1) symmetry is broken and the bosons condense into a superfluid, while at large λ surf the spin SO(3) symmetry is broken and the spins align into a ferromagnet. At small λ surf and large β surf both symmetries are broken. The main plot shows surface magnetization on a sweep in λ surf for β surf = 0. We show mL 3/2 , which is independent of system size in the disordered phase, and grows with system size in the ordered phase. We can clearly see that the SO(3) symmetry is broken as λ surf is increased past a value of approximately 4. All data in this section was taken with β bulk = 0, λ bulk = 5.2.
be proportional to the fraction of the system which makes up the surface, which is 1/L. In Fig. 5 we plot ρ J · L 2 . We see that at large λ surf this quantity is independent of system size, which tells us that the currents are gapped everywhere and the surface is an insulator in the boson degrees of freedom. At small λ surf , ρ J · L 2 ∼ L, which tells us that there is a region whose volume fraction is proportional to 1/L where the bosons are in a superfluid phase. We interpret this as evidence that there is a superfluid at the surface. The surface phase where the U (1) boson symmetry is broken is clearly a superfluid of bosons. In the easy-plane picture the spins can also be thought of as having a U (1) symmetry which is broken when the spins are ordered, and so the phase where the spins order can also be thought of as a superfluid of new particles representing the spins. We can now ask whether these superfluids are trivial superfluids, or the exotic superfluids thought to exist on the surface of a bosonic TI. 6 The exotic properties have to do with the charge and statistics of gapped vortex excitations, and unfortunately we do not have access to these properties in our Monte Carlo. Therefore our surface superfluids cannot tell us whether the binding phase is a bosonic TI.
One predicted feature of the superfluids on the surface of the bosonic TI is that they are connected by a direct transition which is a deconfined critical point. 6 As shown in Figs. 4 and 5, our SO(3) and U (1) symmetries appear to break on the opposite sides of the same point, and so it looks that we also have a direct transition between these phases; however, we would This quantity is constant when the boson U (1) symmetry is preserved and increases linearly with system size when it is broken on the surface. Parameters used are the same as in Fig. 4. We see that there is a phase transition at λ surf ≈ 4.
need larger system sizes to see whether there is indeed a direct transition and whether it is continuous. Though this does not definitively show that the binding phase is topological, but it does provide evidence that we have an unusual field theory on the surface.

C. Surface with Zeeman Field
Another exotic phase predicted to exist on the surface of the bosonic TI is a phase which breaks Z T 2 symmetry and has a quantized Hall reponse. We can try to realize this phase by applying a Zeeman field on the surface of our model to explicitly break the Z T 2 symmetry. We add a term similar to Eq. (12) to our action, but only on the surface of the model. We then expect that there is a surface phase which does not break any U (1) symmetry-i.e., an insulator-but which has the surface Hall conductivity quantized to an odd integer, which is different from the even values expected in the (2+1)D bosonic integer quantum Hall effect. For easy-plane spins, we expect that a finite Zeeman field is needed to destroy the surface superfluids to reach this phase, while for SO(3) spins used here and starting in the spin-ordered phase, we expect to immediately transition to the quantum Hall insulator. 23 In our Monte Carlo study, we actually apply Zeeman field h to both surfaces but take the fields on the two surfaces to have opposite signs. We take λ = 5.2 and β = 0 everywhere (including at the surfaces), which puts the bulk regions deep into the binding or trivial phases respectively, while the surfaces start in the spin-ordered phase at h = 0.
Since the Hall conductivity in the surface system of spins and bosons will be due to correlations between the vorticity of the spins and the boson charges, 7 we will measure these correlations directly before moving on to the more complicated Hall conductivity measurement. In the absence of the Zeeman field, we can use the Z T 2 symmetry, which is reflection of the n variables in the ab-plane, to change the sign of hedgehog and hence the boson charge without changing the spin vorticity, and so such correlations vanish. This can be seen in Fig. 2, where, for example, on the top surface there is one clockwise vortex attached to a positive hedgehog (which is in turn bound to a positively charged boson), and another clockwise vortex is attached to a negative hedgehog. Applying a Zeeman field corresponds to only allowing one type of vortex (↑ or ↓) to pass through the surface. In Fig. 2, this means that only solid lines are allowed to pass through the top surface. We can see that this leads to a correlation between the vorticity of the vortex on the surface and the charge of the boson it is binding nearby.
We can further think about Fig. 2 as depicting a slab of the binding phase. Opposite Zeeman fields on the two surfaces give us a quasi-two-dimensional slab on which vortices are bound to charges with definite relation between the vorticity and charge, and these bound states are proliferated. We have studied such a system in a previous work, 7 and found that its Hall conductivity is quantized to be equal to two. It is reasonable to assume that this conductivity is evenly distributed between the two surfaces, leading to the surface Hall conductivity equal to one on each surface. This intuitive argument reproduces the prediction of Ref. 6.
To measure correlations between vortices and bosons we first define the spin vorticity V µν (R) as where s µ (R) ∈ (−π, π] measures the difference between the spin angles at R +μ and R and brings it to (−π, π]: We can then Fourier transform the vorticity as follows: The prime on the sum indicates we are summing over all sites where k min = (2π/L, 0, 0, 0), and the results are shown in Fig. 6. We see that as soon as the Zeeman field is applied, the vortices and charges become correlated. Unlike the Hall conductivity, we do not expect these correlations to approach any universal value. We do note that the correlations on the two surfaces of the system are approximately equal, which is encouraging as we would expect the Hall conductivity on these surfaces to be equal as well. The differences between the correlations on the two surfaces is likely due to the fact that the surfaces are realized differently on a lattice, but we would expect that universal values such as the Hall conductivity would not have these differences.
In order to measure Hall conductivity, we need to couple both the bosons and spins to external probing gauge fields,  1) and (2) using Heisenberg spins. The horizontal axis is the strength of the Zeeman field applied only near the surfaces (and of opposite sign on the two surfaces). We see that as soon as the Zeeman field is turned on, the correlator takes a non-zero value. The value is non-universal, but is approximately the same on the top and bottom surfaces. Data was taken with β = 0, λ = 5.2 everywhere and the binding phase occupying half of the system, cf. Eq. (17). Since we do not know how to properly couple external gauge fields to this system, we are unable to compute the Hall conductivity. and then use linear response theory to determine the conductivity. When we do this we run into a problem, which has to do with the way the hedgehog currents Q µ were defined. We can see from Eq. (4) that ω µν is a discontinuous function of spins since we required it to be brought to (−π, π]. When including the probing gauge fields, this causes the path integral to be a discontinuous function of the probing fields, which prevents us from taking the derivatives needed for linear response theory, and so we do not know how to calculate conductivity in this system. In the next section we will find a way around this problem, while in the present case we can only appeal to the intuitive argument presented above.

III. REALIZING THE TOPOLOGICAL INSULATOR BY BINDING BOSONS TO HEDGEHOGS OF AN EASY-PLANE
The Heisenberg model discussed in the previous section allowed us to realize a "binding" phase of bosons and hedgehogs. However, in this model we were unable to find definitive evidence that the binding phase is in fact a symmetryprotected bosonic topological insulator, although our indirect arguments are compelling. In this section we introduce a new model, which includes the spin degrees of freedom in a CP 1 representation. This theory has spinor matter fields ("spinons") coupled to a compact gauge field and is a faithful representation of the microscopic spin system with shortrange interactions (i.e., such a lattice "field theory" is "emer-gable" from a local microscopic Hamiltonian). We will see that this formulation allows us to make measurements, such as a Witten effect in the bulk and a quantized Hall conductivity on the surface, which indicate that the binding phase is a bosonic TI.

A. Bulk Phase Diagram
The following action represents the spins in the CP 1 representation: Here the spins are represented by two complex bosonic fields z ↑ ,z ↓ ("spinons"), which satisfy |z ↑ | 2 + |z ↓ | 2 = 1. We can write the z fields as a spinor, z ≡ (z ↑ , z ↓ ) T , and extract the spin n = z † σz, where σ ≡ (σ 1 , σ 2 , σ 3 ) is a vector of Pauli matrices. The spinon fields are minimally coupled to a compact gauge field a µ (R). The last term is a Maxwell-like term for the compact gauge field, which appears after partially integrating out the spinon fields. The variables in the above action live on a cubic lattice, where R gives the position on the lattice and µ, ν are directions. The CP 1 model defined above actually has global SO(3) symmetry, similar to the previous section. In this section we find it convenient to break the SO(3) symmetry down to U (1) explicitly by taking the "easy-plane" limit of the CP 1 model. We align all the spins n in the ab-plane, which corresponds to fixing the magnitude of z ↑ and z ↓ , and allowing only phase fluctuations, i.e., z s ≡ 1 √ 2 e iφs . The CP 1 model in the easyplane limit becomes: Here φ ↑ and φ ↓ are 2π-periodic variables which represent the phases of the spinon fields. We have also replaced the cosine in the Maxwell term by a quadratic "Villain" form, with B µν which are unconstrained, integer-valued dynamical variables residing on the plackets of the lattice. Upon summing over B µν in the partition sum, the third term generates a 2π periodic function of ∇ µ a ν − ∇ ν a µ , and therefore this does not change the universality class of the problem.
Using the Villain form of the Maxwell term is advantageous as it allows us to define the hedgehog current: Note that Q µ (r) resides on the links of the lattice whose sites are labelled by r which, as in the previous section, is interpenetrating with the lattice labelled by R. The above definition is analogous to Eq. (5) with B ρσ ↔ ω ρσ /(2π). The B µν variables have the meaning of "Dirac strings" of the hedgehogs. Thus defined, the hedgehogs in the CP 1 model are actually monopoles of the compact gauge field a µ . We will continue to call them hedgehogs in this work to avoid confusion with a different type of monopole introduced later. We can again study the model described by Eqs. (1) and (23) in Monte Carlo. Equilibration becomes difficult in the regime where K and λ are large, and it is necessary to include composite updates which simultaneously change multiple variables. One such update is to change B µν while also changing J µ so that there is no change in the second term in Eq. (1). Another update is to change both a µ and B µν in such a way as to keep the K term in Eq. (23) small.
We can find phase transitions in this model by looking for singularities in the specific heat, which is defined in the same way as in the previous section. We can identify order in the bosonic degrees of freedom by studying the superfluid stiffness and order in the spins by studying the magnetization. In this easy-plane version of the model, the spin degree of freedom is an XY vector with components (n a , n b ). Since (n a , n b ) = (cos φ spin , sin φ spin ) with φ spin = φ ↑ − φ ↓ , the magnetization is given by: The phase diagram in this model is parameterized by β, K, and λ. As in the previous section, we can make the change of variables in Eq. (10) and find that theJ part of the problem decouples from the spin part. The behaviour of the bosons is the same as in the previous section. When λ is small, the physical bosons are essentially independent of the hedgehogs, and are in the superfluid phase. As λ is increased, they become bound to hedgehogs. The transition happens at λ ≈ 4. The locations of the phase transitions in the spin degrees of freedom are independent of λ, though the nature of the various phases are not. In Fig. 7 we show the phase diagram in the β and K variables, for two cases: λ small and λ large. The phase diagram is consistent with the easy-plane CP 1 model in the literature. 27 Let us first consider the case when λ is small. The bosons will be in a superfluid phase for any β and K. The spin system has the following three phases: i) When β and K are small, the spin degrees of freedom are disordered and the hedgehogs are proliferated. The phase is therefore a conventional paramagnet in the spin degrees of freedom. ii) As K is increased, hedgehogs acquire a large energy cost, and become gapped. This phase was studied in Refs. 23, 25, and 27. It is known as the Coulomb phase 28 because it has an emergent gapless photon and gapped excitations that carry charge 1/2 and interact via a Coulomb interaction. iii) Finally, the phase at large β has a large energy penalty for spin fluctuations, and so the spins order. This phase is a conventional ferromagnet in the spin degrees of freedom.
In the case when λ is large, the spin parts of the Coulomb and ferromagnetic phases do not change, since both of these phases suppress hedgehogs. These phases are now trivial insulators in the boson degrees of freedom. On the other hand, in the paramagnetic phase the hedgehogs are proliferated and the bosons become bound to them. This is the binding phase, which we will argue is a topological phase protected by the U (1) spin × U (1) boson and Z T 2 symmetries.

Symmetries When The Spins Are Represented By An Easy-plane CP 1 Model
When representing our spins with an easy-plane CP 1 model we have explicitly broken the SO(3) symmetry from the previous section down to a U (1) symmetry which corresponds to spin rotations in the easy plane. This complicates our discussion of the discrete symmetries of the model. In the previous section all reflections of n were related to each other by an operation in SO (3), but in the easy-plane case, e.g., reflections in the ab plane (n a , n b , n c ) → (n a , n b , −n c ) are now distinct from reflections in the ac plane (n a , n b , n c ) → (n a , −n b , n c ). The result of this is that there are now two Z T 2 symmetries which protect the topological phase, in the sense that as long as any one of these symmetries is preserved, the topological phase cannot be continuously connected to a trivial phase. The first Z T 2 symmetry is the same as that in the previous section; in the variables of Eq. (23) it reads: This symmetry corresponds to a reflection in the ab plane. The second Z T 2 symmetry is a combination of reflections in the ab and ac planes and acts on spins as (n a , n b , n c ) → (n a , −n b , −n c ). In the CP 1 representation, this symmetry is given by: Note that a symmetry which consists only of a reflection in the ac plane does not protect the topological phase, as under this symmetry we can still apply the Zeeman field in the c direction and connect to a trivial phase as in Sec. II A 1. Note also that we cannot apply fields in the ab plane as this would violate the U (1) spin symmetry. If either of the above Z T 2 symmetries are preserved, the topological phase cannot be connected to a trivial phase. A system with only the first symmetry has a total symmetry of U (1) spin × U (1) boson × Z T 2 , while if we consider the second symmetry the direct product in front of Z T 2 is replaced with a semidirect product. (see Sec. V for further discussion of these symmetries) The Zeeman field introduced in the previous section (and used again below) breaks both of these symmetries.
Note that in the presence of separate U (1) spin and U (1) boson symmetries, we could in principle consider unitary symmetries given by Eqs. (26), (27), by simply omitting the complex conjugation. However, if we imagine introducing a tunnelling between the two U (1) symmetries, which would reduce the total symmetry to U (1) , the discrete symmetry should then treat both spin and boson variables in the same way. We find that this condition is satisfied as long as the above symmetries are understood as antiunitary; hence we will refer to these symmetries as Z T 2 .

B. Observation of a Witten effect
The CP 1 representation allows us to make a bulk measurement which can detect whether our system is a bosonic topological insulator. This measurement, predicted in both the fermionic 29 and bosonic [16][17][18][19] TIs, is called a Witten effect, and is the tendency of an external magnetic monopole in a TI to bind a charge of one-half. In order to justify our claim that the binding phase is a bosonic TI, we now demonstrate that our system exhibits a Witten effect.
The first step in measuring a Witten effect in our Monte Carlo is to add external U (1) gauge fields to the system. These external fields correspond to the U (1) spin and U (1) boson symmetries of the system. We will fix configurations for these fields before performing the simulations, which corresponds to putting our system in some external electromagnetic field. These external gauge fields are distinct from the internal, dynamical gauge field a µ . Similarly, the external monopoles introduced in this section are different from the hedgehogs (which are monopoles of the field a µ ) discussed previously. We will introduce magnetic monopoles in the U (1) spin gauge field and will measure U (1) boson charge. The Witten effect is the statement that the external monopoles of the U (1) spin gauge field will bind half of a charge of the U (1) boson symmetry.
Let us first consider the gauge field coupled to the spin degrees of freedom. Here it is convenient to think in terms of a parton description. This is a description of a spin model by using the easy-plane CP 1 model in which the phases φ ↑ and φ ↓ represent the phases of different types of bosonic "partons". These partons each represent one-half of a physical boson. Each parton carries a unit charge under the internal gauge field a. The physical boson carries unit charge under the external gauge field, which we denote A 1 . The partons carry half a charge under this gauge field, with one parton carrying positive charge and the other negative. To modify Eq. (23) to reflect this, we add ±A 1µ /2 inside the cosines on the first line. Partially integrating out the parton fields then gives compact Maxwell terms in the field combinations a µ + A 1µ /2 and a µ − A 1µ /2, with equal couplings due to the Z T 2 symmetry. We can write each in the Villain form, which introduces two integer-valued placket variables B + µν and B − µν . We can expand and recombine these quadratic terms to get separate terms for the a and A 1 fields, leading to the following action: Here B µν = (B + µν + B − µν )/2, and M µν = B + µν − B − µν . Note that B µν and M µν are not completely independent variables but satisfy the conditions that B µν is 1/2 × integer, M µν is integer, and The variables M µν can be interpreted as the Dirac strings of the external monopoles. We can see that when we introduce external monopoles of odd integer strength, the internal hedgehog variables become half-integer valued-this is a crucial observation for the discussion of the Witten effect. 16 Note that in the case of no external field A 1 and no external monopoles, the B µν variables are integer-valued and Eq. (28) reduces to Eq. (23). Note also that the coupling we wrote for (∇ µ A 1ν − ∇ ν A 1µ − 2πM µν ) 2 [the "Maxwell term" for the external gauge field on the fourth line of Eq. (28)] is special to the preceding spinon-generated argument, while it is expected to be renormalized up by the rest of the universe; in fact, we will assume that the A 1 field is essentially externally controlled and is static. We introduce a monopole into our system by making a specific choice for the external variables A 1µ and M µν . First, we choose a configuration of M µν which will lead to a pair of external monopoles. In our system with periodic boundary conditions, it is not possible to have only a single monopole. We will place external monopoles at coordinates (x, y, z) = (0, 0, 0) and (0, 0, L/2), on the lattice labelled by r. The external monopoles will have opposite charges, with the positively-charged one at the origin. All configurations of external gauge fields will be constant in the τ direction. In order to place external monopoles at these locations, we set M xy (R) = 1 whenever X = −1/2, Y = −1/2, and 1/2 ≤ Z ≤ L/2 − 1/2 (the 1/2's come from the R lattice being displaced from the r lattice by half a lattice spacing). All other M µν are set to zero. By Eq. (29), we must also constrain B xy to be odd half-integers on this Dirac string. Now that we have specified the M µν values which introduce external monopoles, we choose values for the A 1µ variables to minimize the action of the Maxwell term on the fourth line of Eq. (28).
In our simulations we will set A 2µ = 0 everywhere, so that it does not affect the system. It will be used only when computing linear responses.
There are in fact multiple configurations of the variables M µν which give the same configuration of external monopoles. The physics of the system is independent of which configuration of M µν we choose because the various configurations are related by the following gauge transformation: where κ µ (R) is an integer-valued field living on the links of the lattice labelled by R. One can use Eqs. (24) and (28) to check that this transformation does not change the action, including the configurations and energetics of the external monopoles and gauge fields, so our results are independent of the specific choice of the Dirac string M µν .
Having introduced external monopoles into our system, we can begin to see why they should bind half a charge of the bosons. The argument goes as follows. First, when modifying the Dirac string variables M µν to insert external monopoles, we were also forced to modify the variables B µν in such a way as to introduce one-half of a hedgehog at the same locations as the external monopoles. However, we saw in the previous section that hedgehog-boson 'molecules', which have similar bare long-range interactions as just hedgehogs (i.e., carry hedgehog number), are proliferated in the binding phase. Therefore the 1/2-hedgehog which we introduced will be screened by a "cloud" of hedgehog-boson molecules drawn from the rest of the system. This screening is analogous to Debye screening in a plasma. The screening cloud will carry a hedgehog number of one-half, but with opposite sign to the first hedghog, leading to a total hedgehog number of zero. Since in the binding phase hedgehogs are bound to charges, the cloud also carries a boson charge of one-half. Therefore we find that half of a boson charge has bound to the external monopole. We have tested this intuition by direct Monte Carlo simulations.
The above discussion is complicated by two degeneracies in Eq. (28). First, there is a degeneracy between Q τ (0, 0, 0, τ ) = 1/2 and Q τ (0, 0, 0, τ ) = −1/2: e.g., when M xy = 1 on the Dirac string, B xy can be either +1/2 or −1/2 with the same energy. This degeneracy is a result of the symmetry in Eq. (27). In what follows it can be helpful to neglect variations in the τ direction and think about Q τ (0, 0, 0, τ ) as a stationary hedgehog charge at the origin. Because of this degeneracy the statistical mechanics has each of the two states equally probable, which in an infinitely long simulation would lead to zero net hedgehog charge, and no observation of the Witten effect. Second, if we were able to fix the hedgehog number in one of these two states, for example Q τ (0, 0, 0, τ ) = 1/2, we have another degeneracy, between J τ (0, 0, 0, τ ) = 0 and J τ (0, 0, 0, τ ) = 1. This leads to an average boson charge of 1/2 at the location of the monopole. This charge cancels the boson charge from the screening cloud, leading to no observation of the Witten effect.
Despite these degeneracies, we may still observe a Witten effect if the degenerate states are metastable, and the Monte Carlo is stuck in one of the two states. For example, in order to get from Q τ = +1/2 to Q τ = −1/2 one needs to modify all of the B µν on the Dirac string, and the φ ↑,↓ and a µ variables nearby. Such a move would be quite unlikely (and impossible in the limit of an infinite separation between the monopoles). Similarly, to get from J τ = 0 to J τ = 1 one needs to insert a boson loop which passes from one monopole to another, and such a step is highly unlikely with local updates.
The results of our numerics shows that even at the small system sizes that we can access, degeneracy in the J variables is always broken. However, the degeneracy in the Q variables is unbroken, and so we do not observe a Witten effect. Let us consider the Q degeneracy more carefully. One of the two degenerate states has Q = +1/2 bound to a positive external monopole at r = 0, and Q = −1/2 bound to a negative external monopole at r = L/2. The other state as Q = −1/2 bound to a positive external monopole at r = 0, and Q = +1/2 bound to a negative external monopole at r = L/2. Note that to change between the two degenerate states we need to modify variables along a string connecting the two external monopoles, and the probability of this happening reduces exponentially with distance, so in the thermodynamic limit this degeneracy would certainly break, even though on our small systems it does not.
We will break this degeneracy in our system by decreasing the probability that the system will flip between degenerate states. To do this we add the following "biasing" term: where γ is some small real number. We have scanned the system by increasing γ from zero and seen no phase transitions, implying that these small γ do not change the phase we are in. It would be surprising if the above biasing term affected the bulk physical properties of the system, as we are only making a modification to a fraction of the system proportional to L −3 . In addition, though this term breaks the symmetry in Eq. (26), it preserves the symmetry in Eq. (27), and the topological phase is protected as long as either of these symmetries is preserved.
We hope that the above term reduces the probability of switching between degenerate states, so that we can observe a Witten effect. From our numerical results we can see that the above term indeed does this for a large range of γ. Therefore we have broken the problematic degeneracies and removed the obstacle to measuring the Witten effect. We would like to stress that the Witten effect is ordinarily defined for a single monopole, in the thermodynamic limit. The problems we have with degeneracies are artifacts of the fact that we are trying to measure a Witten effect in a finite-size system with two external monopoles. If the action were defined on an infinite system with only one monopole, these problems would not arise as it would take an infinitely long time for the system to change between degenerate configurations.
We observe the Witten effect by measuring the total charge enclosed in a sphere of radius w, centered around the location of the monopole. The precise definition of our measurement is: Note that the sphere discussed above is only a sphere in the x, y and z directions, and we have averaged over the τ direction. We performed simulations with L = 10 and show the results in Fig. 8. At w = 0, we are at the location of the monopole. There is nearly no boson charge bound here. At w = 1, we have already included most of the screening cloud, and therefore measure an enclosed charge close to 1/2, as expected. The fact that w = 1 measures a value close to one half shows that the screening length in the system is quite short. At w = 2, we have included the entire screening cloud, so the charge is even closer to 1/2. As w is further increased, we start to include the screening cloud from the other monopole, which is located at a distance L/2 from the first one. This cancels some of the charge from the first monopole, and so the total charge starts to decrease. When w > L/2, the sphere encompasses the screening clouds from both external monopoles, and so there is a total charge of nearly zero. The values of charge are negative because we set the biasing parameter γ such that at the origin there is a monopole of positive charge, Q = +1/2, and the sign of the charge in the screening cloud is the opposite of the sign of the monopole.
In Fig. 8, we have found that the amount of charge at the site of the monopole (w = 0) is nearly zero. However, this is not universal and in fact depends on the choice of γ in Eq. (31). In our simulations, we find that w = 2 is sufficiently far from the monopole to be unaffected by the change in γ. Figure 9 shows simulations taken with different values of γ. We see that though the amount of charge close to the monopoles (near w = 0 and w = L/2) can be affected by changing γ, the value at w = 2 is always very close to one-half, regardless of what γ is used.
Various other measurements can be made to support our conclusions. Measuring the total charge on each site shows that the half-charge is distributed around the monopole in an approximately spherically symmetric way. We can also use the Witten effect to detect phase transitions out of the bosonic TI. Figure 10 shows the total charge on all the nearestneighbours, as a function of K. We note that the quantized Witten effect disappears at K = 0.6, which is where the phase transition to the Coulomb phase is located, see Fig. 7. (In the Coulomb phase the amount of charge isn't necessarily zero, but it is not quantized and in our simulations we found it to be  Fig. 8. We can compare to Fig. 7, and see that at K = 0.6, when the system transitions from the topological phase to the Coulomb phase, the 1/2-quantization of the enclosed charge abruptly stops. zero.) We also observe the disappearance of the bound charge when the system undergoes a transition to trivial insulator as η is decreased to zero.

C. Surface Phase Diagram
The measurement of the Witten effect is evidence that our binding phase is a bosonic topological insulator. We can now study the exotic physics on the surface of this topological phase. We expect to find the surface phases predicted in Ref. 6.
We begin with some analytical arguments that provide a microscopic derivation of the surface field theory proposed in Ref. 6. We define the surface as in Sec. II B. To uncover the exotic physics, we begin by performing a change of variables from the physical boson currents J µ (r) to new integer-valued variables which satisfy The action for the spins and the new G µ (r) variables is simply For simplicity, from now on we consider situation where the trivial phase region and the SPT phase region are deep in their respective phases: in particular, λ bulk , defined as in Eq. (18), is very large. At first we further simplify the situation by taking λ surf to be very large, in which case we expect the variables G µ (R) to be zero everywhere. We also assume that β is small everywhere, so the spin variables want to be deep in the disordered phase. However, near the two surfaces, the spin configurations must satisfy Focusing on the spins near one surface, say at z R , we can view Q z (x, y, z R − 1, τ ) as simply hedgehog numbers in the corresponding (2+1)D spin system spanned by sites (X, Y, Z = z R − 1/2, T ), and the above conditions correspond to complete suppression of hedgehogs in this spin system. Such a (2+1)D Heisenberg O(3) spin model with hedgehog suppression was studied in Ref. 23 and argued to be described by a non-compact CP 1 field theory (N CCP 1 ). On a simple (2+1)D cubic lattice, the Heisenberg model with complete hedgehog suppression actually has spontaneous magnetic order of spins even when the direct spin-spin interactions are zero. 20,21 However, more generic such models can have a spindisordered phase with a propagating "photon," 20,23 as well as other phases such as coexistence of the magnetic order and the propagating photon. 30 We will see that our findings in the present simulations on the surface of the bosonic TI region are consistent with these earlier results.
Let us now proceed more systematically and, in particular, show how we obtain a generic N CCP 1 model on the surface of the bosonic TI region. For simplicity, we take λ bulk to be very large. For finite λ surf , we need to keep G x , G y , G τ degrees of freedom in the (2+1)D "layer" at z R , while all other G µ (r) are zero. Focusing on the spin variables residing on sites (X, Y, Z = z R − 1/2, T ), the hedgehogs in this (2+1)D system are given precisely by Q z (x, y, z R − 1, τ ), which we will denote simply as Q(x, y, τ ). The structure of the surface theory is subject to constraints (37) Here S matter−gauge represents the first term in Eq. (23) restricted to the surface degrees of freedom. The above is a 3D statistical mechanics model, and a µ , (∇ × a) µν ≡ ∇ µ a ν − ∇ ν a µ , and B µν from Eq. (23) can be now defined as 3-vectors and are denoted by bold-face (e.g., µ-th component of B is 1 2 ǫ µνρ B νρ ). We have suppressed position indices to simplify notation.
This surface theory has spins plus integer-valued "currents" G µ sourced and sinked by the hedgehogs of the spin system. When the "line tension" λ surf for the lines formed by these "currents" is large, we intuitively expect that the hedgehogs of the spin system are linearly confined. It is not immediately clear, however, what happens when λ surf is small. Below we argue that the surface is still qualitatively described by the same "hedgehog-suppressed" field theory, which, however, can be in different regimes and have several different phases.
We can deal with the constraints in the partition sum by changing to new variables which satisfy The action becomes There are now no constraints on the G variables, and we can formally sum these out to obtain a local action which is a function of ∇ × a − 2πB ′ , However, any such action with the compact variables a and divergenceless, integer-valued B ′ can be formally viewed as describing a non-compact gauge field! In the limit of large λ surf , the effective action will have essentially lattice Maxwell form with stiffness K, while for intermediate to small λ surf the gauge field energy will have more complicated but still local form. Thus, the field theory at the surface has the spinon matter fields coupled to a non-compact gauge field. In particular, we expect that the surface can be in the same phases as the (2+1)D easy-plane N CCP 1 model.
We can also study how the surface action is coupled to the external gauge fields introduced in Eq. (28). The minimal coupling between J µ and A ext 2µ , combined with the change of variables in Eq. (33), leads to the following term: It is convenient here to represent Q µ (r) as where we have added a formal zero to the defining Eq. (24). Using this in the preceding equation and integrating the second term by parts, we find both an additional bulk term as well as a surface term which results from taking a derivative of η(r). Focusing again on the (2+1)D layer at z R , we can write the surface contributions as where we defined which is precisely the flux of the non-compact gauge field identified in Eq. (40). When this is combined with Eq. (40), we are left with an effective action for the surface with schematic Lagrangian density This action, which we derived from our lattice model, has precisely the easy-plane N CCP 1 form proposed in Ref. 6. In claiming that this is the correct effective action of the surface, we have neglected bulk terms which in general may also contribute to the surface response properties. We do not have an analytical justification for this choice, though from the Monte Carlo study presented in Sec. III E we find that essentially only the surface terms given above contribute to the measured response properties, and it seems plausible that our argument applies in the limit of a sharp boundary between the topological and trivial phases deep in their respective regimes. Note that the above arguments were based on the assumption that the U (1) and Z T 2 symmetries were preserved in the bulk. If the U (1) symmetry is broken in the bulk, the entire derivation of Eq. (36) based on conserved currents is invalid. On the other hand if only the time reversal is broken in the bulk, the derivation naively holds, but in the matter-gauge sector there is no reason for z ↑ and z ↓ to enter symmetrically-in particular, no reason for them to carry precise +1/2 and −1/2 charges, and the field theory written above is not valid. (In fact, the system will have non-quantized σ xy proportional to the length of the system in the z-direction). Since when the symmetry is broken the bulk ceases to be a topological phase, we of course should not expect exotic physics on the surface in this case.
We can confirm the above arguments, which were made in some simplifying limits, by studying the system in Monte Carlo. We can determine the phase diagram of the surface by looking at singularities in the heat capacity. We can also study the magnetization and current-current correlators as described in the previous section. In this phase diagram we set the bulk parameters so that the system is in the topological phase; specifically, we take β bulk = 0.2, K bulk = 0.2, and λ = 8 (cf. bottom panel of Fig. 7). We then vary the surface parameters and obtain the phase diagram shown in Fig. 11, which is in good agreement with the phase diagram of the N CCP 1 model in the literature. Labels on the phase diagram are taken from Ref. 30, and their relation to labels in Fig. 4 is described below. There are three phases in the diagram. At small β surf the φ ↑,↓ variables are disordered, conserving the U (1) spin symmetry; the U (1) boson symmetry is broken and the corresponding Goldstone mode is precisely the propagating photon in the N CCP 1 theory on the surface, hence the label "Photon Phase" in Fig. 11. At large β surf , K surf the partons z ↑ , z ↓ are condensed, leading to a "Higgs Phase" (in the N CCP 1 language) in which the U (1) spin symmetry is broken but the U (1) boson symmetry is preserved. Finally at large β surf and small K surf both the U (1) spin and U (1) boson symmetries are broken. In the N CCP 1 language, the individual z ↑ and z ↓ are gapped so the gauge field a µ is free to fluctuate, but the "molecular field" Ψ mol ∼ z † ↑ z ↓ , which is precisely the easy-plane spin field, Ψ mol ∼ n a + in b , becomes ordered, hence the label "Molecular Phase" in Fig. 11. We emphasize that the microscopic model we are simulating in (3+1)D has a compact gauge field, and we are detecting the presence or absence of U (1) spin and U (1) boson symmetry breaking on the surface by direct measurements. It is remarkable that the surface phase diagram is captured by the N CCP 1 field theory with non-compact gauge field! All of the phases in Fig. 11 break either a U (1) spin or a U (1) boson symmetry and are therefore superfluids. As in the previous section, without access to the properties of their gapped excitations we cannot directly confirm that they are the predicted surface phases. As in the previous section, our phase diagram contains a direct transition between the superfluid phases, which can be viewed as providing some evidence for the proposed surface physics and is also predicted to be a deconfined critical point. 6

D. Symmetric Surface Phase with Topological Order
Vishwanath and Senthil proposed that it is also possible to have a surface phase which is gapped and breaks no symmetries but has intrinsic topological order. 6 Since this phase is not featured in Fig. 11, we need to add another term to our surface action in order to push the system into this phase. The term we need to add is the following parton "pair hopping term:" 6,16 where we have included proper coupling to the gauge fields. Note particularly that the pair field Ψ pair ∼ z ↑ z ↓ does not carry U (1) spin charge. We can now see what happens to the surface phase diagram (Fig. 11) when we increase t pair trying to induce condensation of Ψ pair . For sufficiently large t pair = 2, we get the phase diagram in Fig. 12. We see that a new phase has opened up at small β surf and large K surf , where, as we will argue, Ψ pair is condensed while the individual z ↑ and z ↓ are gapped. When all the couplings β surf , K surf , and t pair are small, there is nothing which can order the spins or gap the bosons. Therefore we are in the photon phase, which conserves U (1) spin and breaks U (1) boson . To get a pairing phase with no broken symmetries, we need to restore the U (1) boson symmetry without breaking the U (1) spin symmetry, which can be acheived by condensing Ψ pair . Let us first recall how the various terms in the action change the system. The β surf term allows hopping of the partons z ↑,↓ , but even when this hopping is strong the fluctations in the gauge field a µ when K surf is small prevent the partons from condensing. The combination Ψ mol ∼ z † ↑ z ↓ can condense, and this breaks U (1) spin and takes us to the molecular phase. We can see from Fig. 11 that the K surf term on its own does not change the phase of the system if β surf is kept small. However, when it is combined with the β surf term it can prevent fluctations in the a µ field. This gaps the photon, and allows z ↑ and z ↓ to condense. This gives us the Higgs phase, where the U (1) boson symmetry has been restored, but the U (1) spin symmetry has been broken.
With this in mind, we can see why the t pair term brings us into the topological phase. The t pair term is similar to the β surf in that it is also a hopping term, though it hops pairs of partons. Therefore when both t pair and β surf are large the two terms cooperate, which is why the U (1) spin symmetry breaks at lower β surf in Fig. 12 than in Fig. 11. However, when β surf is absent, the t pair only hops pairs of spinons, and so it can condense Ψ pair without condensing the individual z ↑,↓ or Ψ mol . When the t pair term is combined with the K surf term the fluctuations of a µ are gapped. Therefore the phase at large t pair , large K surf , and small β surf can restore U (1) boson without breaking U (1) spin , and this is the phase we are looking for.
We have confirmed that the pairing phase breaks neither U (1) symmetry by direct measurements in the spin and boson sectors. Also, Ψ pair is invariant under the Z T 2 symmetry in Eq. (26), so this symmetry is not broken either, and we indeed do not observe any Hall response on the surface. Therefore we believe that this phase is the fully symmetric gapped phase envisioned by Vishwanath and Senthil, which they argued has intrinsic topological order. Specifically, we expect to have gapped spinon excitations carrying 1/2 of the U (1) spin charge; at the same time, we also have vortices in the Ψ pair field which carry 1/2 of the unit flux of the a µ gauge field and hence 1/2 of the U (1) boson charge; finally, the spinons and the vortices in Ψ pair clearly have mutual statistics of π. Unfortunately, we do not have simple direct measurements to confirm the topological order on the surface. However, the indirect evidence for this scenario is very strong.
Thus, it is suggestive to compare Fig. 12 with the phase diagram obtained in Fig. 1 of Ref. 31. In that work, we studied a (2+1)D model with U (1) × U (1) symmetry and mutual statistics between two different species of bosons. When the mutual statistics is π, we get a phase diagram with the same topology as Fig. 12. The phase diagram contains a topological phase, two phases where one of the U (1) symmetries is broken, and a phase where both symmetries are broken. There is a direct transition between the phases with one broken symmetry, which, if it were continuous, is a candidate for a deconfined critical transition. 32 The surface of our bosonic topological phase is thought to have a similar field theory to that in our previous work, 31,32 and so we expect that the interpretations of the phases and phase transitions are the same in both models. We also remark that in Ref. 7 we presented a microscopic local Hamiltonian which has a topological phase with the same content of excitations. However, that phase also breaks Z T 2 symmetry and in particular has σ xy = 1/2, while the present surface phase respects Z T 2 and has no Hall response.

E. Time-Reversal Breaking and Hall Effect on the Surface
Vishwanath and Senthil 6 predict another exotic phase on the surface of the bosonic topological insulator-a phase which breaks the Z T 2 symmetry and has a Hall conductivity quantized to an odd integer (in units of e 2 h ). We can test this prediction in our Monte Carlo simulations. The first step is to break the Z T 2 symmetry on the surface. By examining Eq. (26), we see that one way to break the symmetry is to replace the parameter β in Eq. (23) by the parameters β ↑ and β ↓ , which appear in the terms containing φ ↑ and φ ↓ respectively. When β ↑ = β ↓ , the Z T 2 symmetry is broken; this is roughly like applying the Zeeman field in Sec. II A 1.
We start with K = 0.4 and λ = 8 everywhere, β bulk = 0.2,  Fig. 11, we see that there is a new phase at small β surf and large K surf . We expect that this surface phase is fully symmetric and has intrinsic topological order. and β ↑ = β ↓ = 0.2 on the surface. This system will have its bulk in the topological phase, and its surface in the photon phase. We break the Z T 2 symmetry on one of the surfaces by increasing β ↑ . We expect that a small increase will not change the properties of the system very much, since z ↑ ∼ e iφ ↑ and z ↓ ∼ e iφ ↓ will still be gapped. 23 However, as β ↑ is further increased, z ↑ ∼ e iφ ↑ "condenses" and vortices in the φ ↑ variables will become gapped. In our simulations we see a singularity in the specific heat measured on the surface, indicating that the system has entered a new phase. We expect that this is the phase that will have Hall conductivity quantized at odd integer. In our simulations we also break time-reversal symmetry in the opposite direction on the other surface by increasing β ↓ . As discussed in Sec. II B, in this setup the top and bottom surface taken together have Hall conductivity adding to two.
We can see that unlike in the Heisenberg model, Eq. (23) is a differentiable function of the probing fields A 1 and A 2 , and so we know how to properly couple the external gauge fields and can use linear response theory to compute the Hall conductivity. If our system has a non-zero Hall conductivity, then we can imagine integrating out the internal degrees of freedom to get the following effective action in terms of the external fields at the surface: where bold face denotes three-component vectors appropriate for the (2+1)D surface, e.g., A 1 = (A 1x , A 1y , A 1τ ), and the above form specifies our convention for σ 12 xy (these units are such that e 2 /h = 1). By taking, e.g., A 1 = (A 1x , 0, 0) and A 2 = (0, A 2y , 0), we have Going to momentum space, we can obtain the Hall conductivity by: where Z is the partition sum, and we also took k = (0, 0, k τ ).
Note that A 1 and A 2 reside on lattices dual to each other, and when defining the Fourier transforms we take the convention to transform in the absolute coordinates of the origins of the links (namely, lattice sites on the dual lattice have absolute coordinates displaced from the direct lattice by half of lattice spacing). Starting from the microscopic model, we can evaluate this conductivity from the current-current correlation functions: where the U (1) spin current on a link R, µ is The measurements are performed at the smallest wave-vector k min = (0, 0, 2π/L), as described in Section II A. Note that the above conductivity measures the response of the U (1) spin currents to applied fields coupled to the bosons (or vice-versa). To study SPTs with single U (1), we can take the usual approach in the literature 6 and "glue" the U (1) spin and U (1) charge by identifying A 1 and A 2 in Eq. (46); the conventional definition of σ xy in the case with single U (1) is then related to the above σ 12 xy via σ xy = 2σ 12 xy . In particular, when the gauge fields are identified the Hall conductivity of a twodimensional system of bosons is quantized to 2 times an integer (in units of e 2 /h), while the Hall conductivity on the surface of a topological phase is an odd integer. Therefore where we present numerical values we show 2σ 12 xy so that the results can be easily compared to the literature values. Figure 13 shows our numerical measurements of the Hall conductivity. The horizontal axis is the strength of the Z T 2 symmetry breaking, which we will loosely call Zeeman field. We see that initially there is no quantized Hall conductivity, until the Zeeman field is strong enough to forbid one species of vortex (realized here by "condensing" the corresponding spinon species). At this point the Hall conductivity increases and reaches the value of approximately 1. Though the value observed is actually slightly less than 1, we believe that a large part of this is a finite-size effect, and indeed as the system size is increased the Hall conductivity gets closer to the expected value. Note that we performed this measurement at precisely z = z R . It is a priori possible that the Hall conductivity could be spread among several values of z near the surface, but this spreading is apparently very small, so that including additional layers does not change the result. In addition to the plot shown, we have performed the same measurement for several different values of K, β and λ, and found that the quantized result is independent of these parameters as long as the bulk stays in the topological phase. This odd integer cannot be observed in a purely two-dimensional system with only short-ranged entanglement, and therefore this observation shows that we are measuring the Hall conductivity on the surface of a bosonic TI.

IV. REALIZING SYMMETRY-ENRICHED TOPOLOGICAL PHASES BY BINDING MULTIPLE HEDGEHOGS TO A BOSON
In all of the above sections, we have studied a system where a single boson is bound to a single hedgehog. In this section we will describe the new physics which arises when our system contains bound states of a boson and multiple hedgehogs. We induce such a binding by making the following modification to Eq. (1): Here d is an integer, and for large λ the action will bind a boson to d hedgehogs, since the λ term is minimized by (Q, J) = (d, 1) × integer. When d = 1 the change of variables in Eq. (10) can no longer be applied. Therefore the phase diagram in this case will be different from the d = 1 case. We can determine the phase diagram by performing Monte Carlo simulations. As an example, the phase diagram for d = 3 is shown in Fig. 14. Phase boundaries were determined from singularities in the specific heat. Note that in the Heisenberg model we can define a maximum of one hedgehog per lattice site, and so this model cannot easily be used to describe the binding of multiple hedgehogs. Therefore all results for d = 1 come from the easy-plane CP 1 model for the spins.  Figure 14 presents the phase diagram in the variables λ and K, with fixed β = 0.1. At small λ and K, there is no energy cost for either hedgehogs or bosons, and they are independent. This leads to a paramagnet of spins and a superfluid of bosons. The U (1) spin symmetry is preserved, and the U (1) boson symmetry is broken. In contrast, at large λ and K both hedgehog and boson currents are forbidden, leading to a Coulomb phase of spins and an insulator of bosons. Here both the U (1) spin and U (1) boson symmetries are preserved, but the spin system has intrinsic topological order. At large K but small λ, the Coulomb phase of the spin system survives and the hedghogs are gapped, but the bosons are condensed breaking the U (1) boson symmetry. On the other hand, at very small K and large λ the hedgehogs are proliferated and bosons are bound to them, and we are in the binding phase; here, neither symmetry is broken, and we will argue shortly that this is a Symmetry Enriched Topological (SET) phase.
Note that similar phase diagram (at fixed small β) for d = 1 contains only such four phases, where the binding phase is the SPT phase discussed in the previous sections. For d = 1 these are the only phases, and due to the change of variables in Eq. (10), the phase boundaries are all straight lines. With d = 1, we have a new phase in the middle of the diagram, and no direct transition from the phase in the lower right to the binding phase. The middle phase can be understood as one in which hedgehog currents are proliferated, but their interactions are still too costly for objects with three hedgehogs and a boson to exist. Therefore such bound states are not proliferated, and individual bosons are also gapped. We expect that this phase preserves the U (1) symmetries from both the spins and bosons and is conventional paramagnet/insulator. The topological phase only arises when K is lowered to the point that it does not penalize significantly objects with a hedgehog current of three and λ is increased to strongly penalize any objects other than the (Q, J) = (3, 1) bound states; at this point these bound states can form and proliferate, and the system enters the topological phase. For other values of d, the phase diagram is expected to have a similar form, with the exception of d = 2, where in our studies of small system sizes the middle phase is not observed and there is a line of phase transitions between the lower left and upper right phases.
Let us now focus on the properties of the binding phase at large λ and small K, which binds d hedgehogs to a boson. This phase is distinct from the topological phase discussed earlier in this work. In particular, it has intrinsic topological order. Condensing bound states of d hedgehogs causes the electric field lines in the phase to fractionalize, i.e., it is possible to have electric field lines of strength 1/d. These fractionalized electric field lines are one of the gapped excitations of the system. The other elementary gapped excitation is a single hedgehog, which binds a boson charge of 1/d. The hedgehog has well-defined statistics as it encircles the electric field lines, and we expect it to acquire a phase of 2π/d when this happens around the elementary fractionalized line. The matter fields z ↑ and z ↓ are confined, but still act as sources and sinks for the electric field lines of integer strength, therefore the line topological excitations in the system are only defined up to an integer, and we can say that the system has Z d topological order. 33,34 In the Appendix we formally demonstrate the above properties by first removing the spinon matter fields and considering a CQED×U (1) boson system in which the monopoles of the compact electrodynamics are bound to bosons. Such CQED×U (1) boson models allow changes of variables similar to those possible in U (1) × U (1) models demonstrating SPT and SET phases of bosons in two dimensions 7 , which allow their properties to be readily determined. After the change of variables has been performed, we can couple the additional spinon matter fields to the CQED sector, and this gives us precisely the CP 1 × U (1) boson model studied in this paper.
The Witten effect and Hall effect measurements can be extended to the cases with multiple hedgehogs. For the Witten effect, the amount of bound charge will be modified, since now for each hedgehog there is a boson charge of 1/d. Therefore the screening cloud will have a charge of 1/(2d). We have studied the cases of d = 2, 3 in Monte Carlo and our results, shown in Fig. 15, confirm this expectation. Recall that when measuring the Witten effect we used a "biasing" term to break degeneracy between positive and negative internal monopoles. This biasing term may also introduces some excess internal monopole or boson charge at the location of the external monopole. This excess is screened by the surrounding system. In the d = 1 case, we found that for most values of γ, such as those shown in Fig. 9, the screening length is quite short and so the Witten effect can still be clearly observed. In the case of d > 1, the screening length seems to be much larger, which can lead to fluctuations of charge larger than the Witten effect we are trying to observe. We have chosen the biasing parameter in Fig. 15 in such a way as to minimize these fluctuations. At other values of γ the Witten effect can still be observed but the observation is less clear due to these fluctuations.
We have also measured the surface Hall effect upon break- ing the Z T 2 symmetry on the surface by applying the Zeeman field as in Sec. III E. Our results for the Hall conductivity are shown in Fig. 16. We find that the surface Hall conductivity is given by 1/d, which is one-half of the value found for a two-dimensional bosonic fractional quantum Hall effect. 7 We can again rationalize this observation by considering a slab of the binding phase as in Fig. 2 with the opposite Zeeman fields on the two surfaces, cf. discussion in the paragraph preceding Eq. (20).

V. DISCUSSION AND CONCLUSIONS
In this work we have realized a three-dimensional bosonic topological insulator in a lattice model which can be studied in Monte Carlo. The model works by binding bosons to point topological defects (hedgehogs). We determined the phase diagram in both the bulk and on the surface of the model. We were able to numerically extract signatures of the topological behavior: In the bulk of the model we observed a Witten effect, while on the surface with broken time-reversal symmetry we found a quantized Hall conductance with values distinct from those possible in a purely (2+1)D system. We also found other surface properties consistent with the bosonic TI, including a direct transition between surface phases with different broken symmetries, and a surface phase which breaks no symmetries and may possess topological order of a kind that would break Z T 2 in a purely (2+1)D system. Finally, we can also realize phases with intrinsic topological order in the bulk by binding multiple topological defects to each boson.
Our model can in principle be used to extract other properties of the bosonic topological insulator. One possible future direction is to determine the properties of the surface phase  transitions, especially the transition between the two different surface superfluids. Another direction would be to find direct evidence of the exotic properties of the surface superfluids and surface topologically ordered phase. Such surface phases have generated much recent interest, as their excitations are expected to have properties not possible in a purely two-dimensional system. 6,35,36 It would also be interesting to investigate the surface physics of the SET phases discussed in Sec. IV.
In an earlier work on the two-dimensional bosonic topological phases, 7 we were also able to reconstruct (starting from Euclidean space-time actions) explicit microscopic Hamiltonians which realized the bosonic integer and fractional quantum Hall effects. We have been unable to do the same in the present three-dimensional case, but this would be a very interesting result. More broadly, the idea of binding bosons to topological defects may continue to yield precise models of bosonic topological phases.
The model in the main text can be thought of as having either U (1) ⋊ Z T 2 or U (1) × Z T 2 symmetry and in both cases we are realizing the same phase. Based on the cohomology classification 3 , the symmetry U (1) ⋊ Z T 2 has a Z 2 2 classification in three dimensions; i.e. there are two base phases with non-trivial topology, but one of the base phases comes from Z T 2 only. There is only one phase that involves U (1) in a non-trivial manner, and it is this phase that we are realizing in our construction and its signature is the Witten effect, which we discuss and observe in Sec. III B. Our construction cannot access the phase that comes from the Z T 2 symmetry only because we require the U (1) symmetry. In the U (1) × Z T 2 case, the cohomology classification gives Z 3 2 , where one of the base phases is again from the Z T 2 alone, while the other two base phases involve U (1) in a non-trivial manner. Of the latter two, again we are realizing the one which has the statistical Witten effect as its signature. The other base phase does not have the Witten effect, its signature is that the monopole is the Kramers doublet under time reversal. 6,10 Our model in principle has more symmetries and could be also deformed to produce this other phase. We do not consider this here but it is a possible direction for future work. We are also not considering beyond cohomology phases, which bring yet another phase due to Z T 2 only in either case. 6,12 . The Appendix contains a more abstract compact quantum electrodynamics (CQED) model where multiple hedgehogs are bound to a boson, which realizes generalizations of SPT and SET phases in more abstract settings where lattice gauge theories are included as microscopic degrees of freedom. Such generalizations have been discussed formally recently, 33,[37][38][39] and similar ideas can be useful for constructing explicit models and analysing physical properties of such phases. As a simple demonstration, in a forthcoming publication we will consider a lattice CQED model where multiple monopoles condense and lead to a novel topological phase of CQED with fractionalized Faraday lines. 34 fractionalized boson particle excitations, as well as non-trivial mutual statistics between the particle and line excitations. 33,34 While the CQED×boson setting may appear artificial, this problem is relevant to the problem of SPT and SET phases of bosons in (3+1)D. Specifically, in the main text we had a model with U (1) × U (1) symmetry, and we represented the first U (1) system as an easy-plane CP 1 model, which has two matter fields ("spinons") coupled to a compact gauge field. We then considered binding of the monopoles in the compact gauge field and physical bosons of the second U (1) symmetry. The crucial difference with the CQED×boson model is the additional matter fields present in the CP 1 ×boson model. We will compare the systems without and with such matter fields and will argue that the matter fields destroy the distinctions of the former model, except when protected by additional discrete symmetries.
Our CQED×boson model is written in terms of compact gauge fields for the CQED part and integer-valued conserved currents for the boson part: The above action can be compared to the action of Eqs. (23) and (50) in the main text. Compared to the main text, the above action is missing the β term which couples the gauge field to "spinons". We also allow the option of binding multiple bosons-here binding c bosons and d monopoles. The boson sector has conserved particles and is coupled to a probing field A ext ρ in the standard way, as in the main text. Unlike the main text, the CQED sector has conserved electric field lines and is coupled to a probing rank-2 field h ext µν (R). In the absence of the binding term, the coupling of the CQED sector to h ext µν is standard; 40 the additional piece in the binding term, while not important for much of the discussed longdistance physics, is the correct form keeping together the "Dirac string" B µν and the probing field h ext µν . Note that in the main text the additional matter fields destroy the conservation of the electric field lines, and therefore the CQED sector cannot be probed by such an external rank-2 field. Variables γ µν realize specific "fluctuating boundary conditions" in the compact gauge field variables; in a representation in terms of an integer-valued electromagnetic field tensor, this corresponds to requiring zero total field for each component. Such details of the boundary conditions are not important for the bulk properties but are nice for a precise mathematical treatment in a system with periodic connectedness assumed here.
From variables B µν (R), we define the monopole four-current Q ρ (r) as in Eq. (24) in the main text. We similarly define the four-vector g ext ρ (r) from h ext µν (R): Note that thus defined Q ρ (r) are conserved currents satisfying ρ ∇ ρ Q ρ (r) = 0, and they also satisfy the condition of zero total current in all directions: Q tot,ρ ≡ r Q ρ (r) = 0.
For the boson sector, we use a representation in terms of integer-valued conserved currents J ρ (r), which satisfy ρ ∇ ρ J ρ (r) = 0. We also require that the total boson current is equal to zero for all directions, J tot,ρ = 0. The primed sum over J ρ (r) in the partition sum signifies all such constraints. The condition of zero total current is again convenient for precise treatment on finite systems [namely, for performing change of variables involving J and Q currents, Eq. (A6) below], while for bulk properties one can ignore these details. The "binding" term parametrized by λ is the key interaction in the action Eq. (A2) and wants to have bound states of d monopoles and c bosons when λ is large [i.e., it is minimized when (Q, J) = (d, c) × integer].
We first separate out the monopoles in the CQED sector as follows: Here B (0) µν (R) is an integer-valued field whose monopolicity gives Q ρ (r) and is treated as a fixed function of Q ρ (r). V µ (R) and N µν are independent integer-valued fields, where the latter appear from careful treatment of the boundary conditions. Schematically, the above arises by dividing all configurations of B µν (R) into classes defined by Q ρ (r) and establishing how to recover all members of a given class from one representative. Furthermore, any Q ρ (r) satisfying ρ ∇ ρ Q ρ (r) = 0 and Q tot,ρ = 0 can be represented as above using some integer-valued B µν (R), so the primed sum over Q ρ (r) can be viewed as signifying these constraints.
A subtle point in the above is the redundancy in V µ (R). One way to address it precisely is as follows. We can argue that from all links of the (3+1)D hyper-cubic lattice, we can select a subset of links such that we can take V µ (R) as independent integer-valued variables on these links, while V µ (R) = 0 on all the other links. We can also argue that the original CQED sector in Eq. (A1) can be equivalently formulated using compact gauge fields a µ (R) that are non-zero only on exactly the same links as the independent V µ (R). We will assume this implicitly in all manipulations below, both for these variables and for the related variables V ′ µ (R),Ṽ µ (R), v µ (R), k µ (R), andã µ (R) that will appear below.
We now operate with the constrained sums over the integervalued currents Q ρ (r) and J ρ (r), taken to be "outside-most" sums in the partition sum. We change to new independent summation variables on each link: 32 Here c and d are the same integers as in the binding term in Eq. (A2), while a and b are new integers such that ad−bc = 1, which makes the above transformation invertible in Z. If c and d are mutually prime, which we will assume throughout, we can always find such a and b, while the arbitrariness a → a + kc, b → b + kd does not affect any physical properties discussed below.
Note that the integration over a µ (R) from 0 to 2π and summation over k µ (R) from 0 to d − 1 effectively corresponds to integration overã µ (R) from 0 to 2π. The same holds for γ µν . We can also express P ρ (r) from Eq. (A7) as P ρ (r) = 1 2 ǫ ρσµν ∇ σBµν and see that B (0) P,µν (R), v µ (R), and n µν only appear in the above combination that givesB µν (R). Furthermore, similarly to Eq. (A5), the summation over constrained P ρ (r) and unconstrained v µ (R) and n µν is equivalent to a summation over unconstrainedB µν . The partition sum becomes