marginparsep has been altered.
topmargin has been altered.
marginparwidth has been altered.
marginparpush has been altered.
The page layout violates the ICML style. Please do not change the page layout,
or include packages like geometry, savetrees, or fullpage, which change it for
you. We’re not able to reliably undo arbitrary changes to the style. Please
remove the offending package(s), or layout-changing commands and try again.

AutoNLU: Detecting, root-causing, and fixing NLU model errors

Anonymous Authors1

###### Abstract

Improving the quality of Natural Language Understanding (NLU) models, and more
specifically, task-oriented semantic parsing models, in production is a
cumbersome task. In this work, we present a system called AutoNLU, which we
designed to scale the NLU quality improvement process. It adds automation to
three key steps: detection, attribution, and correction of model errors, i.e.,
bugs. We detected four times more failed tasks than with random sampling,
finding that even a simple active learning sampling method on an uncalibrated
model is surprisingly effective for this purpose. The AutoNLU tool empowered
linguists to fix ten times more semantic parsing bugs than with prior manual
processes, auto-correcting 65% of all identified bugs.

††footnotetext:  1Anonymous Institution, Anonymous City, Anonymous Region,
Anonymous Country. Correspondence to: Anonymous Author
<EMAIL_ADDRESS>

##  1 Introduction

Natural Language Understanding (NLU) models are widely used by virtual
assistants to complete task-oriented requests. For example, given a request in
natural language like, “Play my running playlist,” an NLU model returns a
structured representation that could be more easily executed downstream.
Typically, NLU models classify a user’s request into a domain, and further
parse the request into intents and slots Mesnil et al. (2013); Gupta et al.
(2018). 111NLU can include other NLP tasks such as disambiguation, coreference
resolution, etc. For the purposes of this paper, we use NLU to refer to task-
oriented semantic parsing.

Developers, or linguists, frequently desire to add new domains, or to expand
the coverage of existing domains to support a wider variety of user requests.
Unfortunately, ensuring the quality of new domains is high or continuously
improving existing domains is non-trivial. Each step of the quality
improvement process, from discovering areas where the model is failing,
understanding why it is failing, and finally, applying a fix, is highly
manual.

We first look at discovering model bugs. Traditional approaches to discovering
model bugs are either too fine or too coarse-grained. The traditional software
engineering quality assurance (QA) approach is too fine-grained. In the QA
approach, the model is hand-tested one-by-one, either with utterances that are
construed to be prototypical of that domain or adversarial examples. The
traditional ML approach is too coarse-grained. In the ML approach, a model is
evaluated on a randomly-sampled validation set and test set; metrics on these
sets may obfuscate the presence of bugs in a relatively small domain. Slicing,
as discussed by Ribeiro et al. (2020), can be a great way of finding and
inspecting bugs that would otherwise be hidden. However, there is still an
opportunity cost in annotating a large validation and test set where bugs may
be relatively rare. Active learning is often used to build a training set in a
cost-effective way Monarch (2021); however, it is not typically framed with
the explicit goal of finding bugs.

Once bugs are detected, developers hope to correct them. However, before they
can correct the bug, they must understand what is causing them. Unfortunately,
this is also non-trivial. Model interpretability tools can help, but
fundamentally bugs often stem from training data Breck et al. (2019), which
can be time-consuming to explore.

Finally, the quality improvement process is slowed down by a lack of tools to
quickly fix bugs. While a machine learning engineer’s instinct may be to try a
new architecture, pre-training, or a parameter sweep, a domain confusion could
be caused by lack of training data or even a mislabel. Linguists often have to
create data patches, i.e., data augmentations or relabels, by hand without
tools to assist them.

These three challenges significantly slow down iteration speed, and we propose
a more scalable approach to finding and fixing NLU bugs. In this work:

  * •

We describe and demonstrate a system tailored to scaling the quality
improvement process for task-oriented NLU models (Sections 2, 3.1).

  * •

We present a fresh take on the objective of active learning for NLU, using it
to continuously detect bugs as opposed to improve sample efficiency. We found
that even a simple uncertainty sampling method is four times more effective
than random sampling at finding utterances where the model will make a mistake
(Sections 3.2, 4.1).

  * •

We describe an approach to automatic error attribution and data augmentation
and show that they are effective ways of fixing NLU bugs in large volumes
(Sections 3.3, 3.4, 4.2).

##  2 AutoNLU

We designed a system called AutoNLU, which ties together and provides
automation around three key steps to quality improvement.

##### Detection

AutoNLU re-frames the goal of active learning to maximizing bug detection:
find all observable examples that may be parsed incorrectly, rank them by
importance (e.g., frequency), and annotate as many as possible within the
available annotation budget. We are particularly interested in correcting
“known unknowns” Monarch (2021). Traditional active learning setups, on the
other hand, aim to increase _sample efficiency_ : annotate the fewest possible
examples to reach acceptable performance on a test set. The flaw of this
framing is that it assumes the test distribution is static. In contrast, in
production settings, the test distribution is dynamic and thus it’s desirable
to make model improvements continuously. Moreover, sample efficiency is often
not a concern by the time the model is put in production – models are
bootstrapped with pre-labeled crowd-sourced or generated data.

##### Error Attribution

Once we have obtained a collection of bugs, AutoNLU attempts to attribute
their root cause. We provide a small set of high-level error classifications
that map directly to actionable next steps that a developer can take to fix
the bug. These classifications can be made more specific and accurate over
time.

##### Correction

Finally, based off the error classification from the previous step, AutoNLU
suggests an appropriate fix. For example, if the error is that there is not
enough training data similar to the bug utterance, AutoNLU generates similar
utterances that could be added to the training set.

###  2.1 Defining an NLU Bug

Using the conventions described by Gupta et al. (2018), given a user request
like, “Play my running playlist,” our NLU model should predict it belongs to
the music domain, and produce a parse like [IN:PLAY_MUSIC Play my
[SL:PLAYLIST_NAME running] [SL:MUSIC_TYPE playlist]],

where IN: is the prefix of the intent, and SL: is the prefix of all the slots
within the intent.222An analogy to an intent is a function or API call.
Similarly, and analogy to a slot is a function argument. If the domain or
intent-level classification of the sentence is incorrect, or if an expected
slot is missing or has the incorrect span, or if an extra slot is produced, we
consider the prediction incorrect and the corresponding input a bug. For
example, for a user request like, “Play my holiday cooking playlist,” if the
NLU model produced a parse like [IN:PLAY_MUSIC Play my holiday cooking
[SL:MUSIC_TYPE playlist]], we would consider the request a bug because the
SL:PLAYLIST_NAME slot is missing.

###  2.2 NLU Model Architecture

Throughout this paper, we loosely use the term NLU model to actually describe
a collection of two sets of models. The first set consists a single domain
classifier. The second set contains many sequence-to-sequence parsing models,
which jointly produce the intent and slot predictions. There is a parser per
domain (or per small grouping of similar domains). The domain classifier
routes the request to the appropriate intent-slot parser.

The domain classifier is a BiLSTM model with an MLP output layer Huang et al.
(2015). The intent-slot parsers have a RoBERTA Liu et al. (2019) encoder and
have either a non-autoregressive (NAR) decoder, similar to that described by
Shrivastava et al. (2021), or an autoregressive BERT decoder Chen et al.
(2019a).

##  3 System Design

The high-level system design of AutoNLU is shown in Figure 1.

Figure 1: This diagram illustrates the flow of production data. The components
highlighted in green are discussed in more detail in the following sections.

###  3.1 User Interface

The user interface (UI) for AutoNLU is a table showing a list of bugs for a
domain developer or linguist to fix. By default, the bugs are ranked by their
frequency. However, the ranking could also take into account other factors
such as recency or uncertainty score. The UI also shows the user the failure
reason and a suggested action to fix the bug. Examples are shown in Figures 2
and 3.

Figure 2: AutoNLU produces results as a pandas DataFrame Wes McKinney (2010);
Reback et al. (2021), viewable from a Jupyter notebook Kluyver et al. (2016) .
It displays confidence scores, the expected and predicted parse, the estimated
contribution of the bug to overall accuracy, the failure reason, and where
possible, an auto-generated correction. It also displays additional debugging
information such as the training set annotation, under cached_seqlogical.
Figure 3: This is a screenshot that illustrates the web app version of the
user interface for AutoNLU. It includes diff highlighting between the expected
and predicted parse. Note that the data shown here is dummy data.

###  3.2 Detection

When designing our algorithm for bug detection, we kept the following
principles in mind:

  * •

It should identify utterances that the NLU model will parse incorrectly.

  * •

It should provide a numerical score that can allow us to rank and compare
unlabeled production data (candidates) for sampling.

  * •

Ideally, it should be feasible to run on hundreds-of-thousands to millions of
examples. Else, we would have to decrease the size of our candidate set.

####  3.2.1 Simple Uncertainty Baseline

The algorithm we use for bug detection is least-confidence uncertainty
sampling, similar to that described by Monarch (2021). This method assigns
high uncertainty scores for candidates where the model prediction had low
confidence.

For all n candidates, we assign an uncertainty score of 1 - intent_confidence.
We then select the top k samples for human grading according to this score.

For simplicity, we do not use domain and slot confidence scores when computing
the uncertainty score. If the semantic frame contained nested intents, we only
use the score from the top-level intent. Combining uncertainty scores across
the semantic frame is an interesting future direction.

####  3.2.2 Intrinsic vs. Extrinsic Sampling Methods

One of the key considerations for the detection algorithm was whether we
should use an intrinsic or extrinsic approach. By intrinsic, we mean using the
model confidence, scores, or logits to decide what to sample. By extrinsic, we
mean training a meta-model that can predict when the base, i.e., NLU model
will make a mistake Koçak et al. (2021).

An advantage of the intrinsic approach is that, unlike extrinsic, it does not
require training and maintaining a separate model. However, the downside is
that it may rely on the base model being well-calibrated. Deep learning models
are known to produce over-confident predictions even when they are incorrect
He et al. (2020).

An advantage of the extrinsic approach is that it may allow you to consider
other features that are available in logs but not at inference (e.g., Did the
user repeat their request? – which could indicate frustration from an
incorrect NLU prediction).

We went with the intrinsic, least-confidence uncertainty method because it was
simple to implement and computationally inexpensive, since we already stored
the logged model prediction and confidence scores from run-time on our
candidates. Despite the lack of calibration, we found we were still able to
get good results as described in Section 4.1.

###  3.3 Error Attribution

Once we’ve sampled bugs, the next step is to identify their root cause. The
key principle we held in mind for error attribution is that the error
categories must be actionable, i.e., map to a concrete next step that can be
applied automatically or by a linguist to fix the bug. This mapping is shown
in Table 1 and the corrections are described in Section 3.4.

We determine the error attribution by sequentially checking the cases below,
in the order provided; the attribution falls through to unknown if none of the
criteria are met. Note that we use the term golden annotation to refer to the
annotation given to the bug utterance.

####  3.3.1 Rule Mismatch

In addition to using the NLU model to parse language requests, we have a small
set of high-precision rules that trigger on must-pass requests. Sometimes the
parse produced by a rule conflicts with the golden annotation provided for the
bug. In this case, either the rule or the golden annotation needs to be
updated so that they are in agreement.

####  3.3.2 Mislabeled

A surprisingly common cause of bugs is a labeling error, either in the
training data or of the golden annotation itself. We assign this attribution
whenever we see an exact match of the bug utterance in the training set, the
annotation in training conflicts with the golden annotation, and the model’s
confidence in its prediction is above some threshold $\lambda$. We set
$\lambda>=0.9$. These signals indicate that the golden annotation may be
incorrect.

A golden annotation can be incorrect for many reasons. First, labeling bugs is
a more cognitively difficult task than labeling training data. Bugs are not
pre-sorted by domain, and by nature of being bugs, may be more linguistically
challenging to annotate than typical training data. Second, the people who
create training data and label golden data are usually not the same. Training
data is created under the instruction of linguists, who are subject-matter
experts. Golden annotation of incoming production data is often done with the
help of third-party vendors. And third, annotation guidelines can be difficult
to keep up-to-date. Often the ontology (label set) is updated faster than
annotation guidelines can keep up.

####  3.3.3 Low Training Data

If a bug has made it past the first two filters, there’s a good sign that
there is a real issue with the model. Often bugs arise because there is not
enough training data similar to the bug in order for the model to generalize.
For the example used earlier, [IN:PLAY_MUSIC Play my holiday cooking
[SL:MUSIC_TYPE playlist]], the model may have not seen “holiday cooking,” or
something similar being used as playlist name.

