Duality between the deconfined quantum-critical point and the bosonic topological transition

Recently significant progress has been made in $(2+1)$-dimensional conformal field theories without supersymmetry. In particular, it was realized that different Lagrangians may be related by hidden dualities, i.e., seemingly different field theories may actually be identical in the infrared limit. Among all the proposed dualities, one has attracted particular interest in the field of strongly-correlated quantum-matter systems: the one relating the easy-plane noncompact CP$^1$ model (NCCP$^1$) and noncompact quantum electrodynamics (QED) with two flavors ($N = 2$) of massless two-component Dirac fermions. The easy-plane NCCP$^1$ model is the field theory of the putative deconfined quantum-critical point separating a planar (XY) antiferromagnet and a dimerized (valence-bond solid) ground state, while $N=2$ noncompact QED is the theory for the transition between a bosonic symmetry-protected topological phase and a trivial Mott insulator. In this work we present strong numerical support for the proposed duality. We realize the $N=2$ noncompact QED at a critical point of an interacting fermion model on the bilayer honeycomb lattice and study it using determinant quantum Monte Carlo (QMC) simulations. Using stochastic series expansion QMC, we study a planar version of the $S=1/2$ $J$-$Q$ spin Hamiltonian (a quantum XY-model with additional multi-spin couplings) and show that it hosts a continuous transition between the XY magnet and the valence-bond solid. The duality between the two systems, following from a mapping of their phase diagrams extending from their respective critical points, is supported by the good agreement between the critical exponents according to the proposed duality relationships.


I. INTRODUCTION
A duality in physics is an equivalence of different mathematical descriptions of a system or a state of matter, established through a mapping by change of variables. The simplest example is the particle-wave duality in quantum mechanics, where the duality transformation is a change of basis by Fourier transformation, and the chosen basis dictates the variables used to describe the system. In classical statistical mechanics, the most famous duality is the Kramers-Wannier duality of the two-dimensional (2D) Ising model [1]. Here the low-and high-temperature expansions of the partition function can be related to each other by identifying a one-to-one correspondence between the terms in the two different series, thus establishing an exact mapping between the ordered and disordered phases and the corresponding collective variables. In this case the critical point is also a self-duality point.
In the 3D Ising model, one can instead find a different model whose high-temperature expansion stands in a direct one-to-one correspondence with the low-temperature expansion of the the Ising gauge model [2,3]. Many other examples of dualities have been established, e.g., the wellknown equivalence between the 3D O(2) Wilson-Fisher fixed point and the 3D Higgs transition with a noncompact U(1) gauge field [4][5][6] [7].
In analogy with the Ising examples mentioned above, it is some times possible to transform a quantum field theory at strong coupling into an equivalent dual theory at weak coupling. The untractable original problem can then be solved by means of perturbative methods applied to the dual theory. Such strong-weak duality (and the more general "S-duality" form) was established in certain supersymmetric Yang-Mills theories [8][9][10][11] and Abelian gauge theory without supersymmetry [12][13][14]. In 1D quantum systems (i.e., in 1 + 1 space-time dimensions) a well known fermion-boson duality is achieved by bosonization of an interacting fermion system through a non-local transformation [15][16][17][18][19]. Usually, in the bosonized formalism interactions can be more easily treated than in the original fermion model. In cases where no formal mapping is known, two Lagrangians that look different in the ultraviolet may still flow (under the renormalization group) to the same theory in the infrared, i.e., these seemingly different field theories actually represent exactly the same low-energy physics. Such a duality goes a step beyond the more familiar concept of universality, by which systems (models or real materials) with the same dimensionality and global symmetries exhibit identical scaling properties at their classical or quantum critical points. Such systems share the same effective critical low-energy field-theory arXiv:1705.10670v1 [cond-mat.str-el] 30 May 2017 description. For example, the critical points of the Bose-Hubbard model and the quantum rotor model are in the same universality class. A duality transformation usually changes the description of the system into a form based on nonlocal objects or defects of the original description. On a practical level, the existence of a dual field theory means that there is a non-trivial choice of which description to use, and one of them (and not necessarily the originally most obvious one) may pose a more tractable setup for calculations.
Even if no strong-weak transformation exists (or is known), a difficult or non-tractable strong-coupling problem can some times be shown to be dual to a different strong-coupling problem that is tractable with some specific computational technique. In particular, the dual problem may be more easily solvable using powerful numerical (lattice) methods. In this paper we will explore such a recently proposed duality between two different (2 + 1)D Lagrangians that respectively involve fermionic and bosonic matter fields coupled with a gauge field [20,21]. Both theories are of great current interest in the context of strongly correlated electrons in two dimensions. Our aim here is to identify a duality between the systems by establishing corresponding lattice models realizing the two low-energy theories.
We follow the recent proposal that (2 + 1)D quantum electrodynamics (QED) with noncompact gauge field and two flavors of Dirac fermions is dual to the critical point of the easy-plane NCCP 1 model (the bosonic QED with two flavors of complex bosons) [20,21]. On the lattice, we realize the former with an interacting fermion model with spin-orbit coupling on the bilayer honeycomb lattice (BH), which hosts a quantum phase transition between a (bosonic) symmetry-protected topological state and a trivial Mott insulator [22][23][24]. It was proposed that this transition is described by N = 2 noncompact QED [25,26]. To realize the low-energy physics of the NCCP 1 theory, in this paper we introduce a planar variant of the spin S = 1/2 J-Q Hamiltonian (a Heisenberg model with additional multi-spin couplings [27]), dubbed the easy-plane J-Q (EPJQ) model, and show that it hosts a deconfined quantum-critical point [28][29][30] separating antiferromagnetic (AFM) and dimerized (valence-bondsolid, VBS) ground states (similar to the case of the J-Q model with full spin SU(2) symmetry [31], but with different universality due to the lowered symmetry). Our numerical (quantum Monte Carlo, QMC) results establishes the critical-point universality and duality between the phase diagrams of the two models. With the EPJQ model being much easier to study on large scales with QMC simulations than the fermionic model, the duality that we establish here allows for detailed studies of the topological transition of the latter, through the analogue of the deconfined quantum-critical point. The phase diagrams and dualities studied are summarized in Fig. 1.
The rest of the paper is organized as follows: In Sec. II we present the details of the two field theories and their putative duality, and in Sec. III we define the lattice mod-els and the proposed mappings relating their phase diagrams to each other. In Sec. IV and Sec. V we present the numerical results for the EPJQ and BH models, respectively. We conclude with a brief summary and discussion of the results in Sec. VI. In Appendix A we present further technical details on the analysis of the critical exponents, and in Appendix B we compare results for different variants of the EPJQ model (with different degrees of spin-anisotropy).

