Abelian Higgs model at four loops, fixed-point collision and deconfined criticality

TL;DR

This paper performs a four-loop analysis of the abelian Higgs model to determine the critical number of scalar fields where the phase transition changes from second to first order, with implications for deconfined quantum criticality.

Contribution

It provides the first four-loop calculation of the abelian Higgs model in $4- ext{epsilon}$ expansion and estimates the critical number of fields in three dimensions, revealing a transition from second to first order.

Findings

Critical number of scalar fields $n_c$ is approximately 12.2 in 3D.

Transition changes from second to first order below $n_c$.

Strong evidence for a smooth interpolation of $n_c$ between 2 and 4 dimensions.

Abstract

The abelian Higgs model is the textbook example for the superconducting transition and the Anderson-Higgs mechanism, and has become pivotal in the description of deconfined quantum criticality. We study the abelian Higgs model with $n$ complex scalar fields at unprecedented four-loop order in the $4-\epsilon$ expansion and find that the annihilation of the critical and bicritical points occurs at a critical number of $n_c \approx 182.95\left(1 - 1.752\epsilon + 0.798 \epsilon^2 + 0.362\epsilon^3\right) + \mathcal{O}\left(\epsilon^4\right)\nonumber$. Consequently, below $n_c$, the transition turns from second to first order. Resummation of the series to extract the result in three-dimensions provides strong evidence for a critical $n_c(d=3)$ which is significantly below the leading-order value, but the estimates for $n_c$ are widely spread. Conjecturing the topology of the renormalization group flow between two and four dimensions, we obtain a smooth interpolation function for $n_c(d)$ and find $n_c(3)\approx 12.2\pm 3.9$ as our best estimate in three dimensions. Finally, we discuss Miransky scaling occurring below $n_c$ and comment on implications for weakly first-order behavior of deconfined quantum transitions. We predict an emergent hierarchy of length scales between deconfined quantum transitions corresponding to different $n$.

Figures

Tables

Table 1. Table 1: Padé and Padé-Borel approximants of the ϵ italic-ϵ \epsilon expansion evaluated at d = 3 𝑑 3 d=3 ( ϵ = 1 italic-ϵ 1 \epsilon=1 ). The Padé-Borel [ 2 / 1 ] delimited-[] 2 1 [2/1] is not availabe due to a singularity in the integral of the Borel sum.

Order		$\left[m/n\right]$	Padé	Padé-Borel
$ϵ$		$\left[0/1\right]$	66.48	89.11
${ϵ}^2$		$\left[0/2\right]$	36.42	66.78
		$\left[1/1\right]$	–37.25	–46.53
${ϵ}^3$		$\left[0/3\right]$	25.27	58.42
		$\left[2/1\right]$	129.81	–
		$\left[1/2\right]$	1.81	–13.65

Table 2. Table 2: Estimates from polynomial interpolation and two-sided Padé approximation for the critical number of components n c subscript 𝑛 𝑐 n_{c} at d = 3 𝑑 3 d=3 dimensions for a = 1 𝑎 1 a=1 . The values which belong to a non-monotonous interpolating function n c ​ ( d ) subscript 𝑛 𝑐 𝑑 n_{c}(d) for d ∈ [ 2 , 4 ] 𝑑 2 4 d\in[2,4] are printed in italics.

polynom. Int.		Two-sided Padé
$P_{1,1}$	11.60	$\text{2}\left[2/1\right]$	26.62
		$\text{2}\left[1/2\right]$	3.77
$P_{1,2}$	2.15	$\text{2}\left[2/2\right]$	18.76
		$\text{2}\left[1/3\right]$	6.56
		$\text{2}\left[3/1\right]$	-1.40
$P_{1,3}$	17.92	$\text{2}\left[2/3\right]$	16.32
		$\text{2}\left[3/2\right]$	14.11
		$\text{2}\left[1/4\right]$	9.80
		$\text{2}\left[4/1\right]$	8.06

Table 3. Table 3: List of all relevant n 𝑛 n -point functions and the associated number of Feynman diagrams evaluated in dependence on the number of loops ( A 𝐴 A : gauge field, ϕ italic-ϕ \phi : scalar field).

Loops	1	2	3	4
	$2$	$20$	$370$	$9,291$
	$4$	$27$	$459$	$11,332$
	$5$	$107$	$3,078$	$106,501$
	$32$	$1,042$	$40,164$	$1,735,706$

Table 4. Table 4: Two-sided Padé estimates from four loops at d = 3 𝑑 3 d=3 and different slopes a 𝑎 a for the expansion in d = 2 + ϵ 𝑑 2 italic-ϵ d=2+\epsilon . In the last column we show the mean value of the two-sided Padés at a specific a 𝑎 a and its standard deviation. Slanted numbers denote non-monotonic curves between d ∈ ( 2 , 4 ) 𝑑 2 4 d\in(2,4) .

$a$	$\text{2}\left[2/3\right]$	$\text{2}\left[3/2\right]$	$\text{2}\left[1/4\right]$	$\text{2}\left[4/1\right]$	$\approx$
0.2	15.85	13.93	2.84	7.93	10.14 $\pm$ 5.92
0.5	16.02	14.00	6.08	7.98	11.02 $\pm$ 4.75
1	16.32	14.11	9.80	8.06	12.07 $\pm$ 3.80
2	16.94	14.32	14.12	8.22	13.40 $\pm$ 3.68
5	19.21	14.97	19.20	8.71	15.52 $\pm$ 4.96
10	24.98	16.03	21.82	9.51	18.08 $\pm$ 6.81

Equations

$$n_c\approx n_{c,0}\left(1-1.752ϵ+0.798{ϵ}^2+0.362{ϵ}^3\right)+\mathcal{O}\left({ϵ}^4\right),$$

$$\mathcal{L}=|D_{\mu }\varphi {|}^2+\frac{1}{4}{F}_{\mu \nu }^2+r|\varphi {|}^2+\lambda (|\varphi {|}^2{)}^2.$$

$${\mathcal{L}}^{\prime }=$$

$$+\frac{Z_A}{4}{F}_{\mu \nu }^2-\frac{1}{2\xi }({\partial }_{\mu }{A}_{\mu }{)}^2.$$

$$\alpha =e_0^2{\mu }^{-ϵ}{Z}_A,\lambda ={\lambda }_0{\mu }^{-ϵ}{Z}_{\varphi }^2{Z}_{{\varphi }^4}^{-1}.$$

$${\beta }_i=\frac{\mathrm{d}{g}_i}{\mathrm{d}\ln b}=ϵg_i+\sum _k{\beta }_i^{(k\ell )},$$

$${\beta }_{\alpha }^{(1\ell )}$$

$${\beta }_{\alpha }^{(2\ell )}$$

$${\beta }_{\alpha }^{(3\ell )}$$

$${\beta }_{\lambda }^{(1\ell )}$$

$${\beta }_{\lambda }^{(2\ell )}$$

$$-(4n+10)\alpha {\lambda }^2+\left(\frac{9}{2}n+\frac{21}{2}\right){\lambda }^3,$$

$${\beta }_{\lambda }^{(3\ell )}$$

$$+\left(-\frac{5}{216}{n}^2+\left[18{\zeta }_3-\frac{989}{8}\right]n-\frac{889}{4}-54{\zeta }_3\right){\alpha }^3\lambda$$

$$+\left(\frac{43}{16}{n}^2+\left[18{\zeta }_3+\frac{1749}{16}\right]n+\frac{1093}{8}+126{\zeta }_3\right){\alpha }^2{\lambda }^2$$

