<figure><img src="image_1.png" /></figure>

# Bootstrapping frustrated magnets: the fate of the chiral $O(N) \times O(2)$ universality class

Marten Reehorst, Slava Rychkov, Benoit Sirois, Balt C van Rees

► To cite this version:

Marten Reehorst, Slava Rychkov, Benoit Sirois, Balt C van Rees. Bootstrapping frustrated magnets: the fate of the chiral $O(N) \times O(2)$ universality class. SciPost Physics, 2025, 18 (2), pp.060. 10.21468/SciPostPhys.18.2.060. hal-04609623

HAL Id: hal-04609623

https://hal.science/hal-04609623v1

Submitted on 26 Aug 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
Bootstrapping frustrated magnets:
The fate of the chiral O(N) × O(2) universality class

Marten Reehorst¹,²*, Slava Rychkov³†, Benoit Sirois³,⁴‡ and Balt C. van Rees¹°

1 CPHT, CNRS, Ecole Polytechnique, Institut Polytechnique de Paris, Palaiseau, France

2 Department of Mathematics, King's College London, Strand, London, WC2R 2LS, UK

3 Institut des Hautes Études Scientifiques, 91440 Bures-sur-Yvette, France

4 Laboratoire de Physique de l'École normale supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université de Paris, F-75005 Paris, France

* marten.reehorst@kcl.ac.uk, † slava@ihes.fr, ‡ benoit.sirois@umontreal.ca,

◦ balt.van-rees@polytechnique.edu

Abstract

We study multiscalar theories with O(N) × O(2) symmetry. These models have a sta-
ble fixed point in d dimensions if N is greater than some critical value N_c(d). Previous
estimates of this critical value from perturbative and non-perturbative renormalization
group methods have produced mutually incompatible results. We use numerical confor-
mal bootstrap methods to constrain N_c(d) for 3 ≤ d < 4. Our results show that N_c > 3.78
for d = 3. This favors the scenario that the physically relevant models with N = 2, 3 in
d = 3 do not have a stable fixed point, indicating a first-order transition. Our result ex-
emplifies how conformal windows can be rigorously constrained with modern numerical
bootstrap algorithms.

<figure><img src="image_3.png" /></figure>

Copyright M. Reehorst *et al.*

This work is licensed under the Creative Commons

Attribution 4.0 International License.

Published by the SciPost Foundation.

Received 2024-07-04

Accepted 2025-02-03

Published 2025-02-19

doi:10.21468/SciPostPhys.18.2.060

<figure><img src="image_4.png" /></figure>

Check for
updates

Contents

1 Introduction 2

2 CFT setup 6

2.1 Ideas about separating $C_+$ not used in this work 7

2.2 Unitarity assumption 8

2.3 Gap assumptions 9

3 Conformal bootstrap analysis 10

3.1 Finding islands for large $N$ 10

3.1.1 $d = 3.8$ and $N = 20$ 11

3.1.2 $d = 3$ and $N = 8$ 13

3.2 $N_c(d)$ from the disappearance of the $C_+$ island 13

3.3 Tracking the navigator minima 14

3.3.1 $d = 3.8$ 14

3.3.2 $d = 3$ 16
<table><tbody><tr><td>3.4</td><td>How does N<sub>c</sub><sup>CB</sup>(3) depend on the choice of gapST<sub>0</sub>?</td><td>17</td></tr><tr><td>3.4.1</td><td>Hessian and derivative information</td><td>18</td></tr><tr><td>3.4.2</td><td>Left endpoint behavior</td><td>19</td></tr><tr><td>3.4.3</td><td>Right endpoint behavior</td><td>20</td></tr><tr><td>4</td><td>Review of prior bootstrap studies</td><td>21</td></tr><tr><td>5</td><td>Comparison to prior work, conclusions and outlook</td><td>22</td></tr><tr><td>A</td><td>Frustrated magnets</td><td>23</td></tr><tr><td>B</td><td>Critical remarks about the focus fixed points</td><td>25</td></tr><tr><td>B.1</td><td>Problems with unitarity</td><td>25</td></tr><tr><td>B.2</td><td>Merger curve and bifurcation theory</td><td>25</td></tr><tr><td>B.3</td><td>Further problems of focus fixed points</td><td>27</td></tr><tr><td>C</td><td>Conformal invariance of multiscalar models</td><td>28</td></tr><tr><td>D</td><td>1/N expansion</td><td>29</td></tr><tr><td>E</td><td>ε-expansion results</td><td>30</td></tr><tr><td>E.1</td><td>N<sub>c</sub>(d)</td><td>30</td></tr><tr><td>E.2</td><td>Scaling dimensions</td><td>31</td></tr><tr><td>F</td><td>Numerical bootstrap details</td><td>33</td></tr><tr><td>F.1</td><td>The navigator method</td><td>33</td></tr><tr><td>F.2</td><td>Details for Figs. 5 and 6</td><td>33</td></tr><tr><td>F.3</td><td>The skydive algorithm</td><td>34</td></tr><tr><td>F.4</td><td>Constructing the navigator function</td><td>35</td></tr><tr><td>F.5</td><td>Minimum-following for Fig. 3</td><td>36</td></tr><tr><td>F.6</td><td>Gap sensitivity</td><td>37</td></tr><tr><td colspan="2">References</td><td>37</td></tr></tbody></table>

# 1 Introduction

It has been known for decades now [1] that a multiscalar theory with $O(N) \times O(2)$ global symmetry describes the critical modes of N-component classical stacked triangular antiferromagnets (STAs) and helimagnets (see App. A for a brief review). The Hamiltonian is [1]:

$$ \mathcal{H} = \partial_{\mu}\varphi \cdot \partial_{\mu}\varphi^{*} + m^{2}\varphi \cdot \varphi^{*} + u(\varphi \cdot \varphi^{*})^{2} + v(\varphi \cdot \varphi)(\varphi^{*} \cdot \varphi^{*}), \quad (1) $$

where $\varphi$ is an N-component complex field.<fnref n="1" /> The basic question is:

For the physical values $N = 2, 3$ and in three dimensions, does the model (1) have a stable fixed point of the Renormalization Group (RG) flow? (2)

If yes, for materials within the basin of attraction of this fixed point, experiments should observe second-order phase transitions. In the opposite case, all phase transitions should be first order.

<footnote n="1">¹In this representation, $O(N)$ acts as $\varphi \to R\varphi$, while $O(2)$ acts as $\varphi \to e^{i\alpha}\varphi$ and $\varphi \leftrightarrow \varphi^*$.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 1: RG flow and fixed point of the model (1) in $d = 4 - \epsilon$ dimensions, for $N > N_c(d)$ (left) and $N < N_c(d)$ (right). Figure from [1].</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 2: Comparison between NPRG, fixed-dimensional RG and the $\epsilon$-expansion. Figure from [10].</figcaption></figure>

Although the model and the question are simple enough, the subject has been the arena of a controversy [2–4], with various methods giving contradictory results, as we now review.

A time-honored way to study the critical properties of multiscalar models is through the $\epsilon$-expansion around $d = 4$ [5]. Applied to our model, the $\epsilon$-expansion predicts the existence of a critical curve $N_c(d)$ above which a stable fixed point exists. This stable fixed point is called chiral and we will denote it as $C_+$ in Fig. 1. At $N = N_c(d)$, the chiral fixed point merges and annihilates with an unstable fixed point $C_-$, called antichiral.<fnref n="2" /> Below this collision, only two fixed points remain, both unstable: the Gaussian $G$ and the $O(2N)$ invariant “Heisenberg” fixed point $\mathcal{H}$, which has $v = 0$ in (1). These RG flows at $N > N_c$ and $N < N_c$ are illustrated in Fig. 1. Perturbation theory allows to compute $N_c(d = 4 - \epsilon)$ as an asymptotic series expansion in $\epsilon$, known up to 6 loops [6].<fnref n="3" /> Using various resummation techniques, [6] concludes that $N_c(3) = 5.96(19)$, well above the physical values. To summarize, the $\epsilon$-expansion predicts first-order phase transitions for $N = 2, 3$.<fnref n="4" />

<footnote n="2">²The “chiral fixed point” refers to chirality present in low-T ground states of STAs, see App. A. The “antichiral fixed point” refers not to opposite chirality of the ordered phase (both chiralities may be present), but to the possibility of annihilation with the chiral fixed point.</footnote>

<footnote n="3">³The very first RG studies of the bifundamental scalar model with $O(M) \times O(N)$ symmetry, at the lowest order in the $\epsilon$-expansion, were performed in mid-1970’s in [7–9].</footnote>

<footnote n="4">⁴The $\epsilon$-expansion also predicts further changes to the phase diagram below $N = 2$, see [1] and footnote 22 for more details. We will not be concerned with these transitions.</footnote>
A competing theoretical view is that of the perturbative RG performed directly in $d = 3$ [11–14] (fixed-dimensional (FD) RG). For this method, $\beta$-functions are expanded in powers of the quartic couplings of the model in the dimension of interest. The $\beta$-functions are then resummed before looking for a fixed point. This may be contrasted to the $\epsilon$-expansion, where fixed points and critical exponents are found in the perturbative regime as series-expansions in $\epsilon$, and resummation is performed at the last step when extrapolating the exponents to the desired dimension. For the simpler $O(N)$ model these two methods agree, however they do not agree very well for $O(N) \times O(2)$. The fixed-dimensional RG predicted the region of the $(d, N)$ plane where a stable chiral fixed point exists to be bounded by an S-shaped curve which is not a single-valued function $N_c(d)$, see Fig. 2. Although the boundary agrees near $d = 4$ with the $\epsilon$-expansion, it deviates and takes a couple of turns as $N$ is lowered. The resulting region where a stable chiral fixed point exists is significantly larger than for the $\epsilon$-expansion. Furthermore, the new stable fixed points found within FD schemes but not in the $\epsilon$-expansion turn out to be “of focus type”, meaning that they have complex correction-to-scaling exponents. This is in conflict with the unitarity of the model, requires unlikely crossings in the operator spectrum, and contradicts basic principles of renormalization group such as the gradient flow property.<fnref n="5" /> See Appendix B for a more detailed discussion.

Unfortunately, experimental data and Monte Carlo simulations do little to clear up the picture. Some Monte Carlo results show clear signs of first-order transitions [20, 21] while some claim to confirm the existence of focus-type fixed points [22].<fnref n="6" /> Experiments are plenti- ful, but the consensus is that there is such a large discrepancy of critical exponents in the cases where a second-order phase transition is believed to happen, that it is impossible to discern between true second-order and weak first-order behavior [2, 3]. We are left with an unfortu- nate situation where the most standard methods fail to give a coherent picture, and thus new perspectives from alternative methods are necessary.

In this respect, calculations have also been performed using Functional, or Nonperturbative RG (NPRG). The results rather confirm the picture of the $\epsilon$-expansion [3, 23], producing a curve $N_c(d)$ in a reasonably good agreement with the $\epsilon$-expansion all the way down to three dimensions, see Fig. 2.

In this paper we will approach the problem using another nonperturbative method - the numerical conformal bootstrap.<fnref n="7" /> Since its inception [26], thanks to many technical and con- ceptual improvements such as [27–38], the method has achieved several important results, notably the high precision determination of the critical exponents of the 3D Ising, $O(2)$ and $O(3)$ universality classes [39, 40]. It has already contributed to the resolution of some puzzles such as the liquid helium heat capacity anomaly [41] and the cubic instability of the Heisen- berg magnets [42]. See [43–46] for reviews and results for other models, including with gauge fields and fermions.

Conformal bootstrap relies on the basic fact that a fixed point of the renormalization group flow is normally described by a unitary conformal field theory (see Appendix C), as well as on some assumptions about gaps in the operator spectrum. This represents a very different set of assumptions as compared to the renormalization group studies. In particular, the confor- mal bootstrap analysis does not suffer from truncation and resummation ambiguities inherent

<footnote n="5">⁵The RG flow of multiscalar models is widely expected to be a gradient flow. This has been long known to be true to three loops [15, 16], and recently has been verified at five loops (and even six loops with some assumptions) [17], using the six-loop beta function results from [18, 19]. Refs. [11–14] resummed each component of the beta-function separately - a procedure which may not have preserved the gradient property.</footnote>

<footnote n="6">⁶It should be noted that the lattice model in [22] is reflection positive, rigorously excluding complex critical exponents (see Appendix B.1). This is in stark conflict with the complex correction-to-scaling exponent $\omega = 0.1^{+0.4}_{-0.05} + i 0.7^{+0.1}_{-0.4}$ they report.</footnote>

<footnote n="7">⁷Previously, conformal bootstrap has already been used to study this problem in Refs. [24, 25]. There are some methodological differences between our work and the approaches of these prior works, whose review is postponed to Section 4.</footnote>
to perturbative and nonperturbative RG techniques mentioned above. It is therefore partic-
ularly interesting to inquire what the conformal bootstrap has to say about the $O(N) \times O(2)$
controversy.

For the purposes of conformal bootstrap analysis, question (2) is reformulated as:

For $N = 2,3$ and in $d = 3$, is there a unitary conformal field theory (CFT) with $O(N) \times O(2)$ symmetry and with one relevant singlet scalar operator in its spectrum?

The latter condition is how the fixed point stability is translated into the CFT language, the only relevant singlet being the mass term in (1). Further conditions on the CFT to allow its identification with the RG fixed point of (1), such as a low-lying bifundamental scalar, will be discussed in Section 2.

Answering (3) directly is out of reach with the current state of the art. Instead we will address the question:

What is the shape of the critical curve defining $N_c(d)$ as we lower the spacetime dimension $d$ from 4 to 3?

Recall that we have perturbative control near $d = 4$ (as well as at large $N$ for any $d$), so the existence of the curve is uncontroversial. But could it be that the $\epsilon$-expansion and NPRG are wrong, and the boundary of the region where the stable fixed point exists is indeed S-shaped as in FD-studies [11–14]? We would like to be open-minded about this, even though another aspect of the FD-studies - focus fixed points - is ruled out as contradicting above-mentioned general principles.

Our strategy for investigating (4) is as follows. After making some minimal assumptions we will find a small allowed island in the space of CFT data that matches all the expectations for the $C_+$ fixed point for large $N$. As we then lower $N$ this island varies in size, but eventually it shrinks and disappears entirely at some numerically determined value $N_c^{CB}(d)$. So if one believes that the $C_+$ fixed point was indeed located inside this island then the conclusion is unavoidably that $N_c^{CB}(d)$ provides a rigorous lower bound for $N_c(d)$.<fnref n="8" />

In Fig. 3 we show $N_c^{CB}(d)$ for $3 \le d \le 3.8$. The main salient features of this plot are as follows:

* The bootstrap curve $N_c^{CB}(d)$ shows the same single-valued monotonic behavior as the $\epsilon$-expansion and NPRG curves. There is no sign of a turnaround similar to the one reported in FD-studies.

* At $d = 3$, we have $N_c^{CB}(3) = 3.78$. This is quite a bit lower than the $\epsilon$-expansion and NPRG predictions, but above the physical values $N = 2, 3$.

The use of a vanishing island is the most natural way to determine the endpoint of a conformal window with the numerical conformal bootstrap. We expect it to be useful also in many other cases.

The paper is structured as follows. In Section 2 we describe our bootstrap setup, detailing in particular in Section 2.3 our choice of gap assumptions in important symmetry channels. We then move on to the detailed description of our bootstrap study and results in Section 3. In Section 4 we review previous bootstrap studies of the same problem [24, 25]. In Section 5 we conclude. In particular, we discuss how our results about the shape of the boundary curve relate to the question (3) about the existence of $N = 2, 3$ CFTs. We also have several appendices. Appendix A is dedicated to introducing frustrated noncollinear magnets, Appendix B critically

<footnote n="8">⁸With more computational resources we would obtain a smaller island and therefore potentially a disappearance for $N$ greater than $N_c^{CB}(d)$. This is why $N_c^{CB}(d)$ is a lower bound.</footnote>
<figure><img src="image_7.png" /><figcaption>Figure 3: This plot, the main result of our paper, reports $N_c^{CB}(d)$ extracted by the conformal bootstrap (blue dots, joined by an interpolating blue curve). We also show the results of NPRG [3, Fig. 9] (red asterisks) and of the $\epsilon$-expansion (green band, whose width is the uncertainty obtained as in Appendix E.1). Conformal bootstrap results used spectrum gap assumptions (7), (9), (10) and the derivative order $\Lambda = 31$.</figcaption></figure>

reviews the controversial notion of focus fixed points, and Appendix C reviews why scale invariance is expected to be enhanced to conformal invariance for our models. Appendix D is devoted to the calculation of the scaling dimension at $O(1/N)$ of a certain crucial operator, and Appendix E discusses results obtained in the $\epsilon$-expansion. Details of the numerical methods are collected in Appendix F.

In the ancillary mathematica notebook *EFM_Spectrum_Along_Our_Curves.nb*, we give the full spectrum and OPE coefficients extracted through the extremal functional method [31] along the $\Lambda = 31$ curves of Figs. 3, 8 and 9. The notebook explains the formatting of this data and provides example plots and tables to help visualize it.

We will be assuming that the reader is reasonably familiar with the conformal bootstrap philosophy and standard technology. We refer to the reviews [43–46] and to the pedagogical lecture notes [47,48]. We will be also using more modern tools such as the navigator function [49] and the skydiving algorithm [50].

## 2 CFT setup

We are interested in the existence of the stable RG fixed point $C_+$ of (1) with $O(N) \times O(2)$ global symmetry. The scale invariance of the fixed point is expected to be enhanced to conformal invariance, as we briefly review in Appendix C. We can therefore use the language of CFTs and the tools of the numerical conformal bootstrap. Instead of focusing only on the physically interesting cases of $d = 3$ and $N = 2, 3$, we will explore the wider range $3 \le d \le 4$ and $N \ge 2$. This includes the regimes $d \to 4$ and $N \to \infty$ which are controlled by the $\epsilon$-expansion and the $1/N$ expansion.

Primary CFT operators are characterized by dimension $\Delta$, spin $\ell$, and by a global symmetry irrep $\mathcal{R}$, which in our case are labeled XY, where X is an irrep of $O(N)$ and Y of $O(2)$. In the rest of this paper we will be using real notation in which the complex vector field $\varphi$ of (1) is repackaged as a real matrix field $\phi = \phi_{ai}$ ($a = 1, \dots, N, i = 1, 2$), in the bifundamental (VV) representation of $O(N) \times O(2)^9$.

<footnote n="9">⁹The complex field $\varphi$ in (1) is $\varphi_a = \phi_{a1} + i\phi_{a2}$.</footnote>
Our study will use crossing symmetry constraints for the correlation functions

$$ \langle \phi \phi \phi \phi \rangle, \quad \langle \phi \phi ss \rangle, \quad \langle ssss \rangle, \qquad (5) $$

