# Introducing sgboost: A Practical Guide and Implementation of sparse-group
boosting in R

Fabian Obster a b and Christian Heumannb CONTACT Fabian Obster Email:
<EMAIL_ADDRESS>aDepartment of Business Administration, University of
the Bundeswehr Munich, Werner-Heisenberg-Weg 39, Bavaria, Germany bDepartment
of Statistics, LMU Munich, Ludwigstr. 33, Bavaria, Germany

###### Abstract

This paper introduces the sgboost package in R, which implements sparse-group
boosting for modeling high-dimensional data with natural groupings in
covariates. Sparse-group boosting offers a flexible approach for both group
and individual variable selection, reducing overfitting and enhancing model
interpretability. The package uses regularization techniques based on the
degrees of freedom of individual and group base-learners, and is designed to
be used in conjunction with the mboost package. Through comparisons with
existing methods and demonstration of its unique functionalities, this paper
provides a practical guide on utilizing sparse-group boosting in R,
accompanied by code examples to facilitate its application in various research
domains. Overall, this paper serves as a valuable resource for researchers and
practitioners seeking to use sparse-group boosting for efficient and
interpretable high-dimensional data analysis.

###### keywords:

group selection; variable selection; R package; Ridge regression; within-group
sparsity

††articletype: Application Note

##  1 Introduction

Regularized regression is used to model high-dimensional data to reduce the
risk of overfitting and to perform variable selection. In many cases,
covariates have natural groupings, such as with gene data or categorical data
often found in survey data. In such cases, one may want to select whole groups
of variables or just individual parts (univariate regression model). ‘sgboost‘
implements the sparse-group boosting in R and other useful functions for
sparse group model interpretation unique to boosting and visualization for
group and individual variable selection. The package is available on CRAN [8].
A formula object defining group base-learners and individual base-learners is
utilized for model fitting. Regularization is based on the degrees of freedom
$(df(\cdot))$ of individual and group base-learners as defined in [10]. For an
individual base-learner, the Regularization parameter $\lambda$, and for a
group base-learners the regularization parameter $\lambda^{(g)}$ is defined,
such that $df(\lambda)=\alpha$ and $df(\lambda^{(g)})=1-\alpha$. The package
should be used with the ’mboost’ package [4]. Sparse-group boosting is an
alternative method to sparse-group lasso, employing boosted Ridge regression
[14]. Although there are many methods of variable selection, most focus on
group selection only, e.g. [7], [5] and [16], or individual variable selection
e.g. [1], [15] and [2]. However, it should be noted, that with group variable
selection with overlapping groups, one could also end up with sparse-group
variable selection. There are not many R packages implementing sparse-group
variable selection. There is ’SGL’ [13] implementing [14] and ’sparsegl’ [6]
with a faster implementation of the sparse-group lasso.
The goal of this paper is to provide a practical guide including the code on
how to use sparse-group boosting in R and get the most out of the method. The
code is presented within this manuscript and can be found also on GitHub
together with the used dataset (https://github.com/FabianObster/sgboost-
introduction).

##  2 Sparse-group boosting workflow with sample data

We first simulate the sample data and corresponding group structure with 40
equal-sized groups. Based on a linear regression model we simulate the
response variable y as part of the data.frame with $n=100$ observations and
$p=200$ predictor variables (each group is formed by 5 predictors).

⬇

beta <- c(

rep(5, 5), c(5, -5, 2, 0, 0), rep(-5, 5),

c(2, -3, 8, 0, 0), rep(0, (200 - 20))

)

X <- matrix(data = rnorm(20000, mean = 0, sd = 1), 100, 200)

df <- data.frame(X) %>%

mutate(y = X %*% beta+rnorm(100, mean = 0, sd = 1)) %>%

mutate_all(function(x){as.numeric(scale(x))})

group_df <- data.frame(

group_name = rep(1:40, each = 5),

variable_name = colnames(df)[1:200]

)

###  2.1 Defining the model

Now we use the group structure to describe the sparse group boosting formula
with the function `create_formula()`. We only need the data.frame() describing
the group structure. It should contain two variables one indicating the name
of the variable in the modeling data (`var_name`), and one indicating the
group it belongs to (`group_name`). Additionally, we need to pass the mixing
parameter `alpha` and the name of the outcome variable.

⬇

sgb_formula <- create_formula(

alpha = 0.4, group_df = group_df, outcome_name = "y",

group_name = "group_name", var_name = "variable_name")

This function returns an R-formula consisting of $p$ model terms defining the
individual base-learners and $G$ group base-learners.

⬇

labels(terms(sgb_formula))[[1]]

## bols(X1, df = 0.4, intercept = FALSE)

labels(terms(sgb_formula))[[201]]

## bols(X1, X2, X3, X4, X5, df = 0.6, intercept = FALSE)

###  2.2 Fitting and tuning the model

`sgboost` is to be used in conjunction with the `mboost` package which
provides many useful functions and methods which can also be used for sparse-
group boosting models.
Now we pass the formula to `mboost()` and use the arguments as seems
appropriate. The main hyperparameters are `nu` and `mstop`. For model tuning,
the function `cvrisk` can be used and plotted. Running the cross-
validation/bootstrap in parallel can speed up the process.

⬇

sgb_model <- mboost(

formula = sgb_formula, data = df,

control = boost_control(nu = 1, mstop = 600)

)

cv_sgb_model <- cvrisk(sgb_model)

mstop(cv_sgb_model)

## 204

plot(cv_sgb_model

Figure 1: Out of sample error depending on the boosting iteration

In this example, the lowest out-of-sample risk is obtained at 204 boosting
iterations, so we only use the first 204 updates for the final model.

###  2.3 Interpreting and plotting a sparse-group boosting model

`sgboost` has useful functions to understand sparse-group boosting models and
reflects that final model estimates of a specific variable in the dataset can
be attributed to group base-learners as well as individual base-learners
depending on the boosting iteration.

####  2.3.1 Variable importance

A good starting point for understanding a sparse-group boosting model is the
variable importance. In the context of boosting, the variable importance can
be defined as the relative contribution of each predictor to the overall
reduction of the loss function (negative log-likelihood). `get_varimp()`
returns the variable importance of each base-learner/predictor selected
throughout the boosting process. In the case of the selection of an individual
variable - call it $x_{1}$ \- as well as the group it belongs to
-$x_{1},x_{2},...x_{p}$ -, both base-learners (predictors) will have an
associated variable importance as defined before. This allows us to
differentiate between the individual contribution of $x_{1}$ as its own
variable and the contribution of the group $x_{1}$ belongs to. It is
impossible to compute the aggregated variable importance of $x_{1}$ as it is
unclear how much $x_{1}$ contributes to the group. However, the aggregated
coefficients can be computed using `get_coef()`, which also returns the
aggregated importance of all groups vs. all individual variables in a separate
data.frame. With `plot_varimp()` one can visualize the variable importance as
a barplot. Since group sizes can be large, the function allows to cut of the
name of a predictor after `max_char_length` characters. One can indicate the
maximum number of predictors to be printed through `n_predictors` or through
the minimal variable importance a predictor has to have through `prop`.
Through both parameters, the number of printed entries can be reduced. Note,
that in this case, the relative importance of groups in the legend is based
only on the plotted variables and not the ones removed. Adding information
about the direction of effect sizes, one could add arrows behind the bars [9].
For groups, one can use the aggregated coefficients from `get_coef()`.

⬇

slice(get_varimp(sgb_model =sgb_model_linear)$varimp,1:5)

    # A tibble: 5  6
      reduction blearner        predictor   selfreq type  relative_
                                                          importance
          <dbl> <chr>           <chr>       <dbl>   <chr>     <dbl>
    1    0.297  bols(X1, X2,... X1, X2, ... 0.206   group     0.301
    2    0.288  bols(X18, in... X18         0.0196  indi...   0.292
    3    0.230  bols(X11, X1... X11, X12... 0.25    group     0.233
    4    0.0414 bols(X7, int... X7          0.0784  indi...   0.0419
    5    0.0392 bols(X6, int... X6          0.0833  indi...   0.0397

⬇

get_varimp(sgb_model = sgb_model_linear)$group_importance

        # A tibble: 2  2
      type       importance
      <chr>           <dbl>
    1 group           0.534
    2 individual      0.466

⬇

plot_varimp(sgb_model = sgb_model_linear, n_predictors = 15)

Figure 2: Variable importance of the sparse-group boosting model for simulated
data

In this example, we see that both individual variables and groups were
selected and contributed to the reduction of the loss function. The most
important predictor is the first group, followed by variable 18, and then
followed by group three. This is in line with what was simulated, as variable
18 has the biggest beta value, and groups one and three are full groups,
meaning all variables within the groups have a non-zero beta coefficient.
Groups two and four have within-group sparsity, therefore they were selected
as individual variables rather than groups.

####  2.3.2 Model coefficients

The resulting coefficients can be retrieved through `get_coef()` In sparse-
group boosting models, a variable in a dataset can be selected as an
individual variable or through a group. Therefore, there can be two associated
effect sizes for the same variable. This function aggregates both and returns
them in a data.frame sorted by the effect size `’effect’`.

⬇

slice(get_coef(sgb_model = sgb_model)$raw,1:5)

    # A tibble: 5  5
      variable effect blearner         predictor               type
      <chr>     <dbl> <chr>            <chr>                   <chr>
    1 X18       0.364 bols(X18, int... X18                     individual
    2 X5        0.250 bols(X1, X2, ... X1, X2, X3, X4, X5      group
    3 X15      -0.249 bols(X11, X12... X11, X12, X13, X14, X15 group
    4 X4        0.234 bols(X1, X2, ... X1, X2, X3, X4, X5      group
    5 X11      -0.228 bols(X11, X12... X11, X12, X13, X14, X15 group

⬇

slice(get_coef(sgb_model = sgb_model)$aggregate,1:5)

    # A tibble: 5  4
     variable effect learner               predictor
      <chr>    <dbl> <chr>                 <chr>
    1 X18      0.364 bols(X18, inte...     X18
    2 X15     -0.272 bols(X11, X12,...;    X11, X12, X13, X14, X15; X15
                     bols(X15, inte...
    3 X5       0.250 bols(X1, X2,...       X1, X2, X3, X4, X5
    4 X4       0.234 bols(X1, X2,...       X1, X2, X3, X4, X5
    5 X13     -0.230 bols(X11, X12,...;    X11, X12, X13, X14, X15; X13
                     bols(X13, inte...

We see that the effect sizes differ between the two perspectives. The variable
X15 for example has a more extreme model coefficient of -0.272 in the
aggregated case compared to the coefficient of -0.249 derived only from the
group base-learner. Consequently, also the ordering differs. X11 has a greater
absolute model coefficient from the group than X13, but in the aggregated
version the absolute model coefficient of X13 exceeds the one of X11.

####  2.3.3 Plotting model coefficients and importance

With `plot_effects()` we can plot the effect sizes of the sparse-group
boosting model in relation to the relative importance to get an overall
picture of the model. Through the parameter `’plot_type’` one can choose the
type of visualization. `’radar’` refers to a radar plot using polar
coordinates. Here the angle is relative to the cumulative relative importance
of predictors and the radius is proportional to the effect size. `’clock’`
does the same as `’radar’` but uses clock coordinates instead of polar
coordinates. `’scatter’` uses the effect size as the y-coordinate and the
cumulative relative importance as the x-axis in a classical Scatter plot.

⬇

plot_effects(sgb_model = sgb_model, n_predictors = 5,

base_size = 10)

plot_effects(sgb_model = sgb_model, n_predictors = 5,

plot_type = "clock", base_size = 10)

plot_effects(sgb_model = sgb_model, n_predictors = 5,

plot_type = "scatter", base_size = 10)

Figure 3: Effect size and

####  2.3.4 Coefficient path

`plot_path` calls `get_coef_path()` to retrieve the aggregated coefficients
from a `mboost` object for each boosting iteration and plots it, indicating if
a coefficient was updated by an individual variable or group.

⬇

plot_path(sgb_model = sgb_model)

Figure 4: Coefficient path of a sparse-group boosting model with simulated
data

In the coefficient path shown in Figure 4, we see the change in model
coefficients. Since the path shows the aggregated model coefficients, the path
of one variable in the dataset may have both colors. This is the case with
variable X1 which was first updated through the group and then also as an
individual variable or with variable X15 in reverse order.

##  3 Real data

In this section, we will fit a sparse-group boosting model with sgboost to a
real dataset. We will use behavioral ecological data and an associated group
structure [11] to explain whether farmers in Chile and Tunisia are planning to
take adaptive measures against climate change in the following years. We will
use a logistic regression model for this binary decision. The data consists of
14 groups and 84 variables for the 801 farmers. Groups include vulnerability
to climate change [12], social, biophysical, and economic assets, as well as
perceptions of the farmers. After loading the data and group structure, we
create the formula with mixing parameter $\alpha=0.3$. Then, we pass the
formula to mboost() with 1000 boosting iterations and a learning rate of 0.3.

⬇

model_df <- readRDS(’model_df.RDS’) %>%

mutate_at(index_df$col_names, factor)

index_df <- readRDS(’index_df.RDS’)

sgb_formula <- create_formula(

group_df = index_df, var_name = ’col_names’,

group_name = ’index’, outcome_name = ’S5.4’

)

model <- mboost(

sgb_formula, data = model_df,

family = Binomial(link = ’logit’),

control = boost_control(mstop = 1000, nu = 0.3)

)

cv_model <- cvrisk(model)

model <- model[mstop(cv_model)]

The model is stopped early after 466 Boosting iterations. We examine the
coefficient path and see that in the early stage, individual base-learners
were dominantly selected like the variable ’S1.8b or ’S8.11 river’ which
indicates whether river irrigation is used. Many of the variables were first
included as individual variables and later also through group base-learners
like ’S8.1b’ or ’S2.5c proximity’ (Proximity to extreme weather events), which
we also saw in the simulated data.

⬇

plot_path(model)

Figure 5: Coefficient path using the ecological dataset

In figure 6, we look at the variable importance with the default values,
plotting all 27 selected predictors of which 8 are groups, the latter having a
relative variable importance of 22 percent. The most important base-learner is
the individual variable ’S1.8b’ indicating whether farming journals are being
used and the most important group is the social asset group followed by the
group consisting of the four considered regions.

⬇

plot_varimp(model)

Figure 6: Variable importance using the ecological dataset

Plotting the effect sizes of all predictors having a relative importance of
greater than 1.5 percent shows the tendency for more important variables to
have greater absolute effect sizes. For readability, we set the number of
printed characters per variable to 6 and use the ’scatter’ version of the
plot.

⬇

plot_effects(

model, plot_type = ’scatter’,

prop = 0.015, max_char_length = 6

)

Figure 7: Coefficient plot using the ecological dataset

##  4 The sparse-group boosting algorithm

General functional gradient descent Algorithm ([3])

  1. 1.

Define base-learners of the structure $h:\mathbb{R}^{n\times p}\to\mathbb{R}$

  2. 2.

Initialize $m=0$ and $\widehat{f}^{[0]}\equiv 0$ or
$\widehat{f}^{[0]}\equiv\overline{y}$

  3. 3.

Set $m=m+1$ and compute the negative gradient $\frac{\partial}{\partial
f}l(y,f)$ and evaluate it at $\widehat{f}^{[m-1]}$. Doing this yields the
pseudo-residuals $u_{1},...,u_{n}$ with

$u_{i}^{[m]}=\frac{\partial}{\partial f}l(y_{i},f)|_{f=\widehat{f}^{[m-1]}},$

for all $i=1,...,n$

  4. 4.

Fit the base-learner $h$ with the response $(u_{1}^{[m]},...,u_{n}^{[m]})$ to
the data. This yields $\widehat{h}^{[m]}$, which is an approximation of the
negative gradient

  5. 5.

Update

$\widehat{f}^{[m]}=\widehat{f}^{[m-1]}+\eta\cdot\widehat{h}^{[m]}$

here $\eta$ can be seen as learning rate with $\eta\in]0,1[$

  6. 6.

Repeat steps 3,4 and 5 until $m=M$

‘sgboost‘ uses the following adjustment to the general functional gradient
descent algorithm. The sparse-group-boosting optimization problem is
formulated as the sum of mean squared error (linear regression) or negative
log-likelihood (generalized linear models e.g. logistic regression) within
each boosting iteration in steps 3 and 4. The base-learner with the smallest
evaluated loss function will be selected.
Linear regression:

  * •

Group base-learner:

$\operatorname*{arg\,min}_{g\leq
G}\min_{\beta^{(g)}\in\mathbb{R}^{p}}\left(\frac{1}{2n}\rVert
u^{[m]}-X^{(g)}\beta^{(g)}\rVert_{2}^{2}+\lambda^{(g)}\rVert\beta^{(g)}\rVert_{2}^{2}\right).$

  * •

Individual base-learner:

$\operatorname*{arg\,min}_{g\leq G,j\leq
p_{g}}\min_{\beta^{(g)}_{j}\in\mathbb{R}}\left(\frac{1}{2n}\rVert
u^{[m]}-X^{(g)}_{j}\beta^{(g)}_{j}\rVert_{2}^{2}+\lambda^{(g)}_{j}\rVert\beta^{(g)}_{j}\rVert_{2}^{2}\right).$

If the group base-learner with the smallest loss as defined above has a
smaller loss than the individual base-learner with the smallest loss then the
group base-learner will be selected in the update and vice versa.
Generalized linear models:

  * •

Group base-learners:

$\operatorname*{arg\,min}_{g\in
G}\min_{\beta^{(g)}\in\mathbb{R}}\left(-\text{loglik}(u^{[m]},X^{(g)},\beta^{(g)})+\lambda^{(g)}\rVert\beta^{(g)}\rVert_{2}^{2}\right).$

  * •

Individual base-learners:

$\operatorname*{arg\,min}_{g\in G,j\in
p_{g}}\min_{\beta^{(g)}_{j}\in\mathbb{R}}\left(-\text{loglik}(u^{[m]},X^{(g)}_{j},\beta^{(g)}_{j})+\lambda^{(g)}_{j}\rVert\beta^{(g)}_{j}\rVert_{2}^{2}\right).$

where

  * •

$df(\lambda^{(g)}_{j})=\alpha$, $df(\lambda^{(g)})=1-\alpha$

  * •

$df(\lambda)=tr(2H_{\lambda}-H_{\lambda}^{T}H_{\lambda})$ and $H_{\lambda}$ is
the Ridge Hat matrix of a base-learner

  * •

$X^{(g)}$ is the submatrix of $X$ with columns corresponding to group $g$ and
$X^{(g)}_{j}$ the $j$-th column of group $g$.

  * •

$\beta^{(g)}$ and $\beta^{(g)}_{j}$ are the corresponding parameter vectors of
the features in group $g$.

  * •

$\alpha$ adjusts the weight between group and individual base-learners.

  * •

$\lambda^{(g)},\lambda^{(g)}_{j}$ Ridge penalty parameters.

  * •

$u^{[m]}$ is the current negative gradient vector from the previous boosting
iteration.

## Funding

This research is funded by dtec.bw – Digitalization and Technology Research
Center of the Bundeswehr. dtec.bw is funded by the European Union –
NextGenerationEU. All statements expressed in this article are the authors’
and do not reflect the official opinions or policies of the authors’ host
affiliations or any of the supporting institutions.

## References

  * [1] K.N. Berk, _Forward and backward stepping in variable selection_ , Journal of Statistical Computation and Simulation 10 (1980), pp. 177–185. Available at https://doi.org/10.1080/00949658008810367, Publisher: Taylor & Francis _eprint: https://doi.org/10.1080/00949658008810367.
  * [2] P. Bühlmann and T. Hothorn, _Boosting Algorithms: Regularization, Prediction and Model Fitting_ , Statistical Science 22 (2007), pp. 477–505. Available at https://projecteuclid.org/journals/statistical-science/volume-22/issue-4/Boosting-Algorithms-Regularization-Prediction-and-Model-Fitting/10.1214/07-STS242.full, Publisher: Institute of Mathematical Statistics.
  * [3] J.H. Friedman, _Greedy function approximation: A gradient boosting machine._ , The Annals of Statistics 29 (2001), pp. 1189–1232. Available at https://projecteuclid.org/journals/annals-of-statistics/volume-29/issue-5/Greedy-function-approximation-A-gradient-boosting-machine/10.1214/aos/1013203451.full.
  * [4] T. Hothorn, P. Buehlmann, T. Kneib, M. Schmid, B. Hofner, F. Otto-Sobotka, F. Scheipl, and A. Mayr, _mboost: Model-Based Boosting_ (2023). Available at https://CRAN.Rproject.org/package=mboost.
  * [5] J. Huang, S. Ma, H. Xie, and C.H. Zhang, _A group bridge approach for variable selection_ , Biometrika 96 (2009), pp. 339–355. Available at https://www.jstor.org/stable/27798828, Publisher: [Oxford University Press, Biometrika Trust].
  * [6] X. Liang, A. Cohen, A.S. Heinsfeld, F. Pestilli, and D.J. McDonald, _sparsegl: An R Package for Estimating Sparse Group Lasso_ (2023). Available at http://arxiv.org/abs/2208.02942, arXiv:2208.02942 [stat].
  * [7] L. Meier, S. Van De Geer, and P. Bühlmann, _The Group Lasso for Logistic Regression_ , Journal of the Royal Statistical Society Series B: Statistical Methodology 70 (2008), pp. 53–71. Available at https://doi.org/10.1111/j.1467-9868.2007.00627.x.
  * [8] F. Obster, _sgboost: Sparse-Group Boosting_ (2024). Available at https://CRAN.R-project.org/package=sgboost.
  * [9] F. Obster, H. Bohle, and P.M. Pechan, _The financial well-being of fruit farmers in Chile and Tunisia depends more on social and geographical factors than on climate change_ , Communications Earth & Environment 5 (2024), pp. 1–12. Available at https://www.nature.com/articles/s43247-023-01128-2, Number: 1 Publisher: Nature Publishing Group.
  * [10] F. Obster and C. Heumann, _Sparse-group boosting – Unbiased group and variable selection_ , Tech. Rep. arXiv:2206.06344, arXiv, 2022. Available at http://arxiv.org/abs/2206.06344, arXiv:2206.06344 [stat] type: article.
  * [11] F. Obster, C. Heumann, H. Bohle, and P. Pechan, _Using interpretable boosting algorithms for modeling environmental and agricultural data_ , Scientific Reports 13 (2023), p. 12767. Available at https://www.nature.com/articles/s41598-023-39918-5, Number: 1 Publisher: Nature Publishing Group.
  * [12] P.M. Pechan, H. Bohle, and F. Obster, _Reducing vulnerability of fruit orchards to climate change_ , Agricultural Systems 210 (2023), p. 103713. Available at https://www.sciencedirect.com/science/article/pii/S0308521X2300118X.
  * [13] N. Simon, J. Friedman, T. Hastie, and a.R. Tibshirani, _SGL: Fit a GLM (or Cox Model) with a Combination of Lasso and Group Lasso Regularization_ (2019). Available at https://CRAN.R-project.org/package=SGL.
  * [14] N. Simon, J. Friedman, T. Hastie, and R. Tibshirani, _A Sparse-Group Lasso_ , Journal of Computational and Graphical Statistics 22 (2013), pp. 231–245. Available at http://www.tandfonline.com/doi/abs/10.1080/10618600.2012.681250.
  * [15] Z. Zhang, _Variable selection with stepwise and best subset approaches_ , Annals of Translational Medicine 4 (2016), p. 136. Available at https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4842399/.
  * [16] N. Zhou and J. Zhu, _Group Variable Selection via a Hierarchical Lasso and Its Oracle Property_ (2010). Available at http://arxiv.org/abs/1006.2871, arXiv:1006.2871 [stat].


# Sub-leading conformal dimensions at the O(4) Wilson-Fisher fixed point

Debasish Banerjee  Saha Institute of Nuclear Physics, Bidhan Nagar, Kolkata,
West Bengal 700064, India  Shailesh Chandrasekharan  Department of Physics,
Box 90305, Duke University, Durham, NC 27708, USA

###### Abstract

In this work we focus on computing the conformal dimensions $D(j_{L},j_{R})$
of local fields that transform in an irreducible representation of
$SU(2)\times SU(2)$ labeled with $(j_{L},j_{R})$ at the $O(4)$ Wilson-Fisher
fixed point using the Monte Carlo method. In the large charge expansion, among
the sectors with a fixed large value of $j={\rm max}(j_{L},j_{R})$, the
leading sector has $|j_{L}-j_{R}|=0$ and the sub-leading one has
$|j_{L}-j_{R}|=1$. Since Monte Carlo calculations at large $j$ become
challenging in the traditional lattice formulation of the $O(4)$ model, a
qubit regularized $O(4)$ lattice model was used recently to compute $D(j,j)$.
Here we extend those calculations to the sub-leading sector. Our Monte Carlo
results up to $j=20$ fit well to the form
$D(j,j-1)-D(j)\sim\lambda_{1/2}/\sqrt{j}+\lambda_{1}/j+\lambda_{3/2}/j^{3/2}$,
consistent with recent predictions of the large charge expansion. Taking into
account systematic effects in our fitting procedures we estimate the two
leading coefficients to be $\lambda_{1/2}=2.08(5)$, $\lambda_{1}=2.2(3)$.

##  I Introduction

There has been a resurrection of interest in conformal field theories in
recent years, especially due to the success of the bootstrap approach in
certain problems Rychkov (2016); Simmons-Duffin (2017). It has also become
clear that conformal field theories simplify in sectors with either large spin
Komargodski and Zhiboedov (2013) or large global charge Hellerman  _et al._
(2015). Due to these developments the field has seen a renaissance with
several new results over the past few years Kaviraj  _et al._ (2015a, b);
Alday (2017); Gopakumar  _et al._ (2017); Dey  _et al._ (2018); Caron-Huot
(2017); Hellerman and Maeda (2017); Hellerman  _et al._ (2017); Jafferis  _et
al._ (2018); Alvarez-Gaume  _et al._ (2019); Orlando  _et al._ (2020). A
recent review of the above progress can be found in Ref.Poland  _et al._
(2019); Dondi  _et al._ (2021).

An interesting quantity in a conformal field theory is the conformal dimension
of local fields that transform according to some representation of the
symmetries of the theory. In this work we focus on global symmetries and study
CFTs that emerge in three-dimensional $O(N)$ models at the Wilson-Fisher fixed
point. Recent work Hellerman  _et al._ (2015) showed that in these theories
the conformal dimension $D(Q)$ of local fields that transform under the
representation with charge $Q$ satisfy a large charge expansion of the form

$\displaystyle
D(Q)=\sqrt{\frac{Q^{3}}{4\pi}}\Bigg{(}c_{3/2}+\frac{4\pi}{Q}c_{1/2}+{\cal
O}\Big{(}1/Q^{2}\Big{)}\Bigg{)}+c_{0},$ (1)

where $c_{3/2}$ and $c_{1/2}$ are low-energy constants that need to be
determined non-perturbatively, while $c_{0}\approx-0.094$ can be determined
analytically Monin (2016); de la Fuente (2018). Non-perturbative Monte Carlo
calculations confirming these predictions for the case of the $O(2)$ model,
where the global charges $Q$ are represented by integers, have also been
performed, and it was discovered that in this case $c_{3/2}=1.195(10)$ and
$c_{1/2}=0.075(10)$ Banerjee  _et al._ (2018). These calculations were later
extended to the $O(4)$ model where the local fields transform in some
representation of the $O(4)$ symmetry. We can classify them according to the
irreducible representations of $SU_{L}(2)\times SU_{R}(2)$ labeled with
charges $(j_{L},j_{R})$ where $j_{L},j_{R}=0,1/2,1,3/2,...$. In this case it
is natural to define the charge as $Q=2j$ where $j={\rm max}(j_{L},j_{R})$. In
the large charge expansion, the leading sector is given by $|j_{L}-j_{R}|=0$
sector, for which Monte Carlo calculations give $c_{3/2}=1.068(4)$ and
$c_{1/2}=0.083(3)$ Banerjee  _et al._ (2019). It was shown that the large
charge expansion predicts the conformal dimensions, even at $Q=1$ within a few
percent.

How well does the large charge expansion do in predicting the conformal
dimensions in the sub-leading sector, where $|j_{L}-j_{R}|=1$? This is the
question that motivates our research in this work. A natural quantity to
measure in this case is ${\widetilde{\Delta}}(j)=D(j,j-1)-D(j,j)$, which has
the expansion of the form

$\displaystyle{\widetilde{\Delta}}(j)$
$\displaystyle=\lambda_{0}+\frac{\lambda_{1/2}}{j^{1/2}}+\frac{\lambda_{1}}{j}+\frac{\lambda_{3/2}}{j^{3/2}}+{\cal
O}\Big{(}\frac{1}{j^{2}}\Big{)},$ (2)

valid for large values of $j$. In this expansion, the coefficients of
fractional powers of $j$ like $\lambda_{1/2}$ and $\lambda_{3/2}$ are Wilson
coefficients that cannot in principle be determined within the large charge
effective field theory Gaumé  _et al._ (2021). Their origin is essentially
classical with quantum corrections. On the other hand the coefficients of
integer powers of $j$ like $\lambda_{0}$ and $\lambda_{1}$ arise from purely
quantum mechanical effects and can in principle be calculable analytically,
similar to $c_{0}$ introduced earlier. It is conjectured that $\lambda_{0}=0$
since it is difficult to imagine a calculation that would distinguish between
$D(j,j)$ and $D(j,j-1)$ in the large $j$ limit Err (2020). However, the spin
of the leading conformal field in the $(j,j)$ sector is different from the
$(j,j-1)$ sector. Such differences were first observed in Alvarez-Gaume  _et
al._ (2017) and later clarified in the context of $O(4)$ in Hellerman  _et
al._ (2019); Antipin  _et al._ (2020). For this reason $\lambda_{0}\neq 0$ may
still be possible. On the other hand $\lambda_{1}$ is most likely non-zero but
remains undetermined until now Com (2021).

In this work we design a new Monte Carlo method to compute
${\widetilde{\Delta}}(j)$ in order to explore the validity of Eq. 2. We
present results in the range $1\leq j\leq 20$. Within this range we study if
the conjecture that $\lambda_{0}=0$ is consistent with our results and
estimate the other three unknown coefficients.

##  II Qubit Regularized O(4) model

In order to construct a Monte Carlo method to compute
${\widetilde{\Delta}}(j)$ it is useful to understand how the $SU_{L}(2)\times
SU_{R}(2)$ symmetry is manifest in our lattice model, which is the same as the
one used in Ref. Banerjee  _et al._ (2019). For this purpose it is helpful to
view our model as strongly coupled lattice Quantum Electrodynamics (QED)
constructed with staggered fermions Cecile and Chandrasekharan (2008). When
gauge fields are integrated out exactly, the microscopic degrees of freedom
are made up of bosons with fermionic constituents. These bosons naturally have
a built-in hardcore interaction, and hence the Hilbert space on each lattice
site is finite dimensional. Such bosonic lattice field models with a finite
dimensional Hilbert space, that reproduce a continuum quantum field theory,
can be referred to as a qubit regularized model of the continuum quantum field
theory Singh and Chandrasekharan (2019). Other examples of qubit regularized
models for studying continuum quantum field theories with $O(N)$ symmetries,
have been constructed recently Singh (2019); Bhattacharya  _et al._ (2020).

The lattice action of our qubit regularized $O(4)$ model can be written using
four Grassmann valued lattice fields
$\psi_{1,k},\psi_{2,k},\bar{\psi}_{1,k},\bar{\psi}_{2,k}$ at each lattice site
$k\equiv({\mathbf{r},\tau})$ on a cubic lattice, where we distinguish between
the two-dimensional spatial coordinate ${\bf r}$ and the Euclidean temporal
coordinate $\tau$. The Euclidean action of our model is given by Banerjee  _et
al._ (2019),

$\displaystyle S(\psi,{\overline{\psi}})\ =\ -\sum_{\langle
k,k^{\prime}\rangle}\mathrm{Tr}(M_{k}M_{k^{\prime}})-\frac{U}{2}\sum_{k}\mathrm{Det}(M_{k}),$
(3)

where $\big{(}M_{k}\big{)}_{a,b}=\psi_{a,k}{\overline{\psi}}_{b,k}$ is a
$2\times 2$ matrix defined at each lattice site $k$. The symbol $\langle
k,k^{\prime}\rangle$ refers to neighboring sites $k$ and $k^{\prime}$. The
partition function is defined as usual through the Grassmann integral

$\displaystyle Z\ =\ \int\prod_{a,k}\ [{\overline{\psi}}_{a,k}\psi_{a,k}]\ \
e^{-S(\psi,{\overline{\psi}})}.$ (4)

It is easy to verify that the action is invariant under the $SU_{L}(2)\times
SU_{R}(2)$ transformations given by $M_{k}\rightarrow L\ M_{k}\ R^{\dagger}$
when $k\in\mbox{even sites}$, and $M_{k}\rightarrow R\ M_{k}\ L^{\dagger}$,
when $k\in\mbox{odd sites}$. Here we assume $L$ and $R$ are $2\times 2$
matrices, each of which are elements of the $SU(2)$ group. This means
$(\psi_{1,k},\psi_{2,k})$ on even sites and
$(-i\overline{\psi}_{2,k},i\overline{\psi}_{1,k})$ on odd sites transform as
$SU_{L}(2)$ doublets, while they are singlets of $SU_{R}(2)$. The same fields
on the opposite parity sites transform as $SU_{R}(2)$ doublets and $SU_{L}(2)$
singlets. When $U=0$, the theory has an additional $U(1)$ symmetry:
$\psi_{a,x}\rightarrow e^{i\theta}\psi_{a,x}$ and
${\overline{\psi}}_{a,x}\rightarrow e^{i\theta}{\overline{\psi}}_{a,x}$ for
the odd sites, and $\psi_{a,x}\rightarrow e^{-i\theta}\psi_{a,x}$ and
${\overline{\psi}}_{a,x}\rightarrow e^{-i\theta}{\overline{\psi}}_{a,x}$ on
the even sites. The $U$-term, therefore, mimics the anomalous axial symmetry
of the action in Quantum Chromodynamics (QCD) Cecile and Chandrasekharan
(2008). For this reason the terms in the partition function that arise due to
a non-zero value of $U$ were referred to as instantons in the earlier work. It
is known that instantons are known to break the anomalous axial symmetry in
QCD. In the more recent view point of qubit regularization, the instantons can
be viewed as simply the local Fock vacuum states Singh (2019).

Figure 1: An illustration of an $O(4)$ worldline configuration in $1+1$
dimensions with $N_{I}=12$ instantons (blue dots). The weight of the
configuration is $U^{12}$. Each loop is an oriented loop representing a red
worldline or a green worldline.

It is possible to perform the Grassmann integrations in Eq. 4 exactly and
rewrite the partition as a sum over worldline configurations of pions and
instantons Cecile and Chandrasekharan (2008). One then obtains $Z\ =\
\sum_{[\ell]}U^{N_{I}}$, where $[\ell]$ is a configuration of worldlines,
which is a collection of closed oriented loops, each of which can be in one of
two colors, red or green. In addition, there are isolated sites which do not
belong to the worldlines and are referred to as instantons. $N_{I}$ is the
total number of instantons in the configuration. An illustration of a
worldline configuration on a two-dimensional lattice is shown in Fig. 1. In an
earlier work, we solved our lattice model using this worldline approach and
showed that at the critical coupling of $U_{c}\approx 1.655394$ we can
reproduce the critical scaling of the Wilson-Fisher fixed point at long
distances Banerjee  _et al._ (2019).

##  III Leading Charge Sectors

Our goal is to compute the conformal dimensions $D(j_{L},j_{R})$ of field
operators that transform in some irreducible representation $(j_{L},j_{R})$ of
the $SU_{L}(2)\times SU_{R}(2)$ group at the critical point. We can accomplish
this by computing the correlation function

$\displaystyle C_{j_{L},j_{R}}$ $\displaystyle=\langle\overline{\cal
O}_{j_{L},j_{R}}{\cal O}_{j_{L},j_{R}}\rangle$
$\displaystyle=\frac{1}{Z}\int\prod_{a,k}[{\overline{\psi}}_{a,k}\psi_{a,k}]\
\ e^{-S(\psi,{\overline{\psi}})}\overline{\cal O}_{j_{L},j_{R}}{\cal
O}_{j_{L},j_{R}},$ (5)

where ${\cal O}_{j_{L},j_{R}}$ and $\overline{\cal O}_{j_{L},j_{R}}$ are the
source and sink terms constructed with Grassmann valued fields that transform
in the irreducible representation $(j_{L},j_{R})$. In our work, the source
will be located at the temporal slice $\tau=0$, while the sink will be at the
temporal slice $\tau=L/2$. Thus, at the critical point, we expect

$\displaystyle C_{j_{L},j_{R}}=A_{j_{L},j_{R}}L^{-2D(j_{L},j_{R})}$ (6)

for sufficiently large values of $L$.

In order to construct ${\cal O}_{j_{L},j_{R}}$ and $\overline{\cal
O}_{j_{L},j_{R}}$ that transform under the irreducible representation
$(j_{L},j_{R})$, let us denote ${|{(j_{L},m_{L})};{(j_{R},m_{R})}\rangle}$ as
the $(2j_{L}+1)(2j_{R}+1)$ dimensional orthonormal basis that spans the
irreducible representation space of $(j_{L},j_{R})$. Here $-j_{L}\leq
m_{L}\leq j_{L}$ and $-j_{R}\leq m_{R}\leq j_{R}$. We can label the fields
that transform according to this irreducible representation as ${\cal
O}_{{|{(j_{L},m_{L})};{(j_{R},m_{R})}\rangle}}$. For the sink terms, it is
natural to choose fields that transform in the conjugate representation
${\langle{(j_{L},m_{L})};{(j_{R},m_{R})}|}$. We can label these sink fields as
$\overline{\cal O}_{{\langle{(j_{L},m_{L})};{(j_{R},m_{R})}|}}$. While we can
choose any of these fields as source and sink terms in the $(j_{L},j_{R})$
representation, we find that choosing $m_{L}=j_{L}$ and $m_{R}=j_{R}$ will be
the most convenient choice for numerical work.

Based on the transformation property of $M_{k}$ it is easy to see that the
four local fermion bilinear lattice fields $i\psi_{1,k}\bar{\psi}_{1,k}$,
$-i\psi_{2,k}\bar{\psi}_{2,k}$, $-i\psi_{1,k}\bar{\psi}_{2,k}$,
$-i\psi_{2,k}\bar{\psi}_{1,k}$ transform under the $(1/2,1/2)$ (vector)
representation of $O(4)$. However, the fields transform differently on even
and odd sites. The exact mapping is given in the table in Table 1.

local field | even site | odd site
---|---|---
${\cal O}_{{|{(1/2,1/2)};{(1/2,1/2)}\rangle}}$ | $-i\psi_{1,k}{\overline{\psi}}_{2,k}$ | $-i\psi_{1,k}{\overline{\psi}}_{2,k}$
${\cal O}_{{|{(1/2,-1/2)};{(1/2,-1/2)}\rangle}}$ | $i\psi_{2,k}{\overline{\psi}}_{1,k}$ | $i\psi_{2,k}{\overline{\psi}}_{1,k}$
${\cal O}_{{|{(1/2,1/2)};{(1/2,-1/2)}\rangle}}$ | $i\psi_{1,k}{\overline{\psi}}_{1,k}$ | $-i\psi_{2,k}{\overline{\psi}}_{2,k}$
${\cal O}_{{|{(1/2,-1/2)};{(1/2,1/2)}\rangle}}$ | $-i\psi_{2,k}{\overline{\psi}}_{2,k}$ | $i\psi_{1,k}{\overline{\psi}}_{1,k}$
|  |
Table 1:  Local fermion bilinear fields transform according to the vector
representation of $O(4)$. This table gives the explicit realization of each
component at the lattice site $k$. As can be seen the realization can depend
on whether the site is even or odd.

Given the source fields we can construct the sink fields through the conjugate
representation. This is given in Table 2.

sink |  | source
---|---|---
$\overline{\cal O}_{{\langle{(1/2,1/2)};{(1/2,1/2)}|}}$ | = | $-{\cal O}_{{|{(1/2,-1/2)};{(1/2,-1/2)}\rangle}}$
$\overline{\cal O}_{{\langle{(1/2,-1/2)};{(1/2,-1/2)}|}}$ | = | $-{\cal O}_{{|{(1/2,1/2)};{(1/2,1/2)}\rangle}}$
$\overline{\cal O}_{{\langle{(1/2,1/2)};{(1/2,-1/2)}|}}$ | = | ${\cal O}_{{|{(1/2,-1/2)};{(1/2,1/2)}\rangle}}$
$\overline{\cal O}_{{\langle{(1/2,-1/2)};{(1/2,1/2)}|}}$ | = | ${\cal O}_{{|{(1/2,1/2)};{(1/2,-1/2)}\rangle}}$
|  |
Table 2:  Relationship between source and sink fields on a local site.

Using the information in Table 1 and Table 2 we can identify each worldline in
Fig. 1 as a vector particle carrying an appropriate charge. For example, we
can identify the sources and sinks of red worldlines as ${\cal
R}_{k}=-i\psi_{1,k}{\overline{\psi}}_{2,k}$ and ${\overline{\cal
R}}_{k}=-i\psi_{2,k}{\overline{\psi}}_{1,k}$ on all sites. On the other hand
the sources and sinks of green worldlines depend on the parity of the sites.
The green source is given by ${\cal G}_{k}=i\psi_{1,k}{\overline{\psi}}_{1,k}$
on even sites, and ${\cal G}_{k}=-i\psi_{2,k}{\overline{\psi}}_{2,k}$ on odd
sites. This is reversed for the sink of green lines. We have ${\overline{\cal
G}}_{k}=-i\psi_{2,k}{\overline{\psi}}_{2,k}$ on even sites, and
${\overline{\cal G}}_{k}=i\psi_{1,k}{\overline{\psi}}_{1,k}$ on odd sites.

Figure 2: The figure shows the arrangement of local sources of red worldlines
we used to create a source in the $(j,j)$ representation. The site labeled $0$
is the origin is left empty. Depending on $j$ the sources are placed on the
sites marked $k=1,2...2j$. The first $20$ sites used in our calculations up to
$j=10$ are shown.

In order to construct sources of more general irreducible representations, we
use tensor product of the vector representation of local fields distributed
over several spatial lattice sites on the time slice $\tau=0$. The same sites
on the time slice $\tau=L/2$ are used to construct the sinks. Let us first
construct sources and sinks that transform in the representation $(j,j)$ which
we refer to as the leading sector. We know we can construct the basis state
${|{(j,j)};{(j,j)}\rangle}$ as a tensor product of $2j$ states of the form
${|{(1/2,1/2)};{(1/2,1/2)}\rangle}$. The source ${\cal
O}_{{|{(1/2,1/2)};{(1/2,1/2)}\rangle}}$ can also be constructed in the same
way, placing $2j$ local sources of red worldlines on nearby sites on a spatial
lattice. In our work we choose these lattice sites $k$ labeled as $1,2,...,2j$
around the origin as shown in Fig. 2. Hence one of the possible sources in the
$(j,j)$ representation is then given by ${\cal
O}_{{|{(j,j)};{(j,j)}\rangle}}={\cal R}_{1}{\cal R}_{2}...{\cal R}_{2j}$,
while the corresponding sink is given by $\overline{\cal
O}_{{\langle{(j,j)};{(j,j)}|}}=\overline{\cal R}_{1}\overline{\cal
R}_{2}...\overline{\cal R}_{2j}$. Here the sites chosen for both source and
sink are the same spatial sites on the $\tau=0,L/2$ time slice. These sources
and sinks can be used to compute the correlation function

$\displaystyle C_{j,j}=\big{\langle}\overline{\cal
O}_{{\langle{(j,j)};{(j,j)}|}}{\cal
O}_{{|{(j,j)};{(j,j)}\rangle}}\big{\rangle}$ (7)

and in order to obtain $D(j,j)$ one can use the relation $C_{j,j}\sim
A_{j,j}L^{-2D(j,j)}$ for large values of $L$. In the actual worm algorithm,
one in fact computes the ratio $R_{j}=C_{j,j}/C_{j-1,j-1}$ and fits it to the
form $(A_{j,j}/A_{j-1,j-1})L^{-2\Delta_{j}}$ to compute
$\Delta_{j}=D(j,j)-D(j-1,j-1)$ for each value of $j$. From these differences
and setting $D(0,0)=1$ one can compute $D(j,j)$. Our results from Banerjee
_et al._ (2019) are tabulated in Table 3 for reference.

$j$ | $D(j,j)$ | $j$ | D(j,j)
---|---|---|---
1/2 | 0.515(3) | 1 | 1.185(4)
3/2 | 1.989(5) | 2 | 2.915(6)
5/2 | 3.945(6) | 3 | 5.069(7)
7/2 | 6.284(8) | 4 | 7.575(9)
9/2 | 8.949(10) | 5 | 10.386(11)
|  |  |
Table 3:  Results for the conformal dimensions $D(j,j)$ up to $j=5$ computed
using worldline Monte-Carlo methods in Banerjee  _et al._ (2019). Figure 3:
Worldline Monte Carlo results for the ratio $\tilde{R}_{j}$ as a function of
$L$ for various values of $j$. The solid lines show the fits given in Table 4.

##  IV Subleading Sector

In this work we extend our earlier results in the leading sector and compute
the conformal dimensions of the subleading sector, $D(j,j-1)$ for $j\geq 1$.
For this we need to construct source and sink operators that transform in the
representation $(j,j-1)$. We know we can construct the state
${|{j,j};{j-1,j-1}\rangle}$ by applying the lowering operator $J_{R}^{-}$ to
the ${|{j,j};{j,j}\rangle}$ and then constructing orthogonal states in the
tensor product space. This procedure naturally leads to $2j-1$ orthonormal
states, which we can label with an additional index $M=1,2...,(2j-1)$.
Translating this to the construction of sources, we now introduce the source
${\cal O}_{\ell}={\cal R}_{1}{\cal R}_{2}...{\cal G}_{\ell}...{\cal R}_{2j}$
where the red source on one lattice site $\ell$ is replaced by a green source
where $\ell=1,2...,2j$. Similarly, we introduce the corresponding sinks as
$\overline{\cal O}_{\ell}=\overline{\cal R}_{1}\overline{\cal
R}_{2}...\overline{\cal G}_{\ell}...\overline{\cal R}_{2j}$. We can then argue
that

$\displaystyle{\cal O}_{{|{(j,j)};{(j,j-1)}\rangle}}$
$\displaystyle=\frac{1}{\sqrt{2j}}\sum_{\ell=1}^{2j}{\cal O}_{\ell}.$ (8)

where the right-hand side is a sum over the $2j$ source terms we introduced
above. Note that there are $2j-1$ sources orthogonal to Eq. 8, which can label
with $M=1,2...,2j-1$ as before. These will naturally transform in the
$(j,j-1)$ representation. Explicitly these sources are given by

$\displaystyle{\cal O}^{M}_{{|{(j,j)};{(j-1,j-1)}\rangle}}$
$\displaystyle=\frac{1}{\sqrt{2j}}\sum_{\ell=1}^{2j}e^{i2\pi(\ell-1)M/(2j)}{\cal
O}_{\ell}.$ (9)

We can similarly define the corresponding $2j-1$ sinks as

$\displaystyle\overline{\cal O}^{M}_{{\langle{(j,j)};{(j-1,j-1)}|}}$
$\displaystyle=\frac{1}{\sqrt{2j}}\sum_{\ell=1}^{2j}e^{-i2\pi(\ell-1)M/(2j)}\overline{\cal
O}_{\ell}.$ (10)

Since all $2j-1$ sources and sinks transform under the same irreducible
representation $(j,j-1)$ any combination of them can be used in Eq. 5 to
extract $D(j,j-1)$. Let us define the correlation matrix

$\displaystyle C^{MM^{\prime}}_{j,j-1}=\big{\langle}\overline{\cal
O}^{M}_{{\langle{(j,j)};{(j-1,j-1)}|}}\ {\cal
O}^{M^{\prime}}_{{|{(j,j)};{(j-1,j-1)}\rangle}}\big{\rangle}$ (11)

where we expect $C^{MM^{\prime}}_{j,j-1}\sim
A^{MM^{\prime}}_{j,j-1}L^{-2D(j,j-1)}$ for large values of $L$. Practically it
is more convenient to compute the average of the trace of this correlation
matrix which can be simplified to the form

$\displaystyle C_{j,j-1}\ =\
\sum_{\ell}\Big{\\{}\frac{1}{2j}\big{\langle}\overline{\cal O}_{\ell}{\cal
O}_{\ell}\big{\rangle}-\frac{1}{2j(2j-1)}\sum_{\ell^{\prime}\neq\ell}\big{\langle}\overline{\cal
O}_{\ell}{\cal O}_{\ell^{\prime}}\big{\rangle}\Big{\\}}.$ (12)

Note that we also expect $C_{j,j-1}\sim A_{j,j-1}L^{-2D(j,j-1)}$.

$j$ |  $L$-range | $A_{j,j-1}/A_{j,j}$ | ${\widetilde{\Delta}}(j)$ |  $\chi^{2}$/DOF
---|---|---|---|---
$1$ | $24-128$ | $5.87(25)$ | $0.813(6)$ | $1.11$
$3/2$ | $24-128$ | $2.50(11)$ | $0.750(6)$ | $0.62$
$2$ | $24-96$ | $2.13(06)$ | $0.722(4)$ | $0.26$
$5/2$ | $32-96$ | $1.75(08)$ | $0.685(6)$ | $1.28$
$3$ | $32-96$ | $1.54(08)$ | $0.659(7)$ | $0.93$
$7/2$ | $32-96$ | $1.35(05)$ | $0.633(5)$ | $0.38$
$4$ | $32-96$ | $1.18(04)$ | $0.607(4)$ | $0.40$
$9/2$ | $40-160$ | $1.05(04)$ | $0.586(5)$ | $0.94$
$5$ | $40-160$ | $0.94(04)$ | $0.566(5)$ | $0.89$
$11/2$ | $48-160$ | $0.90(03)$ | $0.555(4)$ | $0.66$
$6$ | $48-160$ | $0.83(03)$ | $0.541(5)$ | $1.40$
$13/2$ | $64-160$ | $0.75(04)$ | $0.525(7)$ | $1.11$
$7$ | $64-160$ | $0.71(03)$ | $0.513(5)$ | $1.18$
$15/2$ | $64-160$ | $0.69(04)$ | $0.506(6)$ | $1.45$
$8$ | $64-160$ | $0.60(03)$ | $0.486(5)$ | $0.77$
$17/2$ | $64-160$ | $0.61(03)$ | $0.484(5)$ | $0.83$
$9$ | $80-160$ | $0.53(04)$ | $0.467(8)$ | $0.98$
$19/2$ | $80-160$ | $0.53(03)$ | $0.463(7)$ | $0.47$
$10$ | $80-160$ | $0.50(02)$ | $0.454(5)$ | $0.63$
$20$ | $96-256$ | $0.28(01)$ | $0.367(3)$ | $0.65$
|  |  |  |
Table 4:  Results of the fit of $\tilde{R}_{j}$ shown in Fig. 3 to the form
$(A_{j}/A_{j-1})L^{-2\tilde{\Delta}_{j}}$. The range of $L$ values used in the
fit are given in the second column. We observe that as $j$ increases this
range needs to involve larger lattice sizes for a good fit.

As in the leading sector, in the worldline algorithm it much easier to compute
the ratio $\tilde{R}_{j}=C_{j,j-1}/C_{j,j}$. For this one constructs a
worldline Monte Carlo method to generate configurations with $2j$ red sources
and sinks that contribute to $C_{j,j}$. In every configuration of this
ensemble, we then imagine flipping each of the $2j$ sources located at the
sites $\ell=1,2...,2j$ to a green source. The worldline of the green source
then naturally travels through the lattice to a sink at some location
$\ell^{\prime}$. Then we compute the contribution to $\tilde{R}_{j}$ from that
configuration using Eq. 12, which means we add $1/2j$ if $\ell=\ell^{\prime}$,
or subtract the value $1/(2j(2j-1))$ if $\ell\neq\ell^{\prime}$ for every
value of $\ell$. Averaging this contribution over the ensemble of
configurations generated by the worldline algorithm gives us the ratio,
$\tilde{R}_{j}$ which is expected to scale as
$(A_{j}/A_{j-1})L^{-2\tilde{\Delta}_{j}}$ where
$\tilde{\Delta}_{j}=D(j,j-1)-D(j,j)$. Using the values of $D(j,j)$ we compute
$D(j,j-1)$.

$j$ | $\lambda_{0}$ | $\lambda_{1/2}$ | $\lambda_{1}$ | $\lambda_{3/2}$ | $\chi^{2}$
---|---|---|---|---|---
range |  |  |  |  | /DOF
1-20 | 0.07(2) | 1.5(1) | -1.0(2) | 0.3(1) | 0.6
1-20 | 0 | 1.42(1) | 0 | -0.69(1) | 48
1-20 | 0 | 1.96(2) | -1.83(6) | 0.69(5) | 1.8
1-20 | 0.16(1) | 0.98(2) | 0 | -0.34(2) | 2.3
1.5-20 | 0.13(1) | 1.07(2) | 0 | -0.47(3) | 0.3
2-20 | 0 | 2.06(3) | -2.27(14) | 1.14(13) | 1.1
8-20 | 0 | 1.81(2) | 0 | -3.5(2) | 1.1
8-20 | 0 | 2.09(4) | -2.02(13) | 0 | 0.6
|  |  |  |  |
Table 5:  Fits of ${\widetilde{\Delta}}_{j}$ to the functional form given in
Eq. 2. We first consider the whole range of $j$ in the first four rows. While
including all four coefficients as fitting parameters gives an excellent fit,
setting the two purely quantum mechanical terms $\lambda_{0}=\lambda_{1}=0$
makes the fit quite bad. Including even one of them is sufficient to improve
the fit considerably. Setting two of the fitting parameters to zero is only
possible by shrinking the allowed region of $j$ considerably.

##  V Results

We have performed a series of Monte Carlo calculations for $1\leq j\leq 10$ in
increments of $1/2$, and for $j=20$. For runs till $j=10$, the lattice sizes
ranged from $L=32$ up to $L=160$. For each value of $j$ we have always found
that our data fits well to the expected form $C_{j,j-1}/C_{j,j}\sim
L^{-2\tilde{\Delta}_{j}}$ for sufficiently large values of $L$. As examples,
fits for $j=1,5$, $10$, and $20$ are illustrated in Fig. 3. From the figure,
we observe that as $j$ increases the range of lattice sizes where a simple
power law emerges changes to larger $L$ values. The range of $L$ values where
we perform the fits and the obtained fit parameters for various values of $j$
are tabulated in 4. At $j=20$ we have extended our calculations up to lattice
sizes of $L=256$.

Figure 4: The plot of ${\widetilde{\Delta}}_{j}$ as a function of $j$ from
Table 4. The lines shown are fits to the form given in Eq. 2 with parameter
values given in the first (solid), fifth (dashes), sixth (dashes+dots) rows
and seventh (dots) rows of Table 5.

Having obtained the values of $\tilde{\Delta}_{j}$ for various values of $j$
we try to extract the constants $\lambda_{i}$’s in the large charge expansion
based on Eq. 2, assuming we can neglect ${\cal O}(1/j^{2})$ terms. We try
several fits to understand the role of the purely quantum terms $\lambda_{0}$
and $\lambda_{1}$. These fits are shown in Table 5. First, we note that our
data for the entire range of $j$ values fits well to the form Eq. 2, if we
assume all four coefficients are non-zero (first row in Table 5). On the other
hand if we drop both quantum terms the fit becomes quite bad (second row in
Table 5). The presence of either of the two quantum terms is sufficient to
bring down the $\chi^{2}/DOF$ considerably. For example setting
$\lambda_{1}=0$, we can get an excellent fit if we just drop the $j=1$ data
(fifth row in Table 5). On the other hand setting $\lambda=0$ and dropping
both $j=1,1.5$ makes the fit acceptable (sixth row in Table 5). If we drop
both quantum terms (i.e., set $\lambda_{0}=\lambda_{1}=0$) we can only get a
good fit in the range $j=8-20$. We believe this is just an artifact of a small
range of $j$ as can be seen in Fig. 4.

Clearly, our data is consistent with the conjecture that $\lambda_{0}=0$ (see
sixth row in Table 5). Assuming this, we can try to determine the leading two
terms more reliably by fitting our data under the constraint
$\lambda_{3/2}=0$. In this case we can get a good fit only in the smaller
range $j=8-20$ (eighth row in Table 5). As expected, this changes
$\lambda_{1/2}$ and $\lambda_{1}$ slightly. Taking such systematic fitting
effects into account, we estimate that $\lambda_{1/2}=2.08(5)$,
$\lambda_{1}=2.2(3)$. Unfortunately, calculations at higher values of $j$ are
difficult since we need larger lattice sizes, but they can be obtained with
more computing resources and can help confirm the conjecture that
$\lambda_{0}=0$.

##  VI Conclusions

In this work we have constructed a Monte Carlo method to compute the sub-
leading conformal dimensions in the large charge expansion at the $O(4)$
Wilson-Fisher fixed point. We used this method to compute
${\widetilde{\Delta}}_{j}$ for several values of $j$ in the range $1\leq j\leq
20$. While our results are consistent with the general predictions of the
expansion given in Eq. 2, we cannot rule out the possibility that
${\widetilde{\Delta}}_{j}$ approaches a non-zero constant $\lambda_{0}$ in the
large $j$ limit. However, our data is consistent with the conjecture that
$\lambda_{0}=0$. Assuming this to be true, we can estimate the leading two
terms in the expansion to be $\lambda_{1/2}=2.08(5)$, $\lambda_{1}=2.2(3)$.
Since $\lambda_{1}$ must be a calculable number within the large charge
effective field theory, we hope our work will motivate someone to calculate it
in the future. Calculations at large $j$ may also help determine $\lambda_{0}$
reliably.

##  VII Acknowledgments

We thank D. Orlando, S. Reffert for extensive discussions and collaborating
with us in the past. We also thank S. Hellerman for useful discussions about
this work. The material presented here is based upon work supported by the
U.S. Department of Energy, Office of Science, Nuclear Physics program under
Award Number DE-FG02-05ER41368.

## References

  * Rychkov (2016) S. Rychkov,  _EPFL Lectures on Conformal Field Theory in D $>$= 3 Dimensions_, SpringerBriefs in Physics (Springer, 2016) arXiv:1601.05000 [hep-th] .
  * Simmons-Duffin (2017) D. Simmons-Duffin, in  _Theoretical Advanced Study Institute in Elementary Particle Physics: New Frontiers in Fields and Strings_ (2017) pp. 1–74, arXiv:1602.07982 [hep-th] .
  * Komargodski and Zhiboedov (2013) Z. Komargodski and A. Zhiboedov, JHEP 11, 140 (2013), arXiv:1212.4103 [hep-th] .
  * Hellerman  _et al._ (2015) S. Hellerman, D. Orlando, S. Reffert, and M. Watanabe, JHEP 12, 071 (2015), arXiv:1505.01537 [hep-th] .
  * Kaviraj  _et al._ (2015a) A. Kaviraj, K. Sen, and A. Sinha, JHEP 11, 083 (2015a), arXiv:1502.01437 [hep-th] .
  * Kaviraj  _et al._ (2015b) A. Kaviraj, K. Sen, and A. Sinha, JHEP 07, 026 (2015b), arXiv:1504.00772 [hep-th] .
  * Alday (2017) L. F. Alday, Phys. Rev. Lett. 119, 111601 (2017), arXiv:1611.01500 [hep-th] .
  * Gopakumar  _et al._ (2017) R. Gopakumar, A. Kaviraj, K. Sen, and A. Sinha, Phys. Rev. Lett. 118, 081601 (2017), arXiv:1609.00572 [hep-th] .
  * Dey  _et al._ (2018) P. Dey, K. Ghosh, and A. Sinha, JHEP 01, 152 (2018), arXiv:1709.06110 [hep-th] .
  * Caron-Huot (2017) S. Caron-Huot, JHEP 09, 078 (2017), arXiv:1703.00278 [hep-th] .
  * Hellerman and Maeda (2017) S. Hellerman and S. Maeda, JHEP 12, 135 (2017), arXiv:1710.07336 [hep-th] .
  * Hellerman  _et al._ (2017) S. Hellerman, S. Maeda, and M. Watanabe, JHEP 10, 089 (2017), arXiv:1706.05743 [hep-th] .
  * Jafferis  _et al._ (2018) D. Jafferis, B. Mukhametzhanov, and A. Zhiboedov, JHEP 05, 043 (2018), arXiv:1710.11161 [hep-th] .
  * Alvarez-Gaume  _et al._ (2019) L. Alvarez-Gaume, D. Orlando, and S. Reffert, JHEP 12, 142 (2019), arXiv:1909.02571 [hep-th] .
  * Orlando  _et al._ (2020) D. Orlando, S. Reffert, and F. Sannino, Phys. Rev. D 101, 065018 (2020), arXiv:1909.08642 [hep-th] .
  * Poland  _et al._ (2019) D. Poland, S. Rychkov, and A. Vichi, Rev. Mod. Phys. 91, 015002 (2019), arXiv:1805.04405 [hep-th] .
  * Dondi  _et al._ (2021) N. Dondi, I. Kalogerakis, D. Orlando, and S. Reffert, JHEP 05, 035 (2021), arXiv:2102.12488 [hep-th] .
  * Monin (2016) A. Monin, Phys. Rev. D 94, 085013 (2016).
  * de la Fuente (2018) A. de la Fuente, Journal of High Energy Physics 2018, 41 (2018).
  * Banerjee  _et al._ (2018) D. Banerjee, S. Chandrasekharan, and D. Orlando, Phys. Rev. Lett. 120, 061603 (2018).
  * Banerjee  _et al._ (2019) D. Banerjee, S. Chandrasekharan, D. Orlando, and S. Reffert, Phys. Rev. Lett. 123, 051603 (2019), arXiv:1902.09542 [hep-lat] .
  * Gaumé  _et al._ (2021) L. A. Gaumé, D. Orlando, and S. Reffert, Phys. Rept. 933, 2180 (2021), arXiv:2008.03308 [hep-th] .
  * Err (2020) (2020), in Ref. Banerjee  _et al._ (2019) it was predicted that $\lambda_{0}\neq 0$, based on an incorrect assumption. This error was later corrected in Ref. Gaumé  _et al._ (2021).
  * Alvarez-Gaume  _et al._ (2017) L. Alvarez-Gaume, O. Loukas, D. Orlando, and S. Reffert, JHEP 04, 059 (2017), arXiv:1610.04495 [hep-th] .
  * Hellerman  _et al._ (2019) S. Hellerman, N. Kobayashi, S. Maeda, and M. Watanabe, JHEP 10, 038 (2019), arXiv:1705.05825 [hep-th] .
  * Antipin  _et al._ (2020) O. Antipin, J. Bersini, F. Sannino, Z.-W. Wang, and C. Zhang, Phys. Rev. D 102, 045011 (2020), arXiv:2003.13121 [hep-th] .
  * Com (2021) (2021), we thank D. Orlando for helping us understand the discussion presented in this paragraph.
  * Cecile and Chandrasekharan (2008) D. Cecile and S. Chandrasekharan, Phys. Rev. D 77, 014506 (2008), arXiv:0708.0558 [hep-lat] .
  * Singh and Chandrasekharan (2019) H. Singh and S. Chandrasekharan, Phys. Rev. D 100, 054505 (2019), arXiv:1905.13204 [hep-lat] .
  * Singh (2019) H. Singh, arXiv:1911.12353 (2019), arXiv:1911.12353 [hep-lat] .
  * Bhattacharya  _et al._ (2020) T. Bhattacharya, A. J. Buser, S. Chandrasekharan, R. Gupta, and H. Singh, arXiv:2012.02153 (2020), arXiv:2012.02153 [hep-lat] .