We have a simple heuristic to check if a bug should get the low training data
error attribution: we check for an exact match of the bug utterance in the
training data, and if none exists, we assign this error. This heuristic could
over-trigger, and in the future we hope to change it to look for similar
training utterances as opposed to exact matches. For example, if the training
set did not have an exact match but had the utterance, “Play my Christmas
baking playlist,” then we should not assign the low training data attribution,
since this example should have been sufficient for the model to generalize.

####  3.3.4 Unknown

Finally, if none of the above criteria pass, we assign the bug to the unknown
category.

###  3.4 Correction

Once the error attribution is known, the appropriate correction can be
applied. We discuss these potential corrections in further detail.

Table 1: Each error attribution maps to an appropriate correction.
Error Attribution | Appropriate Correction
---|---
Rule Mismatch | Fix Rule
Mislabeled | Fix Annotation Conflicts
Low Training Data | Generate Data
Unknown | Generate Rule

####  3.4.1 Data Augmentation

In the low training data scenario, AutoNLU automatically applies templated
data augmentation with slot-based replacement. For example, given the bug
utterance, “Play my holiday cooking playlist,” it may produce the template
Play my $<$SL:PLAYLIST_NAME$>$ playlist. The template can then be expanded to
new training examples, such as “Play my baking playlist,” “Play my road trip
playlist,” etc.

While we used this templated approach for our experiments to have more control
of the generated data, one could conceivably use language models for text
generation instead.

####  3.4.2 Annotation Conflict Correction

In the case where there is a conflict between the golden annotation of the bug
and the training data, AutoNLU shows the user the conflicting annotations. The
user can then manually update the annotations, or use a labeling library for
semantic frames, which we call data transformers, to update them. The data
transformer library supports operations such as renaming or adding intents and
slots on large batches of data.

####  3.4.3 Rule Generation

In the rule mismatch case, we may need to update an incorrect rule. In the
case where the error attribution is unknown, and it’s not possible to manually
inspect and find the root cause of the bug, our only option may be to generate
a new high-precision rule,

Rule generation is simple to do given the golden annotation and utterance as
inputs. For our previous example, we could create the rule
IN:PLAY_MUSIC$<$Play my $<$SL:PLAYLIST_NAME holiday cooking$>$
$<$SL:MUSIC_TYPE playlist$>>$ which would trigger if the same utterance was
received at runtime. The downside of adding a rule as opposed to training data
for the NLU model is, of course, limited generalization.

##  4 Experiments

Our evaluation of the AutoNLU system focused on its ability to detect bugs and
to auto-correct them. We found that even our simple methods for sampling and
data generation led to significant increases in the number of bugs detected
and fixed.

###  4.1 Simple Uncertainty

Our evaluation of uncertainty sampling was split into two stages, an offline
and an online evaluation. In the offline evaluation, we aimed to understand
the relationship between NLU misclassifications and the uncertainty score. In
the online evaluation, we aimed to understand the effectiveness of uncertainty
sampling at finding failed tasks on production data.

####  4.1.1 Offline Evaluation

For our offline evaluation, we computed the simple uncertainty score of the
NLU model on 20,000 crowd-sourced test examples. As shown in Figure 4, there
was a sharp peak of low uncertainty scores. This is not surprising given that
NLU model’s accuracy on this test set was also high, $>90\%$,.

Figure 4: Distribution of uncertainty scores on 20,000 crowd-sourced test
examples.

Next, we observed how the histogram changes if we only looked at misclassified
samples. This allowed us to get a sense of how correlated the uncertainty
metric is with misclassified points. A perfect histogram would have a sharp
peak at 1.0 and very little presence at 0.

As shown in Figure 5 (left), we found that the model is unfortunately highly
certain (has low uncertainty) about most of the misclassified samples.
However, compared to Figure 5 (right), which only contains correctly
classified samples, there is a much greater tendency for misclassified samples
to have high uncertainty scores than there is for correctly classified samples
to have high uncertainty scores.

Thus, despite the NLU model being poorly calibrated – in fact, only 27% of the
misclassified samples have an uncertainty higher than 0.5 – we found that for
sufficiently small values of $k$ i.e. $k<<n$, sampling the top $k$ utterances
by uncertainty score is likely to yield mostly misclassifications i.e., NLU
bugs.

Figure 5: Distribution of uncertainty scores on 20,000 crowd-sourced test
examples. The histogram on the left only contains misclassified samples, and
the histogram on the right only contains correctly classified samples.

####  4.1.2 Online Evaluation

For the online evaluation, we created a candidate pool of 30 days of
production data from September 2021. The size n of the candidate pool was
approximately one million requests. From this candidate set, we independently
applied two sampling strategies. The first sampling strategy was random
selection. The second sampling strategy was simple uncertainty selection,
which we described in Section 3.2.1. We set $k=10,000$. The number of
utterances we wanted to sample was an order of 100 times smaller than the
candidate pool, i.e., $k<<n$.

For each of these sampled sets, we then counted the number of sampled
utterances that were not part of a successful task. A dialog task is a higher-
level abstraction than an intent. It may contain multiple back-and-forth
interactions between the user and the voice assistant. Refer to Liu et al.
(2018) for more details on tasks and computing task success rate.

We did not grade the task success manually, but rather used an automated
metric for estimating task success that correlates closely with the true task
success rate. This automated metric checks whether the final NLG dialog act
was to inform the user by providing a response (as opposed to replying with an
error message). While this may over-estimate true task success, it serves as a
good proxy.

Returning to our comparison of the sampling strategies, we found that
uncertainty sampling yielded four times more task failures than random
sampling. The raw counts are given in Table 2.

Table 2: For each sampling strategy, we counted the number of sampled
utterances that were part of a task failure.
Sampling Strategy | # of Sampled Rows | # of Task Failures
---|---|---
Random | 10K | 1394
Simple Uncertainty | 10K | 5620

Note that an NLU prediction for a given utterance can be correct, even when
there is task failure; inversely, an NLU prediction for a given utterance can
be incorrect, despite a task success. So while these results suggest that
simple uncertainty will find more NLU bugs than random sampling in an online
setting, it does not guarantee it.

###  4.2 Data Augmentation

For our final experiments, we tested the effectiveness of automated data
augmentation at correcting bugs deemed to have the low training data error
attribution. We compared models trained with three different data augmentation
strategies:

  * •

Baseline An NLU model with no data augmentation applied.

  * •

Data Augmentation (Exact Match) For each bug, we added the exact utterance to
the training data. For example, for the bug, “Please call mom,” we also added,
“Please call mom” to training. Each training utterance was given a weight of 5
(copied 5 times).

  * •

Data Augmentation (Templated) For each bug, we generated 0 to 5 templates,
each of which was expanded up to 10 times. For example, for the bug, “Please
call mom,” we created the template Please call <SL:CONTACT>. This would then
be expanded to “Please call dad,” “Please call Jill,” etc. and added to
training.

We used a set of 1200 unique bug utterances, sampled from production data, as
seeds for data generation. We then evaluated the performance of the augmented
models against the original bug utterances, as well as our validation and test
set (which were randomly sampled from production data), and finally a crowd-
sourced test set. The results of the augmentation are shown in Table 3.

Table 3: We show the improvement that automated data augmentation techniques
have on various test sets.
| Unique Bugs | Production Validation | Production Test | Crowd-sourced Test
---|---|---|---|---
|  ($n=1207$) |  ($n=66253$) |  ($n=32854$) |  ($n=95107$)
Baseline | 5.5% | 96.1% | 93.9% | 90.9%
Exact Match | 95.8% | 98.0% | 95.4% | 90.4%
Templated | 70.6% | 97.5% | 95.0% | 90.0%
$\Delta$ Templated | +65.1% | +1.4% | +1.1% | -0.9%

Both the exact match and templated data augmentation strategies led to
significant improvements on the data sets shown. Using templated data
augmentation, we were able to automatically correct 65%, i.e., 789 out of
1207, of identified bugs. This was an order of magnitude higher than the
number of bugs we are able to correct in the prior three months using manual
data augmentation.

Interestingly, exact match performed even better than templated data
augmentation, which went against our hypothesis: we had expected templated
augmentation to lead to a better model, since the additional variety in the
training data could improve generalization. We believe the reason the
performance of templated data was slightly worse was due to noise in the
extracted templates and slot substitutions, which we discuss further in
Section 4.3.

Finally, another point worth mentioning is that while performance on the data
sets sampled from production traffic increase, performance on the crowd-
sourced data set actually dropped slightly. Digging into this, we found the
reason for the drop was because production data and crowd-sourced data
sometimes had the same or similar utterances annotated slightly differently.
Apart from the labeling discrepancy, there could also be a genuine domain
shift between our production and crowd-sourced data.

###  4.3 Qualitative Feedback

In addition to the evaluations above, we also received qualitative feedback
from linguists, who were the primary users of the AutoNLU user interface. We
share a few highlights of the constructive feedback:

Include humans-in-the-loop while generating data. Although the original goal
of AutoNLU was to automate as much of the quality improvement process as
possible, we found it was beneficial, and in many cases necessary, to have a
human i.e, a developer or linguist, review the generated data prior to landing
it in the training set. Templated data generation can produce nonsensical
utterances and sometimes even incorrect annotations. Until better data
generation strategies are developed, a human is needed to catch and remove
incorrect examples.

Make it easy to update annotations across data sets. Another suggestion we
received was to make it easier for linguists to update incorrect annotations
that are spread across multiple data sets. A change in the annotation of an
utterance in one data set should propagate to others.

Track bug fixes landed over time. Finally, linguists requested the ability to
track bugs that have been fixed in the past. This empowers linguists to
measure the impact of their work, and enables tracking of when a bug recurs
again.

##  5 Related Work

Previous work applying active learning to NLU focuses on improving performance
on a fixed test set with fewer training examples Duong et al. (2018); Sen &
Yilmaz (2020); Chen et al. (2019b); instead, we propose framing active
learning as a problem of continuous bug detection.

Valuable libraries and tools for model interpretability and error analysis
exist, some specifically tailored for NLP Kokhlikyan et al. (2020); Tenney et
al. (2020). However, the debugging information provided by these tools is
quite granular. The error classifications provided by AutoNLU are intended to
be higher-level and map to a concrete next-action that a linguist or developer
can take to resolve a production bug.

On the error correction side, there has been recent work in programmatic
labeling Ratner et al. (2017) and frameworks for adversarial data collection
and generation for NLP Kiela et al. (2021); Morris et al. (2020). However,
these stand-alone libraries are not tied to tools for model and training data
analysis. AutoNLU, in contrast, attempts to bring error attribution and
correction together.

Finally, in industry, there has recently been work by HuggingFace called
AutoNLP and work by DeepOpinion also called AutoNLU. Despite their similarity
in name, these libraries have a different purpose than the AutoNLU described
here: they allow users to create new NLP or NLU models in a no- or low-code
way. The AutoNLU we present here has the entirely different goal of automating
the quality improvement process of already-built NLU models.

##  6 Conclusion

##### Key Takeaways

We presented a system called AutoNLU, which aims to automate three steps of
the NLU quality improvement process: error detection, attribution, and
correction.

We showed that error detection can be done with active learning methods. A
simple least-confidence uncertainty baseline can be surprisingly effective at
finding failed tasks, even with an uncalibrated model, as long as the number
of utterances to be sampled $k$ is significantly smaller than the number of
candidates $n$ (in our experiments, $n$ was 100 times larger than $k$). In an
online test, we found four times as many failed tasks with this method as we
did with random sampling.

In addition, we presented some heuristics for attributing the cause of a model
error and strategies for fixing these errors. We showed that templated data
generation, where entities inside of slots are substituted, is an effective
and scalable method for fixing NLU bugs. However, a human-review step may
still be required to inspect the generated data.

Using this method, we were able to auto-fix nearly 800 bugs from production
traffic, a scale ten-fold of what we were able to achieve when doing data
augmentation by hand. Despite these gains, we acknowledge that there are
limitations of automated data generation, such as generating nonsensical
utterances, that will need to be addressed before fully removing humans from
the loop.

##### Future Directions

We believe there are many interesting avenues to improve AutoNLU. On the
active learning front, we hope to test both new intrinsic and extrinsic
sampling strategies to see if they yield an even higher rate of NLU bugs and
failed tasks. Desai & Aly (2021) showed that an extrinsic, SVM model can
accurately predict when a Transformer model will make a mistake. This SVM
could be trained to incorporate features only available at runtime, such as
whether the user had to repeat their request. Another aspect to think about is
the quality of bugs selected to fix. BatchBALD Kirsch et al. (2019) may help
eliminate sampling of redundant bugs and increase their diversity, for
example.

Error attribution could be improved by providing more fine-grained
classifications that rely less on heuristics on the training data but rather a
more nuanced view of the model’s uncertainty. For example, being able to
distinguish between epistemic uncertainty (uncertainty due to a fundamental
lack of knowledge) and aleatoric uncertainty (uncertainty due to probabilities
being distributed across multiple plausible outcomes) Ülkümen et al. (2016)
would also enable us to better decide if a bug is caused because of a lack of
training data or labeling conflicts.

Finally, on the data augmentation front, we are excited about the potential of
using language models to generate data Brown et al. (2020), which may produce
more natural new training examples than the data generated by a rule-based
template system.

##  7 Acknowledgements

We would like to thank Geoffrey Goh for building the data augmentation library
described in Section 3.4.1. We would also like to thank Rongrong Qiu, George
Guliashvili, TJ Singh, Cherry Wang, Shicong Zhao, and Nancy Li for their work
and guidance on building the interface shown in Figure 3.

We would like to thank Shrey Desai, Michael Marlen, Maria Sumner, John
Bourassa, Julien Odent, Bing Liu, Ousmane Dia, Ted Wild, Jeffrey Kyle Chiu,
Zoe Papakipos, Alexander Zotov, Sonal Gupta, Anuj Kumar, Jinsong Yu, and Luna
Dong for their invaluable early discussions, feedback, and support of this
work.

Finally, we would also like to thank Justin Rill, Brian Moran, Claire Lesage,
Tammy Stark, Caitlin Lohman, Jeremy Kahn, and Safiyyah Saleem for being
patient dogfooders of the prototype and providing an invaluable linguist
perspective on this work.

## References

  * Breck et al. (2019) Breck, E., Polyzotis, N., Roy, S., Whang, S. E., and Zinkevich, M. A.  Data validation for machine learning.  In _MLSys_ , 2019.
  * Brown et al. (2020) Brown, T. B., Mann, B., Ryder, N., Subbiah, M., Kaplan, J., Dhariwal, P., Neelakantan, A., Shyam, P., Sastry, G., Askell, A., Agarwal, S., Herbert-Voss, A., Krueger, G., Henighan, T. J., Child, R., Ramesh, A., Ziegler, D. M., Wu, J., Winter, C., Hesse, C., Chen, M., Sigler, E., Litwin, M., Gray, S., Chess, B., Clark, J., Berner, C., McCandlish, S., Radford, A., Sutskever, I., and Amodei, D.  Language models are few-shot learners.  _ArXiv_ , abs/2005.14165, 2020.
  * Chen et al. (2019a) Chen, Q., Zhuo, Z., and Wang, W.  Bert for joint intent classification and slot filling.  _ArXiv_ , abs/1902.10909, 2019a.
  * Chen et al. (2019b) Chen, X. C., Sagar, A., Kao, J. T., Li, T. Y., Klein, C., Pulman, S. G., Garg, A., and Williams, J.  Active learning for domain classification in a commercial spoken personal assistant.  In _INTERSPEECH_ , 2019b.
  * Desai & Aly (2021) Desai, S. and Aly, A.  Diagnosing transformers in task-oriented semantic parsing.  _ArXiv_ , abs/2105.13496, 2021.
  * Duong et al. (2018) Duong, L., Afshar, H., Estival, D., Pink, G., Cohen, P., and Johnson, M.  Active learning for deep semantic parsing.  In _Proceedings of the 56th Annual Meeting of the Association for Computational Linguistics (Volume 2: Short Papers)_ , pp. 43–48, Melbourne, Australia, July 2018. Association for Computational Linguistics.  doi: 10.18653/v1/P18-2008.  URL https://aclanthology.org/P18-2008.
  * Gupta et al. (2018) Gupta, S., Shah, R., Mohit, M., Kumar, A., and Lewis, M.  Semantic parsing for task oriented dialog using hierarchical representations.  In _EMNLP_ , 2018.
  * He et al. (2020) He, J., Zhang, X., Lei, S., Chen, Z., Chen, F., Alhamadani, A., Xiao, B., and Lu, C.-T.  Towards more accurate uncertainty estimation in text classification.  In _EMNLP_ , 2020.
  * Huang et al. (2015) Huang, Z., Xu, W., and Yu, K.  Bidirectional lstm-crf models for sequence tagging.  _ArXiv_ , abs/1508.01991, 2015.
  * Kiela et al. (2021) Kiela, D., Bartolo, M., Nie, Y., Kaushik, D., Geiger, A., Wu, Z., Vidgen, B., Prasad, G., Singh, A., Ringshia, P., Ma, Z., Thrush, T., Riedel, S., Waseem, Z., Stenetorp, P., Jia, R., Bansal, M., Potts, C., and Williams, A.  Dynabench: Rethinking benchmarking in nlp.  _ArXiv_ , abs/2104.14337, 2021.
  * Kirsch et al. (2019) Kirsch, A., van Amersfoort, J. R., and Gal, Y.  Batchbald: Efficient and diverse batch acquisition for deep bayesian active learning.  In _NeurIPS_ , 2019.
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., Granger, B., Bussonnier, M., Frederic, J., Kelley, K., Hamrick, J., Grout, J., Corlay, S., Ivanov, P., Avila, D., Abdalla, S., Willing, C., and development team, J.  Jupyter notebooks - a publishing format for reproducible computational workflows.  In Loizides, F. and Scmidt, B. (eds.), _Positioning and Power in Academic Publishing: Players, Agents and Agendas_ , pp. 87–90, Netherlands, 2016\. IOS Press.  URL https://eprints.soton.ac.uk/403913/.
  * Kokhlikyan et al. (2020) Kokhlikyan, N., Miglani, V., Martin, M., Wang, E., Alsallakh, B., Reynolds, J., Melnikov, A., Kliushkina, N., Araya, C., Yan, S., and Reblitz-Richardson, O.  Captum: A unified and generic model interpretability library for pytorch.  _ArXiv_ , abs/2009.07896, 2020.
  * Koçak et al. (2021) Koçak, M. A., Ramirez, D., Erkip, E., and Shasha, D.  Safepredict: A meta-algorithm for machine learning that uses refusals to guarantee correctness.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 43:663–678, 2021.
  * Liu et al. (2018) Liu, B., Tür, G., Hakkani-Tür, D. Z., Shah, P., and Heck, L.  Dialogue learning with human teaching and feedback in end-to-end trainable task-oriented dialogue systems.  In _NAACL_ , 2018.
  * Liu et al. (2019) Liu, Y., Ott, M., Goyal, N., Du, J., Joshi, M., Chen, D., Levy, O., Lewis, M., Zettlemoyer, L., and Stoyanov, V.  Roberta: A robustly optimized bert pretraining approach.  _ArXiv_ , abs/1907.11692, 2019.
  * Mesnil et al. (2013) Mesnil, G., He, X., Deng, L., and Bengio, Y.  Investigation of recurrent-neural-network architectures and learning methods for spoken language understanding.  In _INTERSPEECH_ , 2013.
  * Monarch (2021) Monarch, R.  _Human-in-the-Loop Machine Learning: Active learning and annotation for human-centered AI_.  Manning, 2021.  ISBN 9781638351030.  URL https://books.google.com/books?id=bNo2EAAAQBAJ.
  * Morris et al. (2020) Morris, J. X., Lifland, E., Yoo, J. Y., Grigsby, J., Jin, D., and Qi, Y.  Textattack: A framework for adversarial attacks, data augmentation, and adversarial training in nlp.  In _EMNLP_ , 2020.
  * Ratner et al. (2017) Ratner, A. J., Bach, S. H., Ehrenberg, H. R., Fries, J. A., Wu, S., and Ré, C.  Snorkel: Rapid training data creation with weak supervision.  _Proceedings of the VLDB Endowment. International Conference on Very Large Data Bases_ , 11 3:269–282, 2017.
  * Reback et al. (2021) Reback, J., jbrockmendel, McKinney, W., den Bossche, J. V., Augspurger, T., Cloud, P., Hawkins, S., gfyoung, Roeschke, M., Sinhrks, Klein, A., Petersen, T., Tratner, J., She, C., Ayd, W., Hoefler, P., Naveh, S., Garcia, M., Schendel, J., Hayden, A., Saxton, D., Darbyshire, J., Shadrach, R., Gorelli, M. E., Li, F., Jancauskas, V., McMaster, A., Zeitlin, M., Battiston, P., and Seabold, S.  pandas-dev/pandas: Pandas 1.3.3, September 2021.  URL https://doi.org/10.5281/zenodo.5501881.
  * Ribeiro et al. (2020) Ribeiro, M. T., Wu, T. S., Guestrin, C., and Singh, S.  Beyond accuracy: Behavioral testing of nlp models with checklist.  In _ACL_ , 2020.
  * Sen & Yilmaz (2020) Sen, P. and Yilmaz, E.  Uncertainty and traffic-aware active learning for semantic parsing.  In _Proceedings of the First Workshop on Interactive and Executable Semantic Parsing_ , pp. 12–17, Online, November 2020. Association for Computational Linguistics.  doi: 10.18653/v1/2020.intexsempar-1.2.  URL https://aclanthology.org/2020.intexsempar-1.2.
  * Shrivastava et al. (2021) Shrivastava, A., Chuang, P. I.-J., Babu, A., Desai, S., Arora, A., Zotov, A., and Aly, A.  Span pointer networks for non-autoregressive task-oriented semantic parsing.  _ArXiv_ , abs/2104.07275, 2021.
  * Tenney et al. (2020) Tenney, I., Wexler, J., Bastings, J., Bolukbasi, T., Coenen, A., Gehrmann, S., Jiang, E., Pushkarna, M., Radebaugh, C., Reif, E., and Yuan, A.  The language interpretability tool: Extensible, interactive visualizations and analysis for nlp models.  In _EMNLP_ , 2020.
  * Wes McKinney (2010) Wes McKinney.  Data Structures for Statistical Computing in Python.  In Stéfan van der Walt and Jarrod Millman (eds.), _Proceedings of the 9th Python in Science Conference_ , pp. 56 – 61, 2010.  doi: 10.25080/Majora-92bf1922-00a.
  * Ülkümen et al. (2016) Ülkümen, G., Fox, C., and Malle, B.  Two dimensions of subjective uncertainty: Clues from natural language.  _Journal of Experimental Psychology: General_ , 145, 07 2016.  doi: 10.1037/xge0000202.


#  Numerical variational simulations of quantum phase transitions in the sub-
Ohmic spin-boson model with multiple polaron ansatz

Yulong Shen and Nengji Zhou<EMAIL_ADDRESS>Department of Physics,
Hangzhou Normal University, Hangzhou 311121, China

###### Abstract

With extensive variational simulations, dissipative quantum phase transitions
in the sub-Ohmic spin-boson model are numerically studied in a dense limit of
environmental modes. By employing a generalized trial wave function composed
of coherent-state expansions, transition points and critical exponents are
accurately determined for various spectral exponents, demonstrating excellent
agreement with those obtained by other sophisticated numerical techniques.
Besides, the quantum-to-classical correspondence is fully confirmed over the
entire sub-Ohmic range, compared with theoretical predictions of the long-
range Ising model. Mean-field and non-mean-field critical behaviors are found
in the deep and shallow sub-Ohmic regimes, respectively, and distinct physical
mechanisms of them are uncovered.

##  I Introduction

In recent years, there has been significant theoretical and experimental
interest in quantum phase transitions occurring at absolute zero temperature,
which are driven by quantum fluctuations that can be tuned by the external
magnetic field or the coupling constant Vojta (2003); Le Hur (2010); Sachdev
(2011). At low temperatures, they can still be detected if the system is
dominated by quantum fluctuations Sondhi  _et al._ (1997); Paiva  _et al._
(2005); Werlang  _et al._ (2010). Most of the works on quantum phase
transitions have focused on closed quantum systems. However, it is also
important to consider the impact of inherent system-environment coupling on
physical features such as decoherence, dissipation, and entanglement Weiss
(2007). As a result, quantum critical phenomena in dissipation systems have
recently received considerable attention Orth  _et al._ (2010); Banchi  _et
al._ (2014); Maile  _et al._ (2018); Rossini and Vicari (2021). The spin-boson
model (SBM), as schematically depicted in the top panel of Fig. 1, is a
prominent example featuring a two-level system (quantum spin) interacting with
harmonic oscillators Leggett  _et al._ (1987); Le Hur (2008); Breuer  _et al._
(2016). Despite its apparent simplicity, the SBM captures the physics of a
diverse range of systems, from defects in solids to quantum thermodynamics,
physical chemistry, and biological systems Lewis and Raggio (1988); Golding
_et al._ (1992); Chakravarty and Rudnick (1995); Engel  _et al._ (2007);
Porras  _et al._ (2008); Collini and Scholes (2009); Ota  _et al._ (2011);
Uzdin  _et al._ (2015). Recent studies have shown that the SBM and its
variants exhibit rich ground-state and dynamic phase transitions with the
increasing system-environment coupling Guo  _et al._ (2012); Nalbach and
Thorwart (2013); Zhou  _et al._ (2014); Wang  _et al._ (2016); Zhou  _et al._
(2018); Wang  _et al._ (2020); Perroni  _et al._ (2023).

Gapless bosonic baths characterized by a power-law spectral density
$J(\omega)\propto\omega^{s}$ can be divided into three distinct regimes
depending on the value of the spectral exponent $s$: sub-Ohmic ($s<1$), Ohmic
($s=1$), and super-Ohmic ($s>1$) Leggett  _et al._ (1987). Due to the
competition between the tunneling and environment-induced dissipation, the SBM
undergoes a quantum phase transition which separates a non-degenerate
delocalized phase from a doubly degenerate localized phase Kehrein and Mielke
(1996); Le Hur (2008), as illustrated in the bottom panel of Fig. 1. The phase
transition is of second order in the sub-Ohmic regime and takes on the
Kosterlitz-Thouless type in the Ohmic regime. In the super-Ohmic regime,
however, there is no phase transition. In addition, the quantum-classical
mapping predicts that the sub-Ohmic SBM should be equivalent to a classical
Ising spin chain with long-range interactions, and mean-field critical
exponents are expected in the deep sub-Ohmic regime with $0<s\leq 0.5$ Leggett
_et al._ (1987); Bulla  _et al._ (2003); Orth  _et al._ (2010); Guo  _et al._
(2012).

Although an exact analytical solution is currently lacking, various powerful
numerical approaches exist to investigate ground-state properties of the SBM.
Among them, the most famous one is the Numerical Renormalization Group (NRG)
developed by Kenneth Wilson in the $1970$s Wilson (1971). Failure of the
quantum-to-classical correspondence was reported in NRG works Bulla  _et al._
(2003); Vojta  _et al._ (2005), since the results of critical exponents did
not align with the mean-field values. Conversely, the mapping was confirmed by
other numerical approaches, such as Quantum Monte Carlo (QMC) Winter  _et al._
(2009), Density-Matrix Renormalization Group (DMRG) Wong and Chen (2008),
Exact Diagonalization (ED) Alvermann and Fehske (2009), and Variational Matrix
Product States (VMPS) Guo  _et al._ (2012). The controversy was subsequently
addressed for the realization that the earlier NRG results were incorrect due
to bosonic state space truncation Tong and Hou (2012). Besides, the ground-
state phase transition can also be investigated with quantum dissipation
dynamic approaches Zhou and Shao (2008); Yan and Shao (2016); Duan  _et al._
(2017); Wang and Shao (2019).

In the shallow sub-Ohmic regime with $0.5<s<1$, $s$-dependent values of
critical exponents were obtained going beyond the mean-field theory Vojta  _et
al._ (2005); Winter  _et al._ (2009); Blunden-Codd  _et al._ (2017); Bruognolo
_et al._ (2017). It suggests that further examination on the validity of the
quantum-to-classical correspondence is required. However, few research has
been conducted except for that carried out by VMPS Guo  _et al._ (2012).
Besides, numerical results in such a regime show significant variations in the
critical couplings. For instance, the NRG one is greater than others by nearly
$10$ percent at $s=0.9$ Alvermann and Fehske (2009); Winter  _et al._ (2009);
Zhang  _et al._ (2010). Very recently, quantum simulations on the SBM have
been realized in experiments of superconducting quantum circuits and trapped
atomic ion crystals Porras  _et al._ (2008); Leppäkangas  _et al._ (2018);
Magazzù  _et al._ (2018); Lemmer  _et al._ (2018), offering the chance of
experimental studies on quantum phase transitions. However, current
understanding of quantum criticality in the SBM remains limited, particularly
in the context of the shallow sub-Ohmic regime.

The variational technique is another popular method of obtaining approximate
energies and wave functions in quantum mechanical systems He (1999); Sorella
(2005), where the form of the trial wave function is crucial. The original
variational work on the SBM utilized the polaronic unitary transformation
proposed by Silbey and Harris Silbey and Harris (1984). This ansatz was later
adapted into an asymmetric form by Chin et al. to investigate quantum phase
transitions in the sub-Ohmic SBM Chin  _et al._ (2011). Mean-field critical
properties were found in the deep sub-Ohmic regime, consistent with those
obtained from other numerical methods. However, the investigation of cases
within the shallow sub-Ohmic regime was unexplored. With the ansatz modified
by superposing more than one nonorthogonal coherent states, the mean-field
value of the magnetization critical exponent $\beta=0.50(1)$ was confirmed
recently not only for $0<s\leq 0.5$, but also for $0.5<s<1$ He  _et al._
(2018). But, unexpectedly, the $\beta$ values in the latter differ greatly
from VMPS ones which are much less than $0.5$ Guo  _et al._ (2012); Frenzel
and Plenio (2013), thus casting doubt on the ability of the variational method
to treat quantum phase transitions across the entire sub-Ohmic regime.

In this article, extensive variational calculations with thousands of
parameters are carried out to study quantum phase transitions in both deep and
shallow sub-Ohmic SBM based on the generalized trial wave function composed of
coherent-state expansions, which has been proved to be valid in tackling
ground-state and dynamic properties Zhou  _et al._ (2014, 2015a, 2015b, 2016,
2018); Qian  _et al._ (2021, 2022). In the continuum limit, transition points
and critical exponents are accurately determined, and the variational results
are examined in detail, in comparison with those obtained from other numerical
approaches. Moreover, the quantum-to-classical correspondence is confirmed in
the whole sub-Ohmic regime, and the nature of the mean-field and non-mean-
field behaviors as well as the finite-size effect are also analyzed. The rest
of the paper is organized as follows. In Sec. II, the model, variational
approach and scaling analysis are described, and in Sec. III, the numerical
results are presented. Finally, Sec. IV includes the conclusion.

##  II Methodology

###  II.1 Model and Method

Numerical studies of quantum phase transitions are performed in the SBM whose
Hamiltonian can be described by

$\hat{H}=\frac{\varepsilon}{2}\sigma_{z}-\frac{\Delta}{2}\sigma_{x}+\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k}+\frac{\sigma_{z}}{2}\sum_{k}\lambda_{k}(b^{\dagger}_{k}+b_{k}),$
(1)

where $\varepsilon$ ($\Delta$) denotes the energy bias (tunneling amplitude),
$b^{\dagger}_{k}$ ($b_{k}$) represents the bosonic creation (annihilation)
operator of the $k$-th bath mode with the frequency $\omega_{k}$, $\sigma_{x}$
and $\sigma_{z}$ are Pauli spin-$1/2$ operators, and $\lambda_{k}$ signifies
the coupling amplitude between the system and bath. To simplify notation, we
fix the Planck’s constant $\hbar=1$, so model parameters
$\Delta,\varepsilon,\lambda_{k}$, and $\omega_{k}$ are dimensionless. With the
coarse-grained treatment based on the Wilson energy mesh, the values of
$\lambda_{k}$ and $\omega_{k}$ in Eq. (1) can be calculated from the
continuous spectral density function
$J(\omega)=2\alpha\omega_{c}^{1-s}\omega^{s}=\sum_{k}\lambda_{k}^{2}\delta(\omega-\omega_{k})$
Bulla  _et al._ (2005); Vojta  _et al._ (2005); Zhang  _et al._ (2010); Zhou
_et al._ (2014); Blunden-Codd  _et al._ (2017), where $\alpha$ denotes the
dimensionless coupling strength, and $\omega_{c}=1$ represents the high-
frequency cutoff. In order to obtain an accurate description of the ground
state for SBM in a high dense spectrum, the factor $\Lambda=1.05$ is set in
the logarithmic discretization procedure. It is much closer to the continuum
limit $\Lambda\rightarrow 1$, compared to those in previous numerical works
Bulla  _et al._ (2003); Le Hur (2008); Chin  _et al._ (2011); Frenzel and
Plenio (2013); Bera  _et al._ (2014a). Furthermore, the convergency test of
the logarithmic discretization factor is also performed, and the results
demonstrate that a value of $\Lambda=1.05$ is sufficient to obtain reliable
outcomes Qian  _et al._ (2022).

The multiple polaron ansatz Zhou  _et al._ (2014), which is also known as
Davydov multi-D1 ansatz, is used in the studies with the numerical variational
method (NVM),

$\displaystyle|\Psi_{g}\rangle$ $\displaystyle=$
$\displaystyle|\uparrow\rangle\sum_{n=1}^{N}A_{n}\exp\left[\sum_{k=1}^{M}\left(f_{n,k}b_{k}^{{\dagger}}-\mbox{H}.\mbox{c}.\right)\right]|0\rangle_{\textrm{b}}$
(2) $\displaystyle+$
$\displaystyle|\downarrow\rangle\sum_{n=1}^{N}B_{n}\exp\left[\sum_{k=1}^{M}\left(g_{n,k}b_{k}^{{\dagger}}-\mbox{H}.\mbox{c}.\right)\right]|0\rangle_{\textrm{b}},$

where H.c$.\\!$ denotes Hermitian conjugate, $|\uparrow\rangle$
($|\downarrow\rangle$) stands for the spin up (down) state, and
$|0\rangle_{\rm b}$ is the vacuum state of the bosonic bath. The variational
parameters $f_{n,k}$ and $g_{n,k}$ represent displacements of the coherent
states correlated to spin configurations, and $A_{n}$ and $B_{n}$ are weights
of coherent states. The subscripts $n$ and $k$ correspond to the ranks of the
coherent superposition state and effective bath mode, respectively.

The ground state $|\Psi_{g}\rangle$ is obtained by minimizing the energy
expressed as $E=\mathcal{H}/\mathcal{N}$, using the Hamiltonian expectation
$\mathcal{H}=\langle\Psi_{g}|\hat{H}|\Psi_{g}\rangle$ and the norm of the wave
function $\mathcal{N}=\langle\Psi_{g}|\Psi_{g}\rangle$,

$\displaystyle\mathcal{H}=$
$\displaystyle\sum_{m,n}A_{m}B_{n}\Gamma_{m,n}(-\Delta)+\frac{\varepsilon}{2}\sum_{m,n}\left(A_{m}A_{n}F_{m,n}-B_{m}B_{n}G_{m,n}\right)$
$\displaystyle+$
$\displaystyle\sum_{m,n}A_{m}A_{n}F_{m,n}\sum_{k}\left[\omega_{k}f_{m,k}f_{n,k}+\frac{\lambda_{k}}{2}(f_{m,k}+f_{n,k})\right]$
(3) $\displaystyle+$
$\displaystyle\sum_{m,n}B_{m}B_{n}G_{m,n}\sum_{k}\left[\omega_{k}g_{m,k}g_{n,k}-\frac{\lambda_{k}}{2}(g_{m,k}+g_{n,k})\right],$

and

$\mathcal{N}=\sum_{m,n}\left(A_{m}A_{n}F_{m,n}+B_{m}B_{n}G_{m,n}\right),$ (4)

where $F_{m,n},G_{m,n}$, and $\Gamma_{m,n}$ are Debye-Waller factors defined
as

$\displaystyle F_{m,n}$ $\displaystyle=$
$\displaystyle\exp\left[-\frac{1}{2}\sum_{k}(f_{m,k}-f_{n,k})^{2}\right],$
$\displaystyle G_{m,n}$ $\displaystyle=$
$\displaystyle\exp\left[-\frac{1}{2}\sum_{k}(g_{m,k}-g_{n,k})^{2}\right],$ (5)
$\displaystyle\Gamma_{m,n}$ $\displaystyle=$
$\displaystyle\exp\left[-\frac{1}{2}\sum_{k}(f_{m,k}-g_{n,k})^{2}\right].$

A set of self-consistency equations are then deduced by the Lagrange
multiplier method,

$\frac{\partial\mathcal{H}}{\partial
x_{i}}-E\frac{\partial\mathcal{N}}{\partial x_{i}}=0,$ (6)

with respect to the variational parameter $x_{i}(i=1,2,3,\cdots,2NM+2N)$.
Finally, the iterative equations are derived,

$\displaystyle A_{n}^{*}$ $\displaystyle=$
$\displaystyle\frac{\sum_{m}B_{m}\Gamma_{n,m}(-\Delta)+2\sum_{m}^{m\neq
n}A_{m}F_{m,n}(aa_{m,n}-E)}{2(E-aa_{n,n})},$ $\displaystyle B_{n}^{*}$
$\displaystyle=$
$\displaystyle\frac{\sum_{m}A_{m}\Gamma_{m,n}(-\Delta)+2\sum_{m}^{m\neq
n}B_{m}G_{m,n}(bb_{m,n}-E)}{2(E-bb_{n,n})},$ $\displaystyle f_{m,k}^{*}$
$\displaystyle=$ $\displaystyle\frac{2\sum_{n}^{n\neq
m}A_{n}F_{m,n}(\omega_{k}f_{n,k}+\lambda_{k}/2+aa_{m,n}f_{n,k}-Ef_{n,k})}{2A_{m}(E-\omega_{k}-aa_{m,m})}$
(7) $\displaystyle+$
$\displaystyle\frac{\sum_{n}B_{n}\Gamma_{m,n}g_{n,k}(-\Delta)+A_{m}\lambda_{k}}{2A_{m}(E-\omega_{k}-aa_{m,m})},$
$\displaystyle g_{m,k}^{*}$ $\displaystyle=$
$\displaystyle\frac{2\sum_{n}^{n\neq
m}B_{n}G_{m,n}(\omega_{k}g_{n,k}-\lambda_{k}/2+bb_{m,n}g_{n,k}-Eg_{n,k})}{2B_{m}(E-\omega_{k}-bb_{m,m})}$
$\displaystyle+$
$\displaystyle\frac{\sum_{n}A_{n}\Gamma_{n,m}f_{n,k}(-\Delta)-B_{m}\lambda_{k}}{2B_{m}(E-\omega_{k}-bb_{m,m})},$

where $aa_{m,n}$ and $bb_{m,n}$ denote

$\displaystyle aa_{m,n}$ $\displaystyle=$
$\displaystyle\sum_{k}\left[\omega_{k}f_{m,k}f_{n,k}+\frac{\lambda_{k}}{2}(f_{m,k}+f_{n,k})\right]+\frac{\varepsilon}{2},$
$\displaystyle bb_{m,n}$ $\displaystyle=$
$\displaystyle\sum_{k}\left[\omega_{k}g_{m,k}g_{n,k}-\frac{\lambda_{k}}{2}(g_{m,k}+g_{n,k})\right]-\frac{\varepsilon}{2},$
(8)

respectively. Using the relaxation iteration technique, one updates the
variation parameter by $x_{i}^{\prime}=x_{i}+t*(x_{i}^{*}-x_{i})$ where
$x_{i}^{*}$ is defined in Eq. (II.1), and $t=0.1$ is the relaxation factor.
The termination criterion of the iteration procedure is set to be
max$\\{x_{i}^{*}-x_{i}\\}<10^{-13}$ over all variational parameters. More than
one hundred random initial states are taken to reduce statistical noise for
each set of model parameters ($\alpha,\Delta,\Lambda,\varepsilon$).
Furthermore, simulated annealing algorithm is also employed to escape from
metastable states.

With the ground-state wavefunction $|\Psi_{g}\rangle$ at hand, the spin
magnetization $\langle\sigma_{z}\rangle$ and spin coherence
$\langle\sigma_{x}\rangle$ can be measured by

$\displaystyle\langle\sigma_{z}\rangle$ $\displaystyle=$
$\displaystyle\frac{\sum_{m,n}\left(A_{m}A_{n}F_{m,n}-B_{m}B_{n}G_{m,n}\right)}{\sum_{m,n}\left(A_{m}A_{n}F_{m,n}+B_{m}B_{n}G_{m,n}\right)}.$
$\displaystyle\langle\sigma_{x}\rangle$ $\displaystyle=$
$\displaystyle\frac{\sum_{m,n}2A_{m}B_{n}\Gamma_{m,n}}{\sum_{m,n}A_{m}A_{n}F_{m,n}+B_{m}B_{n}G_{m,n}}.$
(9)

The local susceptibility $\chi$ is also investigated for a nonvanishing bias
$\varepsilon$,

$\chi(\alpha)=\left.\frac{\partial\langle\sigma_{z}\rangle}{\partial\varepsilon}=\frac{\delta\langle\sigma_{z}\rangle}{\delta\varepsilon}\right|_{\varepsilon=0},$
(10)

where the variation of the spin magnetization $\delta\langle\sigma_{z}\rangle$
is calculated under $\delta\varepsilon=10^{-5}\omega_{c}$.

Besides spin-related observations, a many-body quantum tomography technique
focused on the bosonic bath is also introduced, which allows direct
characterization of the ground-state wave function. It can be accessed via
standard Wigner tomography Lvovsky and Raymer (2009), which is a method used
to reconstruct the Wigner function of the quantum system. It allows for the
characterization of the complete quantum state. The basic principle is to
perform a series of measurements that directly probe the position and momentum
values of the quantum system. By repeating this measurement and reconstruction
process for different measurement settings, one can build up a comprehensive
picture of the quantum state’s Wigner function in the phase space. This
provides insights into the system’s quantum coherence, entanglement, and non-
classical behavior. In this work, the Wigner function $W^{k}(X,P)$ can be
calculated after all the degrees are traced out except the single bath mode
with the quantum number $k$. Thus, the Wigner function $W^{k}(\beta)$ is
introduced in the well-known definition Bera  _et al._ (2014a),

$\displaystyle W^{k}(\beta)$ $\displaystyle=$
$\displaystyle\frac{1}{\mathcal{N}}\int\frac{d^{2}\eta}{\pi^{2}}e^{-(\eta\beta^{*}-\eta^{*}\beta)}\mbox{Tr}\left[\hat{\rho}e^{\eta
b_{k}^{{\dagger}}-\eta^{*}b_{k}}\right]$ $\displaystyle=$
$\displaystyle\frac{1}{\mathcal{N}}\int\frac{d^{2}\eta}{\pi^{2}}e^{-(\eta\beta^{*}-\eta^{*}\beta)}\langle\Psi_{g}|e^{\eta
b_{k}^{{\dagger}}-\eta^{*}b_{k}}|\Psi_{g}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\mathcal{N}}\int\frac{d^{2}\eta}{\pi^{2}}\exp\left(\eta^{*}\beta-\eta\beta^{*}-\eta\eta^{*}/2\right)$
$\displaystyle\times$
$\displaystyle\sum_{m,n}\left[A_{m}A_{n}F_{mn}\exp\left(\eta
f_{m,l}-\eta^{*}f_{n,l}\right)\right.$ $\displaystyle+$
$\displaystyle\left.B_{m}B_{n}G_{mn}\exp\left(\eta
g_{m,l}-\eta^{*}g_{n,l}\right)\right],$

where $\beta=X+iY$ is a complex number. Taking the Gaussian integral over the
variable $\eta$ in the complex space as well as the variable $Y$, one can
obtain

$\displaystyle W^{k}(X)$ $\displaystyle=$
$\displaystyle\frac{\sqrt{2}}{\sqrt{\pi}\mathcal{N}}\sum_{m,n}A_{m}A_{n}\exp\left[-\frac{1}{2}\sum_{l\neq
k}(f_{m,l}-f_{n,l})^{2}\right]$ $\displaystyle\times$
$\displaystyle\exp\left[-2\left(X-(f_{m,k}+f_{n,k})/2\right)^{2}\right]$
$\displaystyle+$
$\displaystyle\frac{\sqrt{2}}{\sqrt{\pi}\mathcal{N}}\sum_{m,n}B_{m}B_{n}\exp\left[-\frac{1}{2}\sum_{l\neq
k}(g_{m,l}-g_{n,l})^{2}\right]$ $\displaystyle\times$
$\displaystyle\exp\left[-2\left(X-(g_{m,k}+g_{n,k})/2\right)^{2}\right].$

In fact, the function $W^{k}(X)$ represents the probability distribution of
the $k$-th bath mode in the position space, where the position operator is
$\hat{x}=(b_{k}+b_{k}^{\dagger})/\sqrt{2}$.

###  II.2 Implementation of the algorithm

In the following, the algorithm of the numerical variational method in this
paper is implemented with the architecture:

Step $1$: Initialize variational parameters randomly. Specifically,
variational parameters $(A_{n},B_{n})$ are uniformly distributed within an
interval $[-1,1]$, and displacement coefficients $(f_{m,k},g_{m,k})$ of the
initial states obey
$f_{m,k}=-g_{m,k}\sim\lambda_{k}/2(\omega_{k}+0.01\omega_{r})$ with a
uniformly distributed random frequency $\omega_{r}$.

Step $2$: Update variational parameters with the relaxation iteration
technique until the preliminary condition max$\\{x_{i}^{*}-x_{i}\\}<10^{-8}$
over all variational parameters is reached. The trail wavefunction and its
energy are then recorded.

Step $3$: Repeat the steps $1$ and $2$ for more than $100$ times, and
subsequently choose the trail wavefunction with the lowest energy as the
candidate for the ground state.

Step $4$: Carry on the iterative procedure in the candidate until the target
precision $10^{-13}$ is reached. If necessary, the simulated annealing
algorithm is also used where the relaxation factor gradually decreases to
$t=0.001$ in order to improve the energy minimization procedure.

Step $5$: Measure the ground-state properties with the spin magnetization
$\langle\sigma_{z}\rangle$, the spin coherence $\langle\sigma_{x}\rangle$, the
local susceptibility $\chi$, and the Wigner function $W^{k}(X)$ in Eqs.
(9)-(II.1).

The FORTRAN source codes that we developed, which allow for the treatment of
quantum phase transitions in the sub-Ohmic SBM, are provided as the
supplementary material. Moreover, the Central Processing Unit (CPU) time and
cost of the memory for each set of model parameters are also estimated, which
are about $2.5$ months on a single processor and $2$ megabytes of the memory.
The type of CPU is “Intel Xeon Gold $6248$R” with the frequency $3.00$GHz, and
we have $176$ logical processors in total on a Linux-based computing cluster.

###  II.3 Scaling analysis

Since the localized-delocalized phase transition in the sub-Ohmic SBM is of
second order, the order parameter $\langle\sigma_{z}\rangle$ and local
susceptibility $\chi$ should obey

$\langle\sigma_{z}\rangle\sim\tau^{\beta},\qquad\langle\sigma_{z}\rangle\sim|\varepsilon|^{1/\delta},\qquad\chi\sim\tau^{-\gamma},$
(13)

where the reduced coupling $\tau=|\alpha-\alpha_{c}|/\alpha_{c}$ denotes the
distance from the critical coupling, and $\beta,\gamma$, and $\delta$ are
critical exponents. Analogous to a classical phase transition, the correlation
length in imaginary time diverges as a function of the distance,
$\xi\sim\tau^{-\nu}$, where $\nu$ is the correlation-length exponent. Thus,
the transition point can be determined by
$\delta\alpha_{c}(L)=\alpha_{c}(L)-\alpha_{c,\infty}\sim L^{-1/\nu}$ with the
finite size $L$ of the system. Note that the correlation length is given by an
inverse energy scale $\xi\sim 1/\omega^{*}$ where $\omega^{*}$ represents a
characteristic energy scale, above which the critical behavior is observed.
Hence, we choose the energy gap of the lowest frequency mode
$\hbar\omega_{\min}\approx\hbar\omega_{c}/\Lambda^{M}$ as the inverse of the
size $1/L$. For a sufficiently large number of bath modes, e.g., $M=430$, the
finite-size effect is already negligibly small.

According to the quantum-to-classical correspondence, the SBM can be mapped to
a one-dimensional classical Ising chain,

$\hat{H}=-\sum_{i,j}J_{ij}S_{i}S_{j}+\hat{H}_{\rm short-range},$ (14)

where $J_{ij}=J/r^{1+s}$ represents the long-range interaction with the
distance $r$ between two spins, $S_{i}=\pm 1$ denotes the classical Ising
spin, and $\hat{H}_{\rm short-range}$ is an additional generic short-range
interaction irrelevant to the critical behavior. In general, it is believed
that the quantum transition of SBM is in the same universality class as that
of the classical model Winter  _et al._ (2009); Chin  _et al._ (2011); Guo
_et al._ (2012). Further analysis gives the mean-field results:

$\beta=1/2,\qquad\gamma=1,\qquad\delta=3,\qquad\nu=1/s.$ (15)

In the non-mean-field regime ($0.5<s<1$), numerical results of the critical
exponents can be estimated from the analytical prediction by the two-loop
renormalization-group theory as well as the hyperscaling relation Fisher  _et
al._ (1972),

$\displaystyle\gamma$ $\displaystyle=$ $\displaystyle
1+\frac{4}{3}\epsilon+4.368\epsilon^{2}+\mathcal{O}(\epsilon^{3}),$
$\displaystyle 1/\nu$ $\displaystyle=$
$\displaystyle\frac{1}{2}+\frac{1}{3}\epsilon-2.628\epsilon^{2}+\mathcal{O}(\epsilon^{3})$
(16) $\displaystyle\beta$ $\displaystyle=$
$\displaystyle\gamma\frac{1-s}{2s},\qquad\delta=\frac{1+s}{1-s},$

with the expansion parameter $\epsilon=s-1/2$.

##  III Numerical results

With numerical variational calculations, quantum criticality of the sub-Ohmic
SBM in a high dense spectrum is investigated in this section for different
values of the spectral exponent $s=0.1,0.2,\cdots,0.9$ in the weak tunneling
taking the setting of the tunneling amplitude $\Delta=0.1$ as an example.
Theoretically, the number of effective bath modes $M\rightarrow\infty$ is
required for the completeness of the environment. Considering the constraint
available computational resources, a large value of $M=430$ is used in main
results. Similarly, the number of coherent-superposition states $N=8$ is set.
Both of them have been confirmed to be sufficient in accurately describe
ground states through the convergency test Zhou  _et al._ (2014); Qian  _et
al._ (2021, 2022). Thus, more than $8600$ variational parameters are used for
each set of the model parameters in calculations. The energy bias
$\varepsilon=0$ is set in the following unless otherwise noted. Statistical
errors of the critical couplings and exponents are estimated by dividing the
total samples into two subgroups. If the fluctuation in the curve is
comparable with or larger than the statistical error, it will be taken into
account.

###  III.1 Estimation of critical couplings and exponents

In Fig. 2(a), the ground-state spin magnetization $\langle\sigma_{z}\rangle$
in the shallow sub-Ohmic SBM is displayed as a function of the coupling
$\alpha$ for $s=0.6,0.7,0.8$, and $0.9$. Quantum phase transition is
demonstrated from the delocalized phase ($\langle\sigma_{z}\rangle=0$) to the
localized phase ($\langle\sigma_{z}\rangle\neq 0$). The transition point
$\alpha_{c}$ is then obtained, which increases with the spectral exponent $s$.
With $\alpha_{c}$ as input, a power-law behavior of $\langle\sigma_{z}\rangle$
is detected with respect to the shift $\alpha-\alpha_{c}$, as shown in Fig.
2(b). In the inset, the fitting error, as defined by the residual sum of
squares over the degree of freedom ($RSS/DoF$), is carefully examined in a
narrow regime of the coupling $\alpha$, taking the case $s=0.9$ as an example.
It provides a measure of the average error or variability per degree of
freedom. A lower value of $RSS/DoF$ indicates a better fit. By judging the
location of the minimum fitting error, one determines the critical coupling
$\alpha_{c}=0.534(1)$ which is more accurate than before. It agrees well with
numerical results $0.548(2)$ and $0.555$ obtained from QMC and VMPS works
Winter  _et al._ (2009); Guo  _et al._ (2012), respectively, showing the
validity of NVM. The critical exponent $\beta$ can also be measured with Eq.
(13). For each case, the value of $\beta$ is significantly less than $1/2$,
suggesting that the mean-field prediction is violated in the shallow sub-Ohmic
SBM. It is in contrast to the claim obtained from the previous variational
work He  _et al._ (2018).

For comparison, the spin magnetization $\langle\sigma_{z}\rangle$ in the deep
sub-Ohmic SBM is plotted in Fig. 3 for different spectral exponent $s\leq 0.5$
on a log-log scale. The power-law growing curves are almost parallel to each
other, quite different from those in Fig. 2(b). The slopes
$\beta=0.489(9),0.496(7),0.496(5),0.495(7)$, and $0.491(6)$ are then measured,
corresponding to $s=0.1,0.2,0.3,0.4$, and $0.5$, respectively. All of them are
in good agreement with the mean-field prediction in Eq. (15). It confirms the
robustness of the mean-field nature for quantum phase transitions in the deep
sub-Ohmic regime $s<0.5$, as it is reported in literatures Winter  _et al._
(2009); Chin  _et al._ (2011); Guo  _et al._ (2012); Frenzel and Plenio
(2013); He  _et al._ (2018).

In Fig. 4, the influence of the bath-mode number $M$ related to the size of
the system $L$ is investigated, taking $M=350$ for the cases of $s=0.6$ and
$0.7$, and $M=300$ for the cases of $s=0.8$ and $0.9$. The values of the
critical exponent $\beta$ are measured from the slopes at different $s$,
comparable with those in Fig. 2(b). It confirms the finite-size effect induced
by the number $M$ is already immaterial in main results. In addition,
calculations with a lager discretization factor $\Lambda=2$ are also
performed, which was commonly employed in earlier works Bulla  _et al._
(2005); Le Hur (2008); Bera  _et al._ (2014b). Even for the case of $s=0.7$,
the magnetization grows as a power law with a mean-field exponent
$\beta\approx 0.5$ over two decades in the inset of Fig. 4, quite different
from $\beta=0.128(1)$ in Fig. 2(b). The overestimation of $\beta$ is believed
to be caused by the artificial mean-filed nature. It suggests that numerical
works with a large $\Lambda$ may lead to a poor approximation of the ground
state. This is the possible reason of the failure in numerical variational
work He  _et al._ (2018). Hence, the continuum limit $\Lambda\rightarrow 1$
should be taken to obtain an accurate description of quantum transitions.

Subsequently, the extension to the biased SBM is performed for further
studies. In the delocalized phase, the local susceptibility $\chi(\alpha)$
defined in Eq. (10) is plotted in Fig. 5 for different $s$ under a tiny bias
$\delta\varepsilon=10^{-5}\omega_{c}$. With the power-law scaling in Eq. (13),
the measurement of the slope yields the value of the critical exponent ranging
from $\gamma=1.12(1)$ to $2.11(1)$ in the shallow sub-Ohmic regime and
$\gamma\approx 1$ in the deep sub-Ohmic regime, respectively. The latter is in
good agreement with the mean-field prediction in Eq. (15), thereby again
lending support to the reliability of the variational results.

In addition, the response of the external field (bias) at the transition point
$\alpha_{c}$ is investigated in Fig. 6 for different values of $s$. Non-linear
increases of the magnetization $\langle\sigma_{z}\rangle$ are found with
respect to the applied field $\varepsilon$ in both shallow and deep sub-Ohmic
regimes, signifying a second-order phase transition. The critical exponent
$1/\delta$ is then determined from simply fitting the slope of the curve. In
the shallow sub-Ohmic regime, the value of the exponent decreases rapidly with
$s$, and approaches zero in the Ohmic case. The latter indicates that the
magnetization should be always unity over the whole coupling ranging of the
localized phase, consistent with previous reports Le Hur (2008); Nazir  _et
al._ (2012); Zhou  _et al._ (2018). In the deep sub-Ohmic one, however, the
measurements of $1/\delta$ coincide with each other, showing the independence
of $s$.

All the results of the transition point $\alpha_{c}$ for different spectral
exponent $s$ are summarized in Table. 1, in comparison with those in the
literature. Based on the extension of the Silbey-Harris polaron ansatz as well
as the Ginzburg-Landau theory, the theoretical value of the critical coupling
in the deep sub-Ohmic SBM has been solved analytically Chin  _et al._ (2011),

$\alpha_{c}=\frac{\sin(\pi
s)e^{-s/2}}{2\pi(1-s)}\left(\frac{\Delta}{\omega_{c}}\right)^{1-s}.$ (17)

The numerical results are then estimated, as shown in the penultimate column
of the table. However, they are significantly smaller than those obtained from
general numerical methods, i.e., NRG, QMC, VMPS, and DMRG Bulla  _et al._
(2003); Vojta  _et al._ (2005); Winter  _et al._ (2009); Guo  _et al._ (2012);
Wong and Chen (2008). It indicates that the single polaron ansatz employed in
this analytical work seems not sophisticated enough to capture full bath
quantum fluctuations which play a crucial role in quantum phase transitions.
Moreover, the analytical solution of $\alpha_{c}$ is still missing in the
shallow sub-Ohmic SBM. In this work, variational results of the transition
point $\alpha_{c}$ are significantly improved by employing the NVM based on
systematic coherent-state decomposition in the continuum limit. The critical
couplings in both deep and shallow sub-Ohmic regimes are determined
accurately, as shown in the last column. The results accord well with those
obtained by QMC, VMPS, and DMRG, although having a notable deviation from the
NRG findings in the shallow sub-Ohmic regime. Therefore, our numerical method,
NVM, demonstrates remarkable superiority, the same as QMC, VMPS, and DMRG.
Moreover, the physical interpretation of our method is notably clearer and
simpler, adding to its advantage over other approaches.

Finally, the measurements of critical exponents $\beta,\gamma$, and $1/\delta$
are summarized in Table. 2, in comparison with those of recent VMPS work Guo
_et al._ (2012) and theoretical predictions in Eqs. (15) and (II.3). In the
deep sub-Ohmic regime, the exponent remains almost unchanged with respect to
the spectral exponent $s$. Compared to the VMPS ones, our NVM results are
closer to the mean-field predictions $\beta=1/2,\gamma=1$, and $1/\delta=1/3$.
In the shallow sub-Ohmic regime, both $\beta$ and $1/\delta$ monotonically
decrease with $s$, showing the opposite tendency of $\gamma$. All of them are
in good agreement with theoretical results obtained from the two-loop
renormalization-group analysis, as shown in the left columns. It provides a
strong confirmation on the validity of the quantum-to-classical correspondence
in the entire sub-Ohmic regime. In addition, our results are also comparable
with VMPS ones, again pointing to the superiority of NVM in studying quantum
phase transitions.

###  III.2 Other observables related to the spin and bath

For a better understanding of quantum phase transitions, the spin coherence
$\langle\sigma_{x}\rangle$ is also investigated in Fig. 7 as a function of the
coupling $\alpha$ for different values of $s$. Two lines with open triangles
and open circles are demonstrated, corresponding to the deep sub-Ohmic case of
$s=0.1$ and shallow sub-Ohmic case of $s=0.9$, respectively. In the former,
distinct values of the slopes $t_{1}$ and $t_{2}$ are measured from two sides
of the transition point $\alpha_{c}=0.00795$, indicating a discontinuity in
the slope of the first derivative of the ground-state energy
$\langle\sigma_{x}\rangle=\partial E_{g}/\partial\Delta$. Therefore, the
second-order phase transition is clearly evidenced in the deep sub-Ohmic SBM.
In the latter, however, the spin coherence $\langle\sigma_{x}\rangle$
decreases smoothly, and a slight kink occurs at the transition point
$\alpha_{c}=0.534$. Specifically, the slope difference $\delta t=t_{2}-t_{1}$
at $s=0.9$ is approximately three orders of magnitude less than that at
$s=0.1$. It means that the sharp transition is significantly softened,
yielding the non-mean-field critical behavior. Furthermore, an exponential
decay of $\delta t(s)$ is demonstrated in the inset with the slope $7.8$. The
vanishing difference $\delta t=0$ is then inferred in the limit $s\rightarrow
1$, supporting the argument that the spin coherence $\langle\sigma_{x}\rangle$
should be continuous at $\alpha_{c}$ in the presence of an Ohmic bath Bera
_et al._ (2014b).

Besides, ground-state properties of the bosonic bath in the vicinity of the
transition point are analyzed, as depicted by the function $W^{k}(X)$ in Eq.
(II.1). Two cases of $s=0.2$ and $0.7$ are considered for instance,
corresponding to the deep and shallow sub-Ohmic SBM, respectively. In the Fig.
8(a), a significant increase in the mean position $\overline{X}=\int
XW^{k}(X)dX$ is observed from the delocalized phase with $\overline{X}=0$ to
the localized phase with $6.29$, as the coupling strength $\alpha$ is changed
by only a paltry amount of $0.00005$. The shift in the mean position
$\delta\overline{X}$ is much larger than the standard deviation $\Delta
X\approx 0.5$ which represents the quantum fluctuation. It indicates the
robustness and mean-field effect of the phase transition Winter  _et al._
(2009); Chin  _et al._ (2011), which is caused by the low-frequency
nonadiabatic bath modes whose displacements diverge with the same sign as
$\omega\rightarrow 0$. Additionally, two degenerate ground state
$\Psi_{\uparrow}(x)$ and $\Psi_{\downarrow}(x)$ are identified, each
possessing opposite displacements to the other, corresponding to the
magnetization $m=\pm\langle\sigma_{z}\rangle$ in Fig. 1.

In the shallow sub-Ohmic regime, however, the criticality of the bosonic bath
differs significantly from that in the deep one. The curves of Wigner
distribution $W^{k}(X)$ mostly overlap for $\alpha=0.230$ and $0.234$ on two
sides of the transition point, as shown in the subfigure (b). It highlights
that the quantum fluctuation holds a crucial role even for the bath mode with
the lowest frequency $\omega_{\rm min}$, resulting in a non-mean-field scaling
behavior. Furthermore, the spin magnetization $\langle\sigma_{z}\rangle=\pm
0.299>0$ is found in both the two-folding degenerate states, where the mean
position $\overline{X}=\pm 0.257$ is less than the quantum fluctuation $\Delta
X\approx 0.5$. This observation implies that the order parameter is highly
sensitive to the numerical errors in the ground state of the bosonic bath,
providing a great challenge to obtain an accurate description of the quantum
phase transition. For this reason, numerical results in earlier works Bulla
_et al._ (2003); Alvermann and Fehske (2009); Winter  _et al._ (2009); Zhang
_et al._ (2010); Guo  _et al._ (2012), including critical couplings and
exponents, exhibit considerable discrepancies when the spectral exponent is
close to $s=1$.

###  III.3 Finite-size scaling

Subsequently, the correlation-length exponent $\nu$ describing the divergence
of the correlation length $\xi$ is investigated by the finite-size scaling.
Specifically, the dependence of the critical coupling $\alpha_{c}$ on the bath
size $L=1/\omega_{\rm min}$ is investigated, and the results are shown in Fig.
9 where the discretization factor $\Lambda=2.0$ is set for convenience. Power-
law dependence of the transition shift, $\delta\alpha_{c}\sim(1/L)^{1/\nu}$,
are found on the bath size not only in the deep sub-Ohmic regime for
$s=0.2,0.3,0.4$ and $0.5$, but also in the shallow one for $s=0.6,0.7,0.8$ and
$0.9$ (not shown). A direct measurement of the slope gives the value of the
exponent $1/\nu$ which varies with the spectral exponent $s$.

As shown in the inset, the exponent $1/\nu$ exhibits a linear increase when
$s\leq 0.5$, in excellent agreement with the mean-field prediction
$1/\nu_{{}_{MF}}=s$. For $s>0.5$, however, it decreases with $s$ and vanishes
in the Ohmic case $s=1$, showing the divergence of the correlation-length
exponent $\nu$. Moreover, our results closely follow the predictions of the
two-loop renormalization-group method within the error bars, but are much leas
than those from the perturbative calculations $1/\nu_{{}_{PE}}=\sqrt{2(1-s)}$.
It indicates that the formula in the latter holds only when $s>0.9$, which is
much closer to the Ohmic case.

##  IV Conclusion

With large-scale numerical variational calculations, comprehensive studies
have been performed on quantum phase transitions in dissipative quantum
systems, taking the sub-Ohmic spin-boson model with a highly dense spectrum as
an example. In these variational calculations, a generalized trial wave
function composed of coherent-state expansions has been employed to
effectively capture quantum entanglements and quantum fluctuations in the bath
which play a crucial role in determining the quantum phase transition. Thus,
systematical investigations into the quantum criticality of both the spin and
bosonic bath have been carried out, leading to the accurate determination of
the transition point $\alpha_{c}$ as well as the critical exponents
$\beta,\nu,\gamma$, and $1/\delta$ for various spectral exponents
$s=0.1,0.2,\cdots,0.9$ with a fixed tunneling constant of $\Delta=0.1$.

In contrast to previous variational results based on the single polaron
ansatz, our results of the critical couplings demonstrate excellent agreement
with those obtained by the QMC, VMPS, and DMRG methods. Although a minor
discrepancy is observed when comparing with the NRG results in the shallow
sub-Ohmic regime. It thereby lends support to the validity of the variational
calculations in this work. Additionally, the quantum-to-classical
correspondence has been fully confirmed over the entire sub-Ohmic range, for
the consistency between the numerical results of the critical exponents in the
SBM and theoretical predictions of one-dimensional classical Ising spin chain
with long-range interactions which are derived from the two-loop
renormalization-group analysis and the hyperscaling relation. Finally, two
distinct influences of quantum fluctuations have been uncovered on the low-
frequency bath modes in the deep and shallow sub-Ohmic regimes, respectively,
which may be responsible for the corresponding manifestation of mean-field and
non-mean-field critical behaviors.

As the supplementary material, we supply the FORTRAN codes for the numerical
variational method, entitled “init-bath.f90” and “pre-bath.f90” corresponding
to Steps 1-3 and Steps 4-5, respectively, to help further development in this
field.

Acknowledgements: This work was supported by National Natural Science
Foundation of China under Grant Nos. $11875120$.

## References

  * Vojta (2003) M. Vojta, Reports on Progress in Physics 66, 2069 (2003).
  * Le Hur (2010) K. Le Hur, in  _Understanding Quantum Phase Transitions._ , edited by L. D. Carr (CRC Press, Boca Raton, 2010) Chap. 9, pp. 217–240.
  * Sachdev (2011) S. Sachdev,  _Quantum Phase Transitions_ , 2nd ed. (Cambridge University Press, Cambridge, England, 2011).
  * Sondhi  _et al._ (1997) S. L. Sondhi, S. M. Girvin, J. P. Carini, and D. Shahar, Rev. Mod. Phys. 69, 315 (1997).
  * Paiva  _et al._ (2005) T. Paiva, R. T. Scalettar, W. Zheng, R. R. P. Singh, and J. Oitmaa, Phys. Rev. B 72, 085123 (2005).
  * Werlang  _et al._ (2010) T. Werlang, C. Trippe, G. A. P. Ribeiro, and G. Rigolin, Phys. Rev. Lett. 105, 095702 (2010).
  * Weiss (2007) U. Weiss,  _Quantum Dissipative Systems_ , 3rd ed. (World Scientific, Singapore, 2007).
  * Orth  _et al._ (2010) P. P. Orth, D. Roosen, W. Hofstetter, and K. Le Hur, Phys. Rev. B 82, 144423 (2010).
  * Banchi  _et al._ (2014) L. Banchi, P. Giorda, and P. Zanardi, Phys. Rev. E 89, 022102 (2014).
  * Maile  _et al._ (2018) D. Maile, S. Andergassen, W. Belzig, and G. Rastelli, Phys. Rev. B 97, 155427 (2018).
  * Rossini and Vicari (2021) D. Rossini and E. Vicari, Physics Reports 936, 1 (2021).
  * Leggett  _et al._ (1987) A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. 59, 1 (1987).
  * Le Hur (2008) K. Le Hur, Annals of Physics 323, 2208 (2008).
  * Breuer  _et al._ (2016) H. P. Breuer, E. M. Laine, J. Piilo, and B. Vacchini, Rev. Mod. Phys. 88, 021002 (2016).
  * Lewis and Raggio (1988) J. T. Lewis and G. A. Raggio, J. Stat. Phys 50, 1201 (1988).
  * Golding  _et al._ (1992) B. Golding, N. M. Zimmerman, and S. N. Coppersmith, Phys. Rev. Lett. 68, 998 (1992).
  * Chakravarty and Rudnick (1995) S. Chakravarty and J. Rudnick, Phys. Rev. Lett. 75, 501 (1995).
  * Engel  _et al._ (2007) G. S. Engel, T. R. Calhoun, E. L. Read, T. K. Ahn, T. Mancal, Y. C. Cheng, R. E. Blankenship, and G. R. Fleming, Nature 446, 782 (2007).
  * Porras  _et al._ (2008) D. Porras, F. Marquardt, J. von Delft, and J. I. Cirac, Phys. Rev. A 78, 010101(R) (2008).
  * Collini and Scholes (2009) E. Collini and G. D. Scholes, Science 323, 369 (2009).
  * Ota  _et al._ (2011) Y. Ota, S. Iwamoto, N. Kumagai, and Y. Arakawa, Phys. Rev. Lett. 107, 233602 (2011).
  * Uzdin  _et al._ (2015) R. Uzdin, A. Levy, and R. Kosloff, Phys. Rev. X 5, 031044 (2015).
  * Guo  _et al._ (2012) C. Guo, A. Weichselbaum, J. von Delft, and M. Vojta, Phys. Rev. Lett. 108, 160401 (2012).
  * Nalbach and Thorwart (2013) P. Nalbach and M. Thorwart, Phys. Rev. B 87, 014116 (2013).
  * Zhou  _et al._ (2014) N. J. Zhou, L. P. Chen, Y. Zhao, D. Mozyrsky, V. Chernyak, and Y. Zhao, Phys. Rev. B 90, 155135 (2014).
  * Wang  _et al._ (2016) L. Wang, L. P. Chen, N. J. Zhou, and Y. Zhao, J. Chem. Phys. 144, 024101 (2016).
  * Zhou  _et al._ (2018) N. J. Zhou, Y. Y. Zhang, Z. G. Lü, and Y. Zhao, Annalen der Physik 530, 1800120 (2018).
  * Wang  _et al._ (2020) Y. Z. Wang, S. He, L. W. Duan, and Q. H. Chen, Phys. Rev. B 101, 155147 (2020).
  * Perroni  _et al._ (2023) C. A. Perroni, A. De Candia, V. Cataudella, R. Fazio, and G. De Filippis, Phys. Rev. B 107, L100302 (2023).
  * Kehrein and Mielke (1996) S. K. Kehrein and A. Mielke, Physics Letters A 219, 313 (1996).
  * Bulla  _et al._ (2003) R. Bulla, N. H. Tong, and M. Vojta, Phys. Rev. Lett. 91, 170601 (2003).
  * Wilson (1971) K. G. Wilson, Phys. Rev. B 4, 3174 (1971).
  * Vojta  _et al._ (2005) M. Vojta, N. H. Tong, and R. Bulla, Phys. Rev. Lett. 94, 070604 (2005).
  * Winter  _et al._ (2009) A. Winter, H. Rieger, M. Vojta, and R. Bulla, Phys. Rev. Lett. 102, 030601 (2009).
  * Wong and Chen (2008) H. Wong and Z. D. Chen, Phys. Rev. B 77, 174305 (2008).
  * Alvermann and Fehske (2009) A. Alvermann and H. Fehske, Phys. Rev. Lett. 102, 150601 (2009).
  * Tong and Hou (2012) N. H. Tong and Y. H. Hou, Phys. Rev. B 85, 144425 (2012).
  * Zhou and Shao (2008) Y. Zhou and J. Shao, J. Chem. Phys. 128, 034106 (2008).
  * Yan and Shao (2016) Y. Yan and J. Shao, Front. Phys. 11, 110309 (2016).
  * Duan  _et al._ (2017) C. Duan, Z. Tang, J. Cao, and J. Wu, Phys. Rev. B 95, 214308 (2017).
  * Wang and Shao (2019) H. Wang and J. Shao, J. Phys. Chem. A 123, 1882 (2019).
  * Blunden-Codd  _et al._ (2017) Z. Blunden-Codd, S. Bera, B. Bruognolo, N. O. Linden, A. W. Chin, J. von Delft, A. Nazir, and S. Florens, Phys. Rev. B 95, 085104 (2017).
  * Bruognolo  _et al._ (2017) B. Bruognolo, N. O. Linden, F. Schwarz, S. S. B. Lee, K. Stadler, A. Weichselbaum, M. Vojta, F. B. Anders, and J. von Delft, Phys. Rev. B 95, 121115 (2017).
  * Zhang  _et al._ (2010) Y. Y. Zhang, Q. H. Chen, and K. L. Wang, Phys. Rev. B 81, 121105(R) (2010).
  * Leppäkangas  _et al._ (2018) J. Leppäkangas, J. Braumüller, M. Hauck, J. M. Reiner, I. Schwenk, S. Zanker, L. Fritz, A. V. Ustinov, M. Weides, and M. Marthaler, Phys. Rev. A 97, 052321 (2018).
  * Magazzù  _et al._ (2018) L. Magazzù, P. Forn-Díaz, R. Belyansky, J. L. Orgiazzi, M. A. Yurtalan, M. R. Otto, A. Lupascu, C. M. Wilson, and M. Grifoni, Nat. Commun. 9, 1403 (2018).
  * Lemmer  _et al._ (2018) A. Lemmer, C. Cormick, D. Tamascelli, T. Schaetz, S. F. Huelga, and M. B. Plenio, New J. Phys. 20, 073002 (2018).
  * He (1999) J. H. He, International Journal of Non-Linear Mechanics 34, 699 (1999).
  * Sorella (2005) S. Sorella, Phys. Rev. B 71, 241103 (2005).
  * Silbey and Harris (1984) R. Silbey and R. A. Harris, J. Chem. Phys. 80, 2615 (1984).
  * Chin  _et al._ (2011) A. W. Chin, J. Prior, S. F. Huelga, and M. B. Plenio, Phys. Rev. Lett. 107, 160601 (2011).
  * He  _et al._ (2018) S. He, L. W. Duan, and Q. H. Chen, Phys. Rev. B 97, 115157 (2018).
  * Frenzel and Plenio (2013) M. F. Frenzel and M. B. Plenio, New J. Phys. 15, 073046 (2013).
  * Zhou  _et al._ (2015a) N. J. Zhou, L. P. Chen, D. Z. Xu, V. Chernyak, and Y. Zhao, Phys. Rev. B 91, 195129 (2015a).
  * Zhou  _et al._ (2015b) N. J. Zhou, Z. K. Huang, J. F. Zhu, V. Chernyak, and Y. Zhao, J. Chem. Phys. 143, 014113 (2015b).
  * Zhou  _et al._ (2016) N. J. Zhou, L. P. Chen, Z. K. Huang, K. W. Sun, Y. Tanimura, and Y. Zhao, J. Phys. Chem. A 120, 1562 (2016).
  * Qian  _et al._ (2021) X. H. Qian, C. Z. Zeng, and N. J. Zhou, Physica A 580, 126157 (2021).
  * Qian  _et al._ (2022) X. Qian, Z. Sun, and N. Zhou, Phys. Rev. A 105, 012431 (2022).
  * Bulla  _et al._ (2005) R. Bulla, H. J. Lee, N. H. Tong, and M. Vojta, Phys. Rev. B 71, 045122 (2005).
  * Bera  _et al._ (2014a) S. Bera, S. Florens, H. U. Baranger, N. Roch, A. Nazir, and A. W. Chin, Phys. Rev. B 89, 121108 (2014a).
  * Lvovsky and Raymer (2009) A. I. Lvovsky and M. G. Raymer, Rev. Mod. Phys. 81, 299 (2009).
  * Fisher  _et al._ (1972) M. E. Fisher, S. k. Ma, and B. G. Nickel, Phys. Rev. Lett. 29, 917 (1972).
  * Bera  _et al._ (2014b) S. Bera, A. Nazir, A. W. Chin, H. U. Baranger, and S. Florens, Phys. Rev. B 90, 075110 (2014b).
  * Nazir  _et al._ (2012) A. Nazir, D. P. S. McCutcheon, and A. W. Chin, Phys. Rev. B 85, 224301 (2012).

Figure 1:  Schematic of SBM where a quantum spin is immersed in a bosonic bath
(environment) represented by a collection of harmonic oscillators. The sketch
of the quantum phase transition for the unbiased SBM is shown at the bottom,
which separates a non-degenerate delocalized phase from a doubly degenerate
localized one. The system energy difference $E(m)-E(m=0)$ $vs.$ the
magnetization $m$ is plotted for these two phases, and the ground-state values
of the magnetizations $m=\pm\langle\sigma_{z}\rangle$ and $0$ are obtained as
the local minima of the energy landscape, respectively.

Figure 2:  (a) The ground-state magnetization $\langle\sigma_{z}\rangle$ in
the shallow sub-Ohmic SBM is displayed as a function of the coupling $\alpha$
for different values of the spectral exponent $s=0.6,0.7,0.8$, and $0.9$ in
the case of $\Delta/\omega_{c}=0.1$ and $\varepsilon=0$. The discretization
factor $\Lambda=1.05$ and numbers of coherent-superposition states and
effective bath modes, $N=8$ and $M=430$, are set. (b) The log-log plot of the
magnetization $\langle\sigma_{z}\rangle\ vs.\ (\alpha-\alpha_{c})$ where
dashed lines represent power-law fits. Taking $s=0.9$ as an example, inset
illustrates the fitting error $RSS/DoF$ at different input value of the
critical coupling $\alpha_{c}$, and the arrow indicates the location of the
minimization.  Figure 3:  In the deep sub-Ohmic regime $s\leq 0.5$, the
magnetization $\langle\sigma_{z}\rangle$ with respect to the critical coupling
deviation $\alpha-\alpha_{c}$ is plotted with open circles, open triangles,
stars, open diamonds, and pluses, corresponding to $s=0.1,0.2,0.3,0.4$, and
$0.5$, respectively. For clarity, curves are shifted downwards or upwards.
Dashed lines show power-law fits.  Figure 4:  Power-law behaviors of
$\langle\sigma_{z}\rangle$ with respect to the shift $\alpha-\alpha_{c}$ are
displayed for the sizes $M=300$ and $350$ in the cases of $s=0.9,0.8$ and of
$s=0.7,0.6$, respectively. In the inset, the discretization factor
$\Lambda=2.0$ is set, much larger than the previous one $\Lambda=1.05$. Dashed
lines represent power-law fits.

Figure 5:  In the delocalized phase, the local susceptibility $\chi(\alpha)$
is plotted as a function of $\alpha_{c}-\alpha$ for the cases in the shallow
sub-Ohmic regime (a) and deep sub-Ohmic regime (b), respectively. The tiny
bias $\delta\varepsilon/\omega_{c}=10^{-5}$ is used to calculate the
susceptibility
$\chi=\partial\langle\sigma_{z}\rangle/\partial\varepsilon|_{\varepsilon=0}$,
and all power-law decaying curves are fitted by dashed lines.

Figure 6:  The spin magnetization $\langle\sigma_{z}\rangle$ at the critical point $\alpha_{c}$ is plotted in response to an external field $\varepsilon$ for various values of s in the shallow sub-Ohmic SBM (a) and deep sub-Ohmic SBM (b), respectively. Dashed lines show power-law fits, and the curves in (b) are shifted up for clarity.  $s$ |  NRG Bulla  _et al._ (2003); Vojta  _et al._ (2005) |  QMC Winter  _et al._ (2009) |  VMPS Guo  _et al._ (2012) |  DMRG Wong and Chen (2008) |  Single polaron Chin  _et al._ (2011) | This work
---|---|---|---|---|---|---
0.1 | 0.008(1) | 0.0076(3) | $\backslash$ | 0.0074(2) | 0.0065 | 0.00795(1)
0.2 | 0.018(1) | 0.0175(2) | 0.0175367 | 0.0162(2) | 0.0168 | 0.01806(3)
0.3 | 0.035(2) | 0.0350(5) | 0.0346142 | 0.0332(5) | 0.0316 | 0.03470(6)
0.4 | 0.064(2) | 0.0604(8) | 0.0605550 | 0.058(1) | 0.0519 | 0.0604(2)
0.5 | 0.106(2) | 0.098(1) | 0.0990626 | 0.099(1) | 0.0784 | 0.0977(3)
0.6 | 0.170(2) | 0.157(2) | 0.1554073 | 0.155(1) | $\backslash$ | 0.1528(3)
0.7 | 0.261(2) | 0.241(2) | $\backslash$ | 0.238(2) | $\backslash$ | 0.2345(8)
0.8 | 0.392(3) | 0.360(2) | $\backslash$ | 0.359(2) | $\backslash$ | 0.354(1)
0.9 | 0.612(3) | 0.548(2) | 0.5555478 | 0.556(2) | $\backslash$ | 0.534(1)
Table 1:  The critical coupling $\alpha_{c}$ in the sub-Ohmic SBM is shown for different values of the spectral exponents $s$. The last two columns correspond to those obtained by NVM on the basis of single polaron ansatz ($N=1$) and multiple poalron ansatz ($N=8$), respectively. With the plotdigitizing technique, the NRG, QMC, and DMRG results are extracted from the figures in literatures, and the corresponding error bar is estimated as the maximum inaccuracy of the readout.  $s$ |  Theoretical values Fisher  _et al._ (1972) | Our work |  VMPS results Guo  _et al._ (2012)
---|---|---|---
$\beta$ | $\gamma$ | $1/\delta$ | $\beta$ | $\gamma$ | $1/\delta$ | $\beta$ | $\gamma$ | $1/\delta$
0.1 | 1/2 | 1 | 1/3 | 0.489(9) | 1.02(1) | 0.329(2) | $\backslash$ | $\backslash$ | $\backslash$
0.2 | 1/2 | 1 | 1/3 | 0.496(7) | 1.04(2) | 0.320(2) | 0.47(7) | $\backslash$ | 0.39(8)
0.3 | 1/2 | 1 | 1/3 | 0.496(5) | 1.00(2) | 0.322(3) | 0.50(3) | $\backslash$ | 0.34(2)
0.4 | 1/2 | 1 | 1/3 | 0.495(7) | 1.01(2) | 0.324(7) | 0.50(3) | $\backslash$ | 0.33(1)
0.5 | 1/2 | 1 | 1/3 | 0.491(6) | 1.02(2) | 0.323(7) | 0.46(3) | $\backslash$ | 0.30(1)
0.6 | 0.392 | 1.18 | 0.250 | 0.406(4) | 1.12(1) | 0.250(3) | 0.38(1) | $\backslash$ | 0.244(6)
0.7 | 0.309 | 1.44 | 0.176 | 0.302(4) | 1.32(1) | 0.181(2) | 0.292(3) | $\backslash$ | 0.171(3)
0.8 | 0.224 | 1.79 | 0.111 | 0.229(3) | 1.73(1) | 0.117(2) | 0.211(2) | $\backslash$ | 0.109(2)
0.9 | 0.124 | 2.23 | 0.053 | 0.128(1) | 2.11(1) | 0.071(4) | 0.132(1) | $\backslash$ | 0.0523(8)
Table 2:  The critical exponents $\beta,\gamma$, and $\delta$ for different
$s$ are measured from extensive calculations with NVM, compared with those in
the VMPS work and theoretical analysis. The latter are obtained from the
analytical work based on the mean-field ($s<0.5$) and two-loop
renormalization-group theories ($s>0.5$) as well as the hyper-scaling law.
Figure 7: The spin coherences $\langle\sigma_{x}\rangle$ in the deep ($s=0.1$)
and shallow ($s=0.9$) sub-Ohmic cases are displayed as a function of the
coupling $\alpha$ in the right-top panel and left-bottom panel, respectively.
As before, other parameters
$\Delta/\omega_{c}=0.1,\varepsilon=0,\Lambda=1.05,N=8$ and $M=430$ are set.
The transition points $\alpha_{c}=0.00795$ and $0.534$ are marked by the
arrows. In the inset, the difference $\delta t=t_{2}-t_{1}$ between the slopes
on two sides of the transition point is plotted against the spectral exponent
$s$. Dashed lines show linear fits, and solid line in the inset represents the
fit with an exponential function.

Figure 8:  Wigner distribution $W^{k}(X)$ which allows direct characterization
of the ground-state wavefunction relevant to the equilibrium bosonic bath is
shown for the delocalized phase $\alpha<\alpha_{c}$ and localized phase
$\alpha>\alpha_{c}$, taking the cases of $s=0.2$ in (a) and $s=0.7$ in (b) as
examples. All bath modes are traced out except the one with the lowest
frequency $\omega_{k}=\omega_{\rm min}$, and the wigner distribution of
$\Psi_{\uparrow}(X)$ and $\Psi_{\downarrow}(X)$ representing twofold-
degenerate ground states are scaled by a factor $1/\sqrt{2}$ for clarity.
Figure 9: The shift of the transition point
$\delta\alpha_{c}=\alpha_{c}(L)-\alpha_{c,\infty}$ is displayed against the
inverse of the size $1/L$ for different values of the spectral exponent $s$ on
a log-log scale. The critical exponent $1/\nu$ extracted from the finite-size
scaling behavior $\delta\alpha_{c}\sim L^{-1/\nu}$ is shown in the inset as a
function of $s$. For comparison, the theoretical prediction of $1/\nu$ is also
given with the dotted line which is obtained from the mean-field ($s\leq 0.5$)
and two-loop renormalization-group ($s>0.5$) methods. Moreover, solid line in
the inset represents the perturbative prediction $\nu=1/\sqrt{2(1-s)}$ when
$s$ close to $1$.