Multi-critical point and unified description of broken-symmetry phases in spin-$\frac{1}{2}$ anti-ferromagnets on a square lattice

Two-dimensional (2D) spin-$\frac{1}{2}$ antiferromagnets can support a large variety of phases, many of them break spin-rotation and/or lattice symmetries. Spin liquids[1], which break none of these symmetries, have been the focus of much recent research activity. They come in many different types as well. One such type, known as chiral spin liquid that breaks time-reversal symmetry, will be of particular relevance to our discussion below. Needless to say, quantum phase transitions among all these phases are also of strong interest.

Broken symmetry phases are traditionally described in Ginzburg-Landau theory, which is a field theory written in terms of the local order parameter associated with the spontaneously broken symmetry (see, e.g., Ref. [2]). In such descriptions, phases with different broken symmetries are described using different order parameter fields, and direct second-order transitions between them require fine-tuning. Instead, the more generic situations are first-order transitions or intermediate phases where both types of orders co-exist. In Ref. [3], Senthil and et.al. argue that such descriptions miss the possibility of deconfined criticality, which is a critical point separating two different broken symmetries facilitating a direct second-order transition between them. Such novel quantum criticality can only be captured in a field theory that describe both types of broken symmetries on equal footing. Specifically, they argue that such deconfined critical points separate the Neel ordered and valence bond solid (VBS) phases of 2D spin-$\frac{1}{2}$ antiferromagnets, which break spin-rotation and lattice translation symmetry, respectively. In the appropriate field theory, the two symmetry-breaking order parameters are dual to each other and thus afford a unified description. Numerous attempts have been made to identify such deconfined critical points, with inconclusive outcomes thus far (see Refs. [4] and [5] for recent attempts for the Heisenberg and XY symmetry classes, respectively, and references therein).

While it is not our goal to resolve the fate of deconfined criticality, our work is motivated by the line of thoughts that lead to it. To this end, we seek to find a field theory that provides a unified description of relevant phases in this description, and beyond. We find by perturbing a theory of two massless Dirac fermions coupled to a single level-one Chern-Simons (CS) gauge field, we can reach the XY-ordered (we only consider 2D spin-$\frac{1}{2}$ antiferromagnets with easy-plane anisotropy in this paper), VBS, chiral spin liquid, and an Ising Neel state in which the Neel order is along the $z$-direction despite the easy-plane anisotropy (which is possible in the presence of frustration). Within this description, a direct second-order transition between the XY-ordered phase and the VBS or Ising Neel phase must go through this massless point, which requires fine-tuning.

The remainder of the paper is organized as follows. In Sec. II, we introduce the spin-$\frac{1}{2}$ XY model and arrive at the multi-critical point by attaching a flux quantum to each hard-core boson that represents an up spin, and perform a mean-field approximation to smear out the flux. This results in two massless Dirac fermions coupled to a level-one CS field. In Sec. III, we discuss the phases that result when various mass terms are added to perturb this critical point. In Sec. IV, we discuss how the mass terms responsible for spontaneous lattice symmetry breaking are generated by fermion interactions. Section V is devoted to deriving the dual bosonic theory of the multi-critical point, where we also make comparison with the existing theory of de-confined criticality. A brief summary is offered in Sec. VI.

We consider the following spin-$\frac{1}{2}$ Hamiltonian on the square lattice:

	$\displaystyle H$	$\displaystyle=-\sum_{<ij>}(S_{i}^{x}S_{j}^{x}+S_{i}^{y}S_{j}^{y})+\cdots$		(1)
		$\displaystyle=H_{0}+\cdots,$		(2)

where $<ij>$ stands for nearest-neighbors, and the ellipsis represents generic additional couplings that respect the XY rotation symmetry and all lattice symmetries unless noted otherwise. Note that the minus sign means the XY coupling is ferromagnetic instead of anti-ferromagnetic; the two are equivalent under a $\pi$ rotation along the $z$-direction for one of the two sub-lattices. The advantage of considering the ferromagnetic XY coupling is the XY-ordered phase only breaks the $O(2)$ spin rotation symmetry, but none of the lattice symmetries. This makes the discussion of broken symmetries in various phases simpler below. The antiferromagnetic nature of Eq. 1 is thus hidden in the ellipsis, which include $S^{z}$ and further neighbor couplings between XY spins in the same sublattice.

We can map half-spin ladder operators to annihilation and creation operators of the hard-core bosons. Accordingly, the nearest-neighbor XY spin coupling in Eq. 1 becomes nearest-neighbor boson hopping,

$\displaystyle H_{0}=-\frac{1}{2}\sum_{<ij>}(b^{\dagger}_{i}b_{j}+b^{\dagger}_{j}b_{i}),$