$$+\left(-\frac{33}{16}{n}^2+\left[-15{\zeta }_3-\frac{461}{16}\right]n-\frac{185}{4}-33{\zeta }_3\right){\lambda }^4$$

$$+\left(\left[\frac{25}{2}-6{\zeta }_3\right]n+\frac{29}{2}+6{\zeta }_3\right)\alpha {\lambda }^3,$$

$${\beta }_i(\{g_i^{\ast }\})=0\forall i,$$

$${\beta }_i(\{g_i\})=\sum _j{B}_{ij}(g_j-g_j^{\ast })+\mathcal{O}(\{g_i^2\}),$$

$$\det (B_{ij}){|}_{\{g_i^{\ast }\}}=0,$$

$${\beta }_{\alpha }^{(1\ell )}=-\frac{n}{3}{\alpha }^2,{\beta }_{\lambda }^{(1\ell )}$$

$${\lambda }_{\pm }^{\ast }=\frac{3(18+n\pm \sqrt{s})}{2n(n+4)}ϵ+\mathcal{O}({ϵ}^2),$$

$$n_{c,0}=6(15+4\sqrt{15})\approx 182.95.$$

$$n_c\approx n_{c,0}+n_{c,1}ϵ+n_{c,2}{ϵ}^2+n_{c,3}{ϵ}^3+\mathcal{O}({ϵ}^4),$$

$$n_{c,1}=$$

$$n_{c,2}=$$

$$+$$

$$n_{c,3}=$$

$$+$$

$$n_c(d=3)\approx 12.2\pm 3.9.$$

$$\displaystyle\frac{d\lambda}{d\ln b}\Big{|}_{\alpha_{\ast}}=\left(1+\frac{18}{n}\right)\epsilon\lambda-(n+4)\lambda^{2}-\frac{54}{n^{2}}\epsilon^{2}\,.$$

$$\displaystyle\ln\!\left(\!\frac{b_{\mathrm{IR}}}{b_{\mathrm{UV}}}\!\right)\!=\!\frac{-2n}{\epsilon\sqrt{-s}}\arctan\!\left(\frac{2n(4+n)\lambda\!-\!(18+n)\epsilon}{\epsilon\sqrt{-s}}\right)\!\bigg{|}_{\lambda_{\mathrm{UV}}}^{\lambda_{\mathrm{IR}}},$$

$$L_{\mathrm{IR}}(n)\simeq \exp \left(\pi f(n)\right),\mathrm{where}f(n)=\frac{2n}{ϵ\sqrt{-s}},$$

$${\beta }_{\alpha }^{(4\ell )}=\left(-\frac{323}{3888}{n}^3+\left[\frac{451}{54}-\frac{38{\zeta }_3}{9}\right]n^2+\frac{3}{16}n\right){\alpha }^5+\left(\frac{5}{72}{n}^3-\frac{41}{9}{n}^2-\frac{37}{8}n\right){\alpha }^4\lambda +\left(\frac{1}{24}{n}^3+\frac{139}{48}{n}^2+\frac{137}{48}n\right){\alpha }^3{\lambda }^2$$

$$+\left(-\frac{5}{48}{n}^3-\frac{7}{16}{n}^2-\frac{1}{3}n\right){\alpha }^2{\lambda }^3$$

$${\beta }_{\lambda }^{(4\ell )}=\left(\left[40{\zeta }_5-\frac{27{\zeta }_3}{2}-\frac{6145}{48}-\frac{{\pi }^4}{180}\right]n^2-\left[\frac{3157{\zeta }_3}{6}+80{\zeta }_5+\frac{25123}{24}+\frac{143{\pi }^4}{180}\right]n-\frac{12{\pi }^4}{5}-\frac{18503}{16}-768{\zeta }_3-960{\zeta }_5\right){\alpha }^2{\lambda }^3$$

$$+\left(\left[\frac{29{\zeta }_3}{2}-\frac{{\pi }^4}{60}-\frac{377}{96}\right]n^2+\left[-\frac{269{\zeta }_3}{2}+50{\zeta }_5+\frac{7{\pi }^4}{10}-\frac{1403}{32}\right]n+\frac{103{\pi }^4}{60}-\frac{185}{12}+150{\zeta }_5-435{\zeta }_3\right)\alpha {\lambda }^4$$

$$+\left(\left[\frac{1}{81}-\frac{2{\zeta }_3}{9}\right]n^3+\left[28{\zeta }_3-\frac{{\pi }^4}{10}-\frac{55709}{1296}\right]n^2+\left[310{\zeta }_5-\frac{578{\zeta }_3}{3}-\frac{13987}{36}-\frac{19{\pi }^4}{60}\right]n-\frac{12751}{16}-\frac{33{\pi }^4}{20}+504{\zeta }_3-390{\zeta }_5\right){\alpha }^5$$

$$+\left(\left[\frac{19{\zeta }_3}{18}+\frac{67}{2592}\right]n^3+\left[20{\zeta }_5-\frac{77{\zeta }_3}{2}+\frac{{\pi }^4}{5}+\frac{12779}{162}\right]n^2+\left[\frac{431{\zeta }_3}{2}+305{\zeta }_5+\frac{13{\pi }^4}{10}+\frac{209}{12}\right]n+\frac{26{\pi }^4}{5}-\frac{19127}{96}+2095{\zeta }_5-\frac{1191{\zeta }_3}{2}\right){\alpha }^4\lambda$$

Full text

I Introduction

The abelian Higgs (AH) model is one of the most fundamental field theories in both condensed matter and particle physics. It serves as the prime textbook example for the superconducting transition and the Anderson-Higgs mechanism Zinn-Justin (1996); Peskin (2018); Herbut (2007); Altland and Simons (2010)111To honor its many discoverers appropriately it should always be referred to as the Anderson-Englert-Brout-Higgs-Guralnik-Hagen-Kibble mechanism Anderson (1963); Englert and Brout (1964); Higgs (1964); Guralnik et al. (1964), of course.. The AH model features a complex scalar field coupled to a fluctuating U(1) gauge field, and it displays two distinct phases separated by a sharp transition: the symmetric phase and the phase with spontaneously broken symmetry. In the context of superconductors, the symmetric phase is related to the normal metallic state and the spontaneously symmetry broken phase to the superconducting Meissner state. For the three-dimensional AH model with a single complex scalar, it has been established employing duality arguments that the transition is related to the one of the XY model Dasgupta and Halperin (1981); Kiometzis et al. (1994), which is known to be continuous, i.e. exhibiting a critical point. Numerical simulations of lattice versions of both models confirm the conjectured mapping, and also suggest a continuous transition Dasgupta and Halperin (1981); Bartholomew (1983); Nguyen and Sudbø (1999).