where $\phi$ and $s$ are the lowest-lying scalar primaries in the bifundamental (VV) and the total singlet (SS) irreps. In the Lagrangian description we have $s = \phi^2$, but we will not use this notation. The operator $s$ is relevant. All subsequent total singlet scalars have to be irrelevant – this is the RG fixed point stability condition.

Operators in the operator product expansion (OPE) $\phi \times s$ transform as VV. In the Lagrangian description, the VV scalars after $\phi$ are the cubic operators $\phi_{ai}\phi^2$ and $\phi_{aj}\phi_{bi}\phi_{bj}$. One linear combination of these, the derivative of the fixed point potential with respect to $\phi_{ai}$, is a descendant of $\phi$.<fnref n="10" /> The orthogonal linear combination is a primary, denoted $\phi'$.<fnref n="11" /> Made out of three scalar fields in the perturbative notation, this primary is expected to be relevant, while the subsequent primaries in the VV sector, made of 5 or more $\phi$'s, are expected to be irrelevant.

Operators in the $\phi \times \phi$ OPE have $X, Y = S, T, A$, where $T, A$ are the symmetric traceless and antisymmetric 2-index tensor irreps. The important operators in this OPE include the stress tensor $T_{\mu\nu}$, which is a spin-2 SS primary of dimension $\Delta = d$, and the conserved currents of the $O(N)$ and $O(2)$ symmetries, which are spin-1 primaries transforming as AS and SA respectively, of dimension $\Delta = d - 1$. Below we will specify gap assumptions above these operators.

A recurrent theme of our study will be how to distinguish our CFT of interest $C_+$ from the Heisenberg CFT $\mathcal{H}$, which has a larger symmetry $O(2N)$. This is a crucial issue since $\mathcal{H}$ is present below $N_c(d)$, and thus has the potential to pollute our bootstrap analysis. A particularly elegant way to disentangle the two theories in the correlator system Eq. (5) turns out to be the ST scalar channel in the $\phi \times \phi$ OPE, as we now proceed to discuss.

In model (1), the lowest operator in this channel is $O_{ST} = \phi_{ai}\phi_{aj}$-trace. At leading order in the $1/N$ expansion it has dimension close to 2 (in any $d$) for $C_+$, whereas it has dimension $d-2$ in $\mathcal{H}$. Since this difference is fairly big, a gap assumption in this channel might allow us to distinguish between both theories. In practice we have enough computational power to assume the existence of $O_{ST}$ and impose a gap to the subleading primary instead. This operator also has different large-$N$ behavior in $\mathcal{H}$ and $C_+$:

$$ \Delta_{ST'} = \begin{cases} 4 + O(1/N), & \text{in } C_+, \\ d + O(1/N), & \text{in } \mathcal{H}. \end{cases} \qquad (6) $$

This is explained in more detail in Appendix D, where also the $O(1/N)$ terms are given. We see that, for $d < 4$ and for $N$ sufficiently large, $\Delta_{ST'}$ is expected to be larger in $C_+$ than in $\mathcal{H}$. The simplest thing to do would be to put a carefully chosen gap above $O_{ST}$, chosen precisely so that we exclude $\mathcal{H}$ but include $C_+$. It turns out that we can do even better: with a smaller gap in the ST channel, which allows both $\mathcal{H}$ and $C_+$, we are able to isolate both theories into separate islands. This will be enough for all practical purposes because we can track the disappearance of the $C_+$ fixed point without interference from the $\mathcal{H}$ fixed point. The usefulness of the ST channel for this purpose is one of the main discoveries of our work.

To summarize, in Table 1 we list operators treated in our study as isolated. The $\phi$ and $s$ appear as both external and internal (i.e. exchanged) operators, the rest only as internal. As we explain in more detail below, we will numerically explore the three-dimensional parameter space corresponding to the varying scaling dimensions $\Delta_{\phi}$, $\Delta_{SS}$ and $\Delta_{ST}$.

## 2.1 Ideas about separating $C_+$ not used in this work

*Remark 2.1.* Another way to distinguish $C_+$ from $\mathcal{H}$ would be as follows. Primaries of $\mathcal{H}$ transform in irreps of $O(2N)$. Decomposing these irreps under $O(N) \times O(2)$, each primary of $\mathcal{H}$

<footnote n="10">¹⁰Just like $\phi^3$ in the Wilson-Fisher fixed point is a descendant of $\phi$ [51].</footnote>

<footnote n="11">¹¹As usual, we denote by $O', O'', \dots$ the subsequent primaries having the same quantum numbers as $O$.</footnote>
<table><caption>Table 1: Operators treated in our study as isolated; ext.=external, int.=internal.</caption><thead><tr><th>name</th><th>ℓ</th><th>R</th><th>Δ</th><th>Note</th></tr></thead><tbody><tr><td>φ</td><td>0</td><td>VV</td><td></td><td>ext./int.</td></tr><tr><td>s</td><td>0</td><td>SS</td><td></td><td>ext./int.</td></tr><tr><td>T<sub>μν</sub></td><td>2</td><td>SS</td><td>d</td><td>int.</td></tr><tr><td>J<sup>O(N)</sup><sub>μ</sub></td><td>1</td><td>AS</td><td>d − 1</td><td>int.</td></tr><tr><td>J<sup>O(2)</sup><sub>μ</sub></td><td>1</td><td>SA</td><td>d − 1</td><td>int.</td></tr><tr><td>O<sub>ST</sub></td><td>0</td><td>ST</td><td></td><td>int.</td></tr></tbody></table>

gives rise to several primaries in different $O(N) \times O(2)$ irreps having exactly the same scaling dimension. These exact degeneracies are not expected in $C_+$. For example, the conserved currents of $O(2N)$ would give, upon reduction under $O(N) \times O(2)$, conserved spin-1 operators in the AT and TA irreps. On the other hand, conserved currents in these irreps are not expected in $C_+$. We could have imposed small gaps above the unitarity bound in the AT and TA spin-1 channels, with the goal of excluding $\mathcal{H}$ and keeping $C_+$. In this work we will not use these gaps, because instead of excluding $\mathcal{H}$, we will be able to isolate $\mathcal{H}$ and $C_+$ into two separate islands, thanks to the ST channel gap assumption.

**Remark 2.2.** Here is one more idea which we cannot rely on in this work, but which could be useful in future studies.<fnref n="12" /> We know that $C_+$ is stable and $\mathcal{H}$ is unstable, hence $\mathcal{H}$ contains one more relevant SS scalar. This operator can be written as $O = h_{IJKL}W^{IJKL}$ where $W$ is a rank-4 symmetric traceless primary of $O(2N)$ and $h_{IJKL}$ is a tensor which breaks $O(2N)$ to $O(N) \times O(2)$. Can we use the existence of $O$ to distinguish $\mathcal{H}$ from $C_+$? Unfortunately, in our setup this will not work, because the $O(2N)$ selection rules of $\mathcal{H}$ preclude the appearance of $O$ in the OPEs $\phi \times \phi, s \times s$ to which we are sensitive. To be sensitive to $O$, we would have to enlarge the setup by including external fields in $T$ representations of $O(N)$ and/or $O(2)$. This is left for the future.

## 2.2 Unitarity assumption

The unitarity assumption about the CFT is going to play an important role in our numerical conformal bootstrap analysis. Namely we will be assuming, as usual, reality of scaling dimensions $\Delta_i$ and of the OPE coefficients $f_{ijk}$, and the unitarity bounds on $\Delta_i$. Strictly speaking, these constraints are satisfied only for integer $N$ and $d$. Indeed it is known that CFTs for non-integer $N$ [52] and non-integer $d$ [53, 54] are not unitary. However, one can argue [54] that violations of unitarity for $d \ge 2$ are restricted to the sector of high-dimension operators, whose exchanges contribute exponentially little to CFT four-point functions of low-dimension operators that we will be studying. At present these effects are likely below the precision of numerical conformal bootstrap algorithms. In agreement with this intuition, the Ising model for $2 \le d \le 4$ was studied in [55–58] under the unitarity assumption, strictly speaking invalid for non-integer $d$, finding no inconsistency. On the other hand going below $d = 2$ without accounting for the loss of unitarity did lead to suspect results [59].

Analogously, one can argue that violations of unitarity should be negligible in the $O(N)$ model provided that $N$ is sufficiently large. A rule of thumb is to require that the dimensions of the exchanged $O(N)$ representations should be non-negative. E.g. the antisymmetric two-index tensor has dimension $N(N-1)/2$, so if it is exchanged, we should impose $N \ge 1$. This intuition is confirmed by the numerical conformal bootstrap analysis of the $O(N)$ CFT in $3 \le d \le 4$, $1 \le N \le 3$ [60]. There, assuming unitarity, the $O(N)$ CFT was isolated to a

<footnote n="12">¹²We thank Ning Su for discussions concerning this remark.</footnote>
bootstrap island whose position was in excellent agreement with the resummed $\epsilon$-expansion. Furthermore, this island disappeared when $N \to 1^+$. The conclusion is that non-unitarity was negligible for $N \ge 1$ but not for $N < 1$. Indeed, a low-lying state CFT state whose norm becomes zero at $N = 1$, and would have become negative for $N < 1$, was identified in [60].

As already mentioned, the $O(N) \times O(2)$ fixed point $C_+$ is expected to merge and annihilate with $C_-$ at a critical value $N_c(d)$. Our working assumption will be that violations of unitarity are negligible for $N > N_c(d)$, when the fixed points $C_\pm$ have real couplings. For $N < N_c(d)$ the fixed points go into the complex plane and continue their life there as “complex CFTs” [61, 62], having complex scaling dimensions and complex OPE coefficients, which is a much more violent violation of unitarity. Thus, while we expect to find solutions to the unitary bootstrap equations at $N > N_c(d)$, we may hope that these solutions will disappear at $N < N_c(d)$. This will be our method of determining $N_c(d)$.

## 2.3 Gap assumptions

The gap in the symmetry channel XY of spin $\ell$ will be denoted $gapXY_\ell$. This means that all operators in this channel, except for the isolated operators from Table 1, are assumed to satisfy $\Delta \ge gapXY_\ell$. The first group of gap assumptions is:

$$ \text{gapSS}_0 = d, \qquad (7a) $$

$$ \text{gapSS}_2 = d + 1, \qquad (7b) $$

$$ \text{gapAS}_1 = (d-1) + 1, \qquad (7c) $$

$$ \text{gapSA}_1 = (d-1) + 1. \qquad (7d) $$

These are easy to motivate. The $SS_0$ assumption is simply the fixed point stability condition of $C_+$ (see as well Remark 2.2). In the $SS_2$, $AS_1$ and $SA_1$ channels, we have respectively $T_{\mu\nu}$ of dimension $d$, and the conserved currents of dimension $d-1$. At large $N$, the gaps above these operators are equal to 2 at leading order. The increment 1 in (7b), (7c), (7d) was chosen in order to be comfortably below this large $N$ value.

Let us discuss next the $VV_0$ channel, which has two relevant operators $\phi$ and $\phi'$. The operator $\phi$ is treated as isolated. The leading order predictions for the dimension of $\phi'$ in the $\epsilon$-expansion and at large $N$ are:

$$ \Delta_{\phi'} = \begin{cases} 3 \times \frac{d-2}{2} + O(\epsilon), \\ \frac{d-2}{2} + 2 + O(1/N). \end{cases} \qquad (8) $$

We are not aware of results beyond the leading order. To be safe, we will use values for $gapVV_0$ comfortably below the $\epsilon$-expansion result from (8) which is the smaller prediction:

$$ \text{gapVV}_0 = d - 2. \qquad (9) $$

As already mentioned, we will also need a gap assumption in the $ST_0$ channel. The precise value of $gapST_0$ will be informed by the value of $\Delta_{ST'}$ in $\mathcal{H}$ and $C_+$, given in Eqs. (D.3), (D.8) to the first subleading order in $1/N$. In Fig. 4 we plot these values as a function of $d$, for the value $N = N_c^\epsilon(d)$, which is the $\epsilon$-expansion for $N_c(d)$, discussed in Appendix E.1. Given these curves, we will consider

$$ \text{gapST}_0 = \frac{d}{2} + \frac{3}{2}, \qquad (10) $$

a conservative gap assumption, since it lies comfortably below the large-$N$ predictions for $\Delta_{ST'}$ in both $C_+$ and $\mathcal{H}$. In Section 3.4, we will also consider the effect of varying $gapST_0$.

In all other channels, $gapXY_\ell$ will be set to the unitarity bound.
<figure><img src="image_8.png" /><figcaption>Figure 4: Comparison of the gap assumption gapST₀ = d/2 + 3/2 with the dimension of O<sub>ST'</sub> in H and C<sub>+</sub> at the first subleading order in 1/N, evaluated at N = N<sub>c</sub><sup>ε</sup>(d) from the ε-expansion. The C<sub>+</sub> curve has a kink due to crossing between the γ<sub>1</sub> and γ<sub>2</sub> coefficients from (D.7).</figcaption></figure>

**3 Conformal bootstrap analysis**

As we saw in the previous section, our setup leads to a three-dimensional parameter space P consisting of points x = (Δ<sub>ϕ</sub>, Δ<sub>SS</sub>, Δ<sub>ST</sub>), which are respectively the dimensions of the first primary scalars in the VV, SS and ST channels. With our gap assumptions not every point in this space is allowed, and our first order of business will be to get an idea of the allowed regions for various N and d.

We will search for allowed points by the navigator function method [49], which was already
used in several bootstrap studies [58, 60, 63-65]. We recall that the navigator function is con-
structed so that it is positive on disallowed points and negative on allowed points. Therefore,
to find allowed points we can use a local minimization algorithm for the navigator function.
Technical details of the navigator approach are discussed in Appendix F. Some parts of our
study also used the new skydive algorithm [50] for increased efficiency, which solves both
the navigator optimization problem and the semi-definite programming problem at the same-
time.

We will first use the navigator method to find islands in the three-dimensional space P
for fixed N, d. If by varying N or d an island disappears, or equivalently a negative local
minimum turns positive, then this is evidence that the corresponding CFT ceases to exist. We
will determine the exact point at which this happens by including N in our search, so we will
'navigate' towards the vanishing point of the island in the four-dimensional space spanned by
x and N. Our notation below will reflect these two different setups: although the navigator
function N always depends on all five variables N, d, and x = (Δϕ, ΔSS, ΔST) we will suppress
the directions that we hold fixed and write N(x) if we move in three dimensions and N(x, N)
if we move in four.

In all our searches we will hold fixed the gap assumptions of the previous subsection.

**3.1 Finding islands for large N**

We will begin by fixing N to be significantly larger than the expected critical value. In this
region the existence of chiral fixed points is not in doubt, and our first task will be to show
that the conformal bootstrap analysis is able to isolate them. We will carry out the analysis
first near d = 4, namely for d = 3.8 and N = 20. We will then carry out similar analysis
for d = 3 and N = 8. With d and N fixed we will navigate in the remaining three variables
x = (Δ<sub>ϕ</sub>, Δ<sub>SS</sub>, Δ<sub>ST</sub>) and look for islands of allowed points.
<figure><img src="image_9.png" /><figcaption>Figure 5: Top: Minimal bounding boxes containing three isolated allowed regions, at $\Lambda = 31$, in the space of $O(N) \times O(2)$-symmetric unitary$^{14}$ CFTs at $N = 20$, $d = 3.8$ under the gap assumptions (7), (9), (10). Bottom: The solid lines are the projections of these allowed regions to the planes of $(\Delta_{\phi}, \Delta_{SS})$ and $(\Delta_{\phi}, \Delta_{ST})$. The dashed lines indicate the same bounds if weakening (10) to $\text{gapST}_0 = \frac{d}{2}+1$. For reference, we also include the perturbative estimates of the $\epsilon$-expansion for the locations of $\mathcal{H}$ (blue), $C_+$ (red) and $C_-$ (green), see Appendix E.2. An allowed region containing the fixed point $C_-$ exists at this $\Lambda$ despite the actual physical fixed point violating the stability assumption (7a), indicating insufficient sensitivity to this constraint. We checked that all three islands still survive at the higher derivative order $\Lambda = 43$ (see footnote 15).</figcaption></figure>

### 3.1.1 $d = 3.8$ and $N = 20$

We begin our search by exploring $d = 3.8$ and $N = 20$. The former is close to $d = 4$, so we will be able to match our numerical results to the $\epsilon$-expansion, and the latter is moderately above the $\epsilon$-expansion prediction $N_c^\epsilon(3.8) = 17.3997(6)$ (see Appendix E.1) so we expect $C_+$ to exist. The result of our investigations at the derivative order $\Lambda = 31$ (see Appendix F for the numerical details) is shown in Fig. 5. Within the shown ranges we find that almost every point is excluded with the exception of three small allowed regions in close vicinity of the perturbative predictions for the scaling dimensions of $\mathcal{H}$, $C_+$ and $C_-$.<fnref n="13" />

Before we discuss the physical consequences of this result, let us explain the precise meaning of the three boxes shown in Fig. 5. For each region we first found a point where the navigator function was negative, indicating that the point was allowed. By continuity this point should be part of a larger allowed region, called an “island”, which one could in principle delineate by carefully tracing its boundary where the navigator function is exactly zero. To save computational time, we have however chosen not to do so, and instead simply determined the extremal dimensions of the islands along the three axes of the search space $\mathcal{P}$. The islands are thus contained within the boxes shown in Fig. 5, or as we say “boxed in.” See Appendix F for further details.

<footnote n="13">¹³A further allowed region exists for larger $\Delta_{\phi}$, outside the range shown in the plot. This region, known as the “peninsula”, is analogous to the one first identified in the three-correlator study for the 3D Ising model [34].</footnote>

<footnote n="14">¹⁴See the caveats discussed in Section 2.2.</footnote>
<figure><img src="image_10.png" /><figcaption>Figure 6: Top: Minimal bounding boxes containing isolated allowed regions in the space of O(N) × O(2) symmetric unitary CFTs for N = 8, d = 3 under the assumptions (7), (9), (10) at Λ = 31. Bottom: Projections of these allowed regions on the (Δ<sub>φ</sub>, Δ<sub>SS</sub>) and (Δ<sub>φ</sub>, Δ<sub>ST</sub>) planes (solid). The dashed lines indicate the same bounds under the weaker gap assumption gapST<sub>0</sub> = d/2 + 1. For reference, we also include the resummed field theory estimates of [6, 66], along with their uncertainty, for C<sub>+</sub>; and the unresummed large-N expansion for H and C<sub>-</sub> (see the end of Appendix E.2). In this plot we did not find a separate isolated region associated to C<sub>-</sub>.</figcaption></figure>

The most encouraging aspect of Fig. 5 is the clean separation between the different fixed points. To see why, suppose that we had instead found one big allowed region encompassing all three theories. For lower N we would then know that this big region can never disappear, simply because $\mathcal{H}$ always exists, which would make it impossible to determine $N_c$. The possible way out of that impasse could have been to exclude $\mathcal{H}$ from the analysis by making stronger gap assumptions, such as increasing gapST<sub>0</sub>. But this would have been very delicate since the estimated values for this gap in $C_+$ and $\mathcal{H}$ are close (see Fig. 4). Fortunately the separated islands that we found, shown in Fig. 5, instead allow us to proceed differently: we can simply determine $N_c(d)$ by the vanishing of the $C_+$ and $C_-$ islands at lower N, all without worrying about the continued existence of the $\mathcal{H}$ island. This is how we will proceed below.

Finally we should comment on the unexpected appearance of the $C_-$ island. The $C_-$ fixed point is RG-unstable, i.e. it has a second relevant singlet scalar, and in principle it should be ruled out by the gap assumption (7a). That we still see an allowed region around $C_-$ means that our bootstrap setup is not sufficiently sensitive to this particular gap assumption (although it is sensitive to other assumptions, because the allowed region is actually quite small). This may be due to the additional singlet scalar being only slightly relevant, since $N = 20$ is not too far above $N_c^\epsilon(3.8)$. We checked that the $C_-$ island, as the other two, survives also at $\Lambda = 43$.<fnref n="15" />

<footnote n="15">¹⁵Namely, we checked that the $\Lambda = 43$ navigator is negative at the minima of the $\Lambda = 31$ navigator inside each of the three islands.</footnote>
### 3.1.2 $d=3$ and $N=8$

In Fig. 6, we show the results of the same exercise carried out for $d=3$ and $N=8$, which is close to but above the $\epsilon$-expansion value $N_c^\epsilon(3) = 5.9(2)$ obtained as in Appendix E.1.<fnref n="16" /> In this case we see two disconnected islands, one around $\mathcal{H}$, the other around $C_+$, The islands are significantly larger than in Fig. 5. Predictions of the (resummed) $\epsilon$-expansion for the scaling dimensions at the $\mathcal{H}$ and $C_+$ fixed points comfortably fall into the corresponding islands. In this figure, unlike for $d=3.8$, there is no separate disconnected island around the $C_-$ fixed point; instead, the predicted $C_-$ scaling dimensions fall inside the same island as $C_+$ (or close to it). Once again, this shows low sensitivity to the gap assumption (7a). While there is no separate $C_-$ island, the navigator function has two separate minima inside the island containing both fixed points, close to their predicted positions, see Fig. 8 below.

## 3.2 $N_c(d)$ from the disappearance of the $C_+$ island

In Section 3.1 we chose $N$ to lie well above its critical value predicted by the $\epsilon$-expansion. We showed bootstrap evidence that for those values of $N$ the $C_\pm$ fixed points do exist, and have operator dimensions consistent with the $\epsilon$-expansion.

We will now lower $N$ and determine the critical $N_c$ at which the island containing $C_+$ and $C_-$ disappears. We will denote the corresponding conformal bootstrap determination as $N_c^{CB}(d)$. This $N_c^{CB}(d)$ should be seen as a rigorous lower bound on the true $N_c(d)$. It is a lower bound because the island can disappear only faster with improved numerical precision, and it is rigorous because all our gap assumptions were rather conservative.

As mentioned above, for this analysis we consider $N$ as a varying parameter of the navigator function, on par with the CFT data variables $x = (\Delta_\phi, \Delta_{SS}, \Delta_{ST})$. Thus, we will be considering four-dimensional allowed regions where the navigator function

$$ \mathcal{N}(\Delta_{\phi}, \Delta_{SS}, \Delta_{ST}, N) \qquad (11) $$

is negative. The $C_+$ island now lives in the four-dimensional parameter space, and finding $N_c^{CB}(d)$ is equivalent to finding the extremal point of the island in the $N$ direction, i.e. the point where the navigator minimum turns from negative to positive. See Appendix F.5 for technical details. In this way we determined $N_c^{CB}(d)$ in the range $3 \le d \le 3.8$ in steps of 0.1.<fnref n="17" />

The result of this analysis was shown in Fig. 3 on page 6. This plot, which directly addresses question (4), is the main result of our paper. In Fig. 3 we also show for comparison $N_c^\epsilon(d)$ from the $\epsilon$-expansion (Appendix E.1) and $N_c^{NPRG}(d)$ from the NPRG results of [3, Fig. 9]. For $d$ close to 4, the curve $N_c^\epsilon(d)$ should be trustworthy. In this region, our curve $N_c^{CB}(d)$ shows a rather good agreement with $N_c^\epsilon(d)$ to linear order in $\epsilon$, while $N_c^{NPRG}(d)$ shows a small but noticeable negative first-order deviation from $N_c^\epsilon(d)$. In the range $3.5 \le d \le 3.8$ the curve $N_c^{CB}(d)$, which as mentioned is a rigorous lower bound on $N_c$, lies above the NPRG curve and rules it out.

As $d$ approaches 3, the difference between the NPRG and $\epsilon$-expansion predictions decreases, and we have:

$$ N_c^\epsilon(3) = 5.96(19), \qquad N_c^{NPRG}(3) = 5.1. \qquad (12) $$

On the other hand our curve lies quite a bit lower than the other two in this range, and it ends in $d=3$ at

$$ N_c^{CB}(3) = 3.78. \qquad (13) $$

<footnote n="16">¹⁶This number is close to $N_c(3) = 5.96(19)$ obtained in [6] from the $\epsilon$-expansion, which was obtained using a different method from the one described in Appendix E.1 and instead was based on a combination of several resummation techniques.</footnote>

<footnote n="17">¹⁷We note that there is a subtlety in our reasoning which has to do with the fate of the islands as we lower $N$ continuously from large values, where they are isolated, to $N_c$, where they disappear. This will be discussed below in Section 3.4.3.</footnote>
Although this value is significantly lower than $\epsilon$ and NPRG values, we stress again that it is a
rigorous lower bound. Importantly, even this lower bound (13) is strong enough to rule out
the physically relevant values $N = 2, 3$.

A final feature of our curve is that it shows no sign of the turnaround behavior predicted
to happen by the fixed-dimensional RG in Fig. 2 for $d \approx 3.2$.

We emphasize again that we used the gap assumptions laid out in Section 2.3, in particular
(7), (9), as well as (10) for gapST₀. As mentioned, we can afford putting this gap conserva-
tively low, since our method does not require eliminating the H island. It is interesting to
inquire how our N_c(d) depends on the imposed gap assumptions, in particular on gapST₀.
This will be done in Section 3.4 below for d = 3, where we will see that the dependence on
gapST₀ is rather weak in a certain range around (10). On the other hand, we could try to
improve N_c^CB(d) by increasing Λ or by considering more correlators. These improvements will
be left for the future.

## 3.3 Tracking the navigator minima

In Section 3.1 we have found isolated bootstrap islands for $N$ much above the critical value, and for two values of $d$. Being in the region of large $N$ and, for $d = 3.8$, also close to $d = 4$, the position of those islands could be with confidence identified with the Heisenberg and the chiral fixed points of the $O(N) \times O(2)$ model. Then, in Section 3.2 we found that bootstrap islands are disappearing along the curve $N = N_c^{CB}(d)$ shown in Fig. 3. The gap assumptions which went into the determination of the islands are expected to be satisfied by both the Heisenberg and the chiral fixed points. The Heisenberg fixed point exists for any $N \ge 1$ and the island containing it is not expected to disappear as $N$ decreases. Naturally, we associated the disappearing islands with the chiral fixed point.

In this section we will provide further evidence confirming that the identification of the
disappearing islands with the chiral fixed point is correct. For this we will track, for two values
of $d$, the minimum of the navigator function with respect to ($\Delta_{\phi}$, $\Delta_{SS}$, $\Delta_{ST}$), within the $C_+$
island, as a function of $N$. Varying $N$ from large values down to $N_c$, we will show that the
minimum varies continuously. We track minima and not the whole islands, since it would be
much more costly numerically to determine the island extensions in every $N$.

3.3.1 $d = 3.8$

The position of $C_+$ navigator minima for $d = 3.8$ and for $N$ ranging from $N = 30$ down to $N_c^{CB}(3.8)$ is shown in Fig. 7.<fnref n="18" /> The leftmost points correspond to $N = N_c^{CB}(3.8) = 17.1585$ which is where the $C_+$ island disappears for our chosen value of $\Lambda = 31$. In the same plots we show by red curves the $\epsilon$-expansion prediction for the same quantities, computed as in Appendix E.2. These curves end at $N_c^{\epsilon} = 17.3997(6)$ for $d = 3.8$. We see that the navigator minimum varies continuously with $N$, tracking closely the $\epsilon$-expansion curve, except for a few leftmost points, close to $N_c$. We interpret the difference between $N_c^{CB}$ and $N_c^{\epsilon}$, as well as the residual deviations between the bootstrap points and the $\epsilon$-expansion curves as numerical artifacts, expected to disappear as $\Lambda \to \infty$. We have indeed observed that these deviations decrease as $\Lambda$ is increased through $\Lambda = 19, 23, 27$ (not plotted) to $\Lambda = 31$.

These plots also allow us to demonstrate another expected feature of the mechanism be-
hind the disappearance of the chiral fixed point at $N \to N_c(d)^+$. As mentioned, this should
happen via “merger and annihilation”: two fixed points $C_{\pm}$ which exist at $N > N_c(d)$ merge

<footnote n="18">¹⁸The gap assumptions used were the same as in Fig. 5, except for the use of the aggressive gapST₀ = 3.85, roughly halfway between the C⁺ and H large-N curves in Fig. 4 at d = 3.8. As N grows, the C⁺ dimension of O<sub>ST</sub> tends to 4, so this gap should be valid for N ≥ N<sub>c</sub>(3.8).</footnote>
<figure><img src="image_11.png" /><figcaption>Figure 7: Diamonds: $\Delta_{\phi}$, $\Delta_{SS}$, $\Delta_{ST}$ at the minimal navigator point within the $C_+$ island, as a function of $N$, for $d = 3.8$ and the derivative order $\Lambda = 31$. The leftmost points corresponds to $N = N_c^{CB}(3.8) = 17.1585$ which is the critical $N$ value at which the $C_+$ island disappears. Red curves: the $\epsilon$-expansion prediction for the same quantities.</figcaption></figure>

into a single CFT at $N = N_c(d)$. At even lower $N$, the fixed points $C_{\pm}$ would have complex couplings and complex anomalous dimensions. These are “complex CFTs” [67], and they become invisible in our studies since we assume real scaling dimensions.

The merger and annihilation scenario was discussed in various contexts in many works, notably [65, 67–70]. Its telltale signature is a near-marginal singlet scalar operator $\mathcal{O}$ whose scaling dimension shows a square-root behavior near the merger point (see e.g. [67], Eq. (2.4))

$$ \Delta_{\mathcal{O}} \approx d \pm \text{const.} \sqrt{N - N_c(d)} \quad (N - N_c(d) \ll 1), \quad (14) $$

where the sign $\pm$ corresponds to $C_{\pm}$. This behavior can be linked to a square-root singularity in the coupling of $\mathcal{O}$ near the merger.

For us $\mathcal{O} = \text{SS}'$, the second-lowest SS scalar. In our numerical study, we found it hard to extract its dimension reliably using the Extremal Functional Method [31], likely because it is too close to the gap assumption (7a). We leave direct numerical bootstrap verification of (14) to future work.

However, since related to a square-root singularity of a coupling, square-root behavior should be present not just in $\Delta_{\mathcal{O}}$ but in any scaling dimension or any OPE coefficient (see [68], or [62], Fig. 7). And indeed, in Fig. 7 we clearly see such a square-root behavior in $\Delta_{SS}$ and $\Delta_{ST}$ data. According to the $\epsilon$-expansion, the expected square-root singularity in these quantities is quite pronounced, so it is easy for us to detect it numerically. For the third quantity $\Delta_{\phi}$, the $\epsilon$-expansion predicts a tiny square-root singularity not visible on the scale of Fig. 7.

**Remark 3.1.** We note a recent conformal bootstrap attempt to observe the merger and annihilation scenario in a different model [65]. That paper studied the critical and tricritical point of the 3-state Potts model as a function of dimension $d \ge 2$. Merger and annihilation is supposed to occur at some $d = d_c$ between 2 and 3. As a sign of this, a square-root behavior of multiple scaling dimensions seems be setting in, as expected, in their Fig. 2, at least sufficiently below the merger point. However in their study this behavior is replaced closer to the merger point
<figure><img src="image_12.png" /><figcaption>Figure 8: Diamonds and dots: positions of the $C_+$ and $C_-$ navigator minima for $d = 3$, as a function of $N$, at the derivative order $\Lambda = 31$. The leftmost points correspond to $N = N_c^{CB}(3) = 3.78$ where the $C$ island disappears. Red, green, blue curves: the $\epsilon$- and $1/N$-expansion predictions. Red error bars: estimate of $\Delta_{ST}$ in $C_+$ from fixed-dimensional RG (see [66, Table III] and Appendix E.2).</figcaption></figure>

by a linear approach, which is not expected theoretically and is likely a numerical artifact. Our $d = 3.8$ results in Fig. 7 show a clearer picture. Note that in our study we managed to isolate the $C_+$ fixed point to an island, which was not the case for the theories studied in [65].

### 3.3.2 $d=3$

We now move to $d = 3$ and show in Fig. 8 the position of navigator minima for $N$ ranging from $N = 8$, the value in Fig. 6, down to $N_c^{CB}(3) = 3.78$ where the $C$ island disappears.

We observe several different features compared to Fig. 7. First, at $N = 8$ the navigator function now has two well separated minima, close to the respective expected positions of $C_\pm$ within the island. We track both minima which we therefore tentatively call the $C_+$ navigator minimum (diamonds) and the $C_-$ minimum (dots).

As we decrease $N$ we find that the $C_+$ minimum is the first to disappear, around $N = 7$. At this point the local minimum collides with a saddle, leading to the disappearance of both. Precisely at this transition, one Hessian eigenvalue is zero, which we checked by explicitly computing the Hessian. However, the navigator function itself is still negative so we are within the island.<fnref n="19" />

The $C_-$ minimum on the other hand, which we recall now lies within the same island, continues to exist for lower values of $N$. As shown in Fig. 8, it varies continuously with $N$ until it becomes a positive minimum at $N = N_c^{CB}(3) = 3.78$, at which point the $C$ island disappears.

Various colored curves in Fig. 8 show the $\epsilon$-expansion or $1/N$-expansion predictions for the positions of various theories. The blue lines show the position of the $O(2N)$ fixed point. We

<footnote n="19">¹⁹While it is very interesting to note that all local minima of the navigator function discovered up to now seem to have their origin in physical theories it is important to keep in mind that the absence of the $C_+$ minimum below $N < 7$ does not imply this theory does not exist. What is more important is that the conformal data values that such a stable CFT would have are not excluded by our current analysis for $N > 3.78$, but are rigorously excluded for $N < 3.78$.</footnote>
see that the $C_-$ minimum remains well separated from this line for all $N$ of interest. The red and green lines show the positions of the $C_+$ fixed points, which according to the $\epsilon$-expansion annihilate at $N = 5.96(16)$. We see that the $C_-$ minimum tracks closely the green $\epsilon$-expansion curve for $N \gtrsim 6.5$.

Fig. 8 in $d = 3$ paints a somewhat less satisfactory picture than Fig. 7 in $d = 3.8$. Ideally we would have found two persistent local minima, one positive (so excluded, corresponding to $C_-$) and the other negative (corresponding to $C_+$), which upon decreasing $N$ would merge precisely at $N_c^{CB}(3)$ and then remain positive for smaller $N$. Such a scenario would then naturally also produce the desired square-root behavior of the scaling dimensions.

Instead we see in Fig. 8 that the $C_+$ minimum disappears too soon, followed by a rather large region between $N_c^{CB}(3) = 3.78$ and $N \approx 6.5$ without good agreement with the $\epsilon$-expansion. All of this could however be numerical artifacts, and in future work it would be interesting to see how the Fig. 8 evolves when increasing $\Lambda$ and whether the features at small $N$ remain or the bounds converge closer to the ideal scenario.<fnref n="20" />

On the positive side, Fig. 8 does demonstrate a continuous connection between the disappearing island at $N = N_c^{CB}(3)$ and the island at $N = 8$ which could be unambiguously associated with the $C_+$ fixed points. It is therefore not in contradiction with our claim that $N_c^{CB}(3)$ is a rigorous lower bound for the disappearance of $C_+$.

### 3.4 How does $N_c^{CB}(3)$ depend on the choice of gapST$_0$?

In this subsection we will investigate the dependence of $N_c^{CB}(d)$ on gapST$_0$. The best possible scenario would be that $N_c^{CB}(d)$ does not change at all if we vary gapST$_0$ within a certain interval. Of course we cannot increase gapST$_0$ above the actual dimension $\Delta_{ST'}$ of the second scalar in the ST representation since this would exclude the $C_+$ fixed point altogether. We probably cannot decrease gapST$_0$ too much either, because then the allowed regions can become so big that the $C_+$ island merges with either the $\mathcal{H}$ island or the larger peninsula region discussed in footnote 13.

The actual variation of $N_c^{CB}(3)$ with gapST$_0$ is shown in Fig. 9. Most importantly it shows a rather mild dependence as long as we choose gapST$_0$ between about 2.8 and 3.4, confirming the robustness of our estimate Eq. (13) in the sense that it does not depend significantly on this gap assumption. Note that the upper value of about 3.4 is roughly compatible with the first-order large $N$ estimate for $\Delta_{ST'}$ which is shown as the red dashed curve in the figure. We unfortunately see no evidence for a flattening of the curve with increasing $\Lambda$, but at the same time it is clear that our result is far from being converged.

Another way to check for the sensitivity with respect to gapST$_0$ is to extract the operator spectrum from the extremal functional at $N_c(3)$, which we did for the blue points in Fig. 9. In an ideal scenario there is simply no operator with a scaling dimension equal to gapST$_0$, as this would mean that we could vary gapST$_0$ a bit without changing the functional and the bound. It is however often the case that numerical artefacts spoil this scenario, and this appears to be the case here. Indeed, as we show in Fig. 10, the scalar spectrum in the ST sector consistently shows an operator at the gap for each value of gapST$_0$. We also verified that its OPE coefficient is non-vanishing: $\lambda_{\phi\phi O_{ST'}}$ varies smoothly from 0.06 to 0.03 in the plotted window. (For comparison, we obtained $\lambda_{\phi\phi O_{ST}} \approx 0.3$ and $\lambda_{\phi\phi O_{ST''}} \approx 0.005$.) It is the square of these coefficients that multiplies the conformal blocks, whose normalization we took to be that of [44].) In future work it will be important to verify that both Fig. 9 flattens out and that $\lambda_{\phi\phi O_{ST'}}$ decreases upon increasing $\Lambda$.

<footnote n="20">²⁰At this point it is for example unclear whether there is any physical significance to the change of slopes in the conformal data that is observed around $N = 4.2$. This would be interesting to investigate further.</footnote>
<figure><img src="image_13.png" /><figcaption>Figure 9: Lower bound on $N_c$ as found by the numerical conformal bootstrap as a function of the assumed value of gapST$_0$ at $\Lambda = 31$ (in blue). Some additional points are shown for $\Lambda = 27, 29, 33, 35$ in respectively orange, green, red, purple. The dotted vertical line corresponds to the previously chosen value gapST$_0$ = 3, yielding $N_c^{CB}(3) = 3.78$ at $\Lambda = 31$ as stated previously in Eq. (13). The dashed curved line indicates the large $N$ estimate of $\Delta_{ST'}$ in $C_+$ (see Eq. (D.8)). It is natural to expect the $C_+$ island to disappear somewhat to the right this curve.</figcaption></figure>

In the remainder of this subsection we discuss interesting subtleties that arise from studying the endpoints of the $\Lambda = 31$ curve in Fig. 9. These endpoints came about as follows. Fig. 9 was obtained by simply repeating the analysis of Section 3.2 for different values of the gapST$_0$, using again the skydive algorithm of [50], see Appendix F. In practice this procedure turned out to be very delicate for small and large values of gapST$_0$. We frequently observed runaway behavior, forcing us to restart with a slightly different initial point. The endpoints of the $\Lambda = 31$ curve are the last points for which we were able to obtain convergence.

### 3.4.1 Hessian and derivative information

Let us first discuss the behavior of the navigator function $N(\Delta_\phi, \Delta_{SS}, \Delta_{ST}, N)$ in the vicinity of $N_c$ for a generic value of gapST$_0$. Recall that $N_c$ is the lowest point in the $N$ direction where

<figure><img src="image_14.png" /><figcaption>Figure 10: The first three scalar operators in the ST sector in the extremal spectrum as a function of gapST$_0$. (More precisely, this is the extremal spectrum corresponding to the blue points with $\Lambda = 31$ in Fig. 9.) The first operator ($\Delta_{ST}$, blue dots) is isolated, the second one ($\Delta_{ST'}$, orange dots) is precisely at the gap imposed by gapST$_0$ (as indicated by the blue area), and the third operator ($\Delta_{ST''}$, green dots) lies above it.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 11: Eigenvalues of the Hessian H of the navigator function at N_c in the three Δ directions, as a function of gapST_0 for d = 3.</figcaption></figure>

$N \le 0$. Since the navigator function is smooth, this must mean that:

$$ \begin{gathered} N(\Delta_{\phi}, \Delta_{SS}, \Delta_{ST}, N_c) = 0, \\ \frac{\partial}{\partial N} N(\Delta_{\phi}, \Delta_{SS}, \Delta_{ST}, N_c) \le 0, \end{gathered} \qquad (15) $$

where saturation of the inequality on the last line would be atypical and would imply a constraint on the higher derivatives such that the $N$ is positive for $N$ a little bit below $N_c$. Likewise the navigator function must be at a local minimum with respect to the other parameters, so if we define its gradient and Hessian in Δ-space at the critical point $N_c$ as

$$ g_i := \frac{\partial}{\partial \Delta_i} N(\Delta_{\phi}, \Delta_{SS}, \Delta_{ST}, N_c), \quad \Delta_i \in (\Delta_{\phi}, \Delta_{SS}, \Delta_{ST}), \qquad (16) $$

$$ H_{ij} := \frac{\partial^2}{\partial \Delta_i \partial \Delta_j} N(\Delta_{\phi}, \Delta_{SS}, \Delta_{ST}, N_c), $$

then we know that

$$ g_i = 0, \quad \text{and} \quad H > 0. \qquad (17) $$

Note that the second derivative in the N direction is not constrained.

Experimentally the observed behavior agrees with the previous discussion for all the points where our algorithm converges. The algorithm itself guarantees that $N = 0$ and $g_i = 0$ at optimality. Properties of the other derivatives of the navigator are shown in Figs. 11 and 12. In the first we show the eigenvalues of H, which are indeed all positive. In the second figure we include the N direction. We first of all see that the gradient is negative, as expected. We also see that including the fourth direction yields a second-order derivative that is generally positive but turns negative at both endpoints.

### 3.4.2 Left endpoint behavior

At the leftmost point we observe from Fig. 12 that the gradient in the N direction rapidly decreases in magnitude and that the four-dimensional Hessian has a negative eigenvalue. Since we have both gradient and Hessian information, we can construct a quadratic approximation to the real navigator function. This approximation predicts that $N$ becomes negative again for $N$ somewhere below $N_c$. For example, for the leftmost point in Fig. 9 at gapST_0 = 2.82,
<figure><img src="image_16.png" /></figure>

Figure 12: Behavior of the navigator function at $N_c$ in the $N$ direction, as a function of gapST$_0$ for $d = 3$. On the left we plot the gradient. On the right we plot the fourth eigenvalue of the Hessian of the four-dimensional navigator $N(\Delta_\phi, \Delta_{SS}, \Delta_{ST}, N)$. We note that its corresponding eigenvector mostly in the $N$ direction, and the other three eigenvalues agree to within 10% with those shown in Fig. 11. This eigenvalue therefore accurately captures the second-order behavior of the navigator when one includes the $N$ direction. (We could have opted to simply plot the second partial derivative in the $N$ direction, but this always comes out positive and therefore does not convey the right information.)

where we obtained $N_c \approx 3.62$, this quadratic approximation would predict that the navigator becomes negative again at $N \approx 3.13$.

The analysis of gradient and Hessian therefore provides evidence for another allowed region that exists for $N < N_c$ and low values of gapST$_0$. This is not in disagreement with the expectations we formulated above, where we stated that gapST$_0$ must be chosen sufficiently large for the $C_+$ island to (a) exist as an isolated island, and then (b) to disappear. It is however surprising that this additional region is disconnected, since it would have been more natural to see a smooth merging of the $C_+$ island with another allowed region. Such a scenario would correspond to the three-dimensional $H$ becoming singular, but Fig. 11 shows no indication of such behavior near the left endpoint at low gapST$_0$.

In future work it will be important to investigate in detail the existence of a disconnected region below $N_c$. Is this region connected to the $H$ island or the peninsula? And is it due to the unnecessarily low value of gapST$_0$ or does it also persist for higher values? And is it related to the existence of fixed points below $N = 2$ predicted by the $\epsilon$ expansion as discussed in footnote 22?

### 3.4.3 Right endpoint behavior

The behavior at the right endpoint is more surprising. Here the main novel characteristic is the steep drop in the lowest eigenvalue of the three-dimensional $H$ as shown in Fig. 11. Naive extrapolation then predicts that the Hessian $H$ would become singular for a slightly higher value of gapST$_0$. We would then necessarily have a direction in $\Delta$ space in which the navigator decreases, simply because the leading term in that direction is the cubic term. The $C_+$ island is then no longer isolated and instead gets connected to another region which does not disappear for $N$ below $N_c$. We do not know what this additional allowed region is, but the natural candidate would be the aforementioned peninsula that we know exists for all $N$.

Although we discovered this behavior at large gapST$_0$, it must actually also persist at lower gapST$_0$ because lowering a gap can only make the allowed region larger (at fixed $N$).

Let us therefore return to our earlier chosen value of gapST$_0$ = 3. If the above scenario is correct then the behavior is as follows. For $N = 8$ Fig. 6 shows an isolated $C_+$ island. For $N$ slightly above $N_c$ there is likewise an isolated island (by continuity of the navigator function
and the positivity of $H$). But for intermediate values of $N$ there may not be a well-defined
island. Instead, we would predict a merger with the peninsula and then again a re-detachment
of an island.<fnref n="21" />

We emphasize that we observed this issue only for larger values of gapST₀ simply because
that is the only place where we investigated these intermediate values of N. Investigating the
structure of the allowed regions in detail should be a priority for future numerical bootstrap
studies.

4 Review of prior bootstrap studies

Let us review the previous bootstrap studies on the O(N)×O(2) model [24, 25]. Comparison
to our work will be done in the conclusions sections.

Ref. [24] was the first bootstrap work on the O(N)×O(2) model. It was a single-correlator
study, using only the 4pt function $\langle\phi\phi\phi\phi\rangle$. Working in $d=3$ and specializing to $N=3$, they
found a mild kink in the gap maximization plot in the ST$_{ℓ=0}$ channel, as a function of $\Delta_\phi$. They
observed that at the kink, $\Delta_\phi$ as well as the dimensions of scalars in the SS, ST, TS, TT, and
AA representations extracted using the extremal functional method agree reasonably well with
the values of these dimensions predicted for the putative chiral O(3)×O(2) fixed point of focus
type found using the fixed-dimension expansion (for one of two RG schemes). Of course the
focus fixed point, if it existed, would not be unitary since it would have a scaling operator with
complex dimension. Ref. [24] assumed unitarity, and in particular real scaling dimensions for
all operators. Be that as it may, Ref. [24] interpreted the existence of the kink and the above-
mentioned agreement as evidence for the actual existence of a chiral O(3)×O(2) unitary fixed
point.<fnref n="22" />

The more recent study [25] developed the initial findings of [24] in two directions:

1. With the same single-correlator setup as in [24], they carried out gap maximization in the ST<sub>ℓ=0</sub> and TS<sub>ℓ=0</sub> channels<sup>23</sup> for N = 3, 4, 5, 10, 20. They found that for large N (N = 10, 20) when the corresponding unitary CFTs are reliably known to exist, there are kinks in these plots which match closely with the chiral (TS<sub>ℓ=0</sub>) and antichiral (ST<sub>ℓ=0</sub>) fixed points.<sup>24</sup> For smaller N the kinks get milder, and their association with actual unitary CFTs is unclear.

Note the surprising inversion: at large *N* the chiral fixed points are associated with *TS*<sub>*ℓ*=0</sub> gap maximization kinks [25], while for *N* = 3 they are associated with *ST*<sub>*ℓ*=0</sub> gap maximization kinks [24].

2. The second part of the study in [25] proceeds in the three-correlator setup $\langle\phi\phi\phi\phi\rangle$, $\langle\phi\phi ss\rangle$, $\langle ssss\rangle$ (i.e. the same as in our paper), imposing an assumption that the $TS_{\ell=0}$ or $ST_{\ell=0}$ gap saturates the corresponding gap maximization bound from the single-correlator analysis,

<footnote n="21"><sup>21</sup>An island that first grows and then shrinks as one decreases *N* can also be found in the numerical investigations of the O(*N*) models, see for example Fig. 1 of [36].</footnote>

<footnote n="22"><sup>22</sup>It should be mentioned that Ref. [24] also discussed the “collinear” (also called “sinusoidal” [1]) fixed point, which has the quartic coupling $\nu < 0$. According to the state-of-the-art resummed $\epsilon$-expansion results, the collinear fixed point exists only when $N < 1.970(3)$ [6], by far excluding the physically interesting value $N = 4$ when such a fixed point, if it existed, could describe the chiral phase transition in QCD with two massless quark flavors. In this case, there are again fixed-dimension schemes which do find a collinear fixed point for $N = 3, 4$ [66, 71]. In the interest of space, and since our focus in this paper is on the chiral fixed point, we will not discuss this additional controversy here (see also [72]).</footnote>

<footnote n="23"><sup>23</sup>Note that the symmetry group is O(2)×O(N) in [25] while it is O(N)×O(2) for us and [24]. Therefore the order of representations is reversed TS↔ST. We are using our conventions here. In [25], ST is called W and TS is called X.</footnote>

<footnote n="24"><sup>24</sup>This attribution goes back to [73], where O(N) × O(3) CFTs were studied in d = 3, for 5 ≤ N ≤ 20. There, the kinks in the TS<sub>ℓ=0</sub> and ST<sub>ℓ=0</sub> gap maximization plots at large N were first identified with the chiral and antichiral fixed points. Refs. [24, 73] were the first bootstrap studies of multiscalar models with more complicated symmetry than O(N).</footnote>
plus several other reasonable gap assumptions on the spectrum. With these assumptions, they
are able to turn kinks into islands of rather complicated shapes. For large $N$, these islands can
be mapped to the chiral and antichiral CFTs.

Then they consider small $N = 3$, and for $\text{ST}_{\ell=0}$ saturated gap (see the inversion noted above) the island they find agrees with the location of the putative $O(3) \times O(2)$ chiral focus-type fixed point (in the same $\overline{\text{MS}}$ scheme that saw better agreement with the single-correlator results of [24]). Concerning the issue that the putative fixed point is focus-type, they say that “It is unclear to us how sizable nonunitarities could have been missed by our bootstrap results.” As they point out, some previously found non-unitary CFT islands disappeared with increasing the constraining power of the numerics [74], and this may also happen to their island. They conclude by saying “Our results provided support for the existence of these fixed points, but we saw no signs of nonunitarity. Overall, we were unable to provide conclusive answers, but we believe that more dedicated bootstrap work with stronger numerics will be able to reach definitive conclusions in the near future.”

# 5 Comparison to prior work, conclusions and outlook

Let us start by emphasizing the aspects in which our study which were different from prior
work [24,25], reviewed in Section 4. The setup of [24,25], working in $d = 3$ and for a discrete sequence of $N$, could not clearly see if and how the putative $O(3) \times O(2)$ CFT connected to well-established $O(N) \times O(2)$ CFTs at large $N$ and for $d$ close to 4. In contrast, our study investigated continuously the 2-dimensional parameter space ($d, N$). We could therefore see the unitary $O(N) \times O(2)$ CFTs disappear as $N$ approaches some critical curve $N_c(d)$ from above. The second difference between our study and [25] is in how the $\text{ST}_{\ell=0}$ channel of the $\phi \times \phi$ OPE was treated in the multiple correlator setup. They set the first operator in this channel to saturate the single correlator bound, while making no assumptions about subsequent operators. We, on the other hand, treated $\Delta_{\text{ST}}$ as a free parameter which was one of the arguments of the navigator function, while restricting the next operator in this channel to be above a gap $\text{gapST}_0$. We saw in Section 3.4 that with this assumption, $N_c^{\text{CB}}(d = 3)$ shows only mild dependence on $\text{gapST}_0$, while remaining safely above 3. It is possible that the setup of [25] was simply not constraining enough to see a large-enough $N_c^{\text{CB}}(d = 3)$. Finally, our study used $\Lambda = 31$, while [25] used $\Lambda \le 25$.<fnref n="25" />

Our curve $N_c^{\text{CB}}(d)$ of Fig. 3 is, in the range $3 \le d \le 3.84$, a monotonic curve consistent in shape and position with the curves extracted from the NPRG and the $\epsilon$-expansion. We have thus ruled out the S-shaped FD curve in Fig. 2, which has a turnaround point at $d \approx 3.2$ [14].<fnref n="26" />

Let us now come back to the question (3) about the existence of unitary stable CFTs with $O(N) \times O(2)$ symmetry for $N = 2, 3$. Being fully agnostic, there still remains small loopholes which might allow their existence. First, on the basis of our results alone, we cannot rule out a turnaround of the $N_c(d)$ curve at some $d_c < 3$, i.e. somewhat below the range explored here. Indeed, in the MZM fixed-dimensional scheme (see footnote 26), such a turnaround does happen for $d$ slightly below 3 (see the discussion in [14] below Fig. 3), although its precise location has not been determined. It seems implausible from the look of our curve in Fig. 3 that such a turnaround would happen just below $d = 3$. In the future it would be interesting to extend our study to smaller $d$ in order to definitively exclude such a turnaround.

<footnote n="25">²⁵The PyCFTboot parameters $m_{max} = 7$, $n_{max} = 9$ in [25] correspond to a subset of $\Lambda = 2n_{max}+m_{max} = 25$ derivatives.</footnote>

<footnote n="26">²⁶The turnaround point $d \approx 3.2$ corresponds to the $\overline{\text{MS}}$ fixed-dimension scheme, generally considered in the literature as more reliable than another scheme called MZM (massive zero-momentum). The fixed points found in the MZM scheme for $N = 2, 3$ lie outside of the region where the beta functions are Borel summable [11], whereas they are found around the boundary of this region in the $\overline{\text{MS}}$ scheme [14].</footnote>
<figure><img src="image_17.png" /><figcaption>Figure 13: Three scenarios for the shape of the curve separating the region of the (d, N) plane where a unitary O(N) × O(2) CFT exists (white) from the one where it does not (gray). The scenario on the left, where the boundary curve has a turnaround point (•) at a d• > 3, is excluded by our work. The two other scenarios, where the turnaround point is at a d• < 3, or where the allowed region is not connected, are still allowed.</figcaption></figure>

Second, even if there is no turnaround, we cannot a priori exclude an allowed (d, N) region which is disconnected from the region at large N. These two still allowed scenarios (as well as the excluded one) are illustrated in Fig. 13. We are however adamant that the CFTs in the allowed region for integer d and N, being unitary, cannot be of focus type, but should have real correction-to-scaling exponent ω.

**Acknowledgments**

We thank Ning Su for discussions related to Remark 2.2, for sharing with us the skydive code before it was publicly released [50], as well as for the simpleboot framework software [75], and for his invaluable support in using these codes. SR thanks Yu Nakayama, Hugh Osborn, Marco Serone and Andreas Stergiou for discussions related to Appendix C, and Yin-Chen He for discussions related to Appendix A.

**Funding information** This work is supported by the Simons Foundation grants 733758 and 488659 (Simons Bootstrap Collaboration). A part of this work was completed during the Bootstrap 2023 program at the Instituto Principia (São Paulo). BvR acknowledges funding from the European Union (ERC "QFTinAdS", project number 101087025). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. The work of MR was also supported by the UK Research and Innovation (UKRI) under the UK government's Horizon Europe funding Guarantee (grant number EP/X042618/1), and the work of BS was also supported by a *Fonds de Recherche du Québec – Nature et technologies* B2X Doctoral scholarship. Some of the computations presented here were conducted in the Resnick High Performance Computing Center, a facility supported by Resnick Sustainability Institute at the California Institute of Technology.

**A Frustrated magnets**

The theoretical study of noncollinear magnets dates back to the 1970s ([1, 76–80], reviewed e.g. in [81]). It is interesting to compare them to the usual ferromagnets described by the fer-
<figure><img src="image_18.png" /><figcaption>Figure 14: The ordered states of an STA (figure from [1]). The shown states have opposite chirality. There is a continuum of states obtained from the shown ones by a global rotation of all spins.</figcaption></figure>

romagnetic O(N) model (N-component spins at each vertex of a lattice, with nearest-neighbor ferromagnetic interactions). The usual ferromagnets have a disordered state at high T and, in d > 2, an ordered state at low T, which breaks O(N) spontaneously to O(N − 1) (“all spins point in the same direction”). The critical behavior is the standard O(N) universality class, also known as the Wilson-Fisher O(N) fixed point, or the O(N) model CFT.

In contrast to the usual ferromagnets, noncollinear magnets have couplings between spins which involve some frustration. This leads to low-temperature states with a more complicated ordering. The symmetry breaking structure is different and the phase transition, if second order, is expected to be in a universality class distinct from O(N).

Frustration is achieved by introducing some degree of antiferromagnetism, i.e. changing the sign of some nearest-neighbor couplings. On the cubic lattice with only nearest-neighbor couplings, every spin will be exactly anti-aligned with its neighbors at low-T, so there is no frustration and the phase transition is still in the O(N) universality class. True frustration may appear on a lattice involving triangular faces.

One type of frustrated magnets which exhibits critical behavior distinct from O(N) are the so-called *stacked triangular antiferromagnets* (STA). Spins are placed at the vertices of planar triangular 2d lattices stacked on top of one another in the third direction. Within each plane, nearest-neighbor interactions $J_\parallel$ are antiferromagnetic. The sign of interactions between the spins on nearby planes, denoted $J_\perp$, is not important. For continuous spins, possible ordered states at low-T will have the structure in each plane shown in Fig. 14, with spins at nearby sites rotated by 120°.

Frustration can also be obtained directly on the cubic lattice, by adding to the ferromagnetic O(N) model a competing antiferromagnetic next-to-nearest neighbor interaction along a given direction. This construction models so-called *helical magnets*, and together with STAs, they have been referred to in the literature as *noncollinear magnets*.

Both systems near criticality enjoy the same continuum description. The procedure to obtain the Landau-Ginzburg-Wilson (LGW) Hamiltonian from a lattice model is straightforward and quite short (see Appendix A of [1]). This Hamiltonian is written in terms of the coarse-grained magnetization field $\vec{S}(x) \in \mathbb{R}^N$. The interaction, assumed to be a polynomial in $\vec{S}$, is expanded around every minimum of the quadratic part of the Hamiltonian density in momentum space. For the usual ferromagnets described by the O(N) model, there is only one minimum in momentum space, located at $p = 0$. This leads to the well-known $\phi^4$-theory in terms of a single N-component field. In contrast, for noncollinear magnets, one finds two minima at two inequivalent Brillouin zone momenta $\pm Q$. The effective Hamiltonian is written in terms of a complex field $\vec{\phi}(x) = e^{iQx}\vec{S}(x)$. O(N) invariance and momentum conservation shows that the numbers N and $N_*$ of $\vec{\phi}$'s and $\vec{\phi}^*$'s should satisfy the relation $N - N_* = 0 \pmod{6}$.<fnref n="27" /> This implies that the effective Hamiltonian has symmetry O(N) $\times \mathbb{Z}_6$. At the level of the quartic effective Hamiltonian, $\mathbb{Z}_6$ gets enhanced to a continuous O(2). One thus obtains

<footnote n="27">²⁷To derive this relation, one observes that 3Q is in the reciprocal lattice.</footnote>
the Hamiltonian (1) or, by writing $\varphi = a + ib$ where $a, b$ are real $N$ component fields,

$$ \mathcal{H} = \frac{1}{2}((\partial_{\mu}a)^2 + (\partial_{\mu}b)^2) + r_0(a^2 + b^2) + u(a^2 + b^2)^2 + v((a \cdot b)^2 - a^2b^2). \quad (A.1) $$

One may wonder why such anisotropic systems as the STAs and helimagnets are described near criticality by the isotropic effective Hamiltonian (A.1)? Let us consider for example the case of STAs. After performing every step listed above in going from the microscopic to the effective Hamiltonian, one ends up with a kinetic term of the schematic form $\int d^d q (|J_\perp|q_\parallel^2 + |J_\perp|q_\perp^2)(a(q) \cdot a(-q) + b(q) \cdot b(-q))$, where $\parallel, \perp$ refers to directions inside of/perpendicular to the triangular lattice planes [82]. The anisotropy depending on the relative strength of the intra-plane and inter-plane couplings $J_\parallel$ and $J_\perp$ is then scaled away by rescaling momenta.

Let us describe the phase transitions expected in (A.1). As the mass $r_0$ is tuned, (A.1) undergoes a phase transition from the disordered state with $\langle a \rangle = \langle b \rangle = 0$ to a nontrivially ordered state whose structure crucially depends on the sign of the coupling $v$. Namely one finds $a \perp b$ for $v > 0$, while $a \parallel b$ for $v < 0$. The first possibility reproduces the noncollinear chiral behavior of Fig. 14, while the second describes so-called spin-density waves. In mean field theory the sign of the coupling $v$ is arbitrary, while in RG theory the sign of $v$ will be the one corresponding to the stable fixed point. Since noncollinear magnets have low-T ground states like in Fig. 14, we expect that the RG fixed point describing their phase transition (if second-order) should have $v > 0$. Indeed the stable fixed point $C_+$ in Fig. 1 is located at $v > 0$.

# B Critical remarks about the focus fixed points

## B.1 Problems with unitarity

The main reason to doubt focus fixed points is that they contradict unitarity. Our model is unitary, i.e. reflection positive in the Euclidean signature. The critical points, if they exist, should be consistent with unitarity. This would imply that the scaling dimensions of scalar operators must be real. Instead, the focus fixed points have a complex correction-to-scaling exponent $\omega$. This means that there exists an operator $\mathcal{O}$ with a complex scaling dimension $\Delta = d + \omega$, whose two point function is given at large distances by

$$ \langle \mathcal{O}(0)\mathcal{O}(x) \rangle \sim \frac{1}{r^{2\Delta}}. \quad (B.1) $$

The complex conjugate operator $\mathcal{O}^*$ will have the complex conjugate correlator. Assuming conformal invariance (see Appendix C), only operators of equal scaling dimensions have non-vanishing 2pt functions, hence $\langle \mathcal{O}(0)\mathcal{O}^*(x) \rangle$ vanishes at large distances. Now consider a real operator $\Phi = \mathcal{O} + \mathcal{O}^*$. Its 2pt function is given by

$$ \langle \Phi(0)\Phi(x) \rangle \sim \frac{1}{r^{2\Delta}} + \frac{1}{r^{2\Delta^*}} = \frac{2 \cos(2b \log r)}{r^{2a}}, \quad (B.2) $$

where $a = \text{Re}\Delta$, $b = \text{Im}\Delta$. If $b \neq 0$, there are distances at which $\langle \Phi(0)\Phi(x) \rangle < 0$. This contradicts reflection positivity.

## B.2 Merger curve and bifurcation theory

As discussed in the introduction, both the NPRG and the $\epsilon$-expansion predict a monotonic single-valued curve $N = N_c(d)$, while some fixed-dimension (FD) RG schemes predict a more complicated S-shaped curve (see Fig. 2).
As was emphasized by [69], the curve along which two fixed points (for us, $C_+$ and $C_-$) merge and annihilate (the merger curve) corresponds to a saddle-node bifurcation in the space of RG flows. This bifurcation is known to be a stable codimension-1 bifurcation for any number of parameters. In particular, for the 2-dimensional parameter space $(d, N)$, we expect the saddle-node bifurcation to occur, generically, along a smooth 1-dimensional line in the parameter space. Nothing can be predicted about the shape of this line on the basis of bifurcation theory alone. Our merger curve could as well be monotonic or S-shaped.

However, FD not only predicts the S-shape of the merger curve, but also focus fixed points. How can the appearance of focus fixed points be accommodated by bifurcation theory? A hint can be found in [83] which points out the existence on the merger curve of a special point S where both eigenvalues of the stability matrix vanish: $\omega_1 = \omega_2 = 0$, whereas at generic points on this boundary only one eigenvalue vanishes. Using this hint, we identify the point S as a Bogdanov-Takens (BT) bifurcation [84]. It is a codimension-2 bifurcation with the property that the stability matrix (for two dimensional flows) has the form:

$$ \begin{pmatrix} 0 & A \\ 0 & 0 \end{pmatrix}, \quad A \neq 0. \qquad (B.3) $$

By a change of coordinates, the BT bifurcation flow can be brought to the normal form:

$$ \dot{y}_1 = y_2, \qquad (B.4) $$

$$ \dot{y}_2 = \beta_1 + \beta_2 y_1 + y_1^2 + \sigma y_1 y_2, $$

where $\sigma = \pm 1$. In what follows we consider $\sigma = +1$.<fnref n="28" /> Here $y_1, y_2$ should be thought of as coordinates in the $u, v$ space near the point S, while $\beta_1, \beta_2$ are coordinates in the parameter space $(d, N)$ near $(d_S, N_S)$. In the new coordinates, the point S is realized for $\beta_1 = \beta_2 = 0$ and is located at $y_1 = y_2 = 0$.

The phase portraits of RG flows around the BT bifurcation are shown in Fig. 15.<fnref n="29" /> To the right of the point S, at $\beta_1 > \beta_2^2/4$, we have region ①, in which there are no real fixed points. In terms of N and d this is the region delineated by the merger curve, on which S is a distinguished point. Rotating counterclockwise, we cross the part of the merger curve labeled $T_+$ beyond which two fixed points appear, one stable and one unstable. So far this is just a saddle-node bifurcation, but if we go a bit further we cross a new curve $L_+$ on which a stable node-to-spiral transition happens (the two real eigenvalues $\omega_1, \omega_2$ collide and go into the complex plane).<fnref n="30" /> In region ③ we have a saddle and a stable spiral. Rotating further still, we reach curve P, where a (global) saddle homoclinic bifurcation happens. Namely, on P we have a homoclinic orbit which encircles the stable spiral fixed point and joins the saddle point to itself. On the other side of P, in region ④, the homoclinic orbit turns into an unstable cycle which limits the basin of attraction of the stable fixed point (which is still a spiral). This unstable cycle shrinks to a point on line H (the Andronov-Hopf bifurcation), and in region ⑤

$$ \beta_1 = 8\beta_2 + (16 - 2\beta_2)\sqrt{4 - \beta_2} - 32 = \beta_2^2/4 - \beta_2^4/1024 + O(\beta_2^5). $$

The distance between L and T was artificially increased in Fig. 15 for readability. It is significant that region ② (as well as region ⑥), where all exponents are real, is non-empty. Indeed, the direct annihilation of a spiral and a saddle is non-generic. The spiral has to first turn into a node, and only then annihilate. Region ② grows as one gets further away from S, see e.g. [87].

<footnote n="28">²⁸ $\sigma = -1$ is obtainable by flipping the sign of $y_2$ and $t$. Flipping the sign of $t$ would turn stable foci into unstable ones, which is not what we need.</footnote>

<footnote n="29">²⁹ Many of these features were mentioned in [83], however the identification with the BT bifurcation is made here for the first time. Previously, BT bifurcations were found in perturbative RG flows of a scalar-fermion model in $d = 3 - \epsilon$ dimensions with $O(N) \times O(M)$ symmetry, $\phi^6$ interactions, and $N = 1$ supersymmetry, for noninteger values of N and M [86].</footnote>

<footnote n="30">³⁰ Near S, curve L lies very close to curve T, its equation being</footnote>
<figure><img src="image_19.png" /><figcaption>Figure 15: Phase portraits of RG flows in the vicinity of a BT bifurcation. See similar figures in [84] and [85], Fig. 137. Compared to [84], we added the lines $L_±$ which are the loci of the spiral to node transitions. (Note that these are technically not bifurcations but they are physically important.)</figcaption></figure>

we have an unstable saddle and an unstable spiral. The unstable spiral turns into an unstable node along $L_-$.

### B.3 Further problems of focus fixed points

Using the previous subsection, we can point out the following further problems raised by focus fixed points. A basic principle of RG theory is the gradient flow property (see footnote 5), which says that the beta function equations should have the form

$$ \beta^I = G^{IJ} \partial_J A, \qquad (B.5) $$

where A is a scalar function of the couplings and $G^{IJ}$ is a symmetric, positive-definite metric in the coupling space. This property is widely expected, although not yet proven. This property was originally found to hold up to three loops in a general multiscalar model in $d = 4 - \epsilon$ dimensions [15, 16]. In $d = 4$ dimensions, a weaker form of the gradient flow property is proven to all orders of perturbation theory [88], where the metric $G^{IJ}$ may have an antisymmetric part starting from four loops (it is not known if the antisymmetric part actually occurs).31 Later the stronger form of the gradient property, with a symmetric metric, was checked at four loops [89] and, very recently, at five and six loops both in $d = 4$ and $d = 4 - \epsilon$ dimensions [17], using the six-loop beta function computed in [18].

If the gradient flow property holds, the fixed points are identified with the critical points of the potential function A. The stability matrix at a fixed point is then given by:

$$ \Gamma_K^I = \partial_K \beta^I = G^{IJ} \partial_J A. \qquad (B.6) $$

<footnote n="31">31 It is also recognized that for the property to hold, the beta function should be computed in a particular “gauge”, where it is known as the B-function [88]. The B-function fixed the beta-function ambiguity which exists in any multi-field theory, in which a renormalization group step can be accompanied by an infinitesimal rotation in the field space.</footnote>
Although this stability matrix is not necessarily symmetric, it is a product of two symmetric matrices. The eigenvalue problem $\Gamma_K^I v^K = \lambda v^I$ for the stability matrix is equivalent to the generalized eigenvalue problem

$$ (\partial_J \partial_K A) v^K = \lambda G_{JK} v^K. \qquad (B.7) $$

This implies that all eigenvalues $\lambda$ must be real. This excludes focus fixed points.

Another consequence of the gradient flow property is that A changes monotonically along the RG flow:<fnref n="32" />

$$ \dot{A} = -\beta^I \partial_I A = -G^{IJ} \partial_I A \partial_J A < 0. \qquad (B.8) $$

This excludes cyclic RG flows. We have seen in the previous section that the RG flow diagram around a BT bifurcation contains RG cycles in region ④ - a further conflict.

Our final remark concerns the fact that focus fixed points appear when two correction-to-scaling critical exponents $\omega_1, \omega_2$ collide and then go to the complex plane. This happens along the curve L in Fig. 15. Such level crossings between scaling dimensions of operators having the same symmetry are normally not allowed without finetuning. See [58, 90, 91] for recent discussions.

All of the above reasons lead us to believe that the focus type chiral fixed points found in [11–14] are artifacts of uncontrolled resummations. For previous exchanges on the controversy surrounding focus fixed points in the context of the $O(N) \times O(2)$ model, see [4, 10, 92, 93].<fnref n="33" />

## C Conformal invariance of multiscalar models

In this appendix we recall why multiscalar fixed points are believed to be conformally invariant and not just scale invariant. The multiscalar fixed points in any d are local, i.e. they have a local conserved stress tensor $T_{\mu\nu}$ of scaling dimension d. The condition for a local fixed point in $d > 2$ to be scale invariant without being conformal is that the trace of the stress tensor be a divergence of a local vector operator $V_\mu$ known as a virial current [95]

$$ T_{\mu}^{\mu} = \partial_{\mu} V^{\mu}, \qquad (C.1) $$

where $V^\mu \neq J^\mu + \partial_\nu L^{\mu\nu}$ with $J_\mu$ a conserved current ($\partial_\mu J^\mu = 0$) and $L_{\mu\nu} = L_{\nu\mu}$. Here and below, we do not keep track of the terms in $T_\mu^\mu$ which vanish by the equations of motion.

Since scaling dimensions of non-conserved vector operators are not protected, we do not generically expect an interacting fixed point to have a virial current candidate of precisely dimension $d-1$ so that it can appear in the r.h.s. of (C.1). Hence generically scale invariance in presence of interactions implies conformal invariance [96, 97].<fnref n="34" />

For multiscalar fixed points one can also give an argument not relying on the genericity assumption. The argument below is in perturbation theory in $d = 4-\epsilon$ dimensions (to all orders in $\epsilon$). Consider a multiscalar theory of n real scalar fields with a general quartic interaction:

$$ \frac{1}{4!} \lambda_{ijkl} \phi_i \phi_j \phi_k \phi_l =: \lambda^I \mathcal{O}_I. \qquad (C.2) $$

Here we introduced the notation $\lambda^I$ for all couplings. The stress tensor trace along the RG flow has the form

$$ T_{\mu}^{\mu} = \beta^{I} \mathcal{O}_{I} + \partial_{\mu} J_{\nu}^{\mu} + a_{ij} \partial^{2}(\phi_{i} \phi_{j}), \qquad (C.3) $$

<footnote n="32">³²This would hold even if $G^{IJ} \rightarrow G^{IJ} + B^{[IJ]}$ in (B.5), as the antisymmetric part $B^{[IJ]}$ drops out.</footnote>

<footnote n="33">³³It has also been reported to us that these problematic fixed point also have unstable quartic potentials, which would be yet another reason to discard them [94]. We thank Marco Serone for discussions.</footnote>

<footnote n="34">³⁴Some interacting models with scale invariance without conformal invariance do exist [98–100]. They all have an additional noncompact “shift” symmetry which protects the virial current dimension [100].</footnote>
where $\beta^I = d\lambda^I/d\log\mu$ is the beta function, and $J_\nu^\mu$ is a flavor current associated with the $O(n)$ Lie algebra element $\nu$, $J_\nu^\mu = \nu_{ij}\phi_i \overleftrightarrow{\partial_\mu}\phi_j$. The last term $a_{ij}\partial^2(\phi_i\phi_j)$ corresponds to “improvable terms” [95] and we will omit it in the discussion below since by itself it does not preclude the theory from being conformal.

At the fixed point we have $\beta_I = 0$. It looks like the theory may not be conformal if the second term $\partial_\mu J_\nu^\mu$ appears in (C.3) with a $\nu$ corresponding to a broken generator of $O(n)$. We do not know of a completely general argument to all orders in $\epsilon$ that this may not happen. For arguments at low orders in perturbation theory see [95, 101].$^{35}$

An all-order argument can however be given for fixed points with a sufficiently large symmetry group $G \subset O(n)$. Namely, suppose that the representation in which the broken currents transform does not contain a singlet. Then, since the stress tensor is a singlet, the broken currents cannot appear in the r.h.s. of (C.3), and scale invariance implies conformal invariance.

For the $O(N) \times O(2)$ model studied in this paper, the broken currents transform in the bifundamental representation, and the above argument applies.

## D 1/N expansion

We review here the 1/N expansion of model (1) [1, 102–104], and compute the leading 1/N correction to the dimension of ST'. The Lagrangian in the form relevant to the large-N expansion reads [104]

$$ \mathcal{L} = \frac{1}{2}(\partial_{\mu}\phi)^2 + \frac{1}{2}\sigma\phi_{ai}\phi_{ai} + \frac{1}{2}T_{ij}\phi_{ai}\phi_{aj} - \frac{\sigma^2}{16u + 4v} - \frac{1}{8v}\text{Tr}T^2, \quad (D.1) $$

where we introduced two auxiliary fields $\sigma$ and $T_{ij}$. They are scalars belonging to the SS and ST irreps, and they morally replace $\phi^2$ and $\phi_{ai}\phi_{aj} - \frac{1}{2}\phi^2\delta_{ij}$. They have scaling dimensions

$$ \Delta_{\sigma}, \Delta_{T_{ij}} = 2 + O(1/N), \quad (D.2) $$

in any $d$, in contrast with the classical dimension $\Delta = d-2$ of the fields they replace. Scaling dimensions and OPE coefficients may be obtained as expansions in $1/N$, the coefficients of these expansions being computable by standard diagrammatic techniques.

At large $N$, all four fixed points on the left of Fig. 1 exist. They are distinguished by which of the auxiliary fields propagate: at $G$, of course none does; at $H$, only $\sigma$ does, leading to the well known $1/N$ expansion of the $O(N)$ model; at $C_-$, only $T_{ij}$ propagates; and finally, both auxiliary fields propagate at $C_+$. The difference in field content between $C_+$ and $H$ leads to an appreciable difference between the large-$N$ spectrum of both theories, the most obvious differences being in the scalar ST sector. At $C_+$, $O_{ST} = T_{ij}$, followed by two subleading fields $O_{ST'}$ and $O_{ST''}$ with $\Delta = 4+O(1/N)$ resulting from the mixing of $(T^2)_{ij} - \frac{1}{2}\text{Tr}T^2\delta_{ij}$ and $\sigma T_{ij}$. In contrast, $T_{ij}$ is absent at $H$, so that $O_{ST} = \phi_{ai}\phi_{aj} - \frac{1}{2}\phi^2\delta_{ij}$, followed by $O_{ST'} = \sigma O_{ST}$ with $\Delta = d+O(1/N)$.

Although imposing a gap of the order of $d-2$ in the ST channel might have been enough to differentiate between both theories, we have found more constraining power by assuming the existence of $O_{ST}$, and imposing a gap above it. To motivate the value of this gap, it helps to know the $1/N$ correction to the dimension of $O_{ST'}$ in $H$ and $C_+$. In $H$, this operator descends from the subleading symmetric traceless 2-tensor operator of $O(N)$, which has been computed long ago in [105]:

$$ \Delta_{ST'} = \Delta_{ST'}^{(1)} + O(1/N^2), \quad \Delta_{ST'}^{(1)} = d - \frac{2(d^2 - 3d + 4)\Gamma(d-1)}{\Gamma(1 - \frac{d}{2})\Gamma(\frac{d}{2} + 1)\Gamma(\frac{d}{2})^2 N} \quad (H) \quad (D.3) $$

<footnote>$^{35}$S.R. thanks Yu Nakayama, Hugh Osborn and Andy Stergiou for discussions.</footnote>
In $C_+$, this computation requires resolving the mixing problem between $O_{ST'}$ and $O_{ST''}$. Ref. [106] gives a nice and concise description of the procedure to follow to compute anomalous dimensions of composite operators at large-$N$. The example provided in this paper lists and explicitly computes all diagrams contributing to the anomalous dimension of $\sigma^2$ in the $O(N)$ model. Since $C_+$ differs from $\mathcal{H}$ by the addition of a new $T_{ij}$ propagator as well as a new $T\phi\phi$ vertex [102]:

$$ T_{ij}\phi_{ak}\phi_{bl} \rightarrow \delta_{ab}K_{ijkl}, \qquad (D.4) $$
$$ K_{ijkl} = \frac{1}{2}(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk} - \delta_{ij}\delta_{kl}), $$

it should be clear that we can recycle some of the results of [106]. Indeed, the diagrams contributing to the renormalization of $(T^2)_{ij}$—trace and $\sigma T_{ij}$ will have the same geometry as those contributing to $\sigma^2$. Here is an example of such a diagram arising in the computation of the anomalous dimension of $\sigma^2$:

<figure><img src="image_20.png" /></figure>

One can think of this diagram as computing the renormalization of $\sigma^2$ by considering the three point function $\langle\sigma^2\sigma\sigma\rangle$. Diagram (D.5) has the following two nonzero cousins in $C_+$:

<figure><img src="image_21.png" /></figure>

These are determinable from Diagram (D.5) up to $O(2)$ group structure factors, coming from $T_{ij}$ propagators and $T\phi\phi$ vertices, and differences in symmetry factors. There is luckily no need to compute any loop integrals. After the dust settles, the leading order anomalous dimensions of this mixing problem are given by$^{36}$

$$ \begin{aligned} \gamma_1 &= \frac{4(12-22d+13d^2-2d^3)\sin(\frac{\pi d}{2})\Gamma(d-2)}{\pi d \Gamma(\frac{d}{2}-1)\Gamma(\frac{d}{2})} \frac{1}{N}, \\ \gamma_2 &= \frac{2^{d-3}(24-28d+5d^2+d^3)\sin(\frac{\pi d}{2})\Gamma(\frac{d-1}{2})}{\pi^{3/2}\Gamma(\frac{d}{2}+1)} \frac{1}{N}. \end{aligned} \qquad (D.7) $$

The scaling dimensions of $O_{ST'}$ in $C_+$ is then given by

$$ \Delta_{ST'} = \Delta_{ST'}^{(1)} + O(1/N^2), \qquad \Delta_{ST'}^{(1)} = 4 + \min(\gamma_1, \gamma_2) + O(1/N^2). \qquad (D.8) $$

# E $\epsilon$-expansion results

## E.1 $N_c(d)$

In the $\epsilon$-expansion, the location $(u^*, v^*)$ of fixed points and critical exponents in $d-4-\epsilon$ are obtained as asymptotic series expansions in $\epsilon$. To locate the critical value of $N$ where $C_+$

<footnote>$^{36}$We checked this result by comparing it at $O(\epsilon, 1/N)$ to the anomalous dimensions of the corresponding operators in the $\epsilon$-expansion. These do not seem to have been previously computed, but they are easily obtained using the OPE method described e.g. in [54].</footnote>
merges and annihilates with $C_-$, one can ask that the location of both fixed points coincide:

$$u_{-}^{*}(\epsilon, N_c(4-\epsilon)) = u_{+}^{*}(\epsilon, N_c(4-\epsilon)), \qquad (E.1)$$

$$v_{-}^{*}(\epsilon, N_c(4-\epsilon)) = v_{+}^{*}(\epsilon, N_c(4-\epsilon)). \qquad (E.2)$$

Equivalently, one can look for the value of $N$ at which the critical exponent $\omega$, related to the scaling dimension of the subleading SS scalar, is zero:

$$\omega_{\pm}(\epsilon, N_c(4-\epsilon)) = 0. \qquad (E.3)$$

In both cases, one obtains $N_c(4-\epsilon)$ as a series expansion in $\epsilon$. We have ([6], Table I, first line)$^{37}$

$$N_c(4-\epsilon) = 21.798 - 23.431\epsilon + 7.0882\epsilon^2 - 0.0321\epsilon^3 + 4.2650\epsilon^4 - 8.4436\epsilon^5 + O(\epsilon^6). \quad (E.4)$$

This series is divergent, and to obtain the best possible estimate of $N_c(d)$ one would need to resum it. There are many ways to do so, see [6] for investigations on how to do so optimally. Here we used the method of [6], Eqs. (25),(26) which goes back to [107], Eq. (7). Rather than working with the original series, we first input the expected $d=2$ behavior:

$$N_c(4-\epsilon) = 2 + (2-\epsilon)a(\epsilon), \qquad (E.5)$$

and then expand $1/(a(\epsilon))$. The coefficients of the series for $1/(a(\epsilon))$ are better behaved than those of the original series, and we sum it directly, using the contribution of the last term as the error estimate. The so obtained $N_c(d)$ will be denoted $N_c^\epsilon(d)$. This method is simple and gives results in good agreement with other methods, based on Padé and Borel resummation as in [6]. This procedure was also used to compare the $\epsilon$-expansion to the results of the NPRG in [108].

## E.2 Scaling dimensions

In several parts of the paper we referred to the values of the critical exponents evaluated in the $\epsilon$-expansion. We first explain the general methodology used in the literature [6, 66, 102].

The $\epsilon$-expansion is easiest performed for $N > N_c(4) = 4(3 + \sqrt{6})$. In this range of $N$ the fixed point couplings, and therefore the critical exponents, can be expanded as (asymptotic) power series in integer powers of $\epsilon$ with real, $N$-dependent, coefficients. However in several places in this paper we need the critical exponents for $N < N_c(4)$. In this range of $N$ we cannot simply use the series computed for $N > N_c(4)$, since those series have some coefficients proportional to $\sqrt{N - N_c(4)}$, and become complex for $N < N_c(4)$.

At the formal level, what is happening is that the assumption of expanding the couplings in integer powers of $\epsilon$ breaks down for $N = N_c(4) - O(\epsilon)$. One can see for example that for $N = N_c(4-\epsilon)$ given by (E.4) the fixed point couplings have an expansion in half-integer powers. On general grounds, we expect that for $N = N_c(4-\epsilon) + \delta N$ with $\delta N > 0$, the critical exponents should have dependence on $\delta N$ involving positive half-integer powers. We can determine this expansion in $\delta N$ by substituting $N = N_c(4-\epsilon) + \delta N$ into the series derived for $N > N_c(4)$, and reexpand in $\epsilon$. This was proposed below Eq. (3.29) of [102], referred there as a “trick”, although it’s not a trick but a mathematically legitimate way to match two asymptotic expansions with overlapping ranges of validity. As a consistency check, one observes that all terms involving negative half-integer powers of $\delta N$, which could appear a priori when expanding $\sqrt{\delta N + O(\epsilon)}$ in $\epsilon$, cancel order by order in the expansion in $\epsilon$.

<footnote>$^{37}$The first two terms are given exactly by $12+4\sqrt{6}-(12+14\sqrt{2}/3)\epsilon$. This is useful for checking the cancellation of $1/\sqrt{\delta N}$ singularities in the reexpansion of critical exponents discussed in Appendix E.2.</footnote>
One thus obtains critical exponents as a series in $\epsilon$ (with integer powers) with coefficients
real functions of $\delta N$ with at most positive square-root singularities at $\delta N = 0$. Substituting,
for a given $N$, $\delta N = N - N_c^\epsilon(d)$ where $N_c^\epsilon(d)$ is as in Appendix E.1, we get critical exponents
as a series in $\epsilon$ with numerical coefficients. For small $\epsilon$, these series can be used by directly
summing the first few terms. For larger $\epsilon$, such as in $d = 3$, these series need to be Borel or
Padé resummed to be useful.

We now give precise references to sources for each figure.

In Fig. 5, to position $\mathcal{H}, C_\pm$ we used expansion of critical exponents listed in the ancillary
file of [25]. This corresponds to order $\epsilon^3$ for $\Delta_\phi$ and to order $\epsilon^2$ for $\Delta_{SS}$ and $\Delta_{ST}$. In their
conventions:

* $\Delta_{\phi}$ is deltaPhiChiralEps($C_+$), deltaPhiAntiEps($C_-$), deltaPhiONEps($\mathcal{H}$)

* $S \equiv SS$ and $W \equiv ST$

* $\Delta_{SS}$, $\Delta_{ST}$ is deltaScalarsChiralEps[[i]]($C_+$), deltaScalarsAntiEps[[i]]($C_-$), deltaScalarsONEps[[i]]($\mathcal{H}$), where $i = 1,2$ for $\Delta_{SS}$, $\Delta_{ST}$.

As described above, we substitute $N = N_c(4 - \epsilon) + \delta N$, reexpand in $\epsilon$, and subsitute
$\delta N = N - N_c^\epsilon(d)$. Since we are at small $\epsilon = 0.2$, the resulting series is simply summed.

In Fig. 7, we obtained, as described in the previous paragraph, $\Delta_{\phi}$, $\Delta_{SS}$ and $\Delta_{ST}$ at the $C_+$
fixed point as a series in $\epsilon$, to order $\epsilon^2$ for $\Delta_{SS}$ and $\Delta_{ST}$ and to order $\epsilon^3$ for $\Delta_{\phi}$, with coefficients
depending on $\delta N$. We plugged $\epsilon = 0.2$ and plotted the resulting function, too bulky to report
here, as a function of $\delta N = N - N_c^\epsilon(3.8)$.

In Fig. 6, since we are dealing with $\epsilon = 1$, the $\epsilon$-expansion series from [25] are too short
to be useful, and anyway they would have to be resummed. Fortunately for $C_+$ this has been
done elsewhere in the literature, so we borrow the available results. Thus, for $C_+$, we ex-
tract the scaling dimensions of $\phi$ and SS from the $n = 8$ “Final” values of $\eta = 0.042(2)$ and
$\nu = 0.745(11)$ in Table 9 of [6]. These critical exponents are related to the scaling dimen-
sions as $\Delta_{\phi} = \frac{d-2}{2} + \frac{\eta}{2}$ and $\Delta_{SS} = d - \frac{1}{\nu}$. $\Delta_{ST}$ for $C_+$ was obtained from the critical exponent
$y_4 = 1.13(8)$ in the $\overline{MS}$ scheme of Table III of [66]. The two are related as $\Delta_{ST} = 3 - y_4$.
We have not found in the literature equally precise resummed estimates for $C_-$, or for $\mathcal{H}$ (at
$N = 8$ that we need for Fig. 6). For them, we used the large-$N$ results at the highest known
order. For $C_-$, this is $1/N$ order, and the results are tabulated in the ancillary file of [25]. For
$\mathcal{H}$, $\Delta_{\phi}$ is known to order $1/N^3$ [109], and the other two to order $1/N^2$ [110, 111]. The three
are nicely tabulated in the ancillary file of [112].

In Fig. 6, since we are dealing with $\epsilon = 1$, the $\epsilon$-expansion series has to be resummed to be useful. Fortunately for $C_+$ this has been done elsewhere in the literature, so we borrow the available results. Thus, for $C_+$, we extract the scaling dimensions of $\phi$ and SS from the $n = 8$ “Final” values of $\eta = 0.042(2)$ and $\nu = 0.745(11)$ in Table 9 of [6]. These critical exponents are related to the scaling dimensions as $\Delta_{\phi} = \frac{d-2}{2} + \frac{\eta}{2}$ and $\Delta_{SS} = d - \frac{1}{\nu}$. $\Delta_{ST}$ for $C_+$ was obtained from the critical exponent $y_4 = 1.13(8)$ in the $\overline{MS}$ scheme of Table III of [66]. The two are related as $\Delta_{ST} = 3-y_4$. We have not found in the literature equally precise resummed estimates for $C_-$, or for $\mathcal{H}$ ($N=8$ that we need for Fig. 6). For these cases, we instead used the large-$N$ results at the highest known order. For $C_-$, this is $1/N$ order, and the results are tabulated in the ancillary file of [25]. For $\mathcal{H}$, $\Delta_{\phi}$ is known to order $1/N^3$ [109], and the other two to order $1/N^2$ [110, 111]. The three are nicely tabulated in the ancillary file of [112].
# F Numerical bootstrap details

## F.1 The navigator method

The central problem in the numerical conformal bootstrap is classifying parameter space into allowed and disallowed regions. With the navigator function technology [49] this is done by splitting the problem into two steps. First step: look for an allowed point by minimizing the navigator $N(x)$ from a starting point $x_0$. In Appendix F.4 we will show how to construct a navigator function for the class of numerical conformal bootstrap problems that we consider in this paper. We recall that this function should be negative for allowed values of $x$, positive for disallowed values, and zero at the boundaries between these regions. It should also be differentiable.

Once an allowed point $x_*$ is found, the second step is: look for the boundaries of the allowed region around $x_*$ which is the zero set of the navigator. Both steps can be phrased as non-convex optimization problems, and algorithms for solving them were given in [49]. A particular second-step problem concerns finding the maximal extent of the allowed region in a fixed direction, which can be phrased as

$$ \text{maximize } n^T x \text{ over all } x \text{ such that } N(x) \le 0. \qquad (\text{F.1}) $$

Here $n$ is a vector describing the direction in parameter space that we want to maximize, while still remaining in the same isolated allowed region to which $x_*$ belongs.

As in the main text, we generally suppress the dependence of the navigator function on any parameters that are held fixed in the optimization process (like $d$ and $N$ in some cases).

## F.2 Details for Figs. 5 and 6

The above strategy is precisely how the islands in Figs. 5 and 6 were “boxed in.” First we located an allowed point in each isolated allowed region, by minimizing $N(x)$ starting from many initial starting points. Subsequently, within each islands we determined the minimal and maximal allowed values for all three parameters $\Delta_\phi, \Delta_{SS}, \Delta_{ST}$ that the navigator function depends on.

We stress again that the optimizations problems encountered in both steps are non-convex. In the first step, the navigator $N(x)$ has multiple minima. All navigator minima which we observed in our searches, starting from many initial points, fell into one of the following three categories:

1. A negative minumum in one of the isolated allowed regions corresponding to an allowed physical CFT (like $H$ or $C_+$).

2. A negative minimum near the location of a CFT that is (barely) excluded by the bootstrap assumptions, like $C_-$ in Section 3.1.

3. As the parameters $N$ and gapST$_0$ vary, some of the isolated allowed regions may disappear, but the corresponding minima may still remain, although they are now positive.

In this sense, all seen navigator minima were related to either a CFT satisfying the bootstrap assumptions or a CFT that “nearly” obeys the bootstrap assumptions.

In the second step, when finding the boundaries of the islands, non-convexity again plays a role. By extremizing from multiple different points and in additional directions $\vec{n}$ we found that at certain derivative orders $\Lambda$ the isolated allowed regions are not convex. Specifically for some $\Lambda$ and for some gapST$_0$, we detected a narrow “bridge” between the $C_+$ and $H$ islands. However, we ascertained that for the $\Lambda$ and gapST$_0$ used in Fig. 5 and Fig. 6 no such bridge
<table><caption>Table 2: Parameters for skydive.</caption><tr><td>Λ</td><td>31</td></tr><tr><td>x</td><td>30</td></tr><tr><td>order</td><td>120</td></tr><tr><td>spins</td><td>{0,...,40} ∪ {46,47,51,52,55,56,59,60}</td></tr><tr><td>precision</td><td>768</td></tr><tr><td>dualityGapThreshold</td><td>10<sup>-20</sup></td></tr><tr><td>primalErrorThreshold</td><td>10<sup>-15</sup></td></tr><tr><td>dualErrorThreshold</td><td>10<sup>-15</sup></td></tr><tr><td>initialMatrixScalePrimal</td><td>10<sup>20</sup></td></tr><tr><td>initialMatrixScaleDual</td><td>10<sup>20</sup></td></tr><tr><td>feasibleCenteringParameter</td><td>0.3</td></tr><tr><td>infeasibleCenteringParameter</td><td>0.3</td></tr><tr><td>stepLengthReduction</td><td>0.7</td></tr><tr><td>maxComplementarity</td><td>10<sup>100</sup></td></tr><tr><td>dualityGapUpperLimit</td><td>0.2</td></tr><tr><td>externalParamInfinitesimal</td><td>10<sup>-40</sup></td></tr><tr><td>findBoundaryObjThreshold</td><td>10<sup>-20</sup></td></tr><tr><td>betaScanMin</td><td>0.3</td></tr><tr><td>betaScanMax</td><td>1.01</td></tr><tr><td>betaScanStep</td><td>0.1</td></tr><tr><td>stepMinThreshold</td><td>0.1</td></tr><tr><td>stepMaxThreshold</td><td>0.6</td></tr><tr><td>primalDualObjWeight</td><td>0.2</td></tr><tr><td>navigatorWithLogDetX</td><td>False</td></tr><tr><td>navigatorAutomaticShift</td><td>False</td></tr><tr><td>gradientWithLogDetX</td><td>True</td></tr><tr><td>betaClimbing</td><td>1.5</td></tr><tr><td>optimalbeta</td><td>False</td></tr><tr><td>stickToGCP</td><td>False</td></tr></table>

exists. To check this we extremized from different starting points and in additional directions $\vec{n}$, including from starting points in each isolated allowed region in the direction of the other allowed regions.

### E.3 The skydive algorithm

For initial explorations of plots like in Figs. 5 and 6, up to $Λ = 27$, we used Algorithm 1 (for the first step) and Algorithm 2 (for the second step) from [49]. For the final computations at $Λ = 31$ we instead used the skydive algorithm [50]. Both of these methods use a BFGS-like quasi-Newton method. The algorithms from [49] required computing the exact navigator function $N(x)$ at each point $x$ along the minimization path. The skydive algorithm operates instead with an approximate navigator function $N_μ(x)$ where $μ$ is an internal SDPB parameter which can be thought of as describing the accuracy of the approximation. We have $\lim_{μ→0} N_μ(x) = N(x)$. In skydive, $μ$ is lowered gradually as the minimization progresses. As a result the total minimization cost in skydive is comparable to a single evaluation of the navigator function. This algorithm offered a significant speed-up, allowing us to use $Λ = 31$ for most of our results, and even use $Λ = 43$ in some checks.

In skydive, the navigator function and its derivative are sometimes estimated poorly in the early steps when $μ$ is large, which can degrade the performance. To mitigate this, it is
essential to choose an appropriate value for the dualityGapUpperLimit. This value sets an upper limit on $\mu$ for each step, which ensures that the approximation of the navigator function is good enough to guide the search towards the minimum of the true navigator function $N(x)$. Setting this value too small leads to a slow search, while setting it too large often results in the search steering off towards the allowed “peninsula” at large external dimensions, despite being in the region of attraction of another minimum of the true navigator. In the first step problem, we found the value of 0.05 to be effective for searches starting relatively far from the allowed point, while 0.005 is appropriate for searches starting close to the navigator minimum (the same holds for searches of the second step).<fnref n="38" />

The complete workflow was as follows. We used scalar_blocks [113] to generate conformal blocks; simpleboot [75] to setup the semi-definite problems and to manage the computations on the cluster; either SDPB [35, 37] to solve the SDPs (i.e. to compute the navigator function) and BFGS-like Algorithms 1 and 2 from [49] (implemented in simpleboot) to move in the parameter space, or the new skydive to perform the last two steps at once [50]. These codes depend on a number of parameters that must be set by the user, the most crucial being the derivative order $\Lambda$ that dictates the dimension of the space of functionals SDPB will search in. For the most important computations at $\Lambda = 31$ we used skydive, with the parameters given in Table 2.

## F.4 Constructing the navigator function

In order to start applying the navigator method in Section 3.1, we have to first make a choice of the particular navigator function we are going to use. Ref. [49] introduced two different navigator constructions: the GFF navigator and the $\Sigma$ navigator. For the entirety of the paper, we chose to work with the GFF navigator, which has an advantage of a natural normalization and an upper bound: $N(x) \le 1$. Say we are considering a given set of crossing equations

$$ \vec{E}(x) = 0, \qquad (F.2) $$

depending on a set of parameters $x$, usually scaling dimensions or ratios of OPE coefficients. We add a term $\lambda \vec{M}_{add}$ to this crossing equation such that the augmented crossing equation

$$ \vec{E}(x) + \lambda \vec{M}_{add} = 0, \qquad (F.3) $$

always has a solution for some positive $\lambda$. The navigator function is defined [49] as the minimal value of $\lambda$ for which a solution exists. Different terms $\vec{M}_{add}$ lead to different navigator functions. The GFF navigator function is obtained by adding the contribution of every operator below the assumed gaps from a solution where each external field is an independent generalized free field (GFF). This idea is completely general and was already used to study the Ising model [49] and the $O(N)$ model [60].

Here we use it for $O(N) \times O(2)$ symmetric CFTs. We have chosen the additional terms to correspond to the solution to crossing where $\phi$ is an $O(2N)$ GFF, and $s$ is another, independent GFF. The $O(2N)$ invariant GFF solution for $\phi$ also has the symmetry of its $O(N) \times O(2)$ subgroup, and thus solves the $O(N) \times O(2)$ crossing equation. To explicitly write $\vec{M}_{add}$ in this case, we need to know how operators in the $O(2N)$ invariant solution decompose under $O(N) \times O(2)$. The nontrivial decomposition rules (including the weights in the conformal block decomposition) are:

$$ T \rightarrow \frac{1}{2}TT + \frac{1}{2}AA + \frac{1}{2}TS + \frac{1}{N}ST, \qquad (F.4a) $$

$$ A \rightarrow \frac{1}{2}AT + \frac{1}{2}TA + \frac{1}{2}AS + \frac{1}{N}SA. \qquad (F.4b) $$

<footnote n="38">³⁸This choice of course strongly depends on the normalization of the navigator functions that is used. Here we are assuming a navigator function that has been normalized to be bounded by 1.</footnote>
The terms in the r.h.s. are the irreps XY of O(N)×O(2) appearing in the decomposition of irreps T,A of O(2N). The coefficients multiplying them are the relative weights needed to pass from the O(2N) invariant conformal block decomposition of $\langle\phi\phi\phi\phi\rangle$ to the O(N) × O(2) invariant one. These are obtained as follows. Consider the O(2N) generalized free field $\phi_I \equiv \phi_{ai}$. Its 4pt function decomposes into contributions of the different O(2N) irreps $\mathcal{R}$ present in the $\phi_I \times \phi_J$ OPE, weighted by tensor structures:

$$ \langle\phi_I\phi_J\phi_K\phi_L\rangle = \langle\phi_{ai}\phi_{bj}\phi_{ck}\phi_{dl}\rangle = T_{IJKL}^S G_S + T_{IJKL}^T G_T + T_{IJKL}^A G_A; \quad (F.5) $$

$$ \begin{gathered} T_{IJKL}^S = \delta_{IJ}\delta_{KL} = \delta_{ij}\delta_{kl}\delta_{ab}\delta_{cd}, \\ T_{IJKL}^A = \delta_{IK}\delta_{JL} - \delta_{IL}\delta_{JK} = \delta_{ik}\delta_{jl}\delta_{ac}\delta_{bd} - \delta_{il}\delta_{jk}\delta_{ad}\delta_{bc}, \end{gathered} \quad (F.6) $$

$$ T_{IJKL}^T = \delta_{IK}\delta_{JL} + \delta_{IL}\delta_{JK} - \frac{2}{2N}\delta_{IJ}\delta_{KL} = \delta_{ik}\delta_{jl}\delta_{ac}\delta_{bd} + \delta_{il}\delta_{jk}\delta_{ad}\delta_{bc} + \frac{2}{2N}\delta_{ij}\delta_{kl}\delta_{ab}\delta_{cd}. $$

We obtain the branching rules (F.4) by decomposing these tensor structures in terms of products $T_{ijkl}^X T_{abcd}^Y$, where XY is an irrep of O(N) × O(2) appearing in the decomposition of $\mathcal{R}$.

## F.5 Minimum-following for Fig. 3

As explained in Section 3.2, finding the minimal allowed value of N for Fig. 3 can be seen as a second-step optimization problem for the navigator function depending on N as the extra fourth parameter. We used the skydive algorithm at $\Lambda = 31$ to solve this problem.

The initial points for the skydive algorithm were obtained using an auxiliary run at $\Lambda = 19$, using the following different strategy inspired by [60], which we call “minimum-following.” In this strategy we first minimize the navigator function with respect to the parameters $x = (\Delta_\phi, \Delta_{SS}, \Delta_{ST})$. We then track the location of this minimum as we vary N and d.

In this section we write $\mathcal{N}(x, N, d)$ to emphasize that our algorithm also moves us around in N and d.

For example, suppose we start at an allowed value of p where $N \gg N_c(d)$, so at the minimum $x_{min}(N, d)$ we find that the navigator function is negative:

$$ \mathcal{N}(x_{\min}(N, d), N, d) < 0. \quad (F.7) $$

If we now decrease N we can follow this minimum until we reach $N_c$, where $\mathcal{N} = 0$. Then, once we have found this $N_c$ at some value of d, we may follow this extremum $N_c$ as we change d from say d near 4 to d = 3.

When we have minimized the navigator function over its internal parameters x, with N close to $N_c(d)$, we can predict the shape of $N_c(d)$ by solving the following equation

$$ \mathcal{N}(x_{\min} + \delta x_{\min}, N + \delta N, d + \delta d) = 0 \quad (F.8) $$

to first order in $\delta N$ and $\delta d$. Note that the linear term in $\delta x_{\min}$ vanishes because we are at the navigator minimum. For example, if we are at a fixed d (so $\delta d = 0$), we can step towards $N_c(d)$ by taking the step

$$ \delta N = - \frac{\mathcal{N}(x_{\min}, N, d)}{\frac{\partial \mathcal{N}}{\partial N}(x_{\min}, N, d)}, \quad (F.9) $$

which is the Newton step following from the first-order expansion of (F.8). The partial derivative in the denominator can be computed via finite differences or with the SDP gradient formula of Section 4 of [49]. Then, once we have determined $N_c$ at this d, we can take a step $\delta d$. To
remain on the critical curve $N_c(d)$ we should take a corresponding step $\delta N = N'_c(d) \times \delta d$, where from (F.8):

$$N'_c(d) = -\frac{\frac{\partial N}{\partial d}(x_{\min}, N, d)}{\frac{\partial N}{\partial N}(x_{\min}, N, d)} \qquad (F.10)$$

This minimum-following method has advantages and disadvantages which make it complementary to the direct method based on extremizing the 4-variable navigator island in the N direction. By construction, it sticks close to the navigator minima, thus exploring less parameter space. This is both a positive, since the searches are less likely to be attracted by another unimportant navigator minimum, such as the peninsula, but also a downside, since it might miss important minima not connected smoothly enough to previous minima. The direct method proved less robust, as it sometimes makes too large a step and escapes from the region of interest (the $C_+$ island) to disconnected regions like the $H$ island or the peninsula. However, when successful, the direct method is more computationally efficient. Also, the direct method can be parallelized, running separate computation for several d's, while the minimum-following method cannot be easily parallelized as it explores the shape of the curve $N_c(d)$ locally in d.

In the end, both methods proved useful. We first used minimum-following to find the entire curve $N_c(d)$ for $\Lambda = 19$ and some particular set of gap assumptions.<fnref n="39" /> To move to higher $\Lambda$ as in Fig. 3, we used the allowed points of zero navigator obtained from minimum-following to initialize parallel skydive computations of $N_c(d)$ for the 9 values $d = 3.0, 3.1, \dots, 3.8$.

## F.6 Gap sensitivity

We finish this appendix by reporting some tests concerning the impact of gap assumptions (7b)-(7d) on $N_c(d)$. Since we have chosen them very conservatively, $N_c(d)$ should not depend much on them. This is indeed what we observe. More precisely, we see that varying these gaps has a moderate effect on the shapes of the island, but a much smaller effect on the value of $N_c(d)$. E.g. we see that $N_c^{CB}$ in $d = 3.8$ varies by 0.02 when varying the gaps above the stress tensor and the conserved currents between 0.7 and 1.2 (these gaps are set to 1 in (7b)-(7d)).

A similar remark concerns the gap$V_0$ assumption (9). We observed (again in $d = 3.8$) hardly any effect on $N_c^{CB}(d)$ when we varied gap$V_0$ between 1.7 to 2.5. On the other hand the negative navigator minimum disappeared when we set this gap to a much larger value gap$V_0 = 3.2$. This behavior is consistent with the estimates $\Delta_{VV'} \approx 2.7-2.9$ in $d = 3.8$ from (8).

Finally, the dependence of the bounds in $d = 3$ on gap$ST_0$ was discussed around Fig. 9 in the main text.

## References

[1] H. Kawamura, Renormalization-group analysis of chiral transitions, Phys. Rev. B **38**, 4916 (1988), doi:10.1103/PhysRevB.38.4916.

[2] H. Kawamura, Universality of phase transitions of frustrated antiferromagnets, J. Phys.: Condens. Matter **10**, 4707 (1998), doi:10.1088/0953-8984/10/22/004.

<footnote n="39">³⁹The assumptions were slightly different than those presented in Section 2.3: we assumed the existence of an isolated $\phi'$, with then a large gap$V_0 = 3 + 1.875(d-3)$. We also imposed more aggressive gaps of 1.2 above currents instead of the 1 of (7b)-(7d), and a large gap$ST_0 = 3.6 + 0.375(d-3)$. With these more aggressive assumptions, we obtained the first critical curve at relatively low $\Lambda = 19$. Later on when we moved to higher $\Lambda$, we relaxed our assumptions to the ones given in the main text.</footnote>
[3] B. Delamotte, D. Mouhanna and M. Tissier, *Nonperturbative renormalization-group approach to frustrated magnets*, Phys. Rev. B **69**, 134413 (2004), doi:10.1103/PhysRevB.69.134413.

[4] B. Delamotte, M. Dudka, Y. Holovatch and D. Mouhanna, *Relevance of the fixed dimension perturbative approach to frustrated magnets in two and three dimensions*, Phys. Rev. B **82**, 104432 (2010), doi:10.1103/PhysRevB.82.104432.

[5] K. Wilson, *The renormalization group and the $\epsilon$ expansion*, Phys. Rep. **12**, 75 (1974), doi:10.1016/0370-1573(74)90023-4.

[6] M. V. Kompaniets, A. Kudlis and A. I. Sokolov, *Six-loop $\epsilon$ expansion study of three-dimensional $O(n) \times O(m)$ spin models*, Nucl. Phys. B **950**, 114874 (2020), doi:10.1016/j.nuclphysb.2019.114874.

[7] A. Aharony, *Critical behavior of amorphous magnets*, Phys. Rev. B **12**, 1038 (1975), doi:10.1103/PhysRevB.12.1038.

[8] P. Bak and D. Mukamel, *Physical realizations of $n \ge 4$-component vector models. III. Phase transitions in Cr, Eu, MnS$_2$, Ho, Dy, and Tb*, Phys. Rev. B **13**, 5086 (1976), doi:10.1103/PhysRevB.13.5086.

[9] T. Garel and P. Pfeuty, *Commensurability effects on the critical behaviour of systems with helical ordering*, J. Phys. C: Solid State Phys. **9**, L245 (1976), doi:10.1088/0022-3719/9/10/001.

[10] B. Delamotte, Y. Holovatch, D. Ivaneyko, D. Mouhanna and M. Tissier, *Spurious fixed points in frustrated magnets*, (arXiv preprint) doi:10.48550/arXiv.cond-mat/0609285.

[11] A. Pelissetto, P. Rossi and E. Vicari, *Critical behavior of frustrated spin models with non-collinear order*, Phys. Rev. B **63**, 140414 (2001), doi:10.1103/PhysRevB.63.140414.

[12] P. Calabrese, P. Parruccini and A. I. Sokolov, *Chiral phase transitions: Focus driven critical behavior in systems with planar and vector ordering*, Phys. Rev. B **66**, 180403 (2002), doi:10.1103/PhysRevB.66.180403.

[13] P. Calabrese, P. Parruccini and A. I. Sokolov, *Critical thermodynamics of a three-dimensional chiral model for $N > 3$*, Phys. Rev. B **68**, 094415 (2003), doi:10.1103/PhysRevB.68.094415.

[14] P. Calabrese, P. Parruccini, A. Pelissetto and E. Vicari, *Critical behavior of $O(2) \otimes O(N)$ symmetric models*, Phys. Rev. B **70**, 174439 (2004), doi:10.1103/PhysRevB.70.174439.

[15] D. J. Wallace and R. K. P. Zia, *Gradient flow and the renormalisation group*, Phys. Lett. A **48**, 325 (1974), doi:10.1016/0375-9601(74)90449-6.

[16] D. J. Wallace and R. K. P. Zia, *Gradient properties of the renormalisation group equations in multicomponent systems*, Ann. Phys. **92**, 142 (1975), doi:10.1016/0003-4916(75)90267-5.

[17] W. H. Pannell and A. Stergiou, *Gradient properties of perturbative multiscalar RG flows to six loops*, Phys. Lett. B **853**, 138701 (2024), doi:10.1016/j.physletb.2024.138701.

[18] A. Bednyakov and A. Pikelner, *Six-loop beta functions in general scalar theory*, J. High Energy Phys. **04**, 233 (2021), doi:10.1007/JHEP04(2021)233.
[19] M. V. Kompaniets and E. Panzer, *Minimally subtracted six-loop renormalization of O(n)-symmetric φ⁴ theory and critical exponents*, Phys. Rev. D **96**, 036016 (2017), doi:10.1103/PhysRevD.96.036016.

[20] V. T. Ngo and H. T. Diep, *Stacked triangular XY antiferromagnets: End of a controversial issue on the phase transition*, J. Appl. Phys. **103**, 07C712 (2008), doi:10.1063/1.2837281.

[21] V. T. Ngo and H. T. Diep, *Phase transition in Heisenberg stacked triangular antiferromagnets: End of a controversy*, Phys. Rev. E **78**, 031119 (2008), doi:10.1103/PhysRevE.78.031119.

[22] Y. Nagano, K. Uematsu and H. Kawamura, *Monte Carlo study of the critical properties of noncollinear Heisenberg magnets: O(3) × O(2) universality class*, Phys. Rev. B **100**, 224430 (2019), doi:10.1103/PhysRevB.100.224430.

[23] M. Tissier, B. Delamotte and D. Mouhanna, *Frustrated Heisenberg magnets: A nonperturbative approach*, Phys. Rev. Lett. **84**, 5208 (2000), doi:10.1103/PhysRevLett.84.5208.

[24] Y. Nakayama and T. Ohtsuki, *Bootstrapping phase transitions in QCD and frustrated spin systems*, Phys. Rev. D **91**, 021901 (2015), doi:10.1103/PhysRevD.91.021901.

[25] J. Henriksson, S. R. Kousvos and A. Stergiou, *Analytic and numerical bootstrap of CFTs with O(m) × O(n) global symmetry in 3D*, SciPost Phys. **9**, 035 (2020), doi:10.21468/SciPostPhys.9.3.035.

[26] R. Rattazzi, V. S. Rychkov, E. Tonni and A. Vichi, *Bounding scalar operator dimensions in 4D CFT*, J. High Energy Phys. **12**, 031 (2008), doi:10.1088/1126-6708/2008/12/031.

[27] F. Caracciolo and S. Rychkov, *Rigorous limits on the interaction strength in quantum field theory*, Phys. Rev. D **81**, 085037 (2010), doi:10.1103/PhysRevD.81.085037.

[28] R. Rattazzi, S. Rychkov and A. Vichi, *Bounds in 4D conformal field theories with global symmetry*, J. Phys. A: Math. Theor. **44**, 035402 (2010), doi:10.1088/1751-8113/44/3/035402.

[29] D. Poland, D. Simmons-Duffin and A. Vichi, *Carving out the space of 4D CFTs*, J. High Energy Phys. **05**, 110 (2012), doi:10.1007/JHEP05(2012)110.

[30] S. El-Showk, M. F. Paulos, D. Poland, S. Rychkov, D. Simmons-Duffin and A. Vichi, *Solving the 3D Ising model with the conformal bootstrap*, Phys. Rev. D **86**, 025022 (2012), doi:10.1103/PhysRevD.86.025022.

[31] S. El-Showk and M. F. Paulos, *Bootstrapping conformal field theories with the extremal functional method*, Phys. Rev. Lett. **111**, 241601 (2013), doi:10.1103/PhysRevLett.111.241601.

[32] M. Hogervorst and S. Rychkov, *Radial coordinates for conformal blocks*, Phys. Rev. D **87**, 106004 (2013), doi:10.1103/PhysRevD.87.106004.

[33] F. Kos, D. Poland and D. Simmons-Duffin, *Bootstrapping the O(N) vector models*, J. High Energy Phys. **06**, 091 (2014), doi:10.1007/JHEP06(2014)091.

[34] F. Kos, D. Poland and D. Simmons-Duffin, *Bootstrapping mixed correlators in the 3D Ising model*, J. High Energy Phys. **11**, 109 (2014), doi:10.1007/JHEP11(2014)109.

[35] D. Simmons-Duffin, *A semidefinite program solver for the conformal bootstrap*, J. High Energy Phys. **06**, 174 (2015), doi:10.1007/JHEP06(2015)174.
[36] F. Kos, D. Poland, D. Simmons-Duffin and A. Vichi, *Bootstrapping the O(N) archipelago*, J. High Energy Phys. **11**, 106 (2015), doi:10.1007/JHEP11(2015)106.

[37] W. Landry and D. Simmons-Duffin, *Scaling the semidefinite program solver SDPB*, (arXiv preprint) doi:10.48550/arXiv.1909.09745.

[38] R. S. Erramilli, L. V. Iliesiu, P. Kravchuk, W. Landry, D. Poland and D. Simmons-Duffin, *blocks_3d: Software for general 3d conformal blocks*, J. High Energy Phys. **11**, 006 (2021), doi:10.1007/JHEP11(2021)006.

[39] S. El-Showk, M. F. Paulos, D. Poland, S. Rychkov, D. Simmons-Duffin and A. Vichi, *Solving the 3d Ising model with the conformal bootstrap II. c-minimization and precise critical exponents*, J. Stat. Phys. **157**, 869 (2014), doi:10.1007/s10955-014-1042-7.

[40] F. Kos, D. Poland, D. Simmons-Duffin and A. Vichi, *Precision islands in the Ising and O(N) models*, J. High Energy Phys. **08**, 036 (2016), doi:10.1007/JHEP08(2016)036.

[41] S. M. Chester, W. Landry, J. Liu, D. Poland, D. Simmons-Duffin, N. Su and A. Vichi, *Carving out OPE space and precise O(2) model critical exponents*, J. High Energy Phys. **06**, 142 (2020), doi:10.1007/JHEP06(2020)142.

[42] S. M. Chester, W. Landry, J. Liu, D. Poland, D. Simmons-Duffin, N. Su and A. Vichi, *Bootstrapping Heisenberg magnets and their cubic instability*, Phys. Rev. D **104**, 105013 (2021), doi:10.1103/PhysRevD.104.105013.

[43] D. Poland and D. Simmons-Duffin, *The conformal bootstrap*, Nat. Phys. **12**, 535 (2016), doi:10.1038/nphys3761.

[44] D. Poland, S. Rychkov and A. Vichi, *The conformal bootstrap: Theory, numerical techniques, and applications*, Rev. Mod. Phys. **91**, 015002 (2019), doi:10.1103/RevModPhys.91.015002.

[45] D. Poland and D. Simmons-Duffin, *Snowmass white paper: The numerical conformal bootstrap*, (arXiv preprint) doi:10.48550/arXiv.2203.08117.

[46] S. Rychkov and N. Su, *New developments in the numerical conformal bootstrap*, Rev. Mod. Phys. **96**, 045004 (2024), doi:10.1103/RevModPhys.96.045004.

[47] D. Simmons-Duffin, *The conformal bootstrap*, in *New frontiers in fields and strings*, World Scientific, Singapore, ISBN 9789813149434 (2017), doi:10.1142/9789813149441_0001.

[48] S. M. Chester, *Weizmann lectures on the numerical conformal bootstrap*, Phys. Rep. **1045**, 1 (2023), doi:10.1016/j.physrep.2023.10.008.

[49] M. Reehorst, S. Rychkov, D. Simmons-Duffin, B. Sirois, N. Su and B. van Rees, *Navigator function for the conformal bootstrap*, SciPost Phys. **11**, 072 (2021), doi:10.21468/SciPostPhys.11.3.072.

[50] A. Liu, D. Simmons-Duffin, N. Su and B. C. van Rees, *Skydiving to bootstrap islands*, (arXiv preprint) doi:10.48550/arXiv.2307.13046.

[51] S. Rychkov and Z. M. Tan, *The ε-expansion from conformal field theory*, J. Phys. A: Math. Theor. **48**, 29FT01 (2015), doi:10.1088/1751-8113/48/29/29FT01.
[52] J. Maldacena and A. Zhiboedov, *Constraining conformal field theories with a higher spin symmetry*, J. Phys. A: Math. Theor. **46**, 214011 (2013), doi:10.1088/1751-8113/46/21/214011.

[53] M. Hogervorst, S. Rychkov and B. C. van Rees, *Truncated conformal space approach in d dimensions: A cheap alternative to lattice field theory?*, Phys. Rev. D **91**, 025005 (2015), doi:10.1103/PhysRevD.91.025005.

[54] M. Hogervorst, S. Rychkov and B. C. van Rees, *Unitarity violation at the Wilson-Fisher fixed point in 4-ε dimensions*, Phys. Rev. D **93**, 125025 (2016), doi:10.1103/PhysRevD.93.125025.

[55] S. El-Showk, M. Paulos, D. Poland, S. Rychkov, D. Simmons-Duffin and A. Vichi, *Conformal field theories in fractional dimensions*, Phys. Rev. Lett. **112**, 141601 (2014), doi:10.1103/PhysRevLett.112.141601.

[56] A. Cappelli, L. Maffi and S. Okuda, *Critical Ising model in varying dimension by conformal bootstrap*, J. High Energy Phys. **01**, 161 (2019), doi:10.1007/JHEP01(2019)161.

[57] C. Bonanno, A. Cappelli, M. Kompaniets, S. Okuda and K. J. Wiese, *Benchmarking the Ising universality class in 3 ≤ d < 4 dimensions*, SciPost Phys. **14**, 135 (2023), doi:10.21468/SciPostPhys.14.5.135.

[58] J. Henriksson, S. R. Kousvos and M. Reehorst, *Spectrum continuity and level repulsion: The Ising CFT from infinitesimal to finite ε*, J. High Energy Phys. **02**, 218 (2023), doi:10.1007/JHEP02(2023)218.

[59] J. Golden and M. F. Paulos, *No unitary bootstrap for the fractal Ising model*, J. High Energy Phys. **03**, 167 (2015), doi:10.1007/JHEP03(2015)167.

[60] B. Sirois, *Navigating through the O(N) archipelago*, SciPost Phys. **13**, 081 (2022), doi:10.21468/SciPostPhys.13.4.081.

[61] V. Gorbenko, S. Rychkov and B. Zan, *Walking, weak first-order transitions, and complex CFTs*, J. High Energy Phys. **10**, 108 (2018), doi:10.1007/JHEP10(2018)108.

[62] V. Gorbenko, S. Rychkov and B. Zan, *Walking, weak first-order transitions, and complex CFTs II. Two-dimensional Potts model at Q > 4*, SciPost Phys. **5**, 050 (2018), doi:10.21468/SciPostPhys.5.5.050.

[63] M. Reehorst, *Rigorous bounds on irrelevant operators in the 3d Ising model CFT*, J. High Energy Phys. **09**, 177 (2022), doi:10.1007/JHEP09(2022)177.

[64] N. Su, *The hybrid bootstrap*, (arXiv preprint) doi:10.48550/arXiv.2202.07607.

[65] S. M. Chester and N. Su, *Upper critical dimension of the 3-state Potts model*, (arXiv preprint) doi:10.48550/arXiv.2210.09091.

[66] P. Calabrese, A. Pelissetto and E. Vicari, *Multicritical behavior in frustrated spin systems with non-collinear order*, Nucl. Phys. B **709**, 550 (2005), doi:10.1016/j.nuclphysb.2004.12.006.

[67] V. Gorbenko, S. Rychkov and B. Zan, *Walking, weak first-order transitions, and complex CFTs*, J. High Energy Phys. **10**, 108 (2018), doi:10.1007/JHEP10(2018)108.

[68] J. L. Cardy, M. Nauenberg and D. J. Scalapino, *Scaling theory of the Potts-model multicritical point*, Phys. Rev. B **22**, 2560 (1980), doi:10.1103/PhysRevB.22.2560.
[69] S. Gukov, RG flows and bifurcations, Nucl. Phys. B **919**, 583 (2017), doi:10.1016/j.nuclphysb.2017.03.025.

[70] D. B. Kaplan, J.-W. Lee, D. T. Son and M. A. Stephanov, Conformality lost, Phys. Rev. D **80**, 125005 (2009), doi:10.1103/PhysRevD.80.125005.

[71] A. Pelissetto and E. Vicari, Relevance of the axial anomaly at the finite-temperature chiral transition in QCD, Phys. Rev. D **88**, 105018 (2013), doi:10.1103/PhysRevD.88.105018.

[72] S. R. Kousvos and A. Stergiou, CFTs with $U(m) \times U(n)$ global symmetry in 3D and the chiral phase transition of QCD, SciPost Phys. **15**, 075 (2023), doi:10.21468/SciPostPhys.15.2.075.

[73] Y. Nakayama and T. Ohtsuki, Approaching the conformal window of $O(n) \times O(m)$ symmetric Landau-Ginzburg models using the conformal bootstrap, Phys. Rev. D **89**, 126009 (2014), doi:10.1103/PhysRevD.89.126009.

[74] Z. Li and N. Su, Bootstrapping mixed correlators in the five dimensional critical $O(N)$ models, J. High Energy Phys. **04**, 098 (2017), doi:10.1007/JHEP04(2017)098.

[75] N. Su, Simpleboot: A Mathematica framework for bootstrap calculations (2019), https://gitlab.com/bootstrapcollaboration/simpleboot.

[76] P. Bak and D. Mukamel, Physical realizations of $n \ge 4$-component vector models. III. Phase transitions in Cr, Eu, MnS$_2$, Ho, Dy, and Tb, Phys. Rev. B **13**, 5086 (1976), doi:10.1103/PhysRevB.13.5086.

[77] Z. Barak and M. B. Walker, First-order phase transitions in Tb, Dy, and Ho, Phys. Rev. B **25**, 1969 (1982), doi:10.1103/PhysRevB.25.1969.

[78] H. Kawamura and S. Miyashita, Phase transition of the two-dimensional Heisenberg antiferromagnet on the triangular lattice, J. Phys. Soc. Jpn. **53**, 9 (1984), doi:10.1143/JPSJ.53.9.

[79] H. Kawamura, Phase transition of the three-dimensional Heisenberg antiferromagnet on the layered-triangular lattice, J. Phys. Soc. Jpn. **54**, 3220 (1985), doi:10.1143/JPSJ.54.3220.

[80] H. Kawamura, Phase transition of the three-dimensional XY antiferromagnet on the layered-triangular lattice, J. Phys. Soc. Jpn. **55**, 2095 (1986), doi:10.1143/JPSJ.55.2095.

[81] H. Kawamura, Universality of phase transitions of frustrated antiferromagnets, J. Phys.: Condens. Matter **10**, 4707 (1998), doi:10.1088/0953-8984/10/22/004.

[82] H. Kawamura, Renormalization-group approach to the frustrated Heisenberg antiferromagnet on the layered-triangular lattice, J. Phys. Soc. Jpn. **55**, 2157 (1986), doi:10.1143/JPSJ.55.2157.

[83] B. Delamotte, Y. Holovatch, D. Ivaneyko, D. Mouhanna and M. Tissier, Fixed points in frustrated magnets revisited, J. Stat. Mech.: Theory Exp. P03014 (2008), doi:10.1088/1742-5468/2008/03/P03014.

[84] J. Guckenheimer and Y. Kuznetsov, Bogdanov-Takens bifurcation, Scholarpedia **2**, 1854 (2007), doi:10.4249/scholarpedia.1854.
[85] V. I. Arnold, *Geometrical methods in the theory of ordinary differential equations*, Springer, New York, USA, ISBN 9780387966496 (1988), doi:10.1007/978-1-4612-1037-5.

[86] C. B. Jepsen and F. K. Popov, *Homoclinic renormalization group flows, or when relevant operators become irrelevant*, Phys. Rev. Lett. **127**, 141602 (2021), doi:10.1103/PhysRevLett.127.141602.

[87] S. Dutta, *A simple circuit model showing feature-rich Bogdanov-Takens bifurcation* (2011), https://ewh.ieee.org/sb/calcutta/jadavpur/Papers/Circuit_Model_for_Bogdanov-Takens_Bifurcation.pdf.

[88] I. Jack and H. Osborn, *Analogs of the c-theorem for four-dimensional renormalisable field theories*, Nucl. Phys. B **343**, 647 (1990), doi:10.1016/0550-3213(90)90584-Z.

[89] I. Jack and C. Poole, *Scheme invariants in $\phi^4$ theory in four dimensions*, Phys. Rev. D **98**, 065011 (2018), doi:10.1103/PhysRevD.98.065011.

[90] C. Behan, *Conformal manifolds: ODEs from OPEs*, J. High Energy Phys. **03**, 127 (2018), doi:10.1007/JHEP03(2018)127.

[91] G. P. Korchemsky, *On level crossing in conformal field theories*, J. High Energy Phys. **03**, 212 (2016), doi:10.1007/JHEP03(2016)212.

[92] B. Delamotte, Y. Holovatch, D. Ivaneyko, D. Mouhanna and M. Tissier, *Reply to: “Comment on “Spurious fixed points in frustrated magnets”, cond-mat/0609285*, (arXiv preprint) doi:10.48550/arXiv.cond-mat/0610613.

[93] A. Pelissetto and E. Vicari, *Comment on “Spurious fixed points in frustrated magnets”, cond-mat/0609285*, (arXiv preprint) doi:10.48550/arXiv.cond-mat/0610113.

[94] M. Serone, *Private communication*.

[95] J. Polchinski, *Scale and conformal invariance in quantum field theory*, Nucl. Phys. B **303**, 226 (1988), doi:10.1016/0550-3213(88)90179-4.

[96] S. Rychkov, *EPFL lectures on conformal field theory in $D \ge 3$ dimensions*, Springer, Cham, Switzerland, ISBN 9783319436258 (2017), doi:10.1007/978-3-319-43626-5.

[97] S. El-Showk, Y. Nakayama and S. Rychkov, *What Maxwell theory in $d \neq 4$ teaches us about scale and conformal invariance*, Nucl. Phys. B **848**, 578 (2011), doi:10.1016/j.nuclphysb.2011.03.008.

[98] Y. Nakayama, *Interacting scale invariant but nonconformal field theories*, Phys. Rev. D **95**, 065016 (2017), doi:10.1103/PhysRevD.95.065016.

[99] A. Mauri and M. I. Katsnelson, *Scale without conformal invariance in membrane theory*, Nucl. Phys. B **969**, 115482 (2021), doi:10.1016/j.nuclphysb.2021.115482.

[100] A. Gimenez-Grau, Y. Nakayama, and S. Rychkov, *Scale without conformal invariance in dipolar ferromagnets, in preparation*.

[101] Y. Nakayama, *Conformal invariance from scale invariance in nonlinear sigma models*, Phys. Rev. D **102**, 065018 (2020), doi:10.1103/PhysRevD.102.065018.

[102] A. Pelissetto, P. Rossi and E. Vicari, *Large-n critical behavior of $O(n) \times O(m)$ spin models*, Nucl. Phys. B **607**, 605 (2001), doi:10.1016/S0550-3213(01)00223-1.
[103] J. A. Gracey, *Chiral exponents in $O(n) \times O(m)$ spin models at $O(1/N^2)$*, Phys. Rev. B **66**, 134402 (2002), doi:10.1103/PhysRevB.66.134402.

[104] J. A. Gracey, *Critical exponent $\omega$ at $O(1/N)$ in $O(N) \times O(m)$ spin models*, Nucl. Phys. B **644**, 433 (2002), doi:10.1016/S0550-3213(02)00818-0.

[105] S.-k. Ma, *Scaling variables and dimensions*, Phys. Rev. A **10**, 1818 (1974), doi:10.1103/PhysRevA.10.1818.

[106] A. N. Vasil'ev and A. S. Stepanenko, *A method of calculating the critical dimensions of composite operators in the massless nonlinear $\sigma$ model*, Theor. Math. Phys. **95**, 471 (1993), doi:10.1007/BF01015903.

[107] P. Calabrese and P. Parruccini, *Five-loop $\epsilon$ expansion for $O(n) \times O(m)$ spin models*, Nucl. Phys. B **679**, 568 (2004), doi:10.1016/j.nuclphysb.2003.12.002.

[108] B. Delamotte, D. Mouhanna and M. Tissier, *Nonperturbative renormalization-group approach to frustrated magnets*, Phys. Rev. B **69**, 134413 (2004), doi:10.1103/PhysRevB.69.134413.

[109] A. N. Vasil'ev, Y. M. Pis'mak and Y. R. Khonkonen, *$1/n$ expansion: Calculation of the exponent $\eta$ in the order $1/n^3$ by the conformal bootstrap method*, Theor. Math. Phys. **50**, 127 (1982), doi:10.1007/BF01015292.

[110] A. N. Vasil'ev, Y. M. Pis'mak and Y. R. Khonkonen, *$1/N$ expansion: Calculation of the exponents $\eta$ and $\nu$ in the order $1/N^2$ for arbitrary number of dimensions*, Theor. Math. Phys. **47**, 465 (1981), doi:10.1007/BF01019296.

[111] J. A. Gracey, *Crossover exponent in $O(N) \phi^4$ theory at $O(1/N^2)$*, Phys. Rev. E **66**, 027102 (2002), doi:10.1103/PhysRevE.66.027102.

[112] J. Henriksson, *The critical $O(N)$ CFT: Methods and conformal data*, Phys. Rep. **1002**, 1 (2023), doi:10.1016/j.physrep.2022.12.002.

[113] C. Behan, P. Kravchuk, W. Landry, and D. Simmons-Duffin, *scalar_blocks* (2019), https://gitlab.com/bootstrapcollaboration/scalar_blocks.