II. CONTINUUM FIELD THEORIES
The bosonic particle-vortex duality mentioned in the introduction was recently generalized to a model with fermionic matter [32][33][34][35][36][37][38][39][40], in the form of a (2 + 1)D QED Lagrangian with a single flavor of a two component Dirac fermion and noncompact gauge field, i.e., N = 1 QED. This theory is dual to that of a noninteracting Dirac fermion in the infrared limit [41]. Based on this N = 1 duality, Ref. [42] showed that (2 + 1)D QED with noncompact gauge field and N = 2 flavors of Dirac fermions is self-dual. This is also a fermionic version of the selfduality of the easy-plane NCCP 1 model (which can be regarded as N = 2 bosonic QED) [28,29,43]. The selfduality of the N = 2 QED Lagrangian was also verified with different derivations [40,44,45]. Unlike the case of N = 1, there is no equivalent noninteracting description of (2 + 1)D QED with N = 2, however. Because of its self-duality, N = 2 QED hosts an (emergent) O(4) symmetry in the infrared, which factorizes into the two independent SU(2) flavor symmetries on the two sides of the self-dual point.
More recently, based on the previous fermion-boson duality [37], it was argued that N = 2 QED is also dual to the easy-plane NCCP 1 model at the critical point [20,21]. These two field theories can be written as where ψ and z are two-component Dirac fermion and complex boson fields coupled to non-compact U(1) gauge fields, a and b, respectively. The duality maps the variables (m, M ) to (h, r). Moreover, both theories in Eq. (1) are individually self-dual. The putative duality between the two theories implies that the easy-plane NCCP 1 model should also have an emergent O(4) symmetry at its critical point, which is not immediately obvious in Eq. (1b). The corresponding O(4) order parameter is where M b is the monopole operator (gauge flux annihilation operator) of the gauge field b. The proposed duality between the two theories in Eq. (1) leads to very strong predictions for relationships between their properties. For example, the scaling dimension of m in N = 2 QED should be precisely the same FIG. 1: Schematic phase diagrams of (a) the bilayer honeycomb (BH) model, and (b) the easy-plane J-Q 3 model. In all models, four phases joint at the deconfined quantum critical point. Phases of the same color can be mapped to each other among the models. Kane-Mele model, which drives the fermion into a quantum spin Hall state with spin Hall conductance σ sH = ±2 (depending on the sign of the spin-orbit coupling λ). With weak interaction V , the bilayer quantum spin Hall state crosses over to a bosonic symmetry protected topological (BSPT) state. However strong interlayer pairhopping interaction V evantually favors a direct product state of anti-bonding Cooper pairs. In the strong interaction limit (V → ∞), the ground state of the BH model with |0 c ⟩ being the fermion vacuum state. This state has no quantum spin Hall conductance, i.e. σ sH = 0. It could be called a trivial Mott state. As the interaction V varies from 0 to ∞, the quantum spin Hall conductance jumps from ±2 to 0 (which can not happen smoothly), so the BSPT phase and the trivial Mott phase must be separated by quantum phase transitions. It was found numerically that there is indeed a direct continuous transition at V = V c .
The low-energy bosonic fluctuations around the critical point form an O(4) vector N = (Re Σ, Im Σ, Re ∆, Im ∆), where Σ carries spin and ∆ carries charge. The BH model Eq. (??) respects the global SO(4) symmetry that rotates the O(4) vector. However if the symmetry is lowered to U(1) spin × U(1) charge , then the BSPT-Mott transition is unstable towards spontaneous symmetry breaking of the remaining symmetries. For example, the SO(4) symmetry could be explicitly lowered to U(1) spin × U(1) charge by the Hubbard-like interaction is the density of the σ =↑, ↓ fermions. The repulsive U > 0 (or attractive U < 0) interaction drives the spin ⟨Σ⟩ ̸ = 0 (or charge ⟨∆⟩ ̸ = 0) condensation, leading to the spin density wave (SDW) (or superconducting (SC)) phase that breaks the U(1) spin (or U(1) charge ) symmetry spontaneously, as illustrated in the schematic phase diagram Fig. ??(b).
The easy-plane J-Q 3 model (abbreviated as the JQ model hereinafter) is a spin model defined on a square lattice. Starting from the spin-1/2 operator S i on each site i, we define the dimmer operator D ij = S i ·S j across each bond ij, then the model Hamiltonian reads where the −(1/2)S z i S z j term introduces the easy-plane anisotropy and breaks the SU(2) spin symmetry down to U(1) spin explicitly. For small Q, the model essentially reduces to an XXZ model, which has an XY antiferromagnetic (AFM) ground state that further breaks the U(1) spin symmetry spontaneously. When Q is large, the dimmer interaction favors an valance bond solid (VBS) ground state, which breaks the lattice C 4 rotation symmetry. Previous QMC study found a direct continuous AFM-VBS transition in the SU(2) spin symmetric J-Q 3 model. Here we show that the continuous transition persists to the easy-plane case.
Near the critical point, the lattice anisotropy becomes irrelevant, such that the C 4 rotation symmetry can be enlarged to an U(1) latt symmetry at low-energy. At the critical point, the U(1) latt further merge with the U(1) spin symmetry to the emergent O(4) symmetry such that the components of the O(4) vector N = (D x , D y , S x , S y ) all become degenerated, where D µ i ≡ D i,i+μ (for µ = x, y). The AFM-VBS transition is unstable towards the axial Zeeman field −h z i S z i , which will drive the system to the spin-polarized paramagnetic (PM) phase with ⟨S z ⟩ ̸ = 0, as illustrated in the schematic phase diagram Fig. ??(c). Nevertheless, the full U(1) latt × U(1) spin symmetry is preserved in the PM phase.
The N = 2 QED fixed point has an emergent O(4) symmetry. The corresponding O(4) vector reads whereψ = (ψ 1 ,ψ 2 ) denotes the two-component dual fermion, originated from the double monopole operator in the QED theory. In the BH model (or the JQ model), the O(4) vector N corresponds to the real and imaginary part of the spin Σ and the charge ∆ operators (or the dimmer D + = D x + iD y and the XY spin S + = S x + iS y operators), as in the last line of Tab. ??.
In this work, we measure the following exponents at In (c), as was shown in Ref. [21,26,46], when tuning the two masses m and M , the N = 2 QED theory also has two symmetrybreaking (SB) phases and two symmetric (SY) phases, one of which is the BSPT state. In all models, the four phases merge at the deconfined quantum critical point. Phases of the same color can be mapped to each other among the models, following the duality relations proposed in the table on the right. The double arrows in (a) and (b) indicate the quantum phase transitions investigated numerically in this paper.
as the scaling dimension of h in the easy-plane NCCP 1 at its critical point r = 0, while the scaling dimension of M should be the same as that of r. Also, as a consequence of these dualities, i.e., the emergent O(4) symmetry of the two theories, the four components of N should all have the same scaling dimension at the critical point.
Although the duality can be observed and "derived" based on various arguments, it has not been rigorously proven yet. Both Eq. (1a) and Eq. (1b) are strongly interacting conformal field theories, and there is no obvious analytical method that can provide rigorous results for either case. However, both theories can presumably be realized using lattice models, which can be simulated using numerical methods. The goal of this work is to compare the quantitative properties of such lattice models and look for evidence of the proposed duality. As we will show in the later sections, within small error bars of the critical exponents obtained using QMC simulations, our results confirm several predictions of the duality.
Because the duality of the two theories in Eq. (1) was derived based on the assumption of the basic fermionboson duality [37,47], a proof of the former duality indirectly also proves the latter. In principle this result can lead to a number of further dualities between different fermionic and bosonic Lagrangians. Thus, the impact of our work is not limited to the proof of the duality between Eq. (1a) and Eq. (1b), but also provides justification for many other cases.

III. LATTICE MODELS
The easy-plane NCCP 1 model is the field theory that presumably describes the deconfined quantum-critical point between an in-plane (XY) antiferromagnet (AFM) and a valence-bond solid (VBS) [28,29,43]. This tran-sition in the case of full SU(2) symmetry of the Hamiltonian has been realized by the J-Q and related models, and these have been extensively simulated numerically using unbiased QMC techniques [27,31,[48][49][50][51][52][53][54][55][56][57]. Although there are studies that indicate that some version of the J-Q model with an inplane spin symmetry and other U(1) symmetric models should lead to a first order transition [58][59][60][61], in this work we demonstrate that a different model, the EPJQ model, instead leads to a continuous transition in some regions of its parameter space. The r and h terms in Eq. (1) correspond to the distance from the critical point, Q − Q c , and the staggered magnetic field h z (−) i S z i , respectively, in the lattice model. The components of the O(4) vector in Eq. (2) correspond to the two-component easy-plane Néel and twocomponent VBS (dimer) order parameters of the EPJQ model.
The N = 2 QED action has been simulated directly using a lattice QED model [62], and the scaling dimension of M was computed in this way. Also, N = 2 QED with a noncompact U(1) gauge field is the effective theory that describes the transition between the bosonic symmetry protected topological (BSPT) state and a trivial Mott state in 2D [25,26]. This transition was also realized in an interacting fermion model on a bilayer honeycomb (BH) lattice introduced in Refs. [22,23] and simulated [22][23][24]63] with a determinantal QMC method (DQMC) [64,65]. The m and M terms in Eq. (1) correspond to two different interactions in the lattice model, namely, the interlayer pair-hopping V − V c , measured with respect to its critical value, and the Hubbard-like on-site interaction U ; the Hamiltonian will be specified in detail below. The lattice model of Ref. [23] has an exact SO(4) symmetry that precisely corresponds to the proposed emergent symmetry of the N = 2 QED. It should further be noted that the fermions in the BH model do not directly correspond to the Dirac fermions of the N = 2 QED action, because the former are not coupled to any dynamical gauge field. The relation between the two systems instead arises from the correspondence of the gauge invariant fields of N = 2 QED to the lowenergy bosonic excitations of the BH model. Using the BH model and the EPJQ model, the duality between the N = 2 QED and the NCCP 1 field theories can be realized on the lattice. The O(4) vector N can also be conveniently defined in both lattice models, with explicit forms that will be explained below. Thus, the two systems can be investigated and compared via unbiased large-scale QMC simulations-the pursuit and achievement of this work.
In the following, we will first define the microscopic lattice models in detail. In the subsequent sections IV and V we present comparative numerical studies of the models and demonstrate strong support for the duality relations listed in the table in Fig. 1.

A. Bilayer honeycomb model
The BH model is a fermionic model defined on a honeycomb lattice [22][23][24]63]. On each site, we define four flavors of fermions (two layers × two spins); The Hamiltonian is where the band and interaction terms are given by where ij and ij denote nearest-neighbor intra-and inter-layer site pairs, respectively. The band Hamiltonian H band is just two copies of the Kane-Mele model [66], which drives the fermion into a quantum spin Hall state with spin Hall conductance σ sH = ±2 (depending on the sign of the spin-orbit coupling λ). Including a weak interaction V , the bilayer quantum spin Hall state automatically becomes a BSPT state [23,24,67], where only the bosonic O(4) vector N remains gapless (and protected) at the edge, while the fermionic excitations are gapped out (as will be discussed in more detail below). However, a strong interlayer pair-hopping interaction V eventually favors a direct product state of anti-bonding Cooper pairs. In the strong interaction limit (V → ∞), the ground state of the BH model reads with |0 c being the fermion vacuum state. This state has no quantum spin Hall conductance, i.e., σ sH = 0, and, more importantly, it is a direct product of local wave functions, hence dubbed trivial Mott insulator state. It was found numerically that there is a direct continuous transition between the BSPT and the trivial Mott phases at V c /t = 2.82(1) [23,63], where the single-particle excitation gap does not close but the excitation gap associated with the bosonic O(4) vector closes and the quantized spin Hall conductance changes from ±2 to 0. The low-energy bosonic fluctuations around the critical point form an O(4) vector, with N = (ReΣ, ImΣ, Re∆, Im∆) and the components are where Σ carries spin and ∆ carries charge. The BH model Eq. (4) respects the global SO(4) symmetry of the vector N . If the symmetry is lowered to U(1) spin × U(1) charge , then, based on the analysis of N = 2 QED, in principle the mass term Mψσ 3 ψ is allowed; hence the BSPT-Mott transition is unstable towards spontaneous symmetry-breaking of the remaining symmetries. The symmetry of the mass term Mψσ 3 ψ is identical to the following Hubbard-like interaction (both forming a (1, 1) representation of the SO(4)): Here ρ iσ is the density operator (for σ =↑, ↓ spins), which counts the number of σ-spin fermions in both layers on site i with respect to half-filling. The repulsive U > 0 (or attractive U < 0) interaction drives spin Σ = 0 (or charge ∆ = 0) condensation, leading to a spin-density wave (SDW) [24] (or superconducting) phase that breaks the U(1) spin [or U(1) charge ] symmetry spontaneously. This process is illustrated in the schematic phase diagram Fig. 1(a).

B. Easy-plane JQ model
Our EPJQ model is a spin-1/2 system with anisotropic antiferromagnetic couplings which we here define on the simple square lattice of L 2 sites and periodic boundary conditions. It is a "cousin" model of the previously studied SU(2) spin J-Q 3 model [49][50][51]which in turn is an extension of the original J-Q, or J-Q 2 , model [27]. Starting from the spin-1/2 operator S i on each site i, we define the singlet-projection operator on lattice link ij; then the model Hamiltonian reads where the ∆S z i S z j term for ∆ ∈ (0, 1] introduces the easyplane anisotropy that breaks the SU(2) spin symmetry down to U(1) spin explicitly. We have studied two cases: the maximally-planar case ∆ = 1 and the less extreme case ∆ = 1/2. In the latter case, we observe very good scaling behaviors indicating a continuous transition, with rapidly decaying (with the system size L) scaling corrections, while for ∆ = 1 the behavior suggests a first-order transition. Thus, the model may harbor a tricritical point separating first-order and continuous transitions somewhere between ∆ = 1/2 and ∆ = 1. However, we will leave the possible tricritical point to future investigation. As far as the duality is concerned, in this section we discuss our results for ∆ = 1/2, and in Appendix B we present results for ∆ = 1.
We set J + Q = 1 in the simulations and define the control parameter as the ratio For small q, the model essentially reduces to an XXZ model, which has an XYAFM ground state that breaks the U(1) spin symmetry spontaneously. When q is large, the dimer interaction favors a VBS (columnar-dimerized) ground state, which breaks the lattice C 4 rotation symmetry as in the SU(2) spin J-Q 2 and J-Q 3 models [27,31,[48][49][50][51], where previous QMC studies found a direct continuous AFM-VBS transition. Here we demonstrate that the continuous transition persists in EPJQ model, Eq. (11), with ∆ = 1/2. The reason for choosing the Q 3 term (three-dimer interaction) instead of the simpler Q 2 interaction (two-dimer coupling) is that it produces a more robust VBS order when the ratio q is large, thus leading to a smaller critical-point value [as in the SU (2) case [49,50]] with more clearly observable flows to the VBS state on that side of the transition. The XYAFM order in the EPJQ model can be directly probed by the local spin components S x and S y , and we will also study the critical fluctuations in the S z component. We will often not write out the staggered phase factor (−1) xi+yi corresponding to AFM order explicitly (here x i and y i refer to the integer-valued lattice coordinates of site i); in fact in the case of the XY-anisotropy in a model with bipartite interactions, the phase can also simply be transformed away with a sublattice rotation (and then the XYAFM phase maps directly onto hardcore bosons in the superfluid state). In our simulations the staggered phase is absent for the S x and S y components but present for the S z component. The AFM-VBS transition is unstable towards an axial Zeeman field, when which drives the system to the spin-polarized phase with S z = 0, as illustrated in the schematic phase diagram in Fig. 1(b). Here we consider only h = 0.
To study the columnar VBS (dimer) order realized in the EPJQ model, we define where i +x and i +ŷ denote neighbors of site i in the positive x and y-direction, respectively. At the critical point, the proposed self-duality (through the putative duality with N = 2 QED) implies that the C 4 rotation symmetry and the U(1) spin symmetry are enlarged into an emergent O(4) symmetry, such that the components of the O(4) vector (after some proper normalization) should all have the same scaling dimension [21].
C. Duality relations Fig. 1 summarizes the intuitive duality relations among the BH, EPJQ, and QED models; this can also be observed from the similarity of their four-quadrant phase diagrams [68]. To numerically prove the validity of these duality relations, in this work we investigate the following critical behaviors at the BSPT-Mott transition in the BH model: where r ij is the lattice vector separating the sites i, j, ξ is the correlation length of the critical O(4) bosonic modes of the system, and the density ρ iσ and pairing ∆ i operators have been defined in Sec. III A. We also study the following expected critical scaling behavior at the AFM-VBS transition in the EPJQ model: where ξ is the correlation length of the easy-plane spins. If the duality in Eq. (1) is correct, and provided that N = 2 QED is indeed the theory for the BSPT-Mott transition, then the exponents defined above must satisfy the following relationships [21]: 3 − 1 ν xy Here η QED is the anomalous dimension of the fermion massψσ 3 ψ, i.e., the M mass term in our notation in Eq. (1a), which was numerically estimated in the recent lattice QED calculations in Ref. [62].

IV. RESULTS FOR THE EPJQ MODEL
In this section we present our numerical results at the continuous AFM-VBS phase transition of the EPJQ model, obtained using large-scale SSE-QMC [69,70] simulations. Here we discuss only the case ∆ = 1/2 in the Hamiltonian Eq. (11); some results for ∆ = 1 are presented in Appendix B. In the SSE simulations we scaled the inverse temperature as β = 2L, corresponding to the dynamic exponent z = 1 (β ∼ L z ) and staying in the regime where the system is close to its ground state for each L. We consider L up to 44.

A. Crossing-point analysis
The first step is to determine the order of the transition and the position of the critical point (if the transition is continuous). To this end, following the recent example in Ref. [31] for the SU(2) J-Q model, we first analyze crossing points of finite-size Binder cumulants, defined for the AFM order parameter as where M 2 xy is the square of the easy-plane magnetization operator, and M 4 xy is its square. The "phenomenological renormalization" underlying the crossing-point analysis and our technical implementations of it are discussed in Appendix A. Here we show our numerical results and analyze them within the scaling relationships presented in the appendix.
As shown in Fig. 2(a), curves of U (q, L) graphed for different L cross each other at points tending to a value q c . In a finite-size system the deviation of U (q, L) from the asymptotic crossing point depends on L in a way that involves a scaling-correction exponent. For a finite-size pair (L, 2L), the crossing is at [q * c (L), U * c (L)] and at a continuous transition one expects where ν xy JQ is the correlation-length exponent and ω is the smallest subleading exponent (which normally arises from the leading irrelevant field). As shown in Fig. 2(b), an extrapolation with the above form to infinite size gives q c = 0.6197(2) (where the number in parenthesis indicates the one-standard-deviation error in the preceding digits, as obtained using numerical error propagation with normal-distributed noise added to the data points) and 1/ν xy JQ + ω = 4.0(2). The finite-size crossing value of the cumulant itself, U * c (L), should approach its thermodynamic limit U c as This form is used in Fig. 2(c) and delivers ω = 2.3(1). In principle we can now extract ν xy JQ , though the so ob-tained value and the independently determined value of ω in general should be viewed with some skepticism. The exponents are often affected by neglected higher-order scaling corrections and should be regarded as an "effective" critical exponents that flow to their correct values upon increasing the system sizes. Nevertheless, the fits with a single correction exponent are very good and this may indicate that the next-order corrections are small.
The correlation length exponent ν xy JQ can be independently and more reliably obtained from the slope of the Binder cumulant as where U 2 (L) is the derivative of U (q, L) over q evaluated at the crossing point between the L and 2L curves (which we extract by interpolating data close to the crossing point by cubic polynomials). The correlation-length exponent in the thermodynamic limit can be extracted from the expected leading finite-size form 1 ν xy, * JQ (L) As shown in Fig. 2(d), the fit to this form is statistically good if the smallest system sizes are excluded, and an extrapolation then gives ν xy JQ = 0.48 (2). Thus, the combination ν xy JQ + ω based on the independently evaluated two exponents is in remarkably good agreement with the value of the sum extracted directly using Eq. (21) with the data in Fig. 2(b). This serves as a good consistency check and again indicates that the higher-order finite-size scaling corrections should be small (i.e., the following correction exponents beyond ω must either have much larger values or the prefactors must be small, or both). Further support for this scenario can be observed in Fig. 2(d), where the data point for the smallest system size shown has a very large deviation from the good fit to the other points, suggesting a very rapidly decaying correction.
In Fig. 2(a) one may worry about the fact that the Binder cumulant forms a minimum extending to negative values as the system size increases. A negative Binder cumulant often is taken as a sign of a first-order transition. However, it is now understood that also some continuous transitions are associated with a negative Binder cumulant in the neighborhood of the critical point, reflecting non-universal anomalies in the order-parameter distribution. Often the negative peak value grows slowly, e.g., logarithmically, with the system size, instead of the much faster volume proportionality expected at a first-order phase transition. This issue is discussed with examples from classical systems in Ref. [71]. Here we do not see any evidence of a fast divergence of the peak value; thus the transition should still be continuous.

B. Anomalous dimensions
To determine the anomlous dimensions, i.e., the critical correlation-function exponents η xy JQ and η z JQ in Eq. (17), we analyze the system-size dependence of the squares of the easy-plane off-diagonal spin order parameter M 2 xy (L) in Eq. (20) and dimer order parameter where the x and y dimer operators are the appropriate Fourier transforms of Eqs. (14) corresponding to a columnar VBS; In the diagonal S z channel we study the system size dependence of the staggered magnetization All these integrated correlation functions should scale as the correlation functions in Eq. (17) with the distance |r ij | replaced by the system length L. This dimer order parameter should be governed by the same exponent η xy JQ as the off-diagonal spin order parameter if the predicted O(4) symmetry is manifested. In contrast, the diagonal magnetic order parameter is associated with a different (larger) anomalous dimension η z JQ , according to the table in Fig. 1. We have evaluated the order parameters at q = 0.620, consistent with the value of q c determined in previous subsection. Results are shown in Fig. 3 for system sizes up to L = 32 and L = 40 for the dimer and spin order parameters, respectively. In Fig. 3(a,b) we show that M 2 xy and D 2 can be fitted with the same exponent, η xy JQ = 0.13 (3), while the fit to M 2 z in Fig. 3(c) delivers a distinctively different exponent; η z JQ = 0.91(3).

V. RESULTS FOR THE BILAYER HONEYCOMB MODEL
In this section we present our numerical results on the continuous phase transition in the BH model, where an interaction-driven phase transition between a BSPT phase and a trivial Mott insulator is investigated via large-scale DQMC simulations [23,63,72] in the groundstate projector version [73]. Acting with the operator e −ΘH on a trial state (a Slater determinant) with the projection 'time' Θ large enough for converging the finite system to its ground state, we simulated linear system sizes L = 12, 15, 18, 21 and 24, with Θ = 50 for L ≤ 18, Θ = 55 for L = 21, and Θ = 60 for L = 24. The imaginary-time discretization step was ∆τ = 0.05, which is sufficiently small to not lead to any significant deviations of scaling behaviors from the ∆τ = 0 limit.

A. The continuous topological phase transition
We first present simulation results supporting a continuous BSPT-Mott transition (as was also previously discussed in Ref. [23,63]). Figure 4(a) shows the derivative of the kinetic energy density of the BH Hamiltonian Eq. (4) with respect to the control parameter V of the phase transition, Here a broad peak develops close to V c , but there is no sign of a divergence, as would be expected at a firstorder transition. Figure 4(  ground state energy density, which can be conveniently evaluated by invoking the Hellmann-Feynman theorem; For this derivative a first-order transition would lead to a sharp kink developing with increasing L (corresponding to a real or avoided level crossing). The smooth behavior supports a continuous phase transition, though of course a very weak first-order transition would produce a visible singular behavior only for larger system sizes than we consider here.
To determine the phase-transition point V c , we have further calculated the zero-frequency susceptibility of the O(4) vector (where we here take the on-site spin-singlet pairing operator), defined as where the dynamic pair-pair correlation function is defined as where k = 0 and ∆ i defined in Eq. (7). As demonstrated in Fig. 4(c), this quantity exhibits a sharp peak, as expected at a gapless critical point with power-law correlations in both space and time. The divergence is considerably slower than the ∝ space-time-volume behavior expected at a first-order transition. Due to the large computational effort needed for these DQMC simulations, we do not have a sufficient density of points close to V c to carry out a systematic analysis of the drift of the peak position, but the data nevertheless allow us to roughly estimate the convergence to the critical point V c /t = 2.82(1).

B. O(4) gap and BH correlation-length exponent
We have also extracted the excitation gaps, ∆ O(4) , corresponding to the O(4) vectors defined in Eq. (7). According to Refs. [23,63] and Eq. (31), the O(4) gap is obtained from the imaginary-time decay of the dynamical O(4) vector correlation function, and, as discussed in Sec. II and Sec. III A, O(4) bosonic modes are expected to become gapless (with power-law correlation) at the BSPT-Mott critical point. Results for ∆ O(4) as a function of V /t for system sizes L = 6, 9, 12, 15, 18 close to V c are presented in Fig. 5(a). As expected, ∆ O(4) from every system size L exhibits a dip close to V c , with the gap minimum decreasing with L as expected with an emergent gapless point at V c . In Fig. 5 To extract the correlation-length exponent ν BH , we have performed data collapse with the O(4) gap away from the critical point, as shown in Fig. 4(c). Here we focus on the regime V > V c , where we find less scaling corrections than for V < V c and an almost perfect datacollapse according to the expected quantum-critical form  Treating both ν BH and V c as free parameters, the best data collapse delivers ν BH = 0.53 (5) and V c /t = 2.80(1). This value of V c agrees quite well with the result V c = 2.82(1) estimated roughly from the susceptibility peak in Fig. 4, adding further credence to the analysis of the critical point even with the rather limited range of system sizes accessible (as compared with the EPJQ model).
Moreover, since we already determined η z JQ = 0.91(3) in the EPJQ model, we can use the duality relationship in Eq. (18a) to predict that the correlation-length exponent of the BH model should be ν BH ∼ 0.49 (2), which is fully consistent the number determined from the O(4) gap.

C. Anomalous dimensions
Finally, we study the critical equal-time correlations in the BH model. Here we have used V = 2.817 for the longest simulations. This value is within the error bars of the critical value V c = 2.82(1) and, as we will also show, there are no statistically detectable differences between data at V = 2.820 and 2.817 for the quantities studied in this subsection.
Using one of the components of the O(4) order parameter, ∆ † i ∆ j with ∆ i defined in Eq. (7), we again construct a squared order parameter. We can use the susceptibility Eq. (31) to define a corresponding equaltime spatially-integrated correlation function, where the normalization gives the same scaling behaviors as in Eq. (16) with the distance |r ij | replaced by L. The analysis illustrated in Fig. 6(a) indicates a very good power-law scaling, with deviations seen only for the smallest system size (which we exclude from the fit). The fit delivers the exponent η ∆ BH = 0.10(1), which is fully consistent with the EPJQ exponent η xy JQ = 0.10(2) obtained in Sec. IV B. Hence, the duality relation Eq. (18c) is satisfied to within the statistical errors.
In principle the anomalous dimension can also be obtained from the susceptibility in Fig. 5(c). Standard scaling arguments give that the peak height of a generic susceptibility χ should scale as χ peak ∝ L 2−η (when the dynamic exponent z = 1). We find that the peak in χ pair scales approximately as L 2 , i.e., η ∆ is very small, but here there appears to be significant scaling corrections. Moreover, there is large variation of the values for the largest size, L = 18, close to the critical point, and we would need additional points to reliably estimate the peak value. We therefore cannot obtain an independent meaningful estimate for η ∆ BH from these data. We next test the duality relation Eq. (18b). With ν xy JQ = 0.48(3) obtained in Sec. IV A we expect η ρ BH ≈ 1. Further, according to Ref. [62], the exponent η ρ BH should be equal to the anomalous mass dimension η QED of N = 2 QED, for which the value η QED = 1.0(2) was obtained from Monte Carlo simulations. Thus, we already have good consistency following from the predicted duality between η QED and ν xy JQ . Turning to the more direct test with the BH model, in Fig. 6(b) we plot the squared order parameter corresponding to the pair-density, For this quantity, as shown in Fig. 6(b), all five system sizes available give results fully consistent with a powerlaw decay, with no statistically visible scaling corrections. The fit to the five data points gives the anomalous dimension η ρ BH = 1.00(1), which is consistent with both η QED and ν xy JQ , but with a significantly smaller statistical error. Thus, the duality relation in Eq. (18b) is also confirmed to within error bars.
It is remarkable that the BH model actually seems to give better results (smaller error bars) for the anomalous dimensions than the EPJQ model, even though system lengths roughly twice as large were used for the latter. The reason is that the statistical errors on the raw data are smaller in the DQMC simulations. It is still possible that there are scaling corrections present that are not clearly visible with such a small range of system sizes, and there may then be some corrections to the exponents beyond the purely statistical error bars (one standard deviation) reported above. The EPJQ results are important in this regard as they seem to show no significant corrections even with the considerably larger range of system sizes. The lack of significant corrections is also supported by the good agreement between the non-reivial BH and EPJQ exponents, as predicted by the duality conjecture.

VI. DISCUSSION
We have performed detailed numerical tests on the recently proposed duality summarized by the predicted exponent relationships between the BH and EPJQ models in Eq. (18). The relationships were confirmed within statistical precision at the level of a few percent in the values of the critical exponents. The thus confirmed duality of the underlying low-energy quantum field theories, Eq. (1), is of great importance and interest in condensedmatter physics, because it relates two seemingly different quantum phase transitions that have been individually under intense studies during the past several years: the bosonic topological phase transition and the easy-plane deconfined quantum phase transition. The duality was derived using the more basic dualities between field theories that involve only one flavor of matter field, and sometimes also a Chern-Simons term of the dynamical gauge field.
As a consequence of confirming the particular relationship between critical exponents, our study also provides quantitative evidence for the underlying basic dualities for theories with one flavor of matter field [32][33][34][35][36][37]. These basic dualities supported by our work represent a significant step in our understanding of (2 + 1)D conformal field theories. They also form the foundation of a large number of other recently proposed dualities [40,42,44,45,[74][75][76][77][78]. Moreover, they lend support to many other dualities that follow from the same logic and reasoning, such as the duality of Majorana fermions discussed in Refs. [79,80].
To follow up on our results and insights presented here, additional numerical investigations are called for to check other predictions made within these proposed dualities. For example, in Ref. [21] it was proposed that the Gross-Neveu fixed point of the N = 2 QED is dual to the SU(2) version of the NCCP 1 model, and also has an emergent SO(5) symmetry. This symmetry has recently been discussed within SU(2) deconfined quantumcriticality as well, and quite convincing results pointing in this direction were seen in a three-dimensional loop model [57]. Scaling with the same anomalous dimension for both spin and dimer correlators had also been observed already some time ago in the SU(2) J-Q model [81]. Although we have identified the N = 2 QED as the bosonic topological phase transition in our bilayer honeycomb lattice model, we have so far been unable to find the corresponding Gross-Neveu fixed point. Identifying the additional interactions that will be required to drive this transition in the BH model is an important topic for forther research.
Following previous computational studies of deconfined quantum phase transition with SU(2) spin-rotation symmetry in J-Q models [27,31,49], we have here identified a lattice model-the EPJQ model-hosting a continuous phase transition between the U (1) (planar) Néel and VBS states. The fact that this phase transition is continuous is in itself an important discovery, given that U(1) deconfined-quantum criticality had essentially been declared non-existent, due to unexplained hints of firstorder transitions in some other planar models and what seems like definite proofs in other cases [58,60,82]. Here (as further discussed in Appendix B) we have shown that the EPJQ model defined in Eq. (11) can host first-order or continuous transitions, depending on the degree of spin-anisotropy parametrized by the Ising coupling ∆ in Eq. (11). At ∆ = 1/2, we find scaling behaviors with apparently much less influence of scaling corrections than in the SU(2) J-Q model at its deconfined critical point [31], i.e., the leading correction exponent ω is much larger in the EPJQ model. Interestingly, in both cases the correlation-length exponent is unusually small, close to 1/2, while well-studied transitions such as the O(N ) transitioins in three dimensions have exponents close to 2/3. Given its small scaling corrections and likely tricritical point between ∆ = 1/2 and ∆ = 1, the EPJQ model opens doors for future detailed studies on exotic phase transitions beyond the Landau paradigm. To determine the critical point and the critical exponents in an unbiased manner, we adopt the crossingpoint analysis applied and tested for 2D Ising and SU(2) J-Q models in Ref. [31]. Such analysis can be further traced back to Fisher's "phenomenological renormalization", which was first numerically tested with transfermatrix results for the Ising model in Ref. [83]. Ref. [31] presented systematic procedures for a statistically sound application of these techniques with Monte Carlo data. For easy reference we here summarize how we have adapted the method to the EPJQ model studied in this paper. For the BH model, due to the much larger computational cost of the DQMC simulations, we do not have data for enough system sizes to carry out the analysis in this way, and we instead applied other scaling methods in Sec. V.
Considering a generic critical point, with δ = q −q c defined as the distance to the critical point q c -for example, here q can be the control parameter q = Q/(J + Q) that we used for the EPJQ model or it could be T − T c for a finite-temperature transition. For any observable O, the standard finite-size scaling form is where we, for the sake of simplicity, only consider one irrelevant field λ and the corresponding subleading exponent ω. At the critical point, one can Taylor expand the scaling function, If one now takes two system sizes, e.g., L 1 = L and L 2 = rL (r > 1), and trace the crossing points δ * (L) of curves O(δ, L 1 ) and O(δ, L 2 ) versus δ, one finds Now if the quantity O is asymptotically size-independent (dimensionless) at the critical point, for example the Binder cumulant (which we write here with a constant and factor corresponding to a planar vector order parameter), then the corresponding exponent κ = 0, the first term in Eq. (A3) with O = U vanishes, and we obtain the following form for the size-dependent crossing point δ * (L): hence the shift of the finite-size critical point q * c (L) is approaching the asymptotic value q c as L −(1/ν+ω) .
In practice, one can interpolate within a set of points for each system size by a fitted polynomial, e.g., of cubic or quadratic order, and then use these polynomials to find the crossing points. This is illustrated in Fig 7. Error bars can be obtained by repeating the fits multiple times to data with Gaussian-distributed noise added. The scaling behavior of q c predicted by Eq. (A5) can be clearly seen in Fig. 2(b) of the main text, from which the result 1/ν+ω = 4.0(2) for the EPJQ model was obtained.
In addition to obtaining q * c and the exponent combination 1/ν + ω from the crossing points of the cumulant (or, in principle, some other dimensionless quantity), one can also use the value U * c of the quantity at q * c , as well as the derivatives at q * c , to acquire ν and ω independently. We next discuss the derivations underlying these forms.
By inserting δ * (L) into Eq. (A2), one can obtain the value of observable at the the finite-size critical point (or, more precisely, the critical point depending on the two sizes, L and rL) q * c (L). It scales as O * (L) = L −κ/ν (a + bL −ω + · · · ). (A6) Again, for a dimensionless quantity (κ = 0) such as U , the deviation of the value at the crossing point from the value in the thermodynamic limit vanishes with increasing size according to U * c (L) − U c ∝ L −ω , an example of which can be seen in Fig. 2(c) of the main text-in this case the power-law fit gave the value ω = 2.3(1) of the subleading exponent.
Here we can take the difference of the logarithms of the two equations and obtain, or, in other words, one can define a finite-size estimate of the correlation-length exponent ν * (L) from the finite-size crossing point as, It can be seen that Eq. (A11) approaches the thermodynamic limit correlation-length exponent ν at the rate 1/ν * (L) − 1/ν = gL −ω + · · · . This behavior is seen nicely in Fig. 2(d), where the extrapolation to infinite size gave ν xy JQ = 0.48 (2). In principle one can also combine the above method for a dimensionless quantity with some other quantity A, e.g., an order parameter or a long-distance correlation function. Interpolating the data for two system sizes, A(L) and A(rL) at the crossing point of the dimensionless quantity, one can take the logarithm of the ratio and analyze it in a manner similar to the slopeestimate of 1/ν, to yield a series of finite-size estimates for the power-law governing the size dependence of A. This method circumvents the need to know the criticalpoint value exactly. In practice, since q c converges fast, it is also appropriate to just use this value and analyze the size-dependence of the quantity A at this fixed coupling value q = q c , as we did in the main text.

Appendix B: Fully planar EPJQ model
In the main text we discussed the XYAFM-VBS transition at a fixed anisotropy parameter ∆ = 1/2 of the EPJQ model. We here provide some more information on the dependence on ∆.
In the extreme planar limit ∆ → 1 of Eq. (11) we have no S z S z J-term and the Hamiltonian is We have analyzed SSE-QMC results for this model in the same way as we did for ∆ = 1/2 in the main text, using a crossing-point analysis. The results of this analysis show a distinctively different behavior for ∆ = 1, pointing to a first-order transition in this case. Results for the L-dependent quantities based on the XY Binder cumulant are shown in Fig. 8. As an aside, we mention here that for the model in Eq. (B1) the off-diagonal spin correlations can be measured as diagonal correlations upon performing a basis rotation with the z spin components transformed into the x components, and vice versa. Sim- over those for generic ∆ < 1. We therefore have results for larger system sizes in this case. It is clear from Fig. 8 that we can obtain a good estimate of the critical point, but the cumulant itself and the correlation-length exponent do not exhibit the clearcut power-law corrections of the type that we saw for the ∆ = 1/2 model in Sec. IV. The fact that 1/ν * is larger than 3 suggests that the transition may actually be of first order in this case. If so, we would expect the values to eventually tend exactly to 3, and the data are consistent with this behavior.
If the transition indeed is of first order, the order parameter should be non-zero at the transition point, reflecting coexistence of the magnetic and non-magnetic phases. Indeed, as shown in Fig. 9(a), an extrapolation using a trivial 1/L correction, as expected asymptotically for a 2D system breaking a continuous symmetry, indicates a clearly non-zero value in the thermodynamic limit. The extrapolated value only changes slightly if a higher-order (1/L 2 ) correction is added (also expected if the order parameter is non-vanishing). In contrast, as shown in Fig. 9(b) the data for the ∆ = 1/2 model are fully consistent with no magnetization in the thermodynamic limit. Here the polynomial fit is strictly not correct, since a non-trivial power is expected at the critical point (which we confirmed in the main paper), but the extrapolation nevertheless indicates consistency with a vanishing order parameter at the transition in this case.
These results strongly suggest that there is a tricritical point separating a continuous and first-order transition in the EPJQ model somewhere between ∆ = 1/2 and ∆ = 1, which we plan to investigate further in a future study.