The field-theoretical analysis of the AH model including its generalized version with $n$ complex scalars has turned out to be surprisingly subtle. The subtlety originates in the presence of fluctuating complex scalars and other fluctuating massless modes – here, coming from the gauge fields. Such a scenario is generic and also arises in quantum phase transitions of electronic systems. Indeed, in the AH model, the determination of the nature of transition as a function of $n$, i.e. the question of whether it is discontinuous or continuous, is a long-standing problem in the theory of critical phenomena Halperin et al. (1974); Dasgupta and Halperin (1981); Bartholomew (1983); Nguyen and Sudbø (1999); Kolnberger and Folk (1990); Radzihovsky (1995); Kiometzis et al. (1994); Herbut and Tešanović (1996, 1997); Lawrie (1997); Bergerhoff et al. (1996); Bartosch (2013); Fejos and Hatsuda (2017). In mean-field approximation, for example, the $n=1$ transition is found to be discontinuous. This is further supported by renormalization group (RG) calculations at one-loop order and in $4-\epsilon$ dimensions, which only for a very large value of $n>n_{c}\sim 183$ show a critical point Halperin et al. (1974). The next-to-leading order loop expansion, on the other hand, shows a dramatic reduction of the critical $n_{c}$, but due to the large magnitude of this correction, it remains rather inconclusive Herbut and Tešanović (1996, 1997) and has led to the belief that a perturbative RG is not well-suited for the question at hand. Therefore, despite the tremendous recent advances in higher-loop calculations for (quasi)relativistic field theories Kompaniets and Panzer (2017); Batkovich et al. (2016); Kompaniets and Panzer (2016); Gracey et al. (2016); Zerf et al. (2016); Mihaila et al. (2017); Zerf et al. (2017); Gracey (2018, 2019); Ihrig et al. (2018a); Zerf et al. (2018, 2019), further efforts to improve the description of the $n$-component Abelian Higgs model have not been undertaken for many years. In this work, we fill this long-standing void by performing a perturbative RG calculation to four-loop order and in the parameter $\epsilon=4-d$, primarily to provide an improved estimate for the critical number of complex scalars $n_{c}$ above which the transition becomes continuous.

In fact, the $n$-component extension of the AH model is not only a playground for theoretical methods, but has been suggested to be the relevant class of models to describe certain quantum phase transitions beyond the Landau-Ginzburg paradigm Senthil et al. (2004a). Such transitions have been argued to appear in two-dimensional spin models at zero temperature and in three-dimensional classical loop models between two competing ordered phases and are putatively separated by a critical point – the deconfined quantum critical point (DQCP). Explicitly, for $n=2$, the 2+1 dimensional AH model was discussed in this context due to its relation to the NCCP1 model Senthil et al. (2004a, b, 2005); Motrunich and Vishwanath (2004); Kaul and Sandvik (2012). Similar scenarios have been devised for $n=3,4$ Harada et al. (2013); Kaul (2011) and can possibly be generalized to even higher $n$. The excitement about the DQCP scenario stems from its unusual phenomenology Senthil et al. (2004a) which – besides featuring a critical point – includes the emergence of higher symmetry and fractionalized degrees of freedom becoming deconfined at the critical point. Indeed, numerical studies of spin models have found indications for such phenomenology, but have also experienced complications due to apparent drift of the values of the critical exponents. More recent numerical analyses Nahum et al. (2015); Serna and Nahum (2019) argued in favor of a weakly first-order transition instead of a continuous one, but still see emergent symmetry as predicted by the DQCP scenario. Absence of a continuous transition with the scaling dimensions as found in Monte Carlo simulations is further supported by bounds from symmetry-enhanced conformal field theories Nakayama and Ohtsuki (2016); Poland et al. (2019). These combined results suggest that the DQCP scenario is more subtle than initially expected from basic arguments.

The finding that similar scaling properties and violations are present in a variety of models suggests that the numerically observed phenomenology is not uniquely tied to a particular lattice realization and also not exclusively explicable by a fourfold anisotropy of the valence bond solid order parameter Nahum et al. (2015); Serna and Nahum (2019); Lou et al. (2009); Pujari et al. (2013); Block et al. (2013); Powell and Chalker (2008); Charrier et al. (2008); Chen et al. (2009). Rather, it can be expected to be related to a general mechanism which is tied to universal physics described by an effective theory at the phase transition, i.e. the AH model.

In fact, weakly first-order transitions naturally appear in complex conformal field theories Nahum et al. (2015); Gorbenko et al. (2018); Wang et al. (2017); Serna and Nahum (2019) where a RG fixed point is complex-valued and, if imaginary parts are small, the complex fixed point slows down and controls the RG flow of real and unitary gapped physical theories. Such behavior naturally occurs in models where two RG fixed points collide, annihilate and move into the complex plane Kubota and Terao (2001); Kaveh and Herbut (2005); Kaplan et al. (2009); Gies and Jaeckel (2006); Herbut and Janssen (2014); Herbut (2016); Gracey et al. (2018), e.g., the AH model below $n_{c}$.

These observations suggest that the $n_{c}$ of the three-dimensional $n$-component AH model should be above but still near $n=2,3,4$, to be compatible with the numerical findings. Through the explicit higher-loop RG analysis presented in this work, we provide the quantitative background for this scenario. Furthermore, we note that the special point of $n=1$, which exhibits a continuous transition, is not continuously connected to the critical point at large $n>n_{c}$ Nahum et al. (2015).

The rest of the paper is organized as follows: We summarize our key results in Sec. II. Then, in Sec. III, we introduce the model, notation and also the renormalization group procedure. Sec. IV discusses the renormalization group fixed points and the mechanism of annihilation of two fixed points which underlies weakly first-order transitions. The epsilon expansion at four loops of the critical number of complex scalars is also presented there. We discuss resummation of the series in epsilon in Sec. V and interpolation between two and four dimensions in Sec. VI. Finally, we discuss implications of our findings on the DQCP scenario in Sec. VII and conclude in Sec. VIII. Technicalities are given in the appendix.

II Key results

We study the $n$-component AH model at four loops and in $4-\epsilon$ dimensions. An annihilation of the critical and bicritical points occurs at

[TABLE]

with $n_{c,0}\approx 182.95$. For $n<n_{c}$, the phase transition turns from second to weakly first order. Based on conjectures about the topology of the RG flow between two and four dimensions, we obtain a smooth interpolating function $n_{c}(d)$. Our best estimate for the three-dimensional case yields $n_{c}(3)\approx 12.2\pm 3.9$. Slightly below $n_{c}$, Miransky scaling occurs exhibiting a correlation length which is exponentially enhanced upon decreasing the distance to the point of fixed-point collision, e.g., by changing $n$. Due to the exponential dependence, small variations of $n$ are expected to lead to large changes in the correlation length unless the argument in the exponential is an exceptionally flat function of $n$. Our calculations do not provide indications of such behavior and we therefore predict a substantial hierarchy of correlation lengths for different $n$ near $n_{c}$ which may be applicable to $n=2,3,4$ and should be observable in future numerical simulations.

III Model and method

The $n$-component AH model is also known as scalar quantum electrodynamics (QED) and is defined in $d$-dimensional euclidean space(time) by the Lagrangian

[TABLE]

Here, $\phi=(\phi_{1},...,\phi_{n})$ describes the $n$-component complex scalar field with mass term $r$ and quartic interaction $\lambda$. It is minimally coupled to the dynamical non-compact U(1) gauge field $A_{\mu}$ via the covariant derivative $D_{\mu}=\partial_{\mu}-ieA_{\mu}$ with charge $e$ and indices $\mu,\nu$ run from [math] to $d-1$. Common summation convention over repeated indices is implied. The gauge field comes with the field strength tensor $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ and we add a gauge fixing term $\mathcal{L}_{\text{gf}}=-\frac{1}{2\xi}(\partial_{\mu}A_{\mu})^{2}$ where $\xi$ denotes the gauge fixing parameter.

For the case of a single complex scalar field, $n=1$ and in three spatial dimensions, the model is paradigmatically used to describe the superconducting transition Ginzburg and Landau (1950) and also the nematic-to-smectic transition in liquid crystals Halperin et al. (1974). Generally, the scalar mass parameter $r$ can be considered as the tuning parameter of the transition towards, e.g., the superconducting state where the amplitude fluctuations of $\phi$ become massive and – if the charge is finite $e\neq 0$ – the phase fluctuations of $\phi$ can be completely absorbed into the gauge field which becomes massive, too. This is the abelian and simplest version of the Higgs mechanism.

III.1 Renormalization group procedure

The AH model has an upper critical dimension of $d_{c}^{+}=4$ where the charge and the quartic coupling, $g_{i}\in\{g_{1}=\alpha=e^{2},g_{2}=\lambda\}$, become marginal simultaneously, i.e. it is perturbatively renormalizable in $d\leq d_{c}^{+}$. Here, we present a perturbative renormalization group analysis in $4-\epsilon$ dimensions using dimensional regularization (DREG) and the modified minimal subtraction scheme ($\overline{\text{MS}}$) to four-loop order. To that end, we introduce the bare Lagrangian, which is the one from Eq. (2) with fields and couplings replaced by their bare counterparts, i.e. $x\mapsto x_{0}$ and $x\in\{\phi,A_{\mu},e,r,\lambda,\xi\}$. Then, the renormalized Lagrangian reads

[TABLE]

with $D_{\mu}\phi=(\partial_{\mu}-ie\mu^{\epsilon/2}A_{\mu})\phi$ and $\mu$ defines the energy scale parametrizing the renormalization group flow. Here, we have introduced explicit $\mu$ dependencies to rescale the dimensionless couplings in $4-\epsilon$ dimensions as well as the wavefunction renormalizations $Z_{\phi}$ and $Z_{A}$ to relate the bare and the renormalized Lagrangian through $\phi_{0}=\sqrt{Z_{\phi}}\phi$ and $A_{0,\mu}=\sqrt{Z_{A}}A_{\mu}$. Accordingly, the bare and the renormalized bosonic mass terms are related by $r=r_{0}\mu^{-2}Z_{\phi}Z_{\phi^{2}}^{-1}$ and we obtain the following relations between the bare and the renormalized couplings

[TABLE]

For completeness, we note that the flow of the gauge-fixing parameter is encoded in the relation $\xi=\xi_{0}Z_{A}^{-1}$. The calculation of the renormalization constants $Z_{x}$ with $x\in\{\phi,\phi^{2},\phi^{4},A\}$ at four-loop order is performed by using an automated protocol which is described in App. A.

From the renormalization contants, we can construct the renormalization group beta functions which are defined as the logarithmic derivatives of the dimensionless renormalized couplings $\{g_{1}=\alpha=e^{2},g_{2}=\lambda\}$ with respect to $b=\mu^{-1}$. Schematically, the beta functions have the form

[TABLE]

where the index $k$ indicates the loop order. Explicitly, we obtain for the gauge coupling $\alpha$,

[TABLE]

The beta functions for the quartic self-interaction coupling $\lambda$ up to three-loop order are given accordingly as

[TABLE]

where $\zeta_{3}=\zeta(3)$ denotes the Riemann $\zeta$-function. The four loop contributions are presented in App. B. Note that the beta functions are gauge parameter independent as expected. Further, the expressions are in full agreement to the two-loop results Kolnberger and Folk (1990); Herbut (2007). Also, in the limit $\alpha\rightarrow 0$ we recover the beta functions of the purely bosonic $2n$-vector model up to four-loop order Kleinert and Schulte-Frohlinde (1991). For completness, we list the beta function of the gauge fixing parameter in the App. B. Note that we also provide the field anomalous dimensions there.

IV Fixed points and their collision

Generally, fixed points (FP) of the RG flow are determined by the solutions of the beta functions

[TABLE]

where $\{g_{i}\}$ denotes the set of couplings and $\{g_{i}^{\ast}\}$ the set of fixed point coordinates. Linearization of the RG flow near a FP determines its critical behavior,

[TABLE]

where $B_{ij}=\partial\beta_{i}/\partial g_{j}$ defines the stability matrix. The eigenvalues $\theta_{i}$ and eigenvectors $\vec{v}_{i}$ of $B_{ij}$ determine the critical exponents and the corresponding directions. Positive eigenvalues describe relevant and negative ones irrelevant directions. A vanishing eigenvalue means the RG flow in this direction becomes marginal.

When two FP solutions coincide, the flow between them is marginal, implying that the stability matrix has a vanishing eigenvalue, i.e. the determinant vanishes,

[TABLE]

providing an additional condition to the fixed point equations in Eq. (12). We use Eq. (14) and the condition of a vanishing eigenvalue, $\theta_{i}=0$ for some $i$, as a criterion for the appearance of a fixed point collision for both, the abelian Higgs model in $4-\epsilon$ dimensions as well as the non-linear sigma model in $2+\epsilon$ dimensions, see below.

IV.1 AH model and fixed point collision

At one-loop order, the beta functions read

[TABLE]

This set of beta functions features four fixed points (FP) where $\epsilon$, i.e. $\beta_{i}(\alpha^{*},\lambda^{*})=0$ and scaling may emerge. In $d=4-\epsilon$, one FP is the trivial Gaussian fixed point, $\alpha^{\ast}=\lambda^{\ast}=0$, and another one is the Wilson-Fisher fixed point of the scalar O($2n$) model, $\alpha^{\ast}=0,\lambda^{\ast}=\epsilon/(n+4)+\mathcal{O}(\epsilon^{2})$.

The two remaining non-trivial non-Gaussian fixed points have identical FP coordinate of the gauge coupling, $\alpha^{*}=3\epsilon/n+\mathcal{O}(\epsilon^{2})$. There are two corresponding FP solutions for the quartic interaction,

[TABLE]

where $s=n^{2}-180n-540$. We refer to ($\alpha^{\ast}$,$\lambda^{*}_{+}$) as the charged FP and to ($\alpha^{\ast}$,$\lambda^{*}_{-}$) as the bicritical FP . The charged FP is irrelevant in both directions, i.e. represents a stable FP. The bicritical FP is unstable to perturbations in the quartic coupling.

The fixed points have a strong dependence on the number of field components $n$, in particular, the quartic coupling only provides real-valued solutions above a certain critical $n>n_{c}$, when the radicand $s$ of the square root in Eq. (16) is positive. At $n=n_{c}$ the two fixed points collide and, for $n<n_{c}$, drift into the complex plane. The one-loop analysis predicts that this happens for values of

[TABLE]

Consequently, for $n\geq n_{c}$ there is a continuous transition while below it is expected to become first-order Halperin et al. (1974). Note that we neglected an unphysical negative solution $n_{c,-}=6(15-4\sqrt{15})+\mathcal{O}(\epsilon)\approx-2.952+\mathcal{O}(\epsilon)$. At the next order around the negative solution $n_{c,-}$, we find it to become even more negative $n_{c,-}\approx-2.952(1+1.62323\epsilon+\mathcal{O}(\epsilon^{2}))$.

IV.2 Fixed-point collision at four loops

In the previous section, we considered the fixed points of the AH model at one-loop level and found that two fixed points collide and disappear in the complex plane when tuning the number of scalar components $n$. The critical number $n_{c}$ below which the fixed points disappear in the complex plane was computed by analyzing the imaginary part of the quartic coupling’s charged FP coordinate. Here, we compute corrections to this value at higher orders.

We use Eqs. (12) and (14) to calculate the fixed point collision point $n_{c}$ in the parameter $n$ at higher orders in the $\epsilon$ expansion to obtain the series expansion

[TABLE]

with $n_{c,0}$ from Eq. (17) and we find

[TABLE]

where, again, $\zeta_{3}=\zeta(3)$ and $\zeta_{5}=\zeta(5)$ denote the Riemann $\zeta$ function. Numerical evaluation then provides the approximate expression presented in Eq. (1). This series expansion exceeds previous estimates by two orders in $\epsilon$. Since the coefficients in the series of Eq. (18) are still decreasing in magnitude, we provide direct evaluation for $\epsilon=1$, see Fig. 2. In the following, we will explore resummation of the series.

V Resummation of $4-\epsilon$ expansion

The expansions from perturbative RG calculations are asymptotic series. Consequently, obtaining reliable estimates for $n_{c}$ at $d=3$ requires resummation. Whether or not such resummation yields reliable results depends on the underlying model, the order of the expansion and the knowledge about the large-order behavior. Here, we explore three resummation schemes for the series in $4-\epsilon$ dimensions, i.e. 1. Padé approximants, 2. Padé-Borel approximants and 3. Borel resummation. The employed schemes presented here are conventional and we use standard notation, e.g., for the Padé approximants $[m/n]$. Explicit definitions are provided in App. C. We display the available Padé and Padé-Borel approximants order by order in Fig. 3 and in Tab. 1. The Borel resummation scheme including the corresponding set of resummation parameters is introduced in App. D and we merely state the results here. Optimization through variation of the resummation parameters according to the principle of minimal sensitivity and the principle of fastest convergence yields a negative weighted mean value for $n_{c}(\epsilon=1)$ with a huge error of $n_{c}(\epsilon=1)\approx-52\pm 45$. In the present case, we find that the three methods do not allow us to extract a precise numerical result for $n_{c}(\epsilon=1)$ and estimates are scattered over a rather large range.

The main conclusion drawn from this one-sided resummation analysis is therefore merely that $n_{c}(\epsilon\!=\!1)$ lies significantly below the leading-order result of $n_{c,0}\approx 183$.

VI Dimensional interpolation

In view of the high order series for $n_{c}$ in $\epsilon$, we would like to obtain a better estimate on $n_{c}(\epsilon=1)$ than the one from the resummations explored in the previous section. This is possible through an interpolation between two and four dimensions and by additionally including the assumption of $n_{c}(d=2+\epsilon)=a\epsilon$, with unknown coefficient $a>0$. Before we carry out the explicit interpolations between two and four dimensions, we briefly review the arguments leading to that conjecture.

VI.1 Nonlinear $\sigma$ model in $2+\epsilon$ dimensions

The first item is suggested from studying the nonlinear sigma model in $2+\epsilon$ dimensions, e.g., in the $\mathbb{CP}^{n-1}$ formulation Hikami (1979); March-Russell (1992), which can be related to the abelian Higgs model in the following way: Upon additional expansion in large $n$ the $2+\epsilon$ expansion of the nonlinear sigma model and the $4-\epsilon$ expansion of the AH model describe the same fixed point March-Russell (1992). This is exhibited by comparing their critical exponents with the large-$n$ result expanded in the respective dimension which turn out to be identical. In App. E, we provide further evidence for this correspondence by exploiting the four-loop expansion in $4-\epsilon$ dimensions for the correlation length exponent for the first time. Therefore, the two expansions near two and four dimensions in the respective models can be employed to provide for a continuous interpolation of the same renormalization group fixed points in the plane $(n,d)$ for large enough $n>n_{c}(d)$.

In the $\mathbb{CP}^{n-1}$ model $n_{c}(d\rightarrow 2)\rightarrow 0$ Hikami (1979); March-Russell (1992). Moreover in $d=2+\epsilon$, the beta function of the model coupling $t$ features a real valued fixed-point solution $t^{\ast}$ for all $n>0$, which for small $\epsilon$ is $t^{\ast}=\epsilon/n+\mathcal{O}(\epsilon^{2})$. Therefore, a fixed-point collision is expected to be exhibited in the symmetry-allowed higher-derivative terms of the non-linear sigma model. For example, the RG scaling of the canonically least irrelevant higher-derivative terms, i.e. the four-derivative terms, is found to be222The scaling of higher-derivative terms is found to be $y=2-2k+\epsilon+2k(k-1)\,t^{\ast}$ for terms $2k$ derivatives, cf. Ref. Lerner and Wegner (1990). $y_{4}=-2+\epsilon+4\,t^{\ast}$ at one-loop order. A fixed-point collision would now be indicated by a vanishing of $y_{4}$, i.e. when the scaling becomes marginal, cf. Eq. (14), which at that order yields $n_{c}(d=2+\epsilon)=2\epsilon+\mathcal{O}(\epsilon^{2})$.

At this point a word of caution seems in order. Concerning a possible fluctuation-induced fixed-point destabilization due to the RG relevance of terms with an even larger number of derivatives in nonlinear sigma and related models, there has been a extended discussion in the literature, e.g., Refs. Lerner and Wegner (1990); Castilla and Chakravarty (1993); Brézin and Hikami (1997), which we here will not attempt to resolve. Instead, we take on a pragmatic approach by assuming that near two dimensions the first fluctuation-induced relevant direction comes from the least irrelevant canonical terms, which is the four-derivative term. Moreover, since $t^{\ast}$ is of $\mathcal{O}(1)$ at that point and therefore not small, the corrections to the leading-order behavior of $y_{4}$ cannot be expected to be small either. We therefore refrain from fixing the coefficient of the epsilon expansion of $n_{c}(d+\epsilon)$ and just write $n_{c}(d=2+\epsilon)=a\epsilon$, with unknown coefficient $a>0$ as stated above, see also Ref. Nahum et al. (2015).

VI.2 Interpolation

A suitable interpolation between the two critical dimensions can now be constructed by using a polynomial ansatz. To that end, we use both epsilon expansions for $n_{c}$, simultaneously, and set up an interpolating function in the interval $d\in[2,4]$. More specifically, we choose a polynomial interpolation with polynomial $P_{i,j}(d)$ of degree $i+j$, where $i$ ($j$) denotes the highest order of the epsilon expansion in $d=2+\epsilon$ ($d=4-\epsilon$) dimensions, i.e $i=1$ and $j\in\{1,2,3,4\}$. We fix the polynomial coefficients of the first $i+1$ terms with the expansion near the lower critical dimension. The remaining $j+1$ higher-order coefficients are then determined from the requirement that the $j$ lowest derivatives of $P_{i,j}(d)$ at $d=4$ correspond to the $4-\epsilon$-expansion. The resulting polynomials are then by construction $i$-loop exact near the lower critical dimension and ($j$-loop) exact near the upper critical dimension. From this calculation, we obtain interpolating functions for general coefficient $a$, i.e. $P_{1,1}(d=3)\approx 11.35+0.25a$, $P_{1,2}(d=3)\approx 2.023+0.125a$, $P_{1,3}(d=3)\approx 17.86+0.0625a$, exhibiting that the dependence on the coefficient $a$ is rather weak.

In Ref. Nahum et al. (2015) the additional conjecture that $n_{c}(d)$ increases monotonically, i.e. $n_{c}^{\prime}(d)>0$, has been put forward, however, without further justification. Here, we briefly discuss the impact of this assumption: while $P_{1,3}(d)$ increases monotonically for $0<a\lesssim 80$, $P_{1,1}(d)$ requires $10\lesssim a\lesssim 200$ and $P_{1,2}(d)$ is never monotonous for $2<d<4$. We show the results for the polynomial interpolation as a function of $d$ in the left panel of Fig. 4 and compiled their estimates at $d=3$ in Tab. 2.

Following the same reasoning, we can set up two-sided Padé approximants $\textit{2}_{[m/n]}(d)$, cf. the right panel of Fig. 4 and App. C for definitions. Tentatively, setting $a=1$, we observe that increasing the order in $\textit{2}_{[m/n]}(3)$ from $m+n=3$ to $m+n=5$ decreases the range spanned by the estimates for $n_{c}(d=3)$ order by order. This behavior is even more pronounced upon disregarding the values which belong to non-monotonous interpolations, i.e. dropping the values for $\textit{2}[3/1]$ and $\textit{2}[4/1]$.

We use the highest order, $m+n=5$, including the non-monotonous $\textit{2}[4/1]$, to calculate the average and error

[TABLE]

This is our best estimate based on the above reasoning. We have also studied the dependence on the parameter $a$ and find that, within an extended range of $a\in[0.2,5]$, this estimate varies mildly, see Tab. 4 in App. C.

VII Deconfined pseudo criticality

The $n$-component abelian Higgs model has been argued to describe the universal properties of the Néel–valence-bond-solid transition Senthil et al. (2004a, b) and related scenarios. Recent numerical simulations of this transition revealed strong violations of scaling and have argued that this transition appears to be weakly first order with an anomalously large correlation length and drifting critical exponents, see Refs. Nahum et al. (2015); Wang et al. (2017); Serna and Nahum (2019) and references therein.

VII.1 Miransky scaling

The behavior described above may be due to emergent walking behavior that appears when two fixed points have just been annihilated and vanished into the complex plane, e.g., just as it happens in the AH model below $n_{c}$. Here, we briefly recap the underlying reasoning using the beta functions at one-loop order, cf. Eqs. (15). At that order for the gauge coupling the fixed point value $\alpha_{\ast}=3\epsilon/n$ can be acquired for any choice of $n\lessgtr n_{c}$. Further, the gauge coupling is irrelevant and we therefore can replace $\alpha=\alpha_{\ast}$ in the beta function for the quartic coupling $\lambda$, to obtain

[TABLE]

This flow equation can be integrated yielding

[TABLE]

with $s=n^{2}-180n-540$, again. For $n\lesssim n_{c}$, where $s<0$ and $|s|\ll 1$, the flow of $\lambda$ proceeds as follows: $\lambda$ starts at some positive value, goes through a walking regime and eventually diverges towards negative values. Therefore, in the infrared and the ultraviolet, the argument of the $\arctan$ is always large, i.e. we can use $\lim_{x\to\pm\infty}\arctan x=\pm\pi/2$. Arbitrarily choosing $b_{\mathrm{UV}}=1$ and renaming $b_{\mathrm{IR}}=L_{\mathrm{IR}}$, we obtain the exponentially large infrared length scale Nahum et al. (2015); Kaplan et al. (2009); Gies and Jaeckel (2006)

[TABLE]

which is also referred to as Miransky scaling. Note that including the running of the infrared-attractive gauge coupling yields power-law corrections to this behavior Braun et al. (2011). Further, higher-loop contributions will generally modify the function $f(n)$.

Now, approaching $n_{c}$ from below, this correlation length grows exponentially to diverge at $n_{c}$. Therefore, to obtain a large correlation length, one needs to be close to $n_{c}$, which at one-loop order, Eq. (24), is $n_{c,0}\sim 183$. Generically, small changes in $n\lesssim n_{c}$ induce a large hierarchy in length scales, for example at $\epsilon=1$, $L_{\mathrm{IR}}(180)/L_{\mathrm{IR}}(179)\sim 833$ or $L_{\mathrm{IR}}(179)/L_{\mathrm{IR}}(178)\sim 97$ which, however, decreases rapidly when $n$ moves away from $n_{c}$. In case, $n$ is further away from $n_{c}$, the correlation length is, of course, exponentially suppressed.

Including higher-loop contributions into this reasoning will significantly change the quantitative aspects of this behavior, but not its qualitative features. For this work, we refrain from providing a corresponding quantitative analysis, since this would require an appropriate resummation of the underlying higher-loop beta functions. We postpone such an analysis to the future and describe the underlying qualitative implications, instead.

VII.2 Hierarchy of length scales

In numerical simulations of deconfined quantum transitions, similar scaling violations have been observed for various lattice models Lou et al. (2009); Pujari et al. (2013); Block et al. (2013); Powell and Chalker (2008); Charrier et al. (2008); Chen et al. (2009) and, notably, also for a range of values for $n$, i.e. $n=2,3,4$ Harada et al. (2013); Kaul (2011). Here, based on our four-loop calculations and the interpolating functions from Sec. VI.2, we will explore whether such scaling violations for $n=2,3,4$ can be consistently be ascribed to the pseudo-critical or walking behavior related to Miransky scaling, cf. Eq. (24), in the three-dimensional AH model. First, we note that our estimate of $n_{c}(3)\approx 12.2\pm 3.9$, Eq. (22), tentatively supports the scenario that $n_{c}(3)$ is of comparable order as the physically relevant $n=2,3,4$. Uncertainties in this determination are still large and it is therefore conceivable that the true $n_{c}(3)$ is closer to $n=2,3,4$. The following reasoning, however, will not rely on the precise number $n_{c}(3)$, but on the behavior of the interpolating function $n_{c}(d)$ for $d\in[2,4]$ as exhibited in Fig. 4.

The essential feature of the infrared length scale as given by the Miransky scaling, Eq. (24), is that the length scale diverges when approaching $n_{c}$ from below. This is related to the singularity in the argument $f(n)$ appearing in the exponential at $n_{c}$. Sufficiently close to $n_{c}$, the slope of $L_{\mathrm{IR}}(n)$ diverges, too, and small variations of $n$ are expected to lead to large changes in $L_{\mathrm{IR}}(n)$. This behavior is generic unless the function $f(n)$ is exceptionally flat just before it develops the singularity, e.g., like an inverse step function $\Theta^{-1}(n_{c}-n)$. Whereas, we cannot quantitatively determine the function $f(n)$ from our present analysis, our interpolating functions for $n_{c}(d)$, cf. Fig. 4, do not provide indications for such exceptionally flat behavior. We therefore expect that near $n_{c}$ a change of $n\to n\pm 1$ significantly impacts the correlation length due to its exponential behavior. On the other hand, if $n$ is further away from $n_{c}$, it is exponentially suppressed and the pseudo-critical or walking regime where, approximate scaling can be observed, should be small.

We conclude that, generically, a substantial hierarchy of characteristic correlation lengths can be expected for different values of $n$. In particular, if the true $n_{c}(3)$ lies close enough to $n=2,3,4$, we predict that the three cases should exhibit this large hierarchy of length scales, which can be probed in numerical simulations by comparing the correlation lengths for $n=2,3,4$.

VIII Conclusion and outlook

We have delivered a study of the $n$-component AH model at four loops and in $4-\epsilon$ dimensions and analysed the critical number of complex scalars $n_{c}$, below which a fixed-point collision appears and the phase transition turns from second to first order. Based on a series of assumptions on the topology of the RG flow between two and four dimensions, we have obtained a smooth interpolating function $n_{c}(d)$ and our best estimate for the three-dimensional case yields $n_{c}(3)\approx 12.2\pm 3.9$. Slightly below $n_{c}$ weakly first order behavior occurs with a correlation length which is governed by Miransky scaling, i.e. it is exponentially suppressed upon increasing the distance to the point of fixed-point collision.

Due to this exponential dependence, even small variations of $n$ are expected to induce substantial changes in the correlation length unless the argument appearing in the exponential is an exceptionally flat function of $n$. We do not find indications of such behavior and therefore predict a large hierarchy of correlation lengths for different $n$, if $n_{c}(3)$ is near these values. This is relevant to the results from different lattice studies of the corresponding deconfined phase transitions for $n=2,3,4$ and we expect that indications of this hierarchy of length scales can be measured in numerical simulations.

In the future, it will be interesting to extract more quantitative estimates for the merging line $d_{c}(n)$ employing appropriate resummations of the beta functions Folk and Holovatch (1998).

Acknowledgments.

We thank A. Nahum for correspondence and A. Läuchli for discussion. M.M.S. is supported by the Deutsche Forschungsgemeinschaft (DFG), Projektnummer 277146847 – SFB 1238 (project C04). This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 764850 (SAGEX). This work was supported by the NSERC of Canada.

Appendix A Tool chain

We calculate the wave-function and vertex renormalization constants using dimensional regularization (DREG) within the modified minimal subtraction scheme ($\overline{\text{MS}}$). In order to obtain the $Z$ factors of specific Green-functions we apply the Infrared Rearrangement procedure (IRR) introduced in Ref. Misiak and Munz (1995); Chetyrkin et al. (1998). This method allows to extract the overall UV divergence appearing in a Feynman diagram in terms of single scale massive tadpole integrals. The method is based on an exact decomposition of every massless propagator into a sum of massive propagators (of increased power) dropping irrelevant UV finite terms, only. Further it avoids contamination from any IR divergences, which in DREG in general can get regulated in terms of $1/\epsilon$ poles, too. However, the ability to calculate just the overall divergence requires a full subtraction of all existing sub-divergences via an explicit insertion of wave-function and vertex counter terms.

The number of Feynman diagrams for each Green function we have calculated is shown in Tab. 3. To factorize the flavor amplitude on diagram level, we introduced an auxiliary (non-propagating) field $\sigma$ which helps to split the original flavour-non-factorizeable four scalar vertex into a sum of products of two scalar-scalar-$\sigma$ vertices.

To be able to calculate the large number of diagrams we first use QGRAF Nogueira (1993) to generate complete sets of them. Then q2e and exp Seidensticker (1999) are used to map all Feynman diagrams on one-scale massive tadpole integral topologies and to generate diagram source files. Next, FORM Vermaseren (2000); Kuipers et al. (2013); Ruijl et al. (2017) is used to process the diagram source files. It rewrites the amplitudes in terms of massive tadpole integrals with different powers of propagators. Finally it replaces all integrals by their tabulated reduction to a set of known master integrals Czakon (2005). The reduction to master integrals is performed by Crusher Marquard and Seidel (shed) and relies on integration-by-parts identities relating integrals with different propagator powers through a system of coupled equations to each other Chetyrkin and Tkachov (1981). The system of equations can be solved with the Laporta algorithm Laporta (2000) such that all appearing integrals can be written in terms of a linear combination of a finite number of master integrals.

The final renormalization step during which all $Z_{x}$-factors – with $x\in\{\phi,\phi^{2},\phi^{4},A,\phi A\phi\}$ – are obtained is carried out by a FORM program, which ensures the insertion of all possible counter terms in an automated way.

We explicitly computed the vertex renormalization constant $Z_{\phi A\phi}$ in order to check that the Ward identity $Z_{\phi}/Z_{\phi A\phi}=1$ holds and the charge renormalization is solely determined by the wave-function renormalization constant $Z_{A}$ of the photon. The mass square vertex renormalization $Z_{\phi^{2}}$ of the scalar was obtained via a single zero momentum mass square operator insertion into the two-point function of the scalar.

The IRR procedure requires to keep track of the degree of the overall UV divergence via a global divergence power counting variable. Besides the physical renormalization constants already mentioned, one has to insert artificial mass counter terms for the photon and scalar starting at the two-loop level. At the one loop level the IRR introduces artificial divergent mass terms proportional to the unphysical regulator mass square, which can be easily dropped. However, at the two-loop level such terms do get multiplied with negative powers of the regulator mass and thus lose their “mass label” and mix with physical divergences. In order to prevent such a mixing one has to subtract these terms via explicit mass counter term insertions. As a welcome check we verified that on the removal of the divergence power counting variable, each of the two artificial mass counter terms adds up to zero, order by order in the loop expansion.

We also kept the full dependence on the gauge parameter $\xi_{P}=1-\xi$ appearing in the photon propagator. As welcome check we could verify that $\xi_{P}$ fully cancels in all beta function and anomalous dimensions except for $\gamma_{\phi},\gamma_{\phi^{2}},\gamma_{\phi^{4}}$ and $\gamma_{\phi A\phi}$, where it is in principle known that $\xi_{P}$ appears only once at the one loop order Kißler and Kreimer (2017).

We further could verify that upon adopting the $\hat{G}$ scheme, which involves the redefinition of odd zeta values in terms of an epsilon expansion containing even zeta values – in our case just $\hat{\zeta}_{3}=\zeta_{3}-\tfrac{3}{4}\zeta_{4}\epsilon$ – all $\zeta_{4}$ terms vanished in our results for the $Z$-factors. So all $\zeta_{4}$ terms in the $Z$-factors can be fully reconstructed from the existing $\zeta_{3}$ ones following the “no-$\pi$ theorem” Baikov and Chetyrkin (2018); Kotikov and Teber (2019).

Appendix B Four-loop renormalization constants

The beta functions up to three loops are given in Sec. III.1. Here, we provide the missing terms at four loops and give the full expressions for the anomalous dimensions and the flow of the gauge fixing parameter. The four-loop contributions to the beta functions read

[TABLE]

The field anomalous dimensions are defined through the relation $\gamma_{x}=\frac{d\ln Z_{x}}{d\ln b}$, for $x\in\{\phi,\phi^{2},A_{\mu}\}$. Schematically, order by order, they can be written as $\gamma_{x}=\sum_{i=1}^{L}\gamma_{x}^{(i\ell)}$. Up to three loops, the contributions to $\gamma_{\phi^{2}}$ read

[TABLE]

The contributions to $\gamma_{\phi}$ are

[TABLE]

[TABLE]

and for $\gamma_{A}$ we obtain

[TABLE]

The contributions at four loops to these field anomalous dimensions are, again, more lengthy. They read

[TABLE]

Finally, we also give the corresponding expression for the renormalization group flow of the gauge fixing parameter

[TABLE]

Appendix C Padé and Padé-Borel approximants

We explore various resummations for the series

[TABLE]

from Eq. (18). The Padé approximants are defined as

[TABLE]

with $L-1=m+n$. The real-valued coefficients $\{a_{i},b_{i}\}$ are uniquely determined by matching its Taylor series to the original expansion, e.g., Eq. (42), order by order.

Assuming the series is Borel summable, we also compute the corresponding Padé-Borel approximants. The Borel transform is defined as

[TABLE]

Using $B_{n_{c}}$, one analytically continues to the real axis by calculating the Borel sum

[TABLE]

Approximating the Borel transform by a Padé approximant, the Padé-Borel resummed series is obtained.

For the interpolation between the epsilon expansions at lower and upper critical dimension we construct the two-sided Padé approximants. To that end, the plain Padé approximants are generalized making the ansatz

[TABLE]

Herein, the coefficients $\{a_{i}\}$ and $\{b_{j}\}$ are fixed such that the function matches the epsilon expansions near $d=2$ and $d=4$. This leads to the relations

[TABLE]

with $n_{c,k}$ from Eq. (42) and, again, the variable $a>0$ describes the unkown slope of $n_{c}$ close to $d=2$. In the text we provide estimates for $n_{c}$ at $d=3$ by setting $a=1$ and in Fig. 5 we show the $a$ dependence of two-sided Padé evaluated at $d=3$. Additionally, we compiled the estimates of the highest order two-sided Padés for different $a$ in Tab. 4.

Appendix D Borel resummation algorithm

Borel resummation with conformal mapping Zinn-Justin (1996) facilitates high-accuracy determination of critical exponents in 3$d$ $O(n)$ symmetric $\phi^{4}$ theories from the $(4-\epsilon)$ expansion Kompaniets and Panzer (2017). For the present application, we make the assumption that the AH model’s asymptotic behavior is qualitatively the same as the one from the ungauged scalar models. More explicitly,we assume that the epsilon expansion for the critical number of components follows a formal power series with factorially increasing coefficients, i.e.

[TABLE]

for large $k$. The parameter $a$ was determined for the AH model to $a\approx 0.0808$ Itzykson et al. (1977) . The generalized Borel(-Leroy) transform Kleinert and Schulte-Frohlinde (1991) of this series is calculated as in Eq. (44) by canceling the factorial growth of the coefficients

[TABLE]

The coefficients behave like $B_{k}^{b}\sim(-a)^{k}$ and therefore follow a geometric series which can written as rational function with a pole at $\epsilon=-1/a$

[TABLE]

This increases the originally vanishing radius of convergence to in a circle in the complex plane with $|\epsilon|<1/a$. The assumption of Borel summability now implies that the Borel transform can be analytically continued to the positive real axis such that the Borel sum (45) – with an extra factor $t^{b}$ in the integrand – is convergent and leads to the correct estimate for $\epsilon=1$.

In order to incorporate the large-order behavior of the asymptotic series, in the next step the Borel transform is analytically continued with a conformal transformation to the real axis, i.e.

[TABLE]

This maps the cut of the negative real axis $(-\infty,-1/a)$ of the Borel transform to the unit circle in the variable $w$. Expanding the Borel transform in $w$ provides a convergent series in the integration domain of the Borel sum.

As a further improvment for finite order truncated Borel sums, we introduce the parameter $\lambda$ Kompaniets and Panzer (2017); Vladimirov et al. (1979),

[TABLE]

which restores the actual large $\epsilon$ behavior $\sim\epsilon^{\lambda}$ of the Borel transform Kompaniets and Panzer (2017). In the analysis presented here, the parameter $\lambda$ is used to improve the sensitivity of the resummation algorithm, see below. On top of that, a homographic transformation with shifting variable $q$

[TABLE]

can be used to optimize the numerical stability. The resummed series is therefore determined by

[TABLE]

The set of resummation parameters $a$, $b$, $\lambda$ and $q$ are optimized by measuring the sensitivity of the summation to a parameter change. For a function $F$ the variation upon changing one of its parameters $x$ in a range $x\in[x_{0},x_{0}+\Delta]$ can be defined as

[TABLE]

For a function $F(x)$, which is weakly varying within a plateau of length $\Delta$, $\mathcal{S}_{x}$ is small. Therefore, we can quantify the sensitivity as the extension of the plateaus. At the highest available order $\mathcal{O}(\epsilon^{3})$, we find the spreads

[TABLE]

The optimal set of parameters is identified by making use of both, the “principle of minimal sensitity” and “principle of fastest convergence” Kompaniets and Panzer (2017).

The different variations are collected in an error estimate Kompaniets and Panzer (2017); Ihrig et al. (2018b),

[TABLE]

whichs is minimized by the wanted set of parameters. This includes als the sensitivities at different orders $L$ and search for a minimum in the dependence on $L$.

Explicitly, we scan the parameter space in the range $(b,\lambda,q)\in[0,90]\times[0,4]\times[0,0.4]$ in steps of $\delta_{b}=1$, $\delta_{\lambda}=0.02,\delta_{q}=0.02$ and compute $E_{{n_{c}},L}^{b,\lambda,q}$ for each parameter set, see Fig. 6. The gobal minimum $\overline{E}\approx 15.079$ of this scan marks the apparently best set of resummation parameters and is found at $(b,\lambda,q)\leavevmode\nobreak\ =\leavevmode\nobreak\ (87.0.66,0.0)$. However, this minimum is not always sharp and there are other sets of parameters which are almost equally likely. Therefore, we compute a weighted mean of all points which lie below a relative error of $E_{i}/\overline{E}<3$ as in Ref. Ihrig et al. (2018b). From the error estimate at the minimum we compile the error for the resummation of $n_{c}$ at $d=3$ with error as $3\overline{E}$ to $\tilde{n}_{c}(d=3)\approx-52\pm 45$.

Appendix E Relation to the $\mathbb{CP}^{n-1}$ model

The AH model is related to the $\mathbb{CP}^{n-1}$ model with action $S=\frac{1}{2t}\int\mathrm{d}^{d}x\left(\partial_{\mu}\bar{\phi}^{a}_{i}\partial_{\mu}\phi^{a}_{i}+(\bar{\phi}^{a}_{i}\partial_{\mu}\phi^{a}_{i})(\bar{\phi}^{b}_{i}\partial_{\mu}\phi^{b}_{i})\right)$ and constraint $\sum_{a=1}^{n}\bar{\phi}^{a}_{i}\phi^{a}_{j}=\delta_{ij}$ Hikami (1979); March-Russell (1992). The relation can be established by considering the non-linear sigma model (NL$\sigma$M) defined on the Grassmannian manifold of $U(n)/[U(n-p)\times U(p)]$. For $p=1$, the symmetric space is isomorphic to $\mathbb{CP}^{n-1}$ Hikami (1981). The beta function $\beta_{t}$ for general $p$ was computed in $d=2+\epsilon$ to four-loop order Hikami (1981, 1983); Wegner (1989). It features an IR-unstable fixed point $\beta_{t}(t_{c})=0$ reading

[TABLE]

which is real valued for all $n>0$. The associated correlation-length exponent $1/\nu$ is

[TABLE]

A further connection can be drawn between the above NL$\sigma$M with target space $U(n)/[U(n-p)\times U(p)]$ and a $SU(p)\times U(1)$ gauge theory with Lagrangian Hikami (1980); March-Russell (1992)

[TABLE]

Here, $D_{\mu}=\partial_{\mu}-\mathrm{i}eA_{\mu}-\mathrm{i}gB_{\mu}$ is the covariant derivative and $F^{2}$ ($G^{2}$) are kinetic terms for the $U(1)$ ($SU(p)$) gauge fields $A_{\mu}$ ($B_{\mu}$). Eq. (56) represents a generalization of the AH model and its beta functions can be deduced from it upon replacing $n\rightarrow pn$ supplemented by the beta function for the $SU(p)$ gauge field March-Russell (1992). For $p=1$ we recover the AH model beta functions. The large-$n$ expansion for Eq. (56) gives correlation-length exponent Hikami (1980)

[TABLE]

and anomalous dimension

[TABLE]

As a check, we expanded the critical exponents of the AH model to first order in $1/n$ and find order by order agreement in $\epsilon$. Further, the large-$n$ inverse correlation-length exponent expanded around $d=2+\epsilon$ agrees with the exponent from the $\mathbb{CP}^{n-1}$ model to all availabe orders in $\epsilon$. This corroborates the evidence that both theories lie in the same universality class for large $n$.