(3)

and the ground state has half-filling in the absence of a net magnetization along the $z$-direction. We will use the spin and boson representations interchangeably below.

To proceed, we map the hard-core bosons to composite fermions (CFs) attached to a flux quantum by coupling them with pure CS theory in lattice, and then make a mean-field approximation to spread out the flux uniformly that results in a $\pi$ (or half) flux per plaquette [6]. With the gauge choice of Fig. 1, the resultant band Hamiltonian takes the form

$$h_{\bm{k}}=\left(\begin{array}[]{cc}0&\sin k_{x}+i\sin k_{y}\\ \sin k_{x}-i\sin k_{y}&0\end{array}\right),$$

(4)

in which $\bm{k}$ is the lattice momentum. Importantly, we have two Dirac points at $(0,0)$ and $(\pi,0)$ where the two bands meet. In the ground state, the lower band is filled while the upper band is empty, so the chemical potential coincides with the Dirac points. Thus, the low-energy physics of the system at this level of approximation is described by two species of massless Dirac fermions coupled to a singleCS gauge field:

$$\mathcal{L}=i\bar{\Psi}\not{D}\Psi+\mathcal{L}_{\text{CS}}[a]+\cdots,$$

(5)

where,

$$\mathcal{L}_{\text{CS}}[a]=\frac{1}{4\pi}\epsilon^{\mu\nu\lambda}a_{\mu}\partial_{\nu}a_{\lambda}=a\wedge da,$$

(6)

is the level-one CS term, $\Psi=(\psi_{1},\psi_{2})^{T}$ combines the two Dirac fields ¹¹1One for each Dirac point, and they each have two components representing the A and B sublattice. Note a rotation is performed on one of the Dirac points so that the two Dirac fields have the same chirality which allows for a unified description in Eq. 5. where $\psi_{i}$ are two-component Dirac spinors, the slash notation is defined for a general three-vector $b_{\mu}$ as $\not{b}=\gamma^{\mu}b_{\mu}$, where $\gamma^{\mu}$ are two by two Dirac matrices obeying the Clifford algebra $\{\gamma^{\mu},\gamma^{\nu}\}=2\eta^{\mu\nu}$, where $\eta^{\mu\nu}$ is the metric of the Minkowski $2+1$ space-time and $\eta^{\mu\nu}=\text{diag}(+,-,-)$, which will be used for raising and lowering the indices throughout the paper and $\{,\}$ is the anti-commutator. $D_{\mu}=\partial_{\mu}-ia_{\mu}-iA_{\mu}$ includes coupling to both the dynamic field $a_{\mu}$ and background gauge field $A_{\mu}$, and the ellipsis represents the less relevant terms like the Maxwell term of $a$. Equation (5) is the same theory discussed in Ref. 8 in a closely related context. As we demonstrate below, a variety of interesting phases supported by Eq. 1 can be accessed by perturbing Eq. 4 with various mass terms for the Dirac fermions.

The most general mass term that couples the two Dirac points takes the form $\bar{\Psi}M\Psi$, where

$$M=m_{0}\mathbbm{1}+m_{1}\sigma_{1}+m_{2}\sigma_{2}+m_{3}\sigma_{3}=m_{0}\mathbbm{1}+\bm{m}\cdot\bm{\sigma}$$

(7)

is a two-by-two Hermitian matrix. In the following, we discuss how such mass terms can be generated beyond the mean-field Hamiltonian Eq. 4, and what phases they generate once added to the critical theory Eq. 5.

From the perspective of the field theory Eq. (5), the massless point for both of the Dirac fermions (the origin in Fig. 3) is multi-critical and the full mass matrix of Eq. (7) must be tuned to zero. For example a direct second-order transition from the XY phase to the VBS phase must go through this point, while a more generic situation is going through the co-existing region or a direct first-order transition. On the other hand, the masses $\bm{m}$ break lattice symmetries. Thus, unlike $m_{0}$, they are not tuning parameters, but are instead generated from (sufficiently) strong interactions that lead to spontaneous symmetry breaking. Accordingly, we consider the following four-Fermi (Gross-Neveu type) interaction:

		$\displaystyle\mathcal{L}_{\text{int}}$	$\displaystyle=\lambda_{0}[(\overline{\psi}_{1}\psi_{1})^{2}+(\overline{\psi}_{2}{\psi}_{2})^{2}]+\lambda_{1}(\overline{\psi}_{1}{\psi}_{1})(\overline{\psi}_{2}{\psi}_{2})$		(17)
		$\displaystyle+$	$\displaystyle\lambda_{2}(\overline{\psi}_{1}{\psi}_{2})(\overline{\psi}_{2}{\psi}_{1})+\lambda_{3}(\overline{\psi}_{1}\psi_{2})^{2}+\lambda^{*}_{3}(\overline{\psi}_{2}\psi_{1})^{2}.$

It is clear that a positive $\lambda_{1}$ favours $m_{3}$, while a negative $\lambda_{2}$ and any $\lambda_{3}$ favour $m_{1,2}$. We can introduce Hubbard-Stratonovich fields $\Phi$ to decouple these interactions, resulting in a Yukawa type of coupling

$$\mathcal{L}_{Y}=\bar{\Psi}\Phi\Psi,$$

(18)

where

$$\Phi=\phi_{0}\mathbbm{1}+\phi_{1}\sigma_{1}+\phi_{2}\sigma_{2}+\phi_{3}\sigma_{3}.$$

(19)

Obviously $\bm{\phi}=(\phi_{1},\phi_{2},\phi_{3})$ is an order parameter field describing the broken lattice symmetry.

The ordering transition is described by an effective field theory in terms of $\bm{\phi}$ obtained from integrating out $\Psi$. This can be done under the generic situation of $m_{0}\neq 0$. Such a transition, if continuous, takes the system from the XY/CSL phase to a mixed phase where spontaneously broken XY/time-reversal symmetry coexist with spontaneously broken lattice symmetry. A direct continuous transition from the XY phase to the VBS phase, however, again requires fine-tuning $m_{0}$ to zero; in this case, the Dirac fermions are massless and cannot be integrated out perturbatively.

Returning to the generic situation of $m_{0}\neq 0$, to determine the order of transition this effective theory describes at the mean-field level, we are interested in the sign of the prefactor of the $|\bm{\phi}|^{4}$ term. Let us denote this prefactor as $\beta_{4}$. We can calculate it diagrammatically. For notational simplicity, we focus on the $\phi_{3}$ term in Eq. 19 as a representative of $\bm{\phi}$; the conclusions below are general. Also, we assume a uniform $|\bm{\phi}|=m$.

The Feynman rules we use are given in Fig. 4. The diagrams up to two loops that contribute to $\beta_{4}$ are shown in Fig. 5²²2There are other possible four external leg and two-loop, irreducible diagrams with one photon propagator. We find those diagrams to be zero under an appropriate regularization scheme.. The expansion in the number of loops is equivalent to a weak coupling expansion in terms of the coupling constant between fermions and CS gauge field [13], the inverse of the square root of the absolute value of the CS level. We do not show this coupling constant explicitly in our calculations for brevity of the notation. We adopt a renormalized perturbation theory approach, in which we replace the bare mass $m_{0}$ with renormalized mass $m_{r}$ in the free propagators and compensate this by adding mass and field-strength renormalization counter-terms.

Before the quantitative computation of $\beta_{4}$, we can have a qualitative discussion on what to expect. The dimension of the $\beta_{4}$ has an inverse mass dimension, i.e., $[\beta_{4}]=[m^{-1}]$. On the other hand, the only dimensionful free parameter in the theory (and relevant Feynman diagrams that generate $\beta_{4}$) is the mass $m_{r}$. We thus expect $\beta_{4}$ is inversely proportional to the mass. We will calculate this proportionality constant below.

We start with calculating the one-loop contribution $D_{1}$ which is shown in Fig. 5:

	$\displaystyle D_{1}$	$\displaystyle=-\frac{1}{4}\text{tr}\Bigg{[}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{(i(\not{p}-m_{r}))^{4}}{(p^{2}-m_{r}^{2})^{4}+i\epsilon}\Bigg{]}{}$
		$\displaystyle=-\frac{1}{4}\text{tr}\Bigg{[}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{p^{4}-4p^{2}m_{r}\not{p}+6p^{2}m^{2}-4\not{p}m_{r}^{3}+m_{r}^{4}}{(p^{2}-m_{r}^{2})^{4}+i\epsilon}\Bigg{]}{}$
		$\displaystyle=-\frac{2}{4}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{p^{4}+6p^{2}m_{r}^{2}+m_{r}^{4}}{(p^{2}-m_{r}^{2})^{4}+i\epsilon},$		(20)

which, after Wick rotating to Euclidean coordinates $p^{\mu}_{E}=p_{E\mu}=(-ip^{0},\bm{p})$ and $\not{p}=-i\not{p}_{E}$, becomes

$\displaystyle D_{1}$

$\displaystyle=-\frac{i}{2}\frac{1}{|m_{r}|}\int\frac{d^{3}p_{E}}{(2\pi)^{3}}\frac{p_{E}^{4}-6p_{E}^{2}+1}{({(p_{E})}^{2}+1)^{4}}.$

(21)

Using the integral identities given in APPENDIX: Integral identities with $n=4$, $d=3$ we find,

	$\displaystyle D_{1}$	$\displaystyle=-\frac{i}{2}\frac{\Gamma(\frac{1}{2})}{(4\pi)^{3/2}\Gamma(4)}\bigg{(}\frac{1}{|m_{r}|}\bigg{)}\Bigg{[}\frac{15}{4}-\frac{18}{4}+{\frac{3}{4}}\Bigg{]},{}$
		$\displaystyle=0.$		(22)

As a preparation for the calculation of the two-loop diagrams, we first calculate the fermion self-energy, $\Sigma$, up to leading (one-loop) order in gauge coupling,

$$-i\Sigma=-i\Sigma_{2}+D_{5},$$

(23)

where $-i\Sigma_{2}$ is shown in Fig. 6 and given as

$\displaystyle-i\Sigma_{2}$

$\displaystyle=2\pi i\int\frac{d^{3}q}{(2\pi)^{3}}\frac{\gamma_{\mu}(\not{q}-m_{r})\gamma_{\nu}}{q^{2}-m_{r}^{2}+i\epsilon}\frac{\epsilon^{\mu\nu\alpha}(p-q)_{\alpha}}{(p-q)^{2}+i\epsilon},$

(24)

and $D_{5}$ is the corresponding counter term. Equation (24) has a linear UV divergence. We can remove this by applying Pauli-Villars regularization, which is equivalent to the following substitution [14] in Eq. 24:

$$\frac{1}{(p-q)^{2}}\rightarrow\frac{1}{(p-q)^{2}}-\frac{1}{(p-q)^{2}-\Lambda^{2}},$$

(25)

where $\Lambda$ is the cutoff. Next, we use Feynman parameters to bring the denominator of Eq. 24 in a spherically symmetric form by using the identity [14]

$$\frac{1}{(q^{2}-m_{r}^{2})^{n}(p-q)^{2}}=\int_{0}^{1}dx\frac{n(1-x)^{n-1}}{[(q-xp)^{2}-\Delta]^{n+1}},$$

(26)

where $\Delta=-p^{2}x(1-x)+(1-x)m_{r}^{2}$ and $n$ is a positive integer. If we substitute Eq. 25 to Eq. 24, apply Eq. 26 with $n=1$ and change the integration variables as $q-xp\to k$, we have

	$\displaystyle-i\Sigma_{2}=i\int_{0}^{1}dx\int\frac{d^{3}k}{(2\pi)^{2}}\gamma_{\mu}(\not{k}+x\not{p}-m_{r})\gamma_{\nu}\epsilon^{\mu\nu\alpha}{}$
	$\displaystyle\times\big{(}-k+p(1-x)\big{)}_{\alpha}\Bigg{[}\frac{1}{[k^{2}-\Delta]^{2}+i\epsilon}-\frac{1}{[k^{2}-\Delta_{\Lambda}]^{2}+i\epsilon}\Bigg{]},$		(27)

where $\Delta_{\Lambda}=-p^{2}x(1-x)+(1-x)m_{r}^{2}+\Lambda^{2}x$. Next, we use the following identities for 2D gamma matrices:

	$\displaystyle\epsilon^{\mu\nu\alpha}\gamma_{\mu}\not{a}\gamma_{\nu}$	$\displaystyle=2ia^{\alpha},$		(28a)
	$\displaystyle\epsilon^{\mu\nu\alpha}\gamma_{\mu}\gamma_{\nu}$	$\displaystyle=-2i\gamma^{\alpha}.$		(28b)

The self-energy is then given as

	$\displaystyle-i\Sigma_{2}=-\int_{0}^{1}dx\int\frac{d^{3}k}{2\pi^{2}}\Big{(}-k^{2}+(1-x)(p^{2}x+m_{r}\not{p})\Big{)}{}$
	$\displaystyle\times\Bigg{[}\frac{1}{[k^{2}-\Delta]^{2}+i\epsilon}-\frac{1}{[k^{2}-\Delta_{\Lambda}]^{2}+i\epsilon}\Bigg{]},$		(29)

where we have removed the terms that are odd in $k$. Next, we perform a Wick rotation and obtain

	$\displaystyle-i\Sigma_{2}=-\frac{2i}{\pi}\int_{0}^{1}dx\int_{0}^{\infty}dk_{E}k_{E}^{2}\Big{(}k_{E}^{2}+(1-x){}$
	$\displaystyle\times(-p_{E}^{2}x-im_{r}\not{p}_{E})\Big{)}\Bigg{[}\frac{1}{[k_{E}^{2}+\Delta_{E}]^{2}}-\frac{1}{[k_{E}^{2}+\Delta_{E\Lambda}]^{2}}\Bigg{]},$		(30)

where $\Delta_{E}=p_{E}^{2}x(1-x)+(1-x)m_{r}^{2}$ and $\Delta_{E\Lambda}=p_{E}^{2}x(1-x)+(1-x)m_{r}^{2}+\Lambda^{2}x$. We can evaluate the integral over $k_{E}$ using

	$\displaystyle I_{1}$	$\displaystyle=\int_{0}^{\infty}dk_{E}k_{E}^{4}\Bigg{[}\frac{1}{[k_{E}^{2}+\Delta_{E}]^{2}}-\frac{1}{[k_{E}^{2}+\Delta_{E\Lambda}]^{2}}\Bigg{]}{},$
		$\displaystyle=\frac{3\pi}{4}\Big{(}\sqrt{\Delta_{E\Lambda}}-\sqrt{\Delta_{E}}\Big{)},$		(31)

and

	$\displaystyle I_{2}$	$\displaystyle=\int_{0}^{\infty}dk_{E}k_{E}^{2}\Bigg{[}\frac{1}{[k_{E}^{2}+\Delta_{E}]^{2}}-\frac{1}{[k_{E}^{2}+\Delta_{E\Lambda}]^{2}}\Bigg{]}{},$
		$\displaystyle=\frac{\pi}{4}\bigg{(}\frac{1}{\sqrt{\Delta_{E}}}-\frac{1}{\sqrt{\Delta_{E\Lambda}}}\bigg{)}.$		(32)

Note that $\lim_{\Lambda\to\infty}\Delta_{E\Lambda}=\Lambda^{2}x$. Finally, after undoing the Wick rotation we have

	$\displaystyle-i\Sigma_{2}$	$\displaystyle=-\frac{i}{2}\int_{0}^{1}dx\Bigg{[}-3\sqrt{\Delta}+\sqrt{x}\Lambda{}$
		$\displaystyle+\frac{1}{\sqrt{\Delta}}\Big{(}p^{2}x(1-x)+\not{p}m_{r}(1-x)\Big{)}\Bigg{]},$		(33)

where we clearly see the linear divergence. The counter-terms will remove this divergence. We define renormalization conditions as

	$\displaystyle-i\Sigma(\not{p}=-m_{r})$	$\displaystyle=0$		(34a)
	$\displaystyle-i\frac{d\Sigma}{d\not{p}}\Big{|}_{\not{p}=-m_{r}}$	$\displaystyle=0,$		(34b)

which fixes the location of the poles and the residue, thus the physical mass [14]. After substituting Eq. 23 to Eq. 34 we have,

	$\displaystyle D_{5}$	$\displaystyle=i(\not{p}\delta_{2}+\delta m_{r}){}$
		$\displaystyle=i\Big{(}-\frac{1}{2}\text{sgn}(m_{r})\not{p}+\frac{\Lambda}{3}-\frac{3}{2}|m_{r}|\Big{)}.$		(35)

Next, we calculate $D_{2}$ which is shown in Fig. 5 and explicitly given as;

$\displaystyle D_{2}$

$\displaystyle=-\frac{1}{4}\text{tr}\Bigg{[}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{(i(\not{p}-m_{r}))^{5}}{(p^{2}-m_{r}^{2})^{5}+i\epsilon}(-i\Sigma)\Bigg{]},$

(36)

we then substitute Eqs. 23, 24 and 35, to Eq. 36 perform a Wick rotation, and let $p_{E}\to p_{E}m_{r}$, which gives;

	$\displaystyle D_{2}=\frac{i}{8\pi^{2}}\frac{1}{m_{r}}\int_{0}^{1}dx\int dp_{E}\ p^{2}_{E}\Bigg{[}(p_{E}^{4}-10p_{E}^{2}+5)(-p_{E}^{2}){}$
	$\displaystyle\times\Big{(}1+\frac{(1-x)}{\sqrt{\Delta_{0}}}\Big{)}-(5p_{E}^{4}-10p_{E}^{2}+1){}$
	$\displaystyle\times\Big{(}-3\sqrt{\Delta_{0}}+\frac{1-x}{\sqrt{\Delta}_{0}}(-p_{E}^{2}x)+3\Big{)}\Bigg{]}\frac{1}{(p_{E}^{2}+1)^{5}},$		(37)

where $\Delta_{0}=(p_{E}^{2}x+1)(1-x)$. It is easy to see that the sign of $D_{2}$ depends on the combination of the sign of $m_{r}$ and sign of the level of the CS term, and the same is true for all two-loop contributions to $\beta_{4}$. Evaluating this integral yields

$$D_{2}=\frac{i}{64\pi}\frac{1}{m_{r}}.$$

(38)

In preparation for the calculation of $D_{3}$, we first need to calculate the vertex correction $\Gamma_{1}$, which is shown in Fig. 6 and explicitly given as

$\displaystyle\Gamma_{1}=$

$\displaystyle=2\pi\int\frac{d^{3}q}{(2\pi)^{3}}\frac{\gamma_{\mu}(i(\not{q}-m_{r}))^{2}\gamma_{\nu}}{(q^{2}-m_{r}^{2})^{2}+i\epsilon}\frac{\epsilon^{\mu\nu\alpha}(p-q)_{\alpha}}{(p-q)^{2}+i\epsilon}.$

(39)

Here, if we check the superficial degree of divergence of $\Gamma_{1}$ by counting the net order of $q$, we naively find a logarithmic UV divergence. However, this is not the actual case, because the leading term of the integrand is an odd function of $q$. As a result, the naive logarithmically divergent term has zero coefficient, and the integral in Eq. 39 actually converges. We apply Eq. 26 with $n=2$, change the integration variables as $q-xp\to k$, and obtain

	$\displaystyle\Gamma_{1}$	$\displaystyle=-4\pi\int_{0}^{1}dx\int\frac{d^{3}k}{(2\pi)^{3}}\gamma_{\mu}(\not{k}+x\not{p}-m_{r})^{2}\gamma_{\nu}\epsilon^{\mu\nu\alpha}{}$
		$\displaystyle\times\big{(}-k+p(1-x)\big{)}_{\alpha}\Bigg{[}\frac{1-x}{[k^{2}-\Delta]^{3}+i\epsilon}\Bigg{]}.$		(40)

After using the Eqs. (28) and removing the odd terms in $k$, we have

	$\displaystyle\Gamma_{1}$	$\displaystyle=i8\pi\int_{0}^{1}dx\int\frac{d^{3}k}{(2\pi)^{3}}\Bigg{[}\not{p}(1-x)(k^{2}+x^{2}p^{2}+m_{r}^{2}){}$
		$\displaystyle-\not{k}(2xp\cdot k)+2m_{r}(-k^{2}+(1-x)xp^{2})\Bigg{]}{}$
		$\displaystyle\times\Bigg{[}\frac{1-x}{[k^{2}-\Delta]^{3}+i\epsilon}\Bigg{]}.$		(41)

Now it is clear that this integral is not divergent, because the term that would produce logarithmic UV divergence is canceled as a result of the removal of the odd terms. We can further simplify this by making the following substitution:

$$\not{k}(p\cdot k)\to\frac{1}{3}k^{2}\not{p},$$

(42)

which is a result of the symmetry of the integral in $k$. Then, we perform a Wick rotation:

	$\displaystyle\Gamma_{1}=8\pi\int_{0}^{1}dx\int\frac{d^{3}k_{E}}{(2\pi)^{3}}\Bigg{[}-i\not{p}_{E}(1-x)(-k_{E}^{2}-x^{2}p_{E}^{2}{}$
	$\displaystyle+m_{r}^{2})-\frac{2k_{E}^{2}\not{p}_{E}}{3}+2m_{r}(k_{E}^{2}-(1-x)xp_{E}^{2})\Bigg{]}\Bigg{[}\frac{1-x}{[k_{E}^{2}+\Delta_{E}]^{3}}\Bigg{]}.$		(43)

We can evaluate the integral over $k_{E}$ using the integral identities given in APPENDIX: Integral identities. We have,

	$\displaystyle\Gamma_{1}$	$\displaystyle=\frac{1}{4}\int_{0}^{1}dx(1-x)\bigg{[}\frac{3(i\not{p}_{E}(1-5x/3)+2m_{r})}{\sqrt{\Delta_{E}}}{}$
		$\displaystyle+\frac{(m_{r}^{2}-x^{2}p^{2}_{E})(1-x)(-i\not{p}_{E})-2m_{r}x(1-x)p^{2}_{E})}{\Delta_{E}^{3/2}}\bigg{]}.$		(44)

Finally, we can calculate $D_{3}$, which is shown in Fig. 5 and explicitly given as

$\displaystyle D_{3}$

$\displaystyle=-\frac{1}{4}\text{tr}\Bigg{[}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{(i(\not{p}-m_{r}))^{4}}{(p^{2}-m_{r}^{2})^{4}+i\epsilon}\Gamma_{1}\Bigg{]}.$

(45)

As before, we make the Wick rotation, let $p_{E}\to p_{E}m_{r}$, and substitute $\Gamma_{1}$, which gives

	$\displaystyle D_{3}$	$\displaystyle=-\frac{i}{16\pi^{2}}\frac{1}{m_{r}}\int_{0}^{1}dx\int_{0}^{\infty}dp_{E}p_{E}^{2}\Bigg{[}-p_{E}^{2}\Bigg{(}\frac{3-5x}{\sqrt{\Delta_{0}}}{}$
		$\displaystyle-\frac{(1-x)(1-x^{2}p_{E}^{2})}{\Delta_{0}^{3/2}}\Bigg{)}(-4p_{E}^{2}+4){}$
		$\displaystyle+\Bigg{(}\frac{6}{\sqrt{\Delta_{0}}}-\frac{2(1-x)xp_{E}^{2}}{\Delta_{0}^{3/2}}\Bigg{)}(p_{E}^{4}-6p_{E}^{2}+1)\Bigg{]}{}$
		$\displaystyle\times\frac{1-x}{(p_{E}^{2}+1)^{4}}.$		(46)

Evaluating this integral yields

$$D_{3}=0.$$

(47)

Next we calculate $D_{4}$. First, we start with $\Gamma_{2}$, which is shown in Fig. 6 and explicitly given as

$\displaystyle\Gamma_{2}=$

$\displaystyle=2\pi\int\frac{d^{3}q}{(2\pi)^{3}}\frac{\gamma_{\mu}(i(\not{q}-m_{r}))^{3}\gamma_{\nu}}{(q^{2}-m_{r}^{2})^{3}+i\epsilon}\frac{\epsilon^{\mu\nu\alpha}(p-q)_{\alpha}}{(p-q)^{2}+i\epsilon},$

(48)

which is convergent. First, we apply Eq. 26 with $n=3$ and change the integration variables as $q-xp\to k$, and we have

	$\displaystyle\Gamma_{2}=$	$\displaystyle=-i6\pi\int_{0}^{1}dx\int\frac{d^{3}k}{(2\pi)^{3}}\gamma_{\mu}(\not{k}+x\not{p}-m_{r})^{3}\gamma_{\nu}\epsilon^{\mu\nu\alpha}{}$
		$\displaystyle\times\big{(}-k+p(1-x)\big{)}_{\alpha}\Bigg{[}\frac{(1-x)^{2}}{[k^{2}-\Delta]^{4}+i\epsilon}\Bigg{]}.$		(49)

Next, we simplify the gamma matrix terms by using Eqs. (28) and show that

	$\displaystyle\epsilon^{\mu\nu\alpha}\gamma_{\mu}(\not{k}+x\not{p}-m_{r})^{3}\gamma_{\nu}\big{(}-k+p(1-x)\big{)}_{\alpha}=2i(k^{2}{}$
	$\displaystyle+x^{2}p^{2}+2xk\cdot p+3m_{r}^{2})(k+xp)\cdot(-k+p(1-x)){}$
	$\displaystyle+(-\not{k}+\not{p}(1-x))[3(k^{2}+x^{2}p^{2}+2xp\cdot k)+m_{r}^{2}]2im_{r},$		(50)

then, we apply Eq. 42 and a Wick rotation, so $\Gamma_{2}$ becomes

	$\displaystyle\Gamma_{2}$	$\displaystyle=12i\pi\int_{0}^{1}dx\int\frac{d^{3}k_{E}}{(2\pi)^{3}}\Bigg{(}\Big{(}k_{E}^{2}-p_{E}^{2}x(1-x)\Big{)}{}$
		$\displaystyle\times(-k_{E}^{2}-x^{2}p_{E}^{2}+3m_{r}^{2})+\frac{2x}{3}p_{E}^{2}k_{E}^{2}(1-2x){}$
		$\displaystyle-im_{r}2xk^{2}\not{p}_{E}-im_{r}\not{p}_{E}(1-x)\Big{(}3(-k_{E}^{2}-x^{2}p_{E}^{2})+m_{r}^{2}\Big{)}\Bigg{)}{}$
		$\displaystyle\times\Bigg{[}\frac{(1-x)^{2}}{[k_{E}^{2}+\Delta_{E}]^{4}}\Bigg{]}.$		(51)

We can evaluate the integral over $k_{E}$ using the integral identities given in APPENDIX: Integral identities. $\Gamma_{2}$ is now

	$\displaystyle\Gamma_{2}$	$\displaystyle=\frac{3i}{16}\int_{0}^{1}dx(1-x)^{2}\Bigg{[}\frac{-5}{\Delta_{E}^{1/2}}+\frac{1}{\Delta_{E}^{3/2}}\bigg{(}\frac{5x}{3}p_{E}^{2}(1-2x){}$
		$\displaystyle+3m_{r}^{2}+im_{r}\not{p}_{E}(3-5x)\bigg{)}+\frac{1}{\Delta_{E}^{5/2}}\Bigg{(}-p_{E}^{2}x(1-x){}$
		$\displaystyle\times(-x^{2}p_{E}^{2}+3m_{r}^{2})-im_{r}\not{p}_{E}(1-x)[-3x^{2}p_{E}^{2}+m_{r}^{2}]\Bigg{)}\Bigg{]}.$		(52)

Finally, we can calculate $D_{4}$ which is shown in Fig. 5 and explicitly given as

$\displaystyle D_{4}$

$\displaystyle=-\frac{1}{4}\text{tr}\Bigg{[}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{(i(\not{p}-m_{r}))^{3}}{(p^{2}-m_{r}^{2})^{3}+i\epsilon}(\Gamma_{2})\Bigg{]}.$

(53)

As before we make Wick rotation, let $p_{E}\to p_{E}m_{r}$ and substitute $\Gamma_{2}$ which gives

	$\displaystyle D_{4}=-\frac{3i}{64\pi^{2}}\frac{1}{m_{r}}\int_{0}^{1}dx\int_{0}^{\infty}dp_{E}p_{E}^{2}\frac{(1-x)^{2}}{(p_{E}+1)^{3}}{}$
	$\displaystyle\times\Bigg{\{}-p_{E}^{2}(-p_{E}^{2}+3)\Bigg{(}\frac{3-5x}{\Delta_{0}^{3/2}}-(1-x)\frac{-3x^{2}p_{E}^{2}+1}{\Delta_{0}^{5/2}}+{}$
	$\displaystyle(1-3p_{E}^{2})\Bigg{[}\frac{-5}{\Delta_{0}^{1/2}}+\frac{1}{\Delta_{0}^{3/2}}\bigg{(}\frac{5x}{3}p_{E}^{2}(1-2x){}$
	$\displaystyle+3\bigg{)}+\frac{1}{\Delta_{0}^{5/2}}\bigg{(}-p_{E}^{2}x(1-x)(-x^{2}p_{E}^{2}+3)\bigg{)}\Bigg{]}\Bigg{)}\Bigg{\}}.$		(54)

If we evaluate this numerically, we have

$$D_{4}=-\frac{i}{32\pi}\frac{1}{m_{r}}.$$

(55)

Finally, we get $\beta_{4}$ by using Eqs. 37, 46 and 54 which gives,

$$\beta_{4}=-\frac{1}{64\pi}\frac{1}{m_{r}},$$

(56)

which is the main result of this section.

We now discuss three different cases.

(i) $m_{r}>0$. This describes a CSL phase. Since $\beta_{4}>0$, its transition into the phase with VBS and/or Ising Neel order is first order at the mean-field level.

(ii) $m_{r}<0$. This describes an XY phase. Since $\beta_{4}<0$, its transition into the phase with VBS and/or Ising Neel order is second order at the mean-field level.

It should be noted that our evaluation of $\beta_{4}$ is only to the lowest order in gauge coupling (or inverse CS level), which is of order one. We cannot rule out the possibility that higher order correction can reverse the sign of $\beta_{4}$ and thus the conclusions above.

(iii) $m_{r}=0$. This is our multi-critical point, at which we can not integrate out the (massless) Dirac fermions perturbatively as done above. One can, nevertheless, perform a non-perturbative calculation of the effective potential [14] $V_{\rm eff}(\phi_{\rm cl})$ in terms of $\phi_{\rm cl}$, which is the vacuum expectation value of $\phi$ where $V_{\rm eff}(\phi_{\rm cl})$ is minimized. Since the fermion theory is massless and contains no scale, one expects its coupling to $\phi_{\rm cl}$ generates a scale-invariant term $|\phi_{\rm cl}|^{3}$, which is easy to verify by calculating the change of fermion ground-state energy due to $\phi_{\rm cl}$ that plays the role of a mass. The non-analyticity of such a term originates from the masslessness of the Dirac fermion. Its presence signals the non-mean-field behavior of the transition into the phases with broken translation symmetry, even if the theory is analyzed at the mean-field level.

The theory of multi-critical point is also discussed in Ref. [21], which is mainly done by considering a mean-field approach by considering the dual version of the theory. Thus, to have a connection with the literature, we also briefly find a dual version of our theory. In Sec. II, we started with the lattice spin model given in Eq. 2, then we mapped it to hard-core bosons, and then mapped those hard-core bosons to non-relativistic fermions in a lattice with a level-one CS term. Then, we found that the continuum limit of this theory is described by two Dirac fermions coupled to the level-one CS term given in Eq. 5. In this section, we will apply a bosonization transformation to Eq. 5, which, in a sense, close the circle of our mappings.