Duality web on a 3D Euclidean lattice and manifestation of hidden symmetries

We generalize our previous lattice construction of the abelian bosonization duality in 2 + 1 dimensions to the entire web of dualities as well as the Nf = 2 self-duality, via the lattice implementation of a set of modular transformations in the theory space. The microscopic construction provides explicit operator mappings, and allows the manifestation of some hidden symmetries. It also exposes certain caveats and implicit assumptions beneath the usual application of the modular transformations to generate the web of dualities. Finally, we make brief comments on the non-relativistic limit of the dualities.


Introduction
The notion of duality, namely that two seemingly different quantum field theories are in fact descriptions of the same physical system, is a powerful concept with far-reaching consequences. Some dualities relate strongly-coupled theories to weakly-coupled ones, allowing one to study strongly coupled systems perturbatively in weakly-coupled dual descriptions. Dualities between two strongly coupled theories do not endow one with such advantage, but they still shed light on aspects such as hidden global symmetries and relevant deformations that are not obvious from looking at only one side. In three spacetime dimensions, dualities have not only deepened our understanding on the non-perturbative aspects of strongly-coupled theories but are also intimately related to the physics of surface states of strongly interacting topological insulators [1,2], superconductor-insulator transitions [3], the half-filled landau level [4], and superuniversality in quantum Hall transitions [5,6], to name just a few. Historically, the fact that in three spacetime dimensions, a vector and second rank antisymmetric tensor (more generally a 1-form and 2-form respectively) have the same number of components, has served as an important motif for dualities. This allows an operator mapping between matter currents and field strengths of a dual U(1) gauge field. This mapping is a key ingredient in the famous boson-vortex duality [7,8] and in Son's JHEP06(2019)038 fermion-fermion duality [4]. The more recently conjectured dualities involve Chern-Simonsmatter theories: various evidence indicates that there are a variety of dualities between two quantum field theories, one with fermionic matter and another with bosonic matter, coupled to Chern-Simons gauge fields, hence bearing the name "bosonization" dualities [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23]. Bosonization dualities makes earlier ideas on statistical transmutation in the context of the quantum Hall effect [24][25][26][27][28][29][30][31] more precise and extends them to gapless (i.e. critical) systems.
The dualities of interest in this paper are dualities that follow from the simplest bosonization duality. The bosonization duality of interest here states that the following two Lagrangians (written in Euclidean signature) flow to the same IR fixed point: AdA. (1.1) Here, φ and ψ are, respectively, a complex boson field and two-component Dirac spinor, and the double arrow indicates that the two theories are dual to one another. A and b refer to the background electromagnetic field and a fluctuating U(1) gauge field respectively. The duality is supposed to hold in the gapped phases with sgn(r) = sgn(m), and most remarkably, also at the gapless critical point r = 0, m = 0. Lately, it was proposed that the S and T "modular" transformations, which explicitly modifies Lagrangians on the both sides, generate an infinite web of dualities starting from the seed duality eq. (1.1) [18,19]. The S and T transformations were employed in describing the global phase diagram of the quantum Hall effect [32], and are defined as the following [33]: (1.2) S can be understood as promoting the classical background gauge field A to the quantum fluctuating gauge field a. Meanwhile, T can be simply interpreted as "adding a Landau level". S and T transformations generate SL(2, Z) group, hence the name "modular". Remarkably, upon considering time-reversal transformation in addition to SL(2, Z) transformations, one can "derive" the boson-vortex duality and the fermion-fermion duality from this web. The duality web hence bridges the two seemingly disparate classes of threedimensional dualities. While the aforementioned dualities pass many non-trivial checks, directly proving these dualities is extremely challenging. Take the seed bosonization duality, eq. (1.1), for instance. The fermion side is free, but the boson side is strongly coupled, both due to the |φ| 4 interaction and the strong gauge interaction. The duality makes statement about strongly-coupled theory without any special symmetries or structures that allow one to employ conventional analytical tools. Recently, there have been microscopic constructions that demonstrate the bosonization duality, one as an exact UV mapping between two partition functions in the Euclidean lattice [34] and another in a wire array [35,36]. These JHEP06(2019)038 methods provide explicit physical insight into the nature of the operators involved in the duality mapping. In this paper, we will focus on the Euclidean lattice construction. The goal of this paper is to extend the Euclidean lattice construction of the bosonization duality eq. (1.1) to other dualities in the duality web. Some of the dualities in the web contains hidden global symmetries, and our lattice construction will make its origin manifest.
There is another issue we want to make manifest with our approach. It has sometimes been advocated that, if the seed duality eq. (1.1) could be rigorously proven, then the entire duality web generated by the S and T transformations would automatically follow. This cannot be taken for granted. While the T transformation involving background fields only is innocent, additional caveat underlies the S transformation. Both in the seed duality eq. (1.1) and the definition of the S transformation, to properly define the strongly coupled dynamical gauge field, we must regularize the theory in the UV. In this continuum theory this is usually done by suppressing its short-distance fluctuations by a Maxwell term. However, once this is taken into account, S and T do not exactly generate an SL(2, Z) group structure; that the SL(2, Z) is a good approximation in the strongly coupled IR should be viewed as part of the conjectural proposal. While this point is usually disguised and overlooked in the continuum presentation, in our microscopic construction its exposure is inevitable.
The paper is organized as follows. In section 2, we explicitly show how to generate the duality web on the lattice, in particular, the implementation of S and T transformation on the Euclidean lattice. We will observe how our lattice construction captures the hidden quantum time-reversal symmetry in the duality web. We also discuss subtleties arising from adding Maxwell terms in the lattice partition function; these subtleties, as we will allude, are present in the continuum language as well though for slightly different reasons. In section 3, we focus on the lattice UV completion of the N f = 2 duality. We discuss non-relativistic limit in section 4. Concluding remarks appear in section 5 2 Generating the duality web on lattice 2.1 Review of the duality web structure In this section, we briefly show how the combination of modular S, T and the time-reversal transformations realize other familiar dualities -particularly, 1) the "fermionization" duality, 2) Son's fermion-vortex duality, and 3) the boson-vortex duality -as part of the duality web. This section is intended to be a review, and our discussion will use continuum field theory; the lattice construction of the same structure, which is the main topic of this paper, will be discussed in section 2.3.
Let's introduce some notation. We will work in the Euclidean signature. We abbreviate the |φ| 4 theory coupled to a U(1) gauge field B as −L WF [B] (B may be both classical or dynamical, depending on the situation. We reserve lower case roman letters for dynamical gauge fields, and capitals for classical, background gauge fields):

JHEP06(2019)038
On the other hand, L +,m Dirac [A] and L −,m Dirac [A] denote Lagrangians of Dirac fermions coupled to a spin-c connection A, 1 with the sign of the superscript denoting the sign of the i 8π ada term that represents the parity anomaly (which is needed for invariance under large gauge transformations): AdA. (2. 2) The duality eq. (1.1) can therefore be written as: In expressing the dualities, sgn(r) = sgn(m) is always understood; the critical point, which is of crucial interest, occurs when r = 0, m = 0. Given a duality, the application of S or T to both sides generates new dualities. Let's start with an illuminating example. Applying T −1 ST −1 to both sides of eq. (2.3) generates the duality between the following two theories: As for the boson side, expanding the last two Chern-Simons terms gives the following: Note that after the expansion, the gauge field a appears nowhere but as a BF coupling with (B −b). Thus, a can be integrated out, acting as a Lagrange multiplier that enforces B = b. We find L E, boson [B] is equivalent to −L r WF [B], and eq. (2.4) says In three spacetime dimensions, a spin-c connection A is much like a U(1) connection, but with the extra property that, if we simultaneously change a fermion boundary condition and change the holonomy of A across that boundary by π, the theory is left invariant. (In higher dimensions one must use a more involved definition; in three dimensions it can be formulated this way because all oriented three-manifolds admit spin structure.) The distinction of spin-c versus ordinary U(1) connection is useful for keeping track of the consistency of the theory [18]; in this paper we will not use further details beyond the consistency check. In this paper, we always reserve A, a for background / dynamical spin-c connection while B, b for background / dynamical U(1) connection.

JHEP06(2019)038
Compared to (2.3), we have essentially "moved the Chern-Simons term" from the boson side to the fermion side, and this may be viewed as a "fermionization" duality. Sometimes one would write using the fact that −L ±,m Dirac [a] differ by a level 1 Chern-Simons i 4π ada. This theory is commonly referred as "Dirac fermion coupled to U(1) −1/2 Chern-Simons gauge field', due to the half-level Chern-Simons term for the gauge field a.
Note that in (2.6), −L r WF [B] is manifestly time-reversal symmetric, while −L fermion is not. Thus, for the duality statement eq. (2.6) to be valid, −L fermion must have an timereversed version dual to itself in the IR. The time-reversed version of is denoted −L T fermion or "Dirac fermion coupled to U(1) +1/2 Chern-Simons gauge field": It is important to note that the mass has flipped sign. Thus, (2.6) and its time-reversed version gives the IR dualities between three theories: where in this picture time-reversal refers to a left-right flip. Son's fermion-fermion duality can be recovered from the aforementioned duality between −L T fermion and −L fermion . Think of applying T S −1 T transformation to eq. (2.9): note that this is precisely the inverse transformation of the operation we applied on eq. (2.3) to generate the duality eq. (2.6). Thus, applying T S −1 T on −L fermion and −L r WF [B], we recover the original theories in eq. (2.3). On the other hand, action of T S −1 T on −L T fermion does not give a free theory. Instead, the new action −L QED is strongly-interacting: The above theory is known as the properly quantized version of QED 3 . The duality between the free Dirac fermion −L +,m Dirac [A] and the QED 3 that we just derived, is Son's fermionvortex duality. The more familiar "illegal" version of L QED in the condensed-matter literature can be recovered by integrating out b -however, there is level-2 Chern-Simons term for b, so integrating out b cannot be justified on general grounds, and yields a loose, local expression with improperly quantized topological terms.
On the other hand, applying T −1 ST −1 transformation to eq. (2.9), we obtain the time-reversed version of the theories in eq. (2.3) as well as that of −L QED , and hence a time-reversed version of Son's duality.
Finally, we will see how the more familiar boson-vortex duality can be recovered from the duality web. Note that (2.8), which is dual to −L r WF [B], has two equivalent expressions. However, the latter expression may be viewed as applying T S −1 to the fermion side of JHEP06(2019)038  Figure 1. Overview of dualities and their relations to each other. Theories in the dashed box are conjectured to be dual to each other; the modular transformation T −1 ST −1 drawn as a dotted red arrow allows one to generate theories on one dashed box from the theories in the different box. We present dualities marked as lines with both arrowheads as exact lattice dualities in this paper. Boson-vortex duality, marked as a dashed line with both arrowheads, cannot be expressed as exact lattice mapping in our paper. eq. (2.3) (but with −m in place of m), and therefore is dual to applying T S −1 to the boson side of eq. (2.3) (but with −r in place of r): (2.11) is a level-1 Chern-Simons term whose dynamics is completely decoupled from other parts of the Lagrangian, so this term can be integrated out. Upon doing this, we arrive at the Lagrangian for the Abelian Higgs model. Hence, we observed that the boson-vortex duality, −L r WF [B] ↔ −L AbHiggs [B] (note that the Wilson-Fisher theory and the abelian Higgs theory have opposite signs of r), is part of the duality web as well.
We summarized the relation between dualities mentioned in this section in figure 1. In [34], we presented exact lattice partition function mappings that can be viewed as UV demonstration of the 3D bosonization duality eq. (2.3), by implementing Chern-Simons terms with heavy fermions. Following the same ideas, in this paper we realize dualities JHEP06(2019)038 in figure 1 as UV lattice mappings as well. We will show that such UV lattice mappings exhibit non-trivial properties anticipated from the continuum viewpoint as well. While we could not write down the boson-vortex duality as an exact lattice mapping, we will show within our formalism that this duality does exist as an IR duality.
Before moving onto the lattice construction, we make some comments on the inclusion of Maxwell terms to regularize the continuum field theories presented in this section. We first focus on the boson side of the seed bosonization duality, i.e., L boson [B] of eq. (2.3). This theory involves a fluctuating gauge field b coupled to the boson, so it is a strongly coupled quantum field theory. Since the gauge charge e 2 runs strong in the IR in 3 spacetime dimensions, a strongly coupled theory is a natural expectation at low energies. To define such a theory perturbatively in the UV, one may regularize the theory with a UV cutoff Λ U V , and add a Maxwell term − (db) 2 4e 2 to the Lagrangian, with the choice e 2 Λ U V . Nonetheless, for energy scales µ of interest, namely µ e 2 Λ U V , the desired behavior will emerge. Similar issues arise in the application of the S transformation as well. Since S turns the background gauge field into the fluctuating gauge field, for the theory to remain welldefined in the UV after the transformation, Maxwell terms for the gauge fields should be supplemented to suppress UV modes of the now fluctuating gauge fields. For example, L E, boson [A] of eq. (2.4) should be modified as the following: Once more, µ e 2 a , e 2 b Λ UV . However, it turns out that for the now-regularized theory to mimic similar IR behavior as one expects from the naive unregularized picture, one should be able to integrate out a to impose constraint B = b. Note that the Lagrangian is still quadratic in a, and after integrating out a, one obtains Note that the delta function δ(B − b) obtained from integrating out a is now "softened" to the superconducting mass term − e 2 a (B−b) 2 16π 2 (in the London gauge); sending e a → ∞ enables us to recover the delta function. In order for the Lagrangian to truly represent the Wilson- Λ UV should be enforced to make UV of the quantum field theory well-defined, and at the same time, to retain the IR physics of the boson side of eq. (2.4).
However, now there is a caveat. Suppose we want to interpret the boson side of eq. (2.4) using the fermion side instead. For that purpose, we would like to invoke the IR duality eq. (2.3) first before we integrate out a. That is, we would want e 2 a e 2 b so that we can integrate out the b gauge field first to obtain the Dirac Lagrangian, which is contrary to the condition above. Therefore, to "generate eq. (2.6) from eq. (2.3)", on the boson side we assumed e 2 b e 2 a , while on the fermion side we assumed e 2 a e 2 b . We have therefore implicitly made an assumption that changing the relative strengths of e b and e a JHEP06(2019)038 over some finite range of energies does not alter the IR physics significantly -including the physics at the critical point. More compactly, such implicit assumption underlies the seemly innocent rule "SS −1 = 1" once the S and S −1 are regularized by Maxwell terms; in other words, it underlies the recognition of S and T as "generating" the algebra of SL(2, Z) to a good approximation in the IR. This implicit assumption has not been emphasized in the literature, but it will inevitably become manifest in our lattice construction.

Review of the Euclidean lattice construction of the bosonization duality
For the readers who are not familiar with our Euclidean lattice version of the bosonization duality [34] and also for the sake of setting up the notation that carries throughout the rest of the paper, here we briefly review the construction presented in our previous work.
Our construction involves a lattice gauge theory constructed on the 3D cubic lattice. The two main ingredients of our construction are the XY model lattice action (2.14) and the Wilson fermion action (2.15) In the expression above, n is an index that labels each lattice site. Unit vectors that span the three-dimensional cubic lattice are labeled byμ =x,ŷ,ẑ; hence, n +μ refers to a site neighboring n. nµ is an index associated with a link pointing from site n to site n +μ. σ µ refers to 2 × 2 pauli matrices. Matter fields such as two-component Grassmann fields ψ and χ, and the boson phase field θ (a U(1) compact variable) live on lattice sites. Meanwhile, the U(1) gauge fields A and b, which are also compact, reside on the links of the lattice. In the action of eq. (2.14), tuning T induces a phase transition from a trivially gapped phase to a superconducting phase in which A is higgsed. The critical point between two phases is a lattice realization of the O(2) Wilson-Fisher fixed point. Meanwhile, the action in eq. (2.15) corresponds to two-component free Dirac fermions supplemented with fourfermi interaction U that takes a role of a "counterterm" in the exact mapping. When U = 0, one can explicitly compute the "bubble diagram" (the leading contribution to the gauge field self-energy) to see that the level-k Chern-Simons term for the gauge field A is generated; k takes different integer values according to the value of M , as the following: The regime of our interest is 1 < M < 3 -since this regime corresponds to a k = 1 gapped phase, one can understand IR lattice Chern-Simons action as generated from integrating out gapped fermion with proper M value. Thus, the gapped fermion action can be used to implement a Chern-Simons term for the gauge field on the boson side of the duality eq. (1.1). Another interesting point is M = 3, the critical point between k = 1 and k = 0 phase. At M = 3, one can check by expanding near k = (0, 0, 0) that a single gapless Dirac cone emerges in IR; one can employ the lattice action eq. (2.15) near M = 3 to implement the free-fermion side of the duality eq. (1.1), with m M −3. Finally, we stress that four-fermi interaction is irrelevant on the "free-fermion line" U = 0. Turning on small U term may renormalize the critical value of M slightly from 3 by order U , yet qualitative features of free-fermion phase diagram remain valid in small but finite U . (Although M = 1 is also a critical point, it is one across which the Chern-Simons level k jumps by 3, and hence not generic. It is an artifact of the use of cubic lattice, so it does not describe universal physics and we do not focus on it.) Given these two main ingredients, upon choosing proper parameters, one can expect two partition functions defined on the lattice, Z B and Z F , to be the lattice realization of boson and fermion side of the bosonization duality eq. (2.3) respectively: (2.17) where the integration measure is simply defined by .
(2. 18) In particular, on the boson side, 1 < M < 3 so that Z W implements level 1 Chern-Simons, and the coupled theory has some critical T c that r ∼ T − T c . On the fermion side, m ∼ M − 3 of the IR Dirac fermion. We also want both U and U to be small so that we may safely use the fact that they are irrelevant.
Remarkably, one can explicitly integrate out θ and b field in Z B to show that where M and U are related to T , M and U as the following: Here I 0 and I 1 denote the modified Bessel functions -their details do not matter, it only matters that their ratio monotonically increases from 1 to infinite as T increase from 0 to infinite. Therefore, started with 1 < M < 3, there must be some T c around which M 3. Of course, the critical mass might not occur exactly at M = 3 because we must also take U into account. A perturbative calculation in U and U (compared to the lattice scale) shows that there is a parameter regime in which Z W [A; M , U ] can be interpreted JHEP06(2019)038 as a single gapless Dirac fermion, and at the same time Z W [A − b; M, U ] as a heavy Dirac fermion generating a level-1 Chern-Simons term. The critical temperature T c thus occurs at an order of U . In addition to a partition function mapping, one can also show that there is the following relation between correlation functions, for any configuration of the background gauge field A: e iθn 1 χ n 1 · · · e iθn k χ n k e −iθ n 1χ n 1 · · · e −iθ n kχ n k ∝ ψ n 1 · · · ψ n kψn 1 · · ·ψ n k (2.21) where χ and ψ respectively denotes fermion field in Identifying e ±iθn as boson creation/annihilation operator and χ n as monopole operators, this operator mapping fits nicely into continuum picture in which a composite object of a boson and a monopole forms a Dirac fermion.
We make a final remark in relation to the Maxwell terms for the fluctuating gauge field b. In the continuum theory, we mentioned that a Maxwell term must be introduced to regularize the strongly coupled theory. In the lattice presentation, since the theory is already regularized, why do we still include a Maxwell term? First, it can be included, so we should include it for generality. But a better reason is that, in the partition function Z B [A; T, M, U ] of eq. (2.17), the idea that integrating out the heavy Wilson fermion Z W [A− b; M, U ] yields predominantly a Chern-Simons term can only be justified if A − b is a small fluctuation. If it fluctuates strongly the connection to the continuum formulation becomes unclear due to the presence of large irrelevant operators (which in turn imply that we are far from the fixed point described by the continuum theory). To suppress these UV modes, one may add the Maxwell term by hand. For example, one can imagine including the following, "non-compact" Maxwell [7], into the partition function Z B [A; T, M, U ]: where labels the lattice plaquettes, × is the lattice curl, and l is a closed integer field on the plaquettes -"closed" meaning that the sum of the l field coming out of the faces of each lattice cube must vanish. 2 Similar to the continuum case discussed previously, we set g 2 (lattice scale) to suppress UV modes and µ g 2 to preserve non-trivial IR dynamics. Note the similarity to the requirements for e 2 in the bosonization duality in the previous section. Hence, in the lattice construction, the similar issue of including Maxwell terms we observed in the continuum formulation persists as well. However, the origin is very different: the lattice constructions are always UV-regularized, but the Maxwell terms are needed to make the connection to the continuum theory manifest. 2 Since the l field is closed, locally it may also be written as the lattice curl of an integer gauge field on lattice links, l = × m. The integer gauge field m may then be combined with the U(1) gauge field b as b + 2πm which takes real value, hence the name "non-compact". On the other hand, if we relax the closedness condition of the l field, then ZM will be the Villainized version of the "compact" Maxwell term [7], whose usual version has −S = (1/2g 2 ) cos × b in the exponent, without the l field. In the Villainized "compact" Maxwell, the integer l field is interpreted as the Dirac strings, whose end points (i.e. lattice cubes out of which l is not closed) are fluctuating monopoles, so that the total flux × b + 2πl appears non-conserved.

JHEP06(2019)038
The current status [22] of the microscopic construction is therefore that, without a Maxwell term (i.e. g 2 → ∞ on the lattice), we have an exact lattice duality mapping, but the connection to the continuum definition of theories (where e 2 is small in the UV) can only be argued based on the idea that e 2 in the continuum runs towards infinite in the IR anyways. Whilst, with a Maxwell term (i.e. g 2 small on lattice), the connection to the continuum theories is better justified, the exact duality map is absent. 3 Later we will see that the implementation of the S transformation is under the same status, and the issue is related to the caveat that also exists in the continuum theory which we mentioned at the end of section 2.1.

Lattice presentation of duality web and time-reversal symmetry
Now, we will see how the above continuum picture of the duality web is realized on a lattice. The main idea is that the T operation in the continuum is to add a Chern-Simons term i 4π AdA; also the T −1 ST −1 operation can be thought of as promoting the background gauge field B to a dynamical field b and to add a Chern- where 1 < M 1 < 3, and we have set the irrelevant coupling U in eq. (2.17) to zero for convenience. We stress that the equivalence between the two partition functions in the equation above is exact, with M and U given in terms of T , M , U = 0, as defined in eq. (2.20).
In the continuum, one could integrate out the gauge field a on the boson side of eq. (2.4) to see that the theory on the boson side is actually φ 4 theory in disguise, if one implicitly assumes the Maxwell term can be neglected in an S transformation. However, such procedure of integrating out the gauge field cannot be done exactly on the lattice partition function Z B [B; T, M, M 1 ]. An intuitive argument that this must not be possible on the lattice is that, Maxwell and other irrelevant terms are necessarily generated along with the Chern-Simons term from the heavy Wilson fermion.
Let's formulate this in detail. Ignoring Maxwell terms, is the identity, i.e. imposes the "Lagrange multiplier constraint" δ(b − B); with the Maxwell terms included, this delta function is softened to the Higgsing of b − B. We expect this on the lattice. Note that Z B [B; T, M, M 1 ] in (2.23) can be written as 3 Some relevant discussions on this point can be found in appendix C of [37].
by a charge conjugation transformation.

JHEP06(2019)038
where Z DW is the partition function of a "doubled" set of Wilson fermions and is given by Thus, the above is a lattice implementation of (T S −1 T ) (T −1 ST −1 ). We expect this partition function, in appropriate range of parameters, to exhibit the Higgsing of b − B.
Remarkably, Z DW can be simulated using Monte Carlo methods without a sign problem if we set M = M 1 (there is no obstruction to doing so) and if we turn off the gauge fields b and B -it is sign problem free simply because the integrand of eq. (2.24) becomes manifestly positive. If we show Z DW at these parameters exhibits a superconducting phase, then it would Higgs b − B at least for small values of b − B (this is in parallel to our e 2 b e 2 a discussion in the continuum). Indeed, in the numerical simulation -whose details we present in the next subsection -we show a fermion bilinear, with charge 1 with respect to B and charge −1 with respect to b, condenses in an extended range of M , Now we can integrate out b above and find that it exactly maps to Note that in Z F , the role of the heavy fermion that simply generates the Chern-Simons term is played by Z * W ; Z W plays the role of the light Dirac fermion, whose mass parameter drives the theory to the criticality. In Once again, as for the lattice action thus, Z F flows to the free-fermion Z W with a = A. However, Z QED is expected to be interacting, and it is straightforward to recognize this theory as the lattice version of QED 3 that appears in fermion-fermion duality. Finally, we comment on the inclusion of Maxwell terms. Consider the fermion side Z F in eq. (2.23) first. For the heavy Wilson fermion factor Z * W [a − B; M 1 , 0] to be reliably interpreted as the Chern-Simons term, one should include Maxwell term eq. (2.22) Z M [a; g a ] on both sides of the partition function. This Maxwell term plays a very similar role to the one played in the boson side of the lattice bosonization duality in the previous subsection. The coupling g a also follows a similar constraint of µ g 2 a (lattice scale). Note that adding this term does not break the exact partition function mapping in eq. (2.23). Also, the theory eq. (2.24) can be still simulated without a sign problem upon including the Maxwell term. 5 On the other hand, let's consider the boson side Z B in eq. (2.23). We claimed its physics should realize that of Z XY , because the Z DW factor eq. (2.24) Higgses out the field b − B. This Higgsing is reliable only for small fluctuations of (b − B). Therefore, to make it reliable, one needs to include the Maxwell term for the gauge field b on the boson side Z b . For g b to suppress dangerous UV large fluctuations, it should be much smaller than the scale encoded by lattice partition functions and g a but much larger than IR modes. Hence, the hierarchy µ g 2 b g 2 a (lattice scale) should be satisfied to make connection to the boson side of the continuum field theory. Now a caveat similar to that at the end of section 2.1 arises. Upon including Maxwell terms for b, the whole exact mapping from the boson side to the fermion side becomes unavailable, as mentioned at the end of the previous subsection. But for the fermion side to be sensibly interpreted as "Dirac fermion plus Chern-Simons term", we expect g 2 a g 2 b , which is opposite to the range of JHEP06(2019)038 parameters for the boson theory to become Z XY . So we encounter the same situation as in the continuum, where different sides of the "new" duality requires different ranges of parameters, in order to have it "derived" from old dualities. It is also worth noting that, upon the inclusion of the Maxwell terms for a and b, the time reversal picture under eq. (2.25) becomes non-exact.
Once more, since not including Maxwell terms make the non-trivial properties of the dualities manifest, it is natural to guess that these Maxwell terms do not affect dynamics drastically; the partition functions we present without Maxwell terms are expected to capture the same IR physics as the ones with the Maxwell terms and hence with clearer relation to the continuum field theories. Hence, in the remainder of the paper, we will omit Maxwell terms for the gauge fields. However, readers should note that this issue persists throughout any lattice dualities presented in this paper, as a manifestation of the same caveat in the continuum theories.

Numerical confirmation of the implementation of "Lagrange multiplier"
In this subsection, we present more detailed analysis of the lattice action Z DW in eq. (2.24). We first present an analytical argument that Z DW should admit a phase with long-ranged superconducting order that provides a Meissener effect for the gauge field B − b. When b is promoted to being a fluctuating gauge field, this superconductivity higgses out B − b, and Z DW effectively implements δ(B − b) in the deep IR. After the analytical argument, we present strong numerical evidence for the emergence of a superconducting phase in Z DW . Once more, we emphasize Z DW does not represent partition function of field theories that appear in dualities -instead, superconductivity of Z DW [B, b, M, M ] in certain parameter ranges of M gives a strong evidence that Z B flows to a Wilson-Fisher fixed point in IR. Hence, our primary interest in numerics is to see a gapped superconducting phase, not a critical phase. We reserve numerical studies of lattice actions corresponding to actual field theories in dualities and their critical points to future works.
The analytical picture that a superconductor-insulator transition should be present in Z DW is the following. Using techniques similar to those used in the lattice construction of the bosonization duality [34], one can expand the exponentials on the lattice links in the action Z DW as follows: (2.29)

JHEP06(2019)038
Here, (σ µ ) * denotes the complex conjugate (not Hermitian conjugate) of the Pauli matrices chosen for Z W . Just as in the derivation of the lattice bosonization duality, here this expansion allows us to exactly integrate out the gauge field. The idea is that, all terms in this expansion are either proportional to e ±ianµ , or have no dependence on a nµ . Those terms that are proportional to e ±ianµ vanish upon integrating out the gauge field. This gives: Now, let us take a closer look at the second line of the above equation that contains the hopping terms that favor simultaneous hoppings of χ and ψ fermions. For concreteness, we choose the following basis: (2.31) In addition, it is convenient to define the following fermion bilinears (the motivation is that (σ µ ± 1)/2 are projectors to the positive / negative eigenvectors of σ µ ): φ n,0 = 1 √ 2 (ψ n,↑ χ n,↑ + ψ n,↓ χ n,↓ ),φ n,0 = 1 √ 2 (χ n,↑ψn,↑ +χ n,↓ψn,↓ ) φ n,x = 1 √ 2 (ψ n,↑ χ n,↓ + ψ n,↓ χ n,↑ ),φ n,x = 1 √ 2 (χ n,↓ψn,↑ +χ n,↑ψn,↓ ) φ n,y = −i √ 2 (ψ n,↑ χ n,↓ − ψ n,↓ χ n,↑ ),φ n,y = i √ 2 (χ n,↓ψn,↑ −χ n,↑ψn,↓ ) φ n,z = 1 √ 2 (ψ n,↑ χ n,↑ − ψ n,↓ χ n,↓ ),φ n,z = 1 √ 2 (χ n,↑ψn,↑ −χ n,↓ψn,↓ ). (2.32) These fermion bilinears can be understood as hard-core bosons formed out of two fermions, and they enable us to write down eq. (2.30) in a more compact way: Let's focus on the hardcore bosons in the second line. Note that hardcore bosons of the form φ n,µ can only hop across links along theμ direction; they effectively "live" in one JHEP06(2019)038 dimension, and in the absence of any coupling between them, they cannot condense in the thermodynamic limit due to the Mermin-Wagner theorem. 6 By contrast, the φ n,0 boson hops along all three lattice directions, and may condense under appropriate choice of microscopic parameters. Thus, we may neglect φ n,µ , integrate out the fermions, and study the effective theory for φ n,0 . Away from critical points, the effective action for φ n,0 will be local and takes the usual Ginzburg-Landau-Wilson form. When the quadratic term is negative, the φ n,0 boson will acquire long range order, and this corresponds to a superconducting phase. Since φ n,0 is charged under the gauge field B −b, it expels precisely this combination of gauge fields via the Meissner effect in the superconducting phase. In the present argument, the role of the fluctuating gauge field a is to "glue" the two fermions ψ, χ, and when this composite field condenses, a superconducting phase results.
To demonstrate the validity of the above picture explicitly, we sample Z DW directly in eq. (2.24) in a determinant quantum Monte Carlo (DQMC) simulation. When M 1 = M 2 = M , B = b = 0, one can compute the above correlator from the determinant-Monte-Carlo like algorithm. 7 Note that while the theory is strongly interacting due to the presence of the fluctuating gauge field a, the fermions can exactly be integrated out on the lattice for each configuration of a. The resulting fermion determinants come in complex-conjugate pairs, and there is no sign problem -as can be seen from the fact that the integrand of eq. (2.24) is manifestly positive. The details of numerical implementations are explained in the appendix. Based on the analytic picture above, we identify the following correlator that diagnoses off-diagonal long-range order that signals the emergence of superconductivity: φ n 1 ,0 φ n 2 ,0 .
(2.34) 6 Of course, couplings between them are not forbidden by symmetry and will therefore be generated. To the extent that such terms are small, the transition temperature at which long range order develops will be parametrically small. 7 The requirement M1 = M2 is merely to avoid the fermion sign problem. Recall that for energy scales µ M1, M2, the precise values of these masses are unimportant. Consequently, we set them equal in order to avoid the fermion sign problem. We present the simulation results done on the system with size 8 ×8×20, with periodic boundary conditions on all directions. We chose one direction longer than the others to see long distance behavior of correlators more clearly; since the system we are investigating is isotropic, it is natural to expect that the behavior along one direction persists in other directions as well.
We see a clear signature of off-diagonal long-range orders, as shown by the blue curves in figure 2, for low values of M : the distance-correlator plot in linear-linear scale shows that at the long distance the correlator reaches some constant, signaling emergence of long-range order. Repeating the same analysis on 1 < M < 4 we found the transition point M ≈ 1.8 where the correlation function behavior sharply transitions into that of an exponentially decaying function: the red curves in figure 2, especially in linear-log scale, show rapid exponential decay of correlator with respect to distance. In figure 2 and all other graphs presenting the results of Monte-Carlo simulations, error bars are explicitly included. However, in case of determining emergence of superconductivity from the correlator eq. (2.34), correlators can be computed with very small errors, often the size of errors on the graphs being smaller than thickness of two bars that represent the range of the error.
Another interesting quantity that probes long-range order and signals sharp transition at M ≈ 1.8 is χ 2F (in continuum, this can be interpreted as a renormalized pair susceptibility), defined as: The plot of χ 2F for different mass values in figure 3 clearly shows that this number sharply drops at M ≈ 1.8, indicating loss of long-range order at this M value.
We present this Monte-Carlo result as a clear evidence that the theory defined by Z DW exhibits a Higgsing of the gauge field B − b (at least for small B − b) in the IR. That the phase transition occurs at M ≈ 1.8 rather than M ≈ 3 is because interaction in general changes UV parameters; it will take different values if we include a Maxwell term for a, which we will explore in future works.

On the boson-vortex duality
So far, we have observed that including additional heavy Wilson fermions to the euclidean lattice construction of the 3D bosonization duality allows one to realize other dualities in the U(1) duality web as exact partition function mappings. However, we did not present how similar methods would realize the boson-vortex duality. Here, we illustrate the difficulty of realizing the boson-vortex duality as an exact partition function mapping; nevertheless, we provide here some numerical evidence for the claim that topological terms involving heavy Wilson fermions can realize this duality.
Recall that in the continuum derivation of the boson-vortex duality, we used something that we have not otherwise used in the previous subsections: the two equivalent expression of (2.8), which relies on M h ] for some 1 < M h < 3. The two implementations certainly are different in the UV. As the continuum "derivation" of the boson-vortex duality from the bosonization duality relied on the identification above, this procedure does not yield an exact mapping on the lattice for the boson-vortex duality.
Below, we present a lattice implementation of the Abelian Higgs model using the elements we have introduced, though this implementation does not have exact mapping relations to the previous lattice bosonization duality via any lattice S and T transformation. Our strategy is to construct the BF Lagrangian out of level 1 Chern-Simons terms as the following: 8 AdA. (2.39) 8 The reason we have the include the seemly useless A is that, the connection in a level 1 Chern-Simons must be spin-c instead of U(1), because a level 1 Chern-Simons is naturally related to a massive Dirac fermion, which couples to a spin-c connection. In fact, all topological gauge terms consistent with the spin-c property may be built out of linear combination of level 1 Chern-Simons terms. 9 In ZBF [A, b, B; M ], the spin-c connection A is included for heavy fermions to satisfy the spin-charge relation -its role is purely auxiliary and may be assumed to be zero for the analysis.

JHEP06(2019)038
Observe that only the first two heavy fermions of Z BF [A, b, B; M ] contribute to the dynamics, since the latter two are not coupled to dynamical gauge fields. We saw in section 2.4 that by tuning M to a proper value and sending T → ∞, one can enter the phase in which a gauge-invariant fermion bilinear condenses. Meanwhile, when T → 0, Z XY [b; T ] Higgs the gauge field. Since the gauge fluctuation plays a key role in condensation of the gaugeinvariant fermion bilinear, Higgsing out the gauge field will "decondense' the composite objects of two fermions. Thus, the partition function faithfully captures a theory in which a bosonic object condenses (stays trivially gapped) in high temperature (low temperature).
Interpreting the fermion bilinear as a monopole operator is consistent with both the phase diagram of the Abelian Higgs model and the recurrent motif of the paper in which lattice heavy fermion operators are often understood as monopole operators in continuum.
When Peskin [7] first suggested the boson-vortex duality, he demonstrated it with a lattice model. The advantage in our lattice implementation, however, is that it explicitly preserves U(1) nature of the gauge fields -i.e. the gauge field cannot be arbitrarily rescaled, and the level must be quantized. The monopole operators are manifestly available and their condensation is explicit.

Implementing generic T and S transformation on the Euclidean lattice
First, let us explore the general structure of continuum quantum field theories in the duality web. One can generically write any element τ ∈ PSL(2, Z) as: Here, m i can take any integer value. Given a Lagrangian −L[A] with the background U(1) gauge field in the Eulicdean signature, τ , as defined in eq. (2.40) corresponds to the following modification of the Lagrangian: Here, I and J are indices that run from 1 to n + 1. We define a I = (a 1 , a 2 , · · · , a n , A); a 1 , a 2 , · · · , a n corresponds to the fluctuating gauge field introduced by S transformation, and A = a n+1 corresponds to the background gauge field. K IJ is a symmetric matrix with JHEP06(2019)038 integer entry, defined as the following: (2.42) Diagonal elements K II = m I corresponds to T m I transformation in eq. (2.40), and offdiagonal elements correspond to BF terms introduced by S transformations in eq. (2.40).
The following two considerations allow one to write down eq. (2.40) in a slightly different form: • one has a freedom to add any integer multiples of A into the fluctuating gauge field to get an equivalent expression, i.e.
for i = 1, · · · , n give equivalent expression since this is tantamount to adding constant to the integration variable.
• The quantum gauge fields a 1 , a 2 , · · · , a n can be "mixed' by SL(n, Z) transformations. For a n × n matrix S ∈ SL(n, Z) The change of the variable preserves integration measure and compact nature of U(1) gauge fields. Thus, rewriting Lagrangian on the right hand side of eq. (2.41) in terms ofã 1 ,ã 2 , · · · ,ã n is a valid variable change.
Above two variable changes can be accomodated by considering the variable change given by a matrix M ∈ SL(n + 1, Z) with some constraint: The constraint on the last row of M is present to make sure that the definition of the background gauge field is not modified. Defining (M −1 ) T KM −1 =K, eq. (2.41) can be equivalently rewritten as: Note thatK IJ is still real symmetric matrix but may not be a band diagonal matrix as the original K-matrix in eq. (2.42) is.

JHEP06(2019)038
The K-matrix multi-component Chern-Simons term eq. (2.46) primarily consist of two terms: the BF term that couple two different gauge fields and the diagonal Chern-Simons term that simply reduce to the Chern-Simons term for single gauge field. The important point is that, all K-matrix elements have corresponding lattice heavy fermion representations. For example, we know that under the proper choice of the parameter M (for later convenience, we set the four-fermi counterterm U = 0 for implementing PSL(2, Z) transformation on the lattice), AdA . (2.48) A Chern-Simons term with higher levels can be simply generated by placing multiple heavy fermions with the same gauge field coupling. Now, let us consider the heavy fermion representation of the BF coupling. Note that (2.49) Hence, we expect following heavy fermion actions to realize BF coupling in IR: (2.50) Similar to the case of the lattice heavy fermion Chern-Simons action, "Complex-conjugated version" of the fermion action on the left side of the above equations can generate BF coupling with opposite sign. Also, BF-coupling with integer coefficients other than ±1 can be generated by multiplying multiple lattice actions, each of which generate unit-level BF couplings. We conclude this subsection with two remarks. First, there are a number of different ways to implement general PSL(2, Z) transformation on the lattice, primarily due to two factors: there is a freedom of changing variable with respect to the restricted SL(n + 1, Z) transformations to obtain different K matrices, and there are multiple ways of implementing the same BF couplings/Chern-Simons terms with lattice heavy fermions. We argue that these different lattice actions all flow to the same IR fixed point, yet we will see that some non-trivial properties of dualities in the web is more manifest in a certain lattice action. Second, the "free-fermion picture" tells us that M should be tuned to be 1 < |M | < 3 to implement a Chern-Simons term with level ±1, However, in actual implementations,

JHEP06(2019)038
gauge fields coupled to the heavy fermion is quantum gauge fields, and their fluctuations generally renormalize the range of M in which the lattice action is expected to flow to the IR fixed point we want.

A microscopic model for N f = 2 self-duality
In this section, we focus on the self-dual N f = 2 three-dimensional QED [19,38,39]. This duality has found applications in condensed matter physics, most interestingly as an effective description of surface states of (3 + 1)D bosonic topological insulators [38,40] and deconfined quantum criticality [41]. This duality is remarkable for it being self-dual, as we will see soon, and has an emergent SU(2) × SU(2) global symmetry at criticality, as we will discuss in later subsections. The goal of this section is to explicitly construct the operator maps of the duality, and also provide a picture for how the SU(2) × SU(2) symmetry emerges at the criticality.
Consider the duality between the following three theories: Before we discuss it on the lattice, let us review how this duality can be obtained from the duality web that we discussed in the previous section. Let us start from two copies of the Wilson-Fisher bosons: One can consider applying the following series of modifications: Note that X and Y are U(1) background gauge fields, so such substitution is benign. Also, such substitution respects compactness of gauge fields (though the inverse of this substitution does not, because this substitution changes the number of gauge fields).

Attach the BF term
3. Promote B f to the fluctuating gauge field b.

JHEP06(2019)038
The above procedure modifies −L DWF [X, Y ] to −L N b =2 [B, B ]. Meanwhile, we know through the duality eq. (2.9) (i.e. eq. (2.6) and its time-reversed version) that −L DWF [X, Y ] should be dual to the following fermionic theories: Repeating the procedure we applied to L DWF , we find The key step in the manipulation is that after modification, the fluctuating gauge field b only appear in the following BF terms: In both cases, one can explicitly integrate out b to implement δ(a 2 − a 1 + 2B) and δ(ã 1 − a 2 + 2B ). This allows the substitution a 1 → a + B, a 2 → a − B,ã 1 →ã − B and a 2 →ã + B . In this substitution, compactness of the gauge field is preserved, given that a 1,2 ,ã 1,2 satisfy the delta function constraints from b.

Lattice realization of N f = 2 self-duality
We will invoke the above procedure on the lattice partition functions to obtain the lattice Euclidean version of the duality eq. (3.1). Using the exact lattice mapping (2.23) and its time reversal version (referring to eq. (2.27)), the following three partition functions are exactly equivalent: We remind that Z DF is obtained from Z DW F by the exactly mapping ; on the other hand, to obtainZ DF , we first replace a 1 = −ã 1 +b 1 +X, a 2 = −ã 2 +b 2 +Y , and then use the same exact mappings that integrate out b 1,2 . The Wilson fermions with mass M are those that correspond the light Dirac fermions in L DF andL DF in the continuum, while those with mass M corresponds to the monopole operators in the Chern-Simons terms. Note that in the duality between JHEP06(2019)038 Z DF andZ DF , the roles of light fermions and monopole operators (heavy fermions) are switched; moreover, the time reversal operation leaving Z DW F invariant maps Z DF and Z DF to each other. Now, one can do a substitution X → B f + B, Y → B f + B , and include a factor Z BF [A, B + B , B f ; M ] (see eq. (2.39) for definition) that plays the role of the BF term. (Recall that we included the background spin-c connection A to make heavy fermions satisfy spin-charge relations; its role is purely auxiliary. We just set A = 0.) As a final step, turn B f to the fluctuating gauge field b. The string of procedures we described is simply a lattice version of the manipulation we described earlier to obtain N f = 2 selfduality from two copies of the boson-fermion duality described. However, we emphasize that the fact that continuum field theory dualities in this paper are IR dualities implies that promoting the gauge field b from the classical one to the quantum one is subtle, the aforementioned procedure in the lattice version is exact.
Following the manipulation, Z DW F is transformed into the following partition function whose action corresponds to the lattice version of −L N b =2 of eq. (3.1): The same manipulation transforms Z DF andZ DF into where we have defined (3.8)

JHEP06(2019)038
The role of Z QF [A, a 1 , a 2 , B, B ; M ] andZ QF [A, a 1 , a 2 , B, B ; M ] is that, upon tuning M properly and interpreting each heavy fermion as implementing level ±1 Chern-Simons term, we expect The key consistency check for the above proposal for the IR limit is whether delta-functions are implemented properly, presumably due to condensing of gauge-invariant composite objects and the subsequent Higgsing of gauge field coupled to the composite objects. This scenario is analogous to the one we studied in section 2.4; we will see that a similar condensation of composite objects occur in the theoryZ QF and Z QF occur through the sign-free Monte-Carlo simulation; we present the details in section 3.3.

Emergence of SU(2) × SU(2) symmetry
One remarkable feature of the duality eq. (3.1) that lead to its application in deconfined quantum criticality [41] is the emergence of a global SU(2) × SU (2)  Thus, at m 2 = −m 1 , the theory −L N f =2 has a manifest SU(2) ψ symmetry acting on the fermion doublet (ψ 1 , ψ 2 ) of mass m 2 , and similarly the theory −L N f =2 has a manifest SU(2) χ symmetry acting on (χ 1 , χ 2 ) of mass −m 2 . Recall that in the duality between −L N f =2 and −L N f =2 , the light fermion in one theory plays the role of a monopole operator (heavy fermion) in the other, so these two SU(2) symmetries are independent. Thus, we have a global SU(2) ψ × SU(2) χ symmetry in total. This means the boson theory −L N b =2 also has this symmetry, though not manifest. One may also note that the background fields B and B can be understood as the σ 3 gauge fields of the SU(2) ψ and SU(2) χ symmetry respectively.
We recall that the duality between −L N f =2 and −L N f =2 can be regarded as a selfduality when m 1 = m 2 . This means at the m 1 = m 2 = 0 critical point, the action of the SU(2) ψ × SU(2) χ symmetry is self-dual. Indeed, as commented above, when m 2 = −m 1 , the two SU(2)'s act on fermion doublets of mass ±m 2 respectively, which becomes the same only at zero mass.
In any UV completion, there shall be two UV parameters that play the roles of m 1 and m 2 in the IR. Then one might hope the UV completion could be such that, if we tune these two UV parameters to satisfy one relation (corresponding to m 1 = m 2 in the IR) the UV theory would have exact self-duality, and if we tune them to satisfy another relation (corresponding to m 1 = m 2 in the IR) the UV theory would have an exact SU(2) ψ ×SU(2) χ symmetry. Unfortunately, this would be hard to simultaneously realize on the lattice, JHEP06(2019)038 because in the above, we have used the relation −L −,m Dirac [a] = −L +,m Dirac [a] − i 4π ada in the continuum, which becomes non-exact on the lattice as we have seen in section 2.5. Instead, in our construction in section 3, the self-dual property is exact at all scales if we set T 1 = T 2 ; on the other hand, the SU(2) ψ × SU(2) χ symmetry emerges in the IR if we make T 1,2 so that their IR mass ∼ M 1,2 − 3 become opposite. The two properties are compatible at the critical point T 1 = T 2 = T c .

Numerical confirmation of the implementation of "Lagrange multiplier"
In this section, we provide the numerical evidence for condensation of the gauge-invariant four-fermion operators in the theoryZ QF [A, a 1 , a 2 , B, B ; M ]. This four-fermion operator is charged under a 1 − a 2 + 2B ; condensation of this operator indicates that Z QF [A, a 1 , a 2 , B, B ; M ] implements δ(a 1 − a 2 + 2B ), at least for IR modes of the gauge fields. We emphasize thatZ QF does not represent N f = 2 dualities and as in section 2.4, our primary focus in nuemrics is to show emergence of U(1)-breaking gapped phase across some parameters rather than investigating nature of a possible critical point.
First, we observe that the theoryZ QF [A, a 1 , a 2 , B, B ; M ] can be simulated without sign problem when a 1 = a 2 = B = B = 0, similar to what we did for Z DF investigated in section 2.4, since fermion determinants once again come in complex conjugate pairs. To see this, we observe thatZ QF [A, a 1 , a 2 , B, B ; M ] can be written as: φ n,0 = 1 √ 2 (ψ n,↑ χ n,↑ + ψ n,↓ χ n,↓ ),φ n,0 = 1 √ 2 (ψ n,↑χn,↑ +ψ n,↓χn,↓ ) φ n,0 = 1 √ 2 (ψ n,↑ χ n,↑ + ψ n,↓ χ n,↓ ),φ n,0 = 1 √ 2 (ψ n,↑χ n,↑ +ψ n,↓χ n,↓ ). (3.11) One can see that the four fermion composite operator φ n,0 φ n,0 is coupled to the gauge field a 2 − a 1 + 2B . While we do not have analytical intuition why this object should condense from expanding Grassmann variables as in section 2.4, we will present numerical results consistent with long-range order associated with φ n,0 φ n,0 , and thusZ QF implements the delta function we desire. Another important consistency check for whetherZ QF implements the correct delta functions is that φ n,0 and φ n,0 should not condense individually; if φ n,0 and φ n,0 condenses as well,Z QF implements two delta functions instead of one and allows  one to integrate out one more gauge fields in IR, which is not an expected behavior in our N f = 2 construction. The most direct probe of long-range orders associated with φ n,0 φ n,0 and φ n,0 is the spatial dependence of the following correlators: Certain observables in this Monte-Carlo simulation converge much slowly compared to the simulation presented in section 2.4, due to the more complicated dynamics. Particularly, we found that it is very difficult to get reasonable estimate for the C 4 (n 1 , n 2 ). However, "pair susceptibilities" written below can be computed in reasonable amount of time and hence will be used to probe long-range orders: Also, we found that adding a Maxwell term for the gauge field b improves the speed of convergence, though its effects on the numerical results are minor; the effect being rather minor is expected, because the Maxwell term is irrelevant. 10 We computed χ 2F and χ 4F for M = 0.8−2.0, setting the inverse Maxwell term strength g 2 = 100, for L × L × L cubes with L = 8, 9, 10; the computed values are plotted in figure 4. On the right panel of figure 4, one see a strong peak in χ 2F centered around M ≈ 1.8; a much weaker peak, more manifest in larger system sizes, can be observed around M ≈ 1.1.
The computed values of χ 2F do not jump otherwise, suggesting that except at those two putative critical points, φ n,0 remains gapped. Meanwhile, as for computed χ 4F 's plotted on the right panel, there is a gradual increase of χ 4F 's as one tunes from M = 0.8 to M = 1.4, followed by a sharp decrease around M ≈ 1.8. This also suggest the phase JHEP06(2019)038 transition around M ≈ 1.1 and M ≈ 1.8; particularly, the sharp drop at M ≈ 1.8 in χ 4F suggests that the phase transition is associated with condensation/localization of φ n,0 φ n,0 .
The finite-size scaling of χ 2F and χ 4F reveals more information about the nature of the phase. First, let us comment on how χ 4F and χ 2F is expected to scale as one increases the volume. We note that in the thermodynamic limit, one may approximate the sum in eq. (3.13) by integrals. Then, as one approaches thermodynamic limit: The above integral allows one to infer correlation function behaviors as the following: 1. When the correlator decays exponentially or with power law with exponent larger than 3, the integral converges to constant value at thermodynamic limit. Thus, in this case, the corresponding "susceptibility' should scale like 1 V .
2. If the correlator has long-range order, χ 2F and χ 4F should converge to a constant in thermodynamic limit.
3. If the correlator decays with power law with exponent α < 3, by similarly using integral one can argue that the corresponding susceptibility should scale like 1 V α/3 . As for χ 2F , at M = 1.4 − 1.7 and M = 1.9, 2.0, the scenario 1 prevails. Strictly speaking, it is also possible that the correlator decays in power-law with large exponent, yet we conjecture φ n,0 is truly gapped in this regime. Meanwhile, at M = 0.8, we cannot observe well-defined scaling behavior of χ 2F -presumably, one needs to study larger systems to extract the correct scaling behavior at this mass value. For M values near M ≈ 1.1 and M ≈ 1.8, scenario 2 prevails, suggesting second-order critical points around this Mass values.
As for the behavior of χ 4F , at M = 1.9, 2.0, scenario 1 prevails once more, indicating φ n,0 φ n,0 is also gapped in this parameter. Hence, we conjecture that the phase in M > 1.8 is a trivial gapped phase adiabatically connected to M = ∞ without phase transition. At all other values, scenario 2 or 3 prevails, indicating that C 4F is at least power-law decaying in the rest of the M -value. Our particular interest is in M value between 1.1 and 1.8, the range we conjecture for φ n,0 φ n,0 condense. However, it is not entirely clear whether scenario 2 or 3 is valid, due to the size of error bar.
To study finite-size scaling behavior of χ 2F and χ 4F of the phase sandwiched between two critical points at M ≈ 1.1 and M ≈ 1.8 more carefully, we computed χ 2F and χ 4F at M = 1.4, g 2 = 100 for L = 11 and plotted them along with the values for L = 8, 9, 10 in figure 5. One can see clearly on the left panel that χ 2F almost precisely scales like 1/V , confirming our conjecture that scenario 1 prevails for χ 2F at this M value. The right panel confirms that the computed χ 4F values are fundamentally inconsistent with scenario 1. The red dotted line seems to fit well with data when the value from L = 8 is excluded, which seems to indicate scenario 3 prevails. However, we note that error bars is are still large here,and while it is clear that φ n,0 φ n,0 is not gapped, the true nature of scaling behavior remains inconclusive. We believe it is possible that the slow decrease in JHEP06 (2019)  χ 4F as one increases the system size is non-universal finite-size effect that should go away in the thermodynamic limit, and eventually χ 4F converges to a constant. This conjecture should be verified more carefully in future works.

Non-relativistic limit
Soon after the original boson-vortex duality was introduced, non-relativistic generalizations of the duality transformation were developed [8,42]. Such generalizations played an important role in the condensed matter realizations of such dualities, enabling, for instance, a more global understanding of the superfluid-Mott insulator quantum phase transitions [43].
On the other hand, the non-relativistic limit of the bosonization and other dualities discussed in this paper can be related to the traditional practice of flux attachment, which is also an important method in condensed matter physics, especially in the quantum Hall effect [26,28]. Let's start with the non-relativistic limit of the free Dirac theory in the continuum, that appeared in eq. (1.1). The non-relativistic limit occurs at finite m > 0, and the fermions are under a chemical potential η such that η m; the range of energy scale of interest is from some µ ∼ m to µ + ∆µ, where ∆µ ∼ η − m m. The chemical potential makes appearance in the Euclidean signature such that A = (A x , A y , A z ) in the first term of −L +,m Dirac [A] becomes (A x , A y , A z + iη); in the Lorentzian signature, iA z → A t , so η appears in the Dirac Lagrangian in the form A t − η as expected. We emphasize that this shift of A z by iη is not made in the parity anomaly term; this should be taken as part of what we mean by the improperly quantized level-1/2 Chern-Simons term.
One may then ask how the chemical potential η makes its appearance on the boson side of eq. (1.1). There are two ways to interpret η on the boson side. To present the simpler interpretation, let's first shift b in eq. (1.1) by A:

JHEP06(2019)038
and then replace A z with A z + iη (Euclidean signature) or A t with A t − η (Lorentzian signature) as in the fermion theory. This means in the bosonic theory, the boson charge density is subjected to a chemical potential η, as is the fermionic charge density in the fermion theory. Alternatively, we may interpret η in the bosonic theory without shifting b to b . We directly shift A z to A z + iη in the Chern-Simons term in eq. (1.1). While the background dA flux can be assumed exact for simplicity, the dynamical db flux might not be, thereby yielding a term −η dz dxdy(db) xy /(2π), where dxdy(db) xy /(2π) is an integer, the total monopole charge over the x, y-directions. 11 This term can be understood as a chemical potential for the monopole charge density. The two interpretations of η in the bosonic theory are equivalent because the duality states that the boson φ and the db monopole form a composite object, which is the fermion ψ in the fermionic theory, therefore we can attribute the chemical potential either to the boson or to the monopole. The two ways of interpreting the chemical potential is obvious in our lattice construc- , whereÃ is to replace A nz on the nz links with A nz + iη. Then, our lattice construction for the bosonization duality [34] may be easily generalized to the mapping between the two partition functions with chemical potential: Here η 1 may be viewed as the chemical potential for the boson, η 2 that for the monopole (heavy Dirac fermion), and in the dual theory they together contribute to the chemical potential for the light Dirac fermion, which is the composite object of the boson and the heavy Dirac fermion.
In what follows, we make another comment about the non-relativistic limit, in connection to the traditional flux attachment. In the traditional non-relativistic flux attachment, a boson may be mapped to any odd number of fluxes attached to a fermion, or any even number of fluxes attached to a boson. On the other hand, in the relativistic duality, say eq. (2.9), a relativistic boson may only be mapped to either +1 or −1 flux attached to a relativistic fermion, but not the versions with more fluxes. 12 Below, we explain why the relativistic duality is so limited compared to the non-relativistic flux attachement; more exactly, we explain why, as we take the non-relativistic limit of the relativistic theories, extra approximate dualities appear. We note that a similar issue appeared in the context of loop models as well [23,44]. 11 One may wonder why the denominator is 2π though naively it seems to be 4π. This has to do with the fact that the expression i 4π bdb is incomplete and cannot be taken literally when non-exact monopole charge is present. While we will not introduce the complete expression, one way to convince oneself about the coefficient is that db/2π is the current coupled to A, so its density component should be (db)xy/2π. 12 If we carry out further S, T transformations as introduced in section 2.6, we will only get relativistic dualities like "a boson with six fluxes attached is dual to a fermion with five fluxes attached", unlike in the non-relativistic case where a boson with six fluxes attached is just another boson.

JHEP06(2019)038
In non-relativistic flux attachment, one may view the attachment of multiple fluxes as the following sequence of operations. A boson can be realized as a fermion with +1 flux attached. The fermion, in turn, may be realized as another boson with +1 flux attached. This new boson, in turn, may be viewed as another fermion with +1 flux attached, and so on. (Similarly if we replace +1 flux by −1.) Therefore the original boson can be realized as any odd / even number of fluxes attached to a fermion / boson. How is the situation different in the relativistic dualities? The right half of eq. (2.9) states that a boson can indeed be realized as a fermion with +1 flux attached. However, it is crucial that this fermion is described by −L −,−m Dirac which has a −1/2 parity anomaly term. The fermion described by −L −,−m Dirac may only be viewed as the original boson with −1 flux attached (i.e. the time-reversed version of eq. (2.3)), but cannot be viewed as another boson with +1 flux attached (which would be −L +,+m Dirac with +1/2 parity anomaly term), so the sequence of operations stops there as opposed to the non-relativistic case. In other words, the difference compared to the non-relativistic case lies in that, in relativistic theory, −L −,−m Dirac and −L +,+m Dirac are distinct. 13 We expect as we approach the non-relativistic limit in the relativistic theory, the distinction between −L −,−m Dirac and −L +,+m Dirac diminishes and then we may, at better and better approximation, carry out the non-relativistic procedure of attaching multiple fluxes. The diminishing of the distinction can be easily understood. Consider −L +,+m Dirac with finite m > 0 in the non-relativistic limit. The Dirac sea (lower band), along with the parity anomaly term (which, in the lattice UV completion, arises from other "doubler" modes in the Dirac sea), contributes a background AdA Chern-Simons at level −sgn(m)/2 + 1/2 = 0, i.e. their net contribution is trivial. The upper band, on the other hand, contains modes with Berry curvature + 1 2 |m|/E 3 in the (p x , p y ) momentum space, where E = p 2 x + p 2 y + m 2 ; if the upper band has a Fermi surface at chemical potential η m, integrating the Berry curvature over the Fermi sea leads to a Hall conductivilty of (1 − |m|/η)/2. The case of −L −,−m Dirac is similar: the Dirac sea and the parity anomaly term together makes a trivial contribution, while the upper band has Berry curvature − 1 2 |m|/E 3 , and, if at chemical potential −η −|m|, Hall conductivity −(1 − |m|/η)/2. In the non-relativistic limit we are interested in a narrow range of energy scale ∆µ ∼ |η − m| |m| near the bottom of the upper band, around which the difference in the Berry curvature or Hall conductivity between −L +,+m Dirac and −L −,−m Dirac are suppressed by powers of ∆µ/|m|. Therefore, as ∆µ becomes small, one may approximately replace −L −,−m Dirac with −L +,+m Dirac , and thus the attachment of multiple fluxes can be approximately carried out, leading to the sequence of operations in the traditional non-relativistic flux attachment.

Conclusion
In this work, we extended the lattice field theory construction of the U(1) bosonization duality to other dualities in the duality web. First, we observed that S and T transformations can be implemented by adding heavy Wilson fermions and promoting additionally 13 One may wonder what if we trade the −L −,−m Dirac at hand with −L +,−m Dirac at the cost of an extra level-1 Chern-Simons. We have already discussed this manipulation -it leads to the boson-vortex duality.

JHEP06(2019)038
introduced background gauge fields to the quantum fluctuating ones. This allows one to construct the lattice version of the boson-fermion duality and Son's fermion-fermion duality.
The cost of formulating the aforementioned dualities as exact mappings between the lattice field theories is that integrating out gauge fields cannot be done exactly anymore; so one may doubt whether naively implementing modular transformations from heavy Wilson fermions does really make the lattice field theory flow to the fixed point we imagine. Through numerics, we provided evidence that integrating out gauge fields can be valid operations in the IR, evidenced from condensation of properly charged objects. In addition, we could recover non-trivial time-reversal properties of boson-fermion dualities from the lattice construction as well. We pursued a similar idea in constructing the exact lattice mapping of N f = 2 self-dual QED. The numerics shows some evidence that condensation of 4-fermion objects induce the correct delta-function structure needed to implement N f = 2 self-dual QED as well; however, more extensive numerical simulations are required to fully justify this scenario.
Caveats regarding the UV regularization of the S transformation defined in eq. (1.2) are made manifest through the microscopic construction. We particularly emphasize that the seemly innocent relation, SS −1 = 1, does not automatically hold once a Maxwell UV regularization is imposed. For the relation to work in the context of generating the duality web, one must make the additional assumption (which is made implicitly in the literature) that changing the relative strengths between two gauge couplings involved in S and S −1 over a finite range does not alter the IR physics.

JHEP06(2019)038
χ, ψ T and χ T are defined likewise. N (a, B, M 1 ) is a 2n × 2n matrix that characterize the lattice action. Note that at B = b = 0 and M 1 = M 2 , the matrix associated with the part of the action with χ fermions and the one associated with ψ fermions are complex conjugate of each other. Thus, each gauge field configuration effectively has a positive weight | det N (a, 0, M )| 2 , and sum of all these weights define the partition function. Also, at B = b = 0 and M 1 = M 2 , four-point correlation functions as the following: In the above equation, α 1 , α 2 , β 1 , β 2 denote indices for Grassman variables. The message we want to deliver from the above equation is that correlation functions can be computed via Monte-Carlo simulation by sampling gauge field configurations so that each gauge field configuration a has a relative probability distribution | det N (a, 0, M )| 2 and taking average of proper matrix elements, in this case N −1 α 1 α 2 (a, 0, M )N * −1 β 1 β 2 (a, 0, M ). One can easily generalize this to higher-point correlation functions -the matrix elements to be taken average can be readily determined by Wick's theorem. One can similarly proceed forZ QF to see that one can simulateZ QF simply by modifying the weight for each gauge field configuration a to be | det N (a, 0, M )| 4 instead of | det N (a, 0, M )| 2 .
As for pratical implementations, one can employ standard Metropolis algorithm to sample gauge field configurations. For each step of the algorithm, the program proposes a gauge field configuration to be updated from {a} to {a } (For the purposes in this paper, one can safely restrict the update to be local; that is, the gauge field configuration This procedure allows one to sample the gauge field configurations whose probability distribution is proportional to | det N (a, 0, M )| 2 . As for simulatingZ QF , one can simply change the Metropolis weight | det N (a ,0,M )| 2 | det N (a,0,M )| 2 to | det N (a ,0,M )| 4 | det N (a,0,M )| 4 . We briefly comment on some practical issues on the Monte-Carlo simulation. At each step of the simulation, one may store gauge field configuration {a} and inverse matrix N −1 (a, 0, M ). This allows one to readily extract matrix elements needed for computing correlation functions from the inverse matrix N −1 (a, 0, M ). Also, storing the inverse matrix in addition to gauge field configuration allows one to the Metropolis probability | det N (a ,0,M )| 2 | det N (a,0,M )| 2 in O(1) time. To see this, for simplicity, let us assume that we consider an update from a to a , and a and a have a different gauge field value only at one link. One can show that thanks to sparsity of the matrix ∆(a , a, M ), ∆(a , a, M )N −1 (a, 0, M )) is a matrix whose non-zero elements lie on only four columns. Also, by performing column operations, one can show that sixteen non-zero elements of ∆(a − a, M )N −1 (a, 0, M )) that form 4 × 4 block diagonals are the only ones that contribute to non-trivial Metropolis probability; all other non-zero elements can be ignored. Thus, while the Metropolis algorithm involves ratio of determinants of large matrices which are numerically costly to compute, storing the inverse matrix and restricting to local updates allow transition probability to be computed in O(1) time.
While computing the Metropolis probability can be implemented in a numerically inexpensive way, updating the matrix inverse N −1 (a, 0, M ) to N −1 (a , 0, M ) when the proposed change is accepted is difficult to be done efficiently and hence act as the bottleneck of this Monte-Carlo simulation. When restricted to local updates, N −1 (a , 0, M ) can be computed from N −1 (a, 0, M ) using Shermann-Morrison formula in O(V 2 ) complexity, V denoting volume of the system we are simulating; this is the most efficient way to update the matrix inverse from our knowledge. While this is numerically costly, this is much faster than computing N −1 (a , 0, M ) from scratch with complexity O(V 3 ).
Finally, recall that when investigating spatial dependence of correlation functions, we made one direction(say z direction) longer than other, forming a "prism" geometry, and investigated how correlators behave on z axis. We note that for this purposes, it is more efficient to compute the following quantities: C(n z ) = 1 L x L y nx,ny φ (nx,ny,0),0 φ (nx,ny,nz−),0 (A.6) C(n z ) is the "spatially-averaged version" of the correlation function we intend to study. Our action is translationally invariant, so if one can average over all gauge configurations, C(n z ) = φ (nx,ny,0),0 φ (nx,ny,nz−),0 for any n x and n y . However, in the Monte-Carlo simulation, we do random sampling of gauge fields, and choosing specific gauge field configurations break translational invariance. The above quantity takes average spatially within the same configuration and the system as well, allowing one to obtain more fastly convergent results.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.