[2,3,4]Ali Bashashati

1]Faculty of Science, University of British Columbia, 2207 Main Mall,
Vancouver, V6T 1Z4, British Columbia, Canada

2]School of Biomedical Engineering, University of British Columbia, 2222
Health Sciences Mall, Vancouver, V6T 2B9, British Columbia, Canada

3]Department of Pathology and Laboratory Medicine, University of British
Columbia, 2211 Wesbrook Mall, Vancouver, V6T 1Z7, British Columbia, Canada

4]Canada’s Michael Smith Genome Sciences Centre, BC Cancer Research Institute,
570 W 7th Ave, Vancouver, V5Z 4S6, British Columbia, Canada

# Benchmarking Histopathology Foundation Models for Ovarian Cancer Bevacizumab
Treatment Response Prediction from Whole Slide Images

Mayur Mallya  0000-0002-9432-4262<EMAIL_ADDRESS>Ali Khajegili Mirabadi
0000-0002-0489-6583<EMAIL_ADDRESS>Hossein Farahani  0000-0002-9503-1875
<EMAIL_ADDRESS>0000-0002-4212-7224<EMAIL_ADDRESS>[  [  [  [

###### Abstract

Purpose: Bevacizumab is a widely studied targeted therapeutic drug used in
conjunction with standard chemotherapy for the treatment of recurrent ovarian
cancer. While its administration has shown to increase the progression-free
survival (PFS) in patients with advanced stage ovarian cancer, the lack of
identifiable biomarkers for predicting patient response has been a major
roadblock in its effective adoption towards personalized medicine. Methods: In
this work, we leverage the latest histopathology foundation models trained on
large-scale whole slide image (WSI) datasets to extract ovarian tumor tissue
features for predicting bevacizumab response from WSIs. Results: Our extensive
experiments across a combination of different histopathology foundation models
and multiple instance learning (MIL) strategies demonstrate capability of
these large models in predicting bevacizumab response in ovarian cancer
patients with the best models achieving an AUC score of 0.86 and an accuracy
score of 72.5%. Furthermore, our survival models are able to stratify high-
and low-risk cases with statistical significance (p $<$ 0.05) even among the
patients with the aggressive subtype of high-grade serous ovarian carcinoma.
Conclusion: This work highlights the utility of histopathology foundation
models for the task of ovarian bevacizumab response prediction from WSIs. The
high-attention regions of the WSIs highlighted by these models not only aid
the model explainability but also serve as promising imaging biomarkers for
treatment prognosis.

###### keywords:

Ovarian Cancer, Bevacizumab Therapy, Foundation Models, Treatment Response
Prediction, Survival Analysis

##  1 Introduction

Ovarian cancer is a highly lethal gynecologic disease and one of the leading
causes of cancer-related deaths among women. In 2020, a total of 313,959 new
cases of ovarian cancer were recorded globally that resulted in 207,252 new
deaths [14]. Among these epithelial ovarian carcinoma (EOC) is the most common
type accounting for about 90$\%$ of all ovarian malignancies [32]. EOC is
heterogeneous disease with distinct subtypes each with their own clinical and
molecular characteristics that can impact the prognosis and treatment
response. High-grade serous carcinoma (HGSC) is the most common and aggressive
subtype of EOC with a disproportionate share of fatalities.

Standard treatment for a newly diagnosed EOC case involves a surgical
cytoreduction followed by paclitaxel and platinum-based chemotherapy
treatment. However, 70$\%$ of the cases are diagnosed at an advanced stage
when the treatment options are not only limited but also ineffective, leading
to treatment resistance and tumor recurrence and thereby contributing to poor
prognosis and patient outcomes [4]. The limited responsiveness to conventional
chemotherapy at an advanced stage coupled and the resulting high mortality
rate has necessitated the usage of targeted therapeutic agents in the
treatment of EOC.

Bevacizumab (clinically known as Avastin) is an extensively studied targeted
therapeutic agent used in the treatment of recurrent EOC [25, 10]. It plays an
important role in the inhibition of tumor angiogenesis by neutralizing the
vascular endothelial growth factor (VEGF), a key signaling protein responsible
for tumor regrowth. The administration of bevacizumab in conjunction with
chemotherapy has shown to increase the progression-free survival (PFS) in
advanced stage EOC [3, 26]. However, the lack of effective biomarkers for
predicting patient responses remains as a challenge for the personalization of
bevacizumab therapy to this day. Additionally, given the high costs and
potential adverse side-effects, identifying EOC patients with favorable
response to bevacizumab therapy becomes an imperative task.

In this work, we use the publicly available ovarian bevacizumab response
dataset of histopathological whole slide images (WSI) to predict the treatment
effectiveness for the patient [34]. Analysis of hematoxylin and eosin (H&E)
stained WSIs is a cost-effective and routine practice in clinical pathology,
offering valuable insights into the tumor microenvironment and morphology. The
WSIs play an integral role in the clinical diagnosis and management of ovarian
cancer. Advancements in deep learning (DL) over the last decade have enabled
the computational analysis of WSIs, which are gigapixel resolution and often
lack detailed annotations from pathologists. Weakly-supervised learning
strategies such as multiple instance learning (MIL) have been prominently
applied and have shown great success across a wide variety of WSI analysis
tasks such as grading [31], subtyping [30, 23], survival analysis [40, 20],
etc. By breaking down the WSI into a bag of multiple patches, MIL methods can
automatically identify the most discriminative patches and aggregate their
features to generate slide-level features which can be used for downstream
tasks.

Due to the small-size nature of medical image datasets, including those in
histopathology, prior works on ovarian bevacizumab response prediction have
either relied on extensive WSI pre-processing strategies for efficient patch
selection from the high-resolution WSIs [35] or leveraged more informative
molecular counterparts for this task [36, 37]. A popular approach for dealing
with small-size datasets is using transfer learning strategies where models
pre-trained on natural image datasets such as ImageNet are used as feature
extractors for histopathology images [1]. Despite looking significantly
different from histopathology images, the ImageNet pre-trained models provide
a strong backbone network which can be fine-tuned for task-specific
applications. However, the past 2 years have seen the dramatic rise of
histopathology foundation models that are trained on massive amounts of tissue
data in a self-supervised fashion, thereby providing strong domain-specific
feature extractors for histopathology. These models have shown to outperform
the models pre-trained on natural images across a wide variety of primary
histopathology analysis tasks such as tumor detection, grading, and subtyping
across pan-cancer data [5, 8, 15, 17, 38]. However, to the best of our
knowledge, their efficacy on secondary tasks such as treatment response
prediction, mutation prediction, etc. from WSIs is not fully explored.

Our contribution in this work is 3 folds. 1) We provide the first
comprehensive benchmark for evaluating the performance of histopathology
foundation models across multiple MIL frameworks for the task of bevacizumab
treatment response prediction. Our study highlights the superior performance
of histopathology foundation models in comparison to the models pre-trained on
natural image datasets emphasizing the utility of domain-specific encoders for
a relatively unexplored secondary histopathology analysis task of treatment
response prediction. 2) Our models are able to identify the patients with
favorable response to bevacizumab therapy with an AUC score of 0.86 and an
accuracy of 72.5% when no effective clinical biomarkers exist for predicting
the treatment effectiveness for this task. Furthermore, our survival analysis
experiments demonstrate statistically significant risk stratification between
high- and low- risk cases even among the patients with the highly aggressive
HGSC subtype. 3) Our models identify high-attention tumorous areas in the WSIs
providing a promising approach for the identification of prognostic imaging
biomarkers for bevacizumab treatment response in ovarian cancer patients.

##  2 Materials and Methods

###  2.1 Histopathology Foundation Models

Self-supervised learning has fueled the recent development of foundation
models by leveraging vast amounts of unlabeled data which are commonly found
in digital pathology. Trained on hundreds of thousands of tissue patches,
typically spanning across multiple cancer types, these models are able to
learn powerful task-agnostic tissue representations. CTransPath [38] and Lunit
[17] foundation models trained on the entire TCGA cohort [39] have shown
improvements across a variety of primary analysis tasks such as tumor
subtyping, mitosis detection, and cell segmentation. Huang et al. fine-tuned
the CLIP [27] model using the histopathology text-image pairs from Twitter to
produce PLIP [15] which improved the tumor detection and tissue grading and
subtyping tasks. UNI [5] and Virchow [33] are trained on some of the largest
private data cohorts with 100,000 and 1,000,000 WSIs respectively. Among
these, only Phikon [8] and Virchow [33] explored their performance on a
secondary analysis task of mutation prediction while only Phikon explored the
task of survival prediction. However, their tests on secondary analysis tasks
were limited and didn’t produce unanimous conclusions on survival prediction
task.

Table 1: Summary of the histopathology foundation models.

Model | Public | Training data | Cohort size | Model size
---|---|---|---|---
CTransPath [38] | ✓ | TCGA + PAIP | 32K | 28M
Lunit-Dino [17] | ✓ | TCGA + TULIP | 37K | 22M
Phikon [8] | ✓ | TCGA | 6K | 86M
PLIP [15] | ✓ | OpenPath | – | 86M
UNI [5] | ✓ | Mass-100K | 100K | 307M
Virchow [33] | X | MSKCC | 1.5M | 632M

00footnotetext: In addition to TCGA, the publicly available datasets include
PAIP [18] and OpenPath [15].

###  2.2 Multiple Instance Learning

The weakly-supervised learning nature of MIL models suits perfectly for the
analysis of local regions-of-interest in gigapixel WSIs. The MIL-based deep
models have seen tremendous success in computational pathology in the recent
years [9]. Ilse et al. [16] proposed the first learnable attention-based
aggregation strategy in MIL models (ABMIL) that outperformed the traditional
pooling methods. Subsequent works such as CLAM [23] and VarMIL [29] built on
top of this method by improving latent representations of the patch and slide
features. The onset of self-attention has enabled the integration of
transformer-based aggregation strategies in MIL models such as TransMIL [30].
The consistent success of MIL models have made it a default strategy for the
analysis of WSIs and in this work we use the aforementioned MIL methods that
have been used extensively in the literature.

###  2.3 Dataset

Ovarian bevacizumab response dataset [34] is a publicly available dataset
provided by the Cancer Imaging Archive (TCIA). The dataset consists of 286
hematoxylin and eosin (H&E) stained whole section WSIs from 78 patients
scanned at $20\times$ magnification from the tissue bank of the Tri-Service
General Hospital and the National Defense Medical Center, Taipei, Taiwan. The
dataset also includes the clinical information of the patients such as the PFS
along with the treatment effectiveness of the bevacizumab treatment. The
ground-truth binary treatment response was identified based on the pre- and
post-treatment CA-125 concentrations, with the bevacizumab treatment being
effective for 160 patient slides and ineffective for the remaining 126 slides
in the cohort. Fig. 1a. shows the patient distribution across the different
ovarian cancer subtypes in the dataset and Fig. 1b. shows the slide-level
distribution along with the subtype-wise binary effectiveness label splits.

For training the model, we divide the dataset into 3 folds with training (70%)
and validation (15%) splits, and we use the same held-out testing split (15%)
to evaluate the models across all folds. Our test set has an equal class
distribution with 37 slides each in the effective and ineffective treatment
classes (total test set size = 74 slides). The data splits are done at the
patient level so that all slides from a patient belong to the same split.
Additionally, we conduct experiments exclusively on the serous subtype
(includes peritoneal serous papillary carcinoma and papillary serous
carcinoma) which is the majority subtype in the dataset as shown in Fig. 1.

Figure 1: a) Patient distribution across the different subtypes in the ovarian
bevacizumab response dataset. b) Slide distribution across the different
subtypes showing the label distribution for each subtype. The ovarian cancer
subtypes include clear cell carcinoma (CC), endometrioid carcinoma (EC),
mucinous carcinoma (MC), peritoneal serous papillary carcinoma (PsPC),
papillary serous carcinoma (PsC), unclassifed carcinoma (UC). The slide labels
include “Effective” corresponding to the patients with a favorable response to
bevacizumab treatment and “Not effective” corresponding the non-responders of
the bevacizumab treatment.

###  2.4 Problem Formulation

Figure 2: Whole slide image (WSI) analysis pipeline depicting the steps
involved in the processing of WSIs that includes pre-processing (tissue area
masking, patch sampling, color normalization), patch-level feature extraction
from foundation models, aggregating patch-level features to produce a slide-
level representation using MIL model, and prediction of treatment response.
Note that we only train the MIL model along with the MLP layers while training
our models with the cross-entropy loss ($L_{CE}$) for a binary classification
task of treatment effectiveness prediction or the Cox partial likelihood loss
($L_{COX}$) for the time-to-event regression task of survival prediction.

We formulate the problem of treatment response prediction from WSIs using two
approaches. First is the binary classification approach where the model
predicts the treatment effectiveness from WSI whether the treatment was
effective ($y=1$) or ineffective ($y=0$). Second is the survival prediction
problem where the model uses the WSIs to predict the hazard score of the
patient relative to other patients in the cohort.

We denote the WSI as $W$ and the corresponding binary treatment label, time-
to-event, and censor status as $y$, $t$, and $e$ respectively. Our dataset can
then be represented as
$\\{(W_{1},y_{1},t_{1},e_{1}),(W_{2},y_{2},t_{2},e_{2}),\cdots(W_{n},y_{n},t_{n},e_{n})\\}$,
where $n$ is the total number of WSIs in the dataset. Each WSI can be treated
as a bag of tissue patches denoted by $p$ which can be represented as
$W=\\{p_{1},p_{2}\cdots p_{k}\\}$, where $k$ is the number of patches
extracted from each WSI. As part of the pre-processing step, each patch is
passed through the color normalization module where a reference patch is used
to normalize the stains across different patches.

Each patch $p_{i}$ is passed through a pre-trained histopathology foundation
model, denoted by $\textbf{E}(\cdot)$ to produce the corresponding feature
representation $f_{i}$,

$f_{i}=\textbf{E}(p_{i};\theta^{*}_{E})$ (1)

where $\theta_{E}$ are the parameters of the pre-trained foundation model and
the asterisk indicates the parameters are frozen. Each WSI can now be
represented as a bag of patch features as, $W=\\{f_{1},f_{2}\cdots f_{k}\\}$.
This bag of features is passed to the MIL model denoted by $\textbf{M}(\cdot)$
in order to aggregate the patch-level features of a WSI $W_{i}$ to form a
slide-level representation $s_{i}$ in a learnable manner. This can be
represented as,

$s_{i}=\textbf{M}(W_{i};\theta_{M})=\textbf{M}(\\{p_{i1},p_{i2}\cdots
p_{ik}\\};\theta_{M})$ (2)

where $\theta_{M}$ denotes the trainable parameters of the MIL model
$\textbf{M}(\cdot)$. The slide-level representation $s_{i}$ is used for
predicting the treatment response for the patient corresponding to the WSI
$W_{i}$. The slide-level representation is passed to the MLP layers denoted by
$\textbf{MLP}(\cdot)$ to predict the binary treatment response $\hat{y_{i}}$
as follows,

$\hat{y_{i}}=\textbf{MLP}(s_{i};\theta_{MLP})$ (3)

where $\theta_{MLP}$ denotes the trainable parameters of the MLP layers. In
the same way, the MLP layers can be used to predict the logarithmic hazard
score $h_{i}$ from the slide-level representation as follows,

$h_{i}=\text{exp}(\textbf{MLP}(s_{i};\theta_{MLP}))$ (4)

Given the model predictions $\hat{y_{i}}$ and $h_{i}$, we use the ground-truth
treatment response $y_{i}$ along with the time-to-event $t_{i}$ and censor
status $e_{i}$ corresponding to WSI $W_{i}$ to train the model parameters. For
training the classification model, we use the binary cross-entropy loss
function $L_{CE}(\cdot)$ and for the survival prediction model, we use the cox
negative partial log-likelihood loss function $L_{COX}$ [6], which can be
calculated as follows,

$L_{CE}=-\frac{1}{N}\sum_{i=1}^{N}y_{i}\cdot\text{log}(\hat{y_{i}})+(1-y_{i})\cdot\text{log}(1-\hat{y_{i}})$
(5)

$L_{COX}=-\sum_{i=1|\forall
e_{i}=1}^{N}\left(\text{log}(h_{i})-\text{log}\sum_{j\in
R(t_{i})}h_{j}\right)$ (6)

Eq. 5 shows the binary cross entropy loss function where $N$ denotes the
number of WSIs in the training set and Eq. 6 denotes the cox negative partial
log-likelihood function where $R(t_{i})=\\{j:t_{j}\geq t_{i}\\}$ denotes the
risk set at time $t_{i}$. In both cases, we minimize the loss functions while
updating the parameters of MIL model $\theta_{M}$ and MLP layers
$\theta_{MLP}$ while the parameters of the encoder model $\theta_{E}$ are
frozen.

###  2.5 Implementation

In our experiments, we use a patch size of $224\times 224$ at the original
$20\times$ magnification and sample $k=300$ patches from each WSI. We use the
Macenko stain normalization technique [24, 2] to alleviate the staining
disparities across different WSIs by performing the normalization across all
extracted patches. Further, for training the models, we use the binary cross-
entropy loss (Eq. 5) for the classification setting and the cox partial log-
likelihood loss (Eq. 6) for the survival prediction setting. For all our
experiments, we use Adam optimizer [19] with a learning rate of $1\times
10^{-3}$ and a weight decay of $1\times 10^{-2}$. All models are trained for a
maximum of $50$ epochs and the model at the best validation performance epoch
during the training process is used for inference. To ensure the statistical
stability of models, we train all the models with 10 different random seeds
across 3 folds and report the 3-fold average metric with the best-performing
seed. We use the accuracy and AUC scores as the metrics to evaluate the
performance of the models as our test set has an equal class distribution. To
evaluate the survival models, we use the commonly used concordance index
(c-index) metric [11] in addition to evaluating the risk stratification using
the Kaplan-Meier (KM) curves with the log-rank test. For all our experiments,
we use the NVIDIA GeForce RTX 3090 and RTX A6000 GPUs to train the models.

##  3 Results and Discussion

In our experiments, we evaluate the performance of the latest histopathology
foundation models on a secondary histopathology image analysis task of ovarian
cancer bevacizumab response prediction. We compare the performance of the
histopathology foundation models with that of the traditionally used encoders
pre-trained on natural image datasets. For the histopathology foundation
models, we use the openly accessible models namely Phikon [8], PLIP [15], UNI
[5], Lunit-Dino [17], and CTransPath [38]. For the natural image encoders, we
use the convolutional backbone models such as ResNet50 [12], DenseNet121 [13],
and ConvNeXt [22] along with the transformer backbone models such as ViT [7]
and Swin [21]. Additionally, we also use the KimiaNet [28] models which are
trained on histopathology images in a supervised fashion unlike the
histopathology foundation models that are trained in a self-supervised manner.
We evaluate the performance of all encoders across multiple MIL models and
average the results to produce a robust benchmark for this task.

Table 2: Results of binary classification of bevacizumab treatment
effectiveness prediction from histopathology images. The columns represent the
MIL models $\textbf{M}(\cdot)$ and the rows represent the pre-trained encoders
$\textbf{E}(\cdot)$. The values represent the percentage accuracy scores of
prediction averaged across a 3-fold experiment on a held-out test set. Note
that encoders in the top half of the table are trained on natural image
datasets while those at the bottom are trained on histopathology datasets.

Encoder | ABMIL | TransMIL | VarMIL | CLAM-SB | CLAM-MB | Average
---|---|---|---|---|---|---
ResNet50 | 50.00 $\pm$ 0.0 | 50.00 $\pm$ 0.0 | 51.35 $\pm$ 1.9 | 50.00 $\pm$ 0.0 | 50.90 $\pm$ 1.2 | 50.45
DenseNet121 | 50.45 $\pm$ 0.6 | 53.15 $\pm$ 2.7 | 53.15 $\pm$ 4.4 | 57.21 $\pm$ 3.1 | 55.86 $\pm$ 3.5 | 53.96
ConvNeXt | 62.16 $\pm$ 3.8 | 55.41 $\pm$ 3.9 | 65.77 $\pm$ 2.7 | 63.96 $\pm$ 5.5 | 65.77 $\pm$ 3.3 | 62.61
ViT | 63.06 $\pm$ 5.5 | 53.60 $\pm$ 2.7 | 65.32 $\pm$ 2.7 | 67.57 $\pm$ 1.1 | 65.32 $\pm$ 10.0 | 62.97
Swin (tiny) | 65.32 $\pm$ 8.4 | 52.70 $\pm$ 5.0 | 65.77 $\pm$ 2.3 | 66.22 $\pm$ 5.0 | 70.72 $\pm$ 4.9 | 64.14
Swin | 64.86 $\pm$ 3.3 | 57.21 $\pm$ 8.9 | 64.86 $\pm$ 3.3 | 68.47 $\pm$ 2.7 | 68.02 $\pm$ 3.3 | 64.68
KimiaNet | 50.90 $\pm$ 1.2 | 52.70 $\pm$ 1.1 | 62.16 $\pm$ 7.7 | 53.15 $\pm$ 0.6 | 58.56 $\pm$ 6.7 | 55.49
KimiaNet (OV) | 57.21 $\pm$ 5.2 | 54.05 $\pm$ 3.9 | 57.66 $\pm$ 6.0 | 54.50 $\pm$ 4.1 | 59.01 $\pm$ 2.3 | 56.48
Phikon | 63.96 $\pm$ 7.3 | 67.12 $\pm$ 2.3 | 67.12 $\pm$ 11.2 | 66.67 $\pm$ 3.8 | 68.02 $\pm$ 4.7 | 66.57
PLIP | 65.32 $\pm$ 2.3 | 50.00 $\pm$ 0.0 | 66.22 $\pm$ 1.9 | 69.82 $\pm$ 2.7 | 71.17 $\pm$ 6.0 | 64.50
UNI | 72.07 $\pm$ 3.0 | 54.95 $\pm$ 8.0 | 69.37 $\pm$ 9.1 | 65.32 $\pm$ 7.2 | 65.77 $\pm$ 3.2 | 65.49
Lunit-Dino | 71.17 $\pm$ 3.1 | 70.72 $\pm$ 4.9 | 68.47 $\pm$ 0.1 | 69.82 $\pm$ 2.5 | 65.32 $\pm$ 4.5 | 69.10
CTransPath | 69.37 $\pm$ 3.5 | 63.06 $\pm$ 0.1 | 72.07 $\pm$ 5.2 | 72.52 $\pm$ 6.0 | 71.62 $\pm$ 3.9 | 69.72

00footnotetext: KimiaNet (OV) refers to the pre-trained KimiaNet model fine-
tuned on an interal ovarian dataset for subtype classification in a supervised
fashion.

Table 2 presents the results of binary classification of bevacizumab treatment
response prediction on a held-out test set with the values averaged across a
3-fold cross-validation. We divide the table into 3 parts – encoders pre-
trained on natural images at the top, encoders pre-trained on histopathology
datasets in a supervised manner in the middle, and the histopathology
foundation model encoders pre-trained in a self-supervised manner at the
bottom of the table. The histopathology foundation models achieve the best
prediction performance across all MIL frameworks with CTransPath achieving a
72.5% accuracy with Clam-SB model. Furthermore, the average column showing the
average performance of an encoder across MIL models shows that histopathology
foundation models achieve superior performance compared to other encoders
irrespective of the choice of the MIL framework. Figure 3a. shows the ROC
curves for the best model from each encoder along with the corresponding AUC
scores and similar to Table 2, the histopathology foundation models have the
highest AUC scores with 4 of the 5 models (except Phikon with an AUC score of
0.7) achieving an AUC score of 0.8 and above.

Figure 3: a) Receiver Operating Characteristic (ROC) curves showing the
treatment response prediction performance along with the AUC score for
different encoders. b) Top 5 high-attention patches sorted in the descending
order of the attention values (left to right) for the best-performing models
from our experiments namely, 1) CTransPath with Clam-SB, 2) CTransPath with
VarMIL, 3) Lunit-Dino with ABMIL.

Table 3 presents the results of bevacizumab treatment survival prediction on a
held-out test set averaged across 3-fold cross-validation experiments. Similar
to the results of the binary classification in Table 2, the histopathology
foundation models outperform the models pre-trained on natural image datasets
on survival prediction task across different MIL models. The histopathology
foundation models on average achieve a higher c-index score as compared to the
other models with CTransPath and UNI achieving an average c-index of 0.64.
Among the natural image encoders, the performance of Swin is on par with the
histopathology foundation models with an average c-index of 0.63 while the
other natural image encoders are significantly worse compared to the best
models. Furthermore, we present the KM survival plots for the best performing
models in Table 3 in Figs. 4 and 5 showing the stratification of high- and
low-risk cases based on the predicted hazard score. Fig. 4 shows the KM plot
corresponding to all cases in the test set while Fig. 5 shows the KM plot
corresponding to only the serous cases of the test set. We use the median PFS
value of this cohort as the threshold for differentiating the high- and low-
risk cases. In the case of Fig. 4 that includes all subtypes, we observe
statistically significant risk stratification on the log-rank test ($p<0.05$)
across all the top-5 best performing models while in the case of Fig. 5 that
includes only the serous subtype, only 2 of the 5 best performing models
achieve statistically significant stratification while the remaining 3 models
achieve a p-value close to the 0.05 cut-off.

Table 3: Results of survival prediction from histopathology images. The
columns represent the MIL models $\textbf{M}(\cdot)$ and the rows represent
the pre-trained encoders $\textbf{E}(\cdot)$. The values represent the c-index
scores of averaged across a 3-fold experiment on a held-out test set.

Encoder | ABMIL | TransMIL | VarMIL | CLAM-SB | CLAM-MB | Average
---|---|---|---|---|---|---
ResNet50 | 0.59 $\pm$ 0.01 | 0.56 $\pm$ 0.02 | 0.57 $\pm$ 0.02 | 0.59 $\pm$ 0.01 | 0.58 $\pm$ 0.01 | 0.58
DenseNet121 | 0.54 $\pm$ 0.04 | 0.55 $\pm$ 0.02 | 0.57 $\pm$ 0.03 | 0.57 $\pm$ 0.01 | 0.57 $\pm$ 0.01 | 0.56
ViT | 0.55 $\pm$ 0.02 | 0.60 $\pm$ 0.02 | 0.54 $\pm$ 0.00 | 0.59 $\pm$ 0.03 | 0.59 $\pm$ 0.03 | 0.57
Swin | 0.62 $\pm$ 0.03 | 0.67 $\pm$ 0.02 | 0.59 $\pm$ 0.02 | 0.63 $\pm$ 0.01 | 0.63 $\pm$ 0.01 | 0.63
KimiaNet (OV) | 0.53 $\pm$ 0.01 | 0.58 $\pm$ 0.01 | 0.56 $\pm$ 0.02 | 0.49 $\pm$ 0.04 | 0.49 $\pm$ 0.04 | 0.53
Phikon | 0.59 $\pm$ 0.01 | 0.60 $\pm$ 0.02 | 0.59 $\pm$ 0.02 | 0.60 $\pm$ 0.03 | 0.60 $\pm$ 0.03 | 0.60
PLIP | 0.62 $\pm$ 0.02 | 0.63 $\pm$ 0.02 | 0.63 $\pm$ 0.02 | 0.63 $\pm$ 0.01 | 0.63 $\pm$ 0.01 | 0.63
Lunit-Dino | 0.63 $\pm$ 0.00 | 0.62 $\pm$ 0.01 | 0.60 $\pm$ 0.01 | 0.59 $\pm$ 0.01 | 0.59 $\pm$ 0.01 | 0.60
CTransPath | 0.66 $\pm$ 0.03 | 0.62 $\pm$ 0.02 | 0.65 $\pm$ 0.03 | 0.63 $\pm$ 0.01 | 0.63 $\pm$ 0.01 | 0.64
UNI | 0.60 $\pm$ 0.02 | 0.65 $\pm$ 0.03 | 0.63 $\pm$ 0.02 | 0.67 $\pm$ 0.03 | 0.67 $\pm$ 0.03 | 0.64

00footnotetext: We do not include ConvNext, Swin (tiny), and KimiaNet due to
the lack of model convergence due to the NaN values in the loss function.

Interestingly, both KimiaNet and KimiaNet-OV (KimiaNet fine-tuned on an
internal ovarian dataset for subtype classification in supervised fashion)
features fail to perform well on this task showing that supervised pre-
training doesn’t generalize for tasks outside of its expertise while
emphasizing the need for task-agnostic self-supervised histopathology pre-
training. Moreover, the KimiaNet models use a convolutional backbone similar
to that of ResNet50, DenseNet121, and ConvNeXt while all other encoders in
Table 2 use a transformer-based architecture, and the performance disparity
due to the choice of model architecture underscores the advantage of
transformer-based architecture over convolutional backbones for histopathology
analysis task.

Figure 4: Kaplan-Meier plots stratifying the high- and low-risk cases in the
test set across the 5 best performing models from Table 3. a) CTransPath with
ABMIL b) CTransPath with VarMIL c) Swin with TransMIL d) UNI with Clam-SB e)
UNI with Clam-MB.  Figure 5: Kaplan-Meier plots stratifying the high- and low-
risk serous cases in the test set across the 5 best performing models from
Table 3. a) CTransPath with ABMIL b) CTransPath with VarMIL c) Swin with
TransMIL d) UNI with Clam-SB e) UNI with Clam-MB.

Model explainability is an important aspect of medical image analysis and the
MIL models in our study use the attention mechanism to rank the patches in the
order of their informativeness to the task. Figure 3b. shows the top-5 high
attention tumor patches from the 3 best performing models from our experiments
in Table 2. We observe that despite the similar performance, these models
assigned different attention to the patches without any common patches among
the top-5 high-attention patches across these models. Nevertheless, this
approach provides a promising direction for the identification of novel
imaging biomarkers for ovarian cancer bevacizumab therapy and we encourage the
future works to validate this approach on other ovarian cohorts to confirm the
findings. Exploring biological insights from the high-attention patches from
different models can also be a valuable research direction for future works.

##  4 Conclusion

Lack of effective biomarkers for predicting the treatment response of ovarian
bevacizumab therapy has been a long-standing hurdle in the personalization of
its treatment. Additionally, the high costs and potential toxicity make it
imperative to identify patient response to ovarian bevacizumab treatment. In
this work, we perform the first comprehensive benchmarking study evaluating
the performance of newly developed histopathology foundation models (along
with traditionally used models trained on natural images) on a relatively
unexplored secondary histopathology analysis task of ovarian bevacizumab
treatment response prediction. In addition to achieving a prediction AUC of
0.86 and an accuracy of 72.5%, our models significantly stratify the low-risk
cases from the high-risk ones. Furthermore, the models can identify
informative tumor regions in the WSIs corresponding to the treatment response
prediction thereby serving as a promising direction for the identification of
prognostic imaging biomarkers for bevacizumab treatment of ovarian cancer.

Acknowledgements The authors acknowledge the Canada’s Michael Smith Genome
Sciences Centre for hosting the computational clusters used in this research.

## Declarations

Funding

This research was funded by the Natural Sciences and Engineering Research
Council of Canada (NSERC), Canadian Institutes of Health Research (CIHR),
Health Research BC, and the BC Cancer Foundation.

Competing interests

The authors declare no competing interests.

Ethics approval and consent to participate

This research did not require ethics approval and consent to participate as
the dataset used in this study is publicly available.

Consent for publication

Not applicable.

Data availability

The ovarian bevacizumab response prediction dataset used in this study is
publicly available and can be downloaded from the cancer imaging archive
(TCIA) database at https://www.cancerimagingarchive.net/collection/ovarian-
bevacizumab-response/

Author contribution

The project was conceptualized and administered by HF and AB. AKM developed
the codebase for running the experiments. MM ran the experiments and wrote the
manuscript. All authors reviewed and approved the submitted version.

## References

  *   * Aitazaz et al [2023] Aitazaz T, Tubaishat A, Al-Obeidat F, et al (2023) Transfer learning for histopathology images: an empirical study. Neural Computing and Applications 35(11):7963–7974
  * Boschman et al [2022] Boschman J, Farahani H, Darbandsari A, et al (2022) The utility of color normalization for ai-based diagnosis of hematoxylin and eosin-stained pathology images. The Journal of Pathology 256(1):15–24
  * Burger et al [2011] Burger RA, Brady MF, Bookman MA, et al (2011) Incorporation of bevacizumab in the primary treatment of ovarian cancer. New England Journal of Medicine 365(26):2473–2483
  * Burges and Schmalfeldt [2011] Burges A, Schmalfeldt B (2011) Ovarian cancer: diagnosis and treatment. Deutsches Ärzteblatt International 108(38):635
  * Chen et al [2024] Chen RJ, Ding T, Lu MY, et al (2024) Towards a general-purpose foundation model for computational pathology. Nature Medicine 30(3):850–862
  * Cox [1972] Cox DR (1972) Regression models and life-tables. Journal of the Royal Statistical Society: Series B (Methodological) 34(2):187–202
  * Dosovitskiy et al [2020] Dosovitskiy A, Beyer L, Kolesnikov A, et al (2020) An image is worth 16x16 words: Transformers for image recognition at scale. arXiv preprint arXiv:201011929
  * Filiot et al [2023] Filiot A, Ghermi R, Olivier A, et al (2023) Scaling self-supervised learning for histopathology with masked image modeling. medRxiv pp 2023–07
  * Gadermayr and Tschuchnig [2024] Gadermayr M, Tschuchnig M (2024) Multiple instance learning for digital pathology: A review of the state-of-the-art, limitations & future potential. Computerized Medical Imaging and Graphics p 102337
  * Garcia et al [2020] Garcia J, Hurwitz HI, Sandler AB, et al (2020) Bevacizumab (avastin®) in cancer treatment: A review of 15 years of clinical experience and future outlook. Cancer treatment reviews 86:102017
  * Harrell Jr et al [1996] Harrell Jr FE, Lee KL, Mark DB (1996) Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors. Statistics in medicine 15(4):361–387
  * He et al [2016] He K, Zhang X, Ren S, et al (2016) Deep residual learning for image recognition. In: Proceedings of the IEEE conference on computer vision and pattern recognition, pp 770–778
  * Huang et al [2017] Huang G, Liu Z, Van Der Maaten L, et al (2017) Densely connected convolutional networks. In: Proceedings of the IEEE conference on computer vision and pattern recognition, pp 4700–4708
  * Huang et al [2022] Huang J, Chan WC, Ngai CH, et al (2022) Worldwide burden, risk factors, and temporal trends of ovarian cancer: a global study. Cancers 14(9):2230
  * Huang et al [2023] Huang Z, Bianchi F, Yuksekgonul M, et al (2023) A visual–language foundation model for pathology image analysis using medical twitter. Nature medicine 29(9):2307–2316
  * Ilse et al [2018] Ilse M, Tomczak J, Welling M (2018) Attention-based deep multiple instance learning. In: International conference on machine learning, PMLR, pp 2127–2136
  * Kang et al [2023] Kang M, Song H, Park S, et al (2023) Benchmarking self-supervised learning on diverse pathology datasets. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp 3344–3354
  * Kim et al [2021] Kim YJ, Jang H, Lee K, et al (2021) Paip 2019: Liver cancer segmentation challenge. Medical image analysis 67:101854
  * Kingma and Ba [2014] Kingma DP, Ba J (2014) Adam: A method for stochastic optimization. arXiv preprint arXiv:14126980
  * Liu et al [2024] Liu P, Ji L, Ye F, et al (2024) Advmil: Adversarial multiple instance learning for the survival analysis on whole-slide images. Medical Image Analysis 91:103020
  * Liu et al [2021] Liu Z, Lin Y, Cao Y, et al (2021) Swin transformer: Hierarchical vision transformer using shifted windows. In: Proceedings of the IEEE/CVF international conference on computer vision, pp 10012–10022
  * Liu et al [2022] Liu Z, Mao H, Wu CY, et al (2022) A convnet for the 2020s. In: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pp 11976–11986
  * Lu et al [2021] Lu MY, Williamson DF, Chen TY, et al (2021) Data-efficient and weakly supervised computational pathology on whole-slide images. Nature biomedical engineering 5(6):555–570
  * Macenko et al [2009] Macenko M, Niethammer M, Marron JS, et al (2009) A method for normalizing histology slides for quantitative analysis. In: 2009 IEEE international symposium on biomedical imaging: from nano to macro, IEEE, pp 1107–1110
  * Monk et al [2013] Monk BJ, Pujade-Lauraine E, Burger R (2013) Integrating bevacizumab into the management of epithelial ovarian cancer: the controversy of front-line versus recurrent disease. Annals of Oncology 24:x53–x58
  * Perren et al [2011] Perren TJ, Swart AM, Pfisterer J, et al (2011) A phase 3 trial of bevacizumab in ovarian cancer. New England Journal of Medicine 365(26):2484–2496
  * Radford et al [2021] Radford A, Kim JW, Hallacy C, et al (2021) Learning transferable visual models from natural language supervision. In: International conference on machine learning, PMLR, pp 8748–8763
  * Riasatian et al [2021] Riasatian A, Babaie M, Maleki D, et al (2021) Fine-tuning and training of densenet for histopathology image representation using tcga diagnostic slides. Medical image analysis 70:102032
  * Schirris et al [2022] Schirris Y, Gavves E, Nederlof I, et al (2022) Deepsmile: Contrastive self-supervised pre-training benefits msi and hrd classification directly from h&e whole-slide images in colorectal and breast cancer. Medical image analysis 79:102464
  * Shao et al [2021] Shao Z, Bian H, Chen Y, et al (2021) Transmil: Transformer based correlated multiple instance learning for whole slide image classification. Advances in neural information processing systems 34:2136–2147
  * Su et al [2022] Su Z, Tavolara TE, Carreno-Galeano G, et al (2022) Attention2majority: Weak multiple instance learning for regenerative kidney grading on whole slide images. Medical Image Analysis 79:102462
  * Torre et al [2018] Torre LA, Trabert B, DeSantis CE, et al (2018) Ovarian cancer statistics, 2018. CA: a cancer journal for clinicians 68(4):284–296
  * Vorontsov et al [2023] Vorontsov E, Bozkurt A, Casson A, et al (2023) Virchow: A million-slide digital pathology foundation model. arXiv preprint arXiv:230907778
  * Wang et al [2022a] Wang CW, Chang CC, Khalil MA, et al (2022a) Histopathological whole slide image dataset for classification of treatment effectiveness to ovarian cancer. Scientific Data 9(1):25
  * Wang et al [2022b] Wang CW, Chang CC, Lee YC, et al (2022b) Weakly supervised deep learning for prediction of treatment effectiveness on ovarian cancer from histopathology images. Computerized Medical Imaging and Graphics 99:102093
  * Wang et al [2022c] Wang CW, Lee YC, Chang CC, et al (2022c) A weakly supervised deep learning method for guiding ovarian cancer treatment and identifying an effective biomarker. Cancers 14(7):1651
  * Wang et al [2023] Wang CW, Lee YC, Lin YJ, et al (2023) Ensemble biomarkers for guiding anti-angiogenesis therapy for ovarian cancer using deep learning. Clinical and Translational Medicine 13(1)
  * Wang et al [2022d] Wang X, Yang S, Zhang J, et al (2022d) Transformer-based unsupervised contrastive learning for histopathological image classification. Medical image analysis 81:102559
  * Weinstein et al [2013] Weinstein JN, Collisson EA, Mills GB, et al (2013) The cancer genome atlas pan-cancer analysis project. Nature genetics 45(10):1113–1120
  * Yao et al [2020] Yao J, Zhu X, Jonnagaddala J, et al (2020) Whole slide images based cancer survival prediction using attention guided deep multiple instance learning networks. Medical Image Analysis 65:101789


Study of anomalous gauge-Higgs couplings using $Z$ boson polarization at LHC

Kumar Raoa, Saurabh D. Rindanib and Priyanka Sarmaha
aPhysics Department, Indian Institute of Technology Bombay,
Powai, Mumbai 400076, India
bTheoretical Physics Division, Physical Research Laboratory,
Navrangpura, Ahmedabad 380009, India

Abstract

We estimate model independent bounds that could be obtained on the anomalous
$ZZH$ vertex using polarization parameters of the $Z$ boson produced in the
Higgstrahlung process at the LHC. We calculate the eight independent
polarization parameters from the spin density matrix elements of the $Z$,
which can probe underlying new physics contributions to $ZH$ production. By
using the approach that connects these polarization observables to the
coefficients in the angular distribution of the decay products of the $Z$, we
estimate the limits on the anomalous $ZZH$ coupling that can be obtained at
the 14 TeV LHC.

##  1 Introduction

In the absence of evidence so far of any definitive beyond the Standard Model
(SM) physics at the Large Hadron Collider (LHC), it becomes important to probe
with high precision the properties of the 125 GeV Higgs ($H$) at the planned
high luminosity phase of the LHC (HL-LHC). This requires precise measurements
of the couplings of the Higgs to electroweak gauge bosons
($V=W^{\pm},Z,\gamma$), its Yukawa couplings to the fermions as well as its
self-couplings. Of these, the $VVH$ couplings, whose form is fixed by the
$SU(2)_{L}\times U(1)_{Y}$ gauge structure of the SM have a particular
importance. Although the present scenario indicates that the couplings of the
Higgs boson are in good agreement with the SM predictions, one would need more
accurate measurements to further constrain the couplings or to see a small
deviation from the SM predictions which could be a hint towards some
underlying new physics. This will require one to go beyond usual observables
like cross sections and differential rates which will be possible with higher
statistics at the HL-LHC.

A large amount of work has been carried out on probing the structure of the
$VVH$ couplings at the LHC and at planned $e^{+}e^{-}$ colliders [1, 2, 3, 4,
5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]. These studies have probed the
most general tensorial form of the $VVH$ coupling by using a variety of
observables involving kinematic distributions of the $Z$ and the charged
leptons from $Z$ decay. Study of Higgs-gauge coupling in the effective field
theory framework at the LHC has been studied in [18, 19] and at future
$e^{+}e^{-}$ colliders in [20, 21, 12].

In this paper, we propose studying the $ZZH$ coupling by making use of the
spin observables of the $Z$ boson. We study the $ZZH$ coupling using the
associated production of the $Z$ with the Higgs at the LHC. The formalism used
connects various angular asymmetries of the decay products of the $Z$ to its
eight independent polarization parameters extracted from the $Z$ production
spin density matrix[23, 22]. With the help of these parameters, we estimate
limits on the anomalous couplings. $Z$ polarization has been studied in the
context of new physics at the LHC [22, 18] and at an $e^{-}e^{+}$ collider
[24, 25, 26]. Analogously, polarization of the $W$ boson produced in
association with the Higgs at the LHC has been studied in [27, 18].

The main significance of our work is that we use completely analytical
expressions for the matrix element for the production and decay of polarized
$Z$ at the partonic level. Hence the angular asymmetries that we calculate
involve no numerical calculations at the partonic level. Only the integrations
over parton distributions have to be done numerically.

Our analytical approach has some overlap with that employed in [18]. However,
we have estimated limits on anomalous couplings that would be obtained at the
HL-LHC, which has not been done in [18].

We consider the process $pp\rightarrow ZHX$, where the vertex
$Z_{\mu}(k_{1})\rightarrow Z_{\nu}(k_{2})H$ has the Lorentz structure

$\Gamma^{V}_{\mu\nu}=\frac{g}{\cos\theta_{W}}m_{Z}\left[a_{Z}g_{\mu\nu}+\frac{b_{Z}}{m_{Z}^{2}}\left(k_{1\nu}k_{2\mu}-g_{\mu\nu}k_{1}.k_{2}\right)+\frac{\tilde{b}_{Z}}{m_{Z}^{2}}\epsilon_{\mu\nu\alpha\beta}k_{1}^{\alpha}k_{2}^{\beta}\right]$
(1)

where $g$ is the $SU(2)_{L}$ coupling and $\theta_{W}$ is the weak mixing
angle. $a_{Z}$ and $b_{Z}$ are invariant under CP, while $\tilde{b}_{Z}$
corresponds to CP violating term in the Lagrangian. In the SM, at tree level,
the coupling $a_{Z}=1$, whereas the other two couplings $b_{Z}$,
$\tilde{b}_{Z}$ vanish. These vanishing couplings are the anomalous couplings
which could arise from loop corrections in the SM or in any extension of SM
with some new particles or interactions. However, we are not concerned with
the predictions of any specific model here and derive the helicity amplitudes
for the process of our interest in a model-independent way using the general
form of the $ZZH$ vertex in Eqn. (1).

The current experimental bound on the $ZZH$ anomalous couplings is obtained by
the CMS collaboration[29, 28]. Although the current data are consistent with
the SM predictions, the constraints are still weak enough to allow for beyond
the SM contributions to the vertex. The $68\%$ confidence level (CL) upper
bounds on the $ZZH$ couplings, assuming them to be real, in our notation
translate to $|{\rm Re}~{}b_{Z}|<0.058$ and $|{\rm
Re}~{}\tilde{b}_{Z}|<0.078$. These limits are obtained from measurements of
ratios of the cross section contributions arising from the different $ZZH$
couplings. Ref[30] obtains possible bounds on the anomalous $ZZH$ coupling at
CLIC. For example, the $95\%$ CL limits obtained are $-0.118<b_{Z}<0.041$ and
$-0.096<\tilde{b}_{Z}<0.096$ at 3 TeV centre of mass energy (c.m.) and 1000
$\text{fb}^{-1}$ integrated luminosity, neglecting systematic uncertainties.
The possibility of a future Large Hadron electron Collider (LHeC) to probe
anomalous $ZZH$ couplings has been studied in [31], where weak limits are
found, viz., $-0.21<b_{Z}<0.43$ and $-0.32<\tilde{b}_{Z}<0.32$ for an electron
beam energy of 60 GeV and mild improvement for a beam energy of 140 GeV, with
proton beam energy of 7 TeV in either case.

##  2 $Z$ Polarization as a Probe

We consider the process $pp\rightarrow ZHX$ at the LHC, which at the partonic
level proceeds via the process

$q(p_{1})+\bar{q}(p_{2})\rightarrow Z^{\alpha}(p)+H(k)$ (2)

through $s$-channel $Z$ exchange. Here $q$ stands for both up type and down
type quarks of any generation, in the massless limit of the initial particles,
with the $ZZH$ vertex given in Eqn.(1). We first compute the helicity
amplitudes for this process considering the following representations for the
transverse and longitudinal polarization vectors of the $Z$:

$\epsilon^{\mu}(p,\pm)=\frac{1}{\sqrt{2}}(0,\mp\cos\theta,-i,\pm\sin\theta),$
(3)
$\epsilon^{\mu}(p,0)=\frac{1}{m_{Z}}(|\vec{p}_{Z}|,E_{Z}\sin\theta,0,E_{Z}\cos\theta),$
(4)

where $E_{Z}$ and $\vec{p}_{Z}$ are the energy and momentum of the $Z$
respectively, with $\theta$ being the polar angle made by the $Z$ with respect
to the quark momentum taken to be along the positive $z$ axis.

The non-zero helicity amplitudes in the limit of massless initial states and
assuming the SM value $a_{Z}=1$ are

$\displaystyle M(-,+,+)$ $\displaystyle=$
$\displaystyle\frac{g^{2}m_{Z}\sqrt{\hat{s}}(c_{V}+c_{A})}{2\sqrt{2}\cos^{2}\theta_{W}(\hat{s}-m_{Z}^{2})}\left[1-\frac{\sqrt{\hat{s}}}{m_{Z}^{2}}(E_{Z}b_{Z}+i\tilde{b}_{Z}|\vec{p}_{Z}|)\right]$
$\displaystyle\times(1-\cos\theta)$ $\displaystyle M(-,+,-)$ $\displaystyle=$
$\displaystyle\frac{g^{2}m_{Z}\sqrt{\hat{s}}(c_{V}+c_{A})}{2\sqrt{2}\cos^{2}\theta_{W}(\hat{s}-m_{Z}^{2})}\left[1-\frac{\sqrt{\hat{s}}}{m_{Z}^{2}}(E_{Z}b_{Z}-i\tilde{b}_{Z}|\vec{p}_{Z}|)\right]$
$\displaystyle\times(1+\cos\theta)$ $\displaystyle M(-,+,0)$ $\displaystyle=$
$\displaystyle\frac{g^{2}\sqrt{\hat{s}}(c_{V}+c_{A})}{2\cos^{2}\theta_{W}(\hat{s}-m_{Z}^{2})}\left[E_{Z}-\sqrt{\hat{s}}b_{Z}\right]\sin\theta$
(7) $\displaystyle M(+,-,+)$ $\displaystyle=$
$\displaystyle\frac{-g^{2}m_{Z}\sqrt{\hat{s}}(c_{V}-c_{A})}{2\sqrt{2}\cos^{2}\theta_{W}(\hat{s}-m_{Z}^{2})}\left[1-\frac{\sqrt{\hat{s}}}{m_{Z}^{2}}(E_{Z}b_{Z}+i\tilde{b}_{Z}|\vec{p}_{Z}|)\right]$
$\displaystyle\times(1+\cos\theta)$ $\displaystyle M(+,-,-)$ $\displaystyle=$
$\displaystyle\frac{-g^{2}m_{Z}\sqrt{\hat{s}}(c_{V}-c_{A})}{2\sqrt{2}\cos^{2}\theta_{W}(\hat{s}-m_{Z}^{2})}\left[1-\frac{\sqrt{\hat{s}}}{m_{Z}^{2}}(E_{Z}b_{Z}-i\tilde{b}_{Z}|\vec{p}_{Z}|)\right]$
$\displaystyle\times(1-\cos\theta)$ $\displaystyle M(+,-,0)$ $\displaystyle=$
$\displaystyle\frac{g^{2}\sqrt{\hat{s}}(c_{V}-c_{A})}{2\cos^{2}\theta_{W}(\hat{s}-m_{Z}^{2})}\left[E_{Z}-\sqrt{\hat{s}}b_{Z}\right]\sin\theta$
(10)

Here the first two entries in $M$ denote the signs of the helicities of the
quark and antiquark respectively and the third entry is the $Z$ helicity.
$\sqrt{\hat{s}}$ is the partonic c.m. energy, and $c_{V}$ and $c_{A}$ are the
respective vector and axial vector couplings of the relevant quark to the $Z$.

The $a_{Z}$ dependence can be easily recovered by multiplying the helicity
amplitude expressions by $a_{Z}$, and then replacing $b_{Z}$ and
$\tilde{b}_{Z}$ by $b_{Z}/a_{Z}$ and $\tilde{b}_{Z}/a_{Z}$, respectively.

We evaluate the elements of the spin-density matrix for $Z$ production, which
can be expressed in terms of the helicity amplitudes as follows

$\rho(i,j)=\overline{\sum}_{\lambda,\lambda^{{}^{\prime}}}M(\lambda,\lambda^{{}^{\prime}},i)M^{\ast}(\lambda,\lambda^{{}^{\prime}},j)$
(11)

the average being over the initial helicities $\lambda$,
$\lambda^{{}^{\prime}}$ of the quark and antiquark respectively and also over
the initial color states. The $Z$ helicity indices $i,j$ can take values
$\pm,0$ and with $i=j$ corresponding to the diagonal elements of Eqn.(11)
which are the squared matrix elements for $Z$ production with definite
polarization. It is known that a complete information of the state of
polarization is encoded in all the density matrix elements. So to attain
maximum possible information, it is necessary to study the full density matrix
description, which also includes the off diagonal elements. The density matrix
elements, for $q\bar{q}\to ZH$, derived from the helicity amplitudes are given
by

$\displaystyle\rho(\pm,\pm)$ $\displaystyle=$
$\displaystyle\frac{g^{4}m^{2}_{Z}s}{8\cos^{4}\theta_{W}(\hat{s}-m_{Z}^{2})^{2}}\left[(c_{V}+c_{A})^{2}(1\mp\cos\theta)^{2}\right.$
(12)
$\displaystyle\hskip-22.76228pt\left.+(c_{V}-c_{A})^{2}(1\pm\cos\theta)^{2}\right]\left[1-2(\text{Re}~{}b_{Z}\mp\beta_{Z}\text{Im}~{}\tilde{b}_{Z})\frac{E_{Z}\sqrt{\hat{s}}}{m^{2}_{Z}}\right.$
$\displaystyle\hskip-22.76228pt\left.+\frac{E_{Z}^{2}\hat{s}}{m^{4}_{Z}}|b_{Z}|^{2}\mp\frac{2E_{Z}P_{Z}\hat{s}}{m^{4}_{Z}}(\text{Im}~{}\tilde{b}_{Z}~{}{\rm
Re}~{}b_{Z}-\text{Im}~{}b_{Z}~{}{\rm Re}~{}\tilde{b}_{Z})\right.$
$\displaystyle\hskip-22.76228pt\left.+\frac{P_{Z}^{2}\hat{s}}{m^{4}_{Z}}|\tilde{b}_{Z}|^{2}\right]$
$\displaystyle\rho(0,0)$ $\displaystyle=$
$\displaystyle\frac{g^{4}E^{2}_{Z}s}{2\cos^{4}\theta_{W}(s-m_{Z}^{2})^{2}}\sin^{2}\theta\,(c_{V}^{2}+c_{A}^{2})\left[1-2{\rm
Re}~{}b_{Z}\frac{\sqrt{s}}{E_{Z}}\right.$ (13)
$\displaystyle\hskip-22.76228pt\left.+\frac{\hat{s}}{E^{2}_{Z}}|b_{Z}|^{2}\right]$
$\displaystyle\rho(\pm,\mp)$ $\displaystyle=$
$\displaystyle\frac{g^{4}m^{2}_{Z}s}{4\cos^{4}\theta_{W}(\hat{s}-m_{Z}^{2})^{2}}\sin^{2}\theta\,(c_{V}^{2}+c_{A}^{2})$
(14) $\displaystyle\times\left[1-2({\rm Re}~{}b_{Z}\pm i\beta_{Z}{\rm
Re}~{}\tilde{b}_{Z})\frac{E_{Z}\sqrt{s}}{m^{2}_{Z}}+\frac{E_{Z}^{2}\hat{s}}{m^{4}_{Z}}|b_{Z}|^{2}\pm
i\frac{2E_{Z}P_{Z}\hat{s}}{m^{4}_{Z}}\right.$
$\displaystyle\hskip-22.76228pt\left.(\text{Im}~{}\tilde{b}_{Z}~{}\text{Im}~{}b_{Z}+{\rm
Re}~{}b_{Z}~{}{\rm
Re}~{}\tilde{b}_{Z})-\frac{2P_{Z}^{2}\hat{s}}{m^{4}_{Z}}|\tilde{b}_{Z}|^{2}\right]$
$\displaystyle\rho(\pm,0)$ $\displaystyle=$
$\displaystyle\frac{g^{4}m_{Z}E_{Z}s}{4\sqrt{2}\cos^{4}\theta_{W}(\hat{s}-m_{Z}^{2})^{2}}\sin\theta$
(15)
$\displaystyle\hskip-22.76228pt\times\left[(c_{V}+c_{A})^{2}(1\mp\cos\theta)-(c_{V}-c_{A})^{2}(1\pm\cos\theta)\right]$
$\displaystyle\hskip-22.76228pt\times\left[1-{\rm
Re}~{}b_{Z}\sqrt{\hat{s}}\frac{(E^{2}_{Z}+m^{2}_{Z})}{E_{Z}m^{2}_{Z}}-i\sqrt{\hat{s}}\frac{E_{Z}}{m^{2}_{Z}}\left({\rm
Im}~{}b_{Z}~{}\beta^{2}_{Z}\pm\tilde{b}_{Z}\beta_{Z}\right)\right.$
$\displaystyle\hskip-22.76228pt\left.\mp\frac{\hat{s}}{m^{2}_{Z}}|b_{Z}|^{2}\pm\frac{\hat{s}P_{Z}}{m^{2}_{Z}E_{Z}}({\rm
Im}~{}b_{Z}+i{\rm Re}~{}b_{Z})({\rm
Re}~{}\tilde{b}_{Z}+i\text{Im}~{}\tilde{b}_{Z})\right]$

where $\beta_{Z}=|\vec{p}_{Z}|/E_{Z}$ is the velocity of the $Z$ in the c.m
frame. The analytical manipulation software FORM [35] has been used to verify
these expressions. We have kept the finite $Z$ width in our numerical
calculations later.

The full density matrix Eqn.(11) on integrating over an appropriate kinematic
range, can be parametrized in terms of the 3 components of the vector
polarization $\vec{P}$ and 5 components of the tensor polarization $T$ of the
$Z$ boson[36]. Defining this as $\sigma(i,j)$ we have

$\sigma(i,j)\equiv\sigma\;\left(\begin{array}[]{ccc}\frac{1}{3}+\frac{P_{z}}{2}+\frac{T_{zz}}{\sqrt{6}}&\frac{P_{x}-iP_{y}}{2\sqrt{2}}+\frac{T_{xz}-iT_{yz}}{\sqrt{3}}&\frac{T_{xx}-T_{yy}-2iT_{xy}}{\sqrt{6}}\\\
\frac{P_{x}+iP_{y}}{2\sqrt{2}}+\frac{T_{xz}+iT_{yz}}{\sqrt{3}}&\frac{1}{3}-\frac{2T_{zz}}{\sqrt{6}}&\frac{P_{x}-iP_{y}}{2\sqrt{2}}-\frac{T_{xz}-iT_{yz}}{\sqrt{3}}\\\
\frac{T_{xx}-T_{yy}+2iT_{xy}}{\sqrt{6}}&\frac{P_{x}+iP_{y}}{2\sqrt{2}}-\frac{T_{xz}+iT_{yz}}{\sqrt{3}}&\frac{1}{3}-\frac{P_{z}}{2}+\frac{T_{zz}}{\sqrt{6}}\end{array}\right)$
(16)

where $\sigma$ is the production cross section,

$\sigma=\sigma(+,+)+\sigma(-,-)+\sigma(0,0).$ (17)

The eight independent vector and tensor polarization observables of the $Z$
can then be constructed using appropriate linear combinations of the
integrated density matrix elements of Eqn.(16):

$\displaystyle P_{x}$ $\displaystyle=$
$\displaystyle\frac{\\{\sigma(+,0)+\sigma(0,+)\\}+\\{\sigma(0,-)+\sigma(-,0)\\}}{\sqrt{2}\sigma}$
(18) $\displaystyle P_{y}$ $\displaystyle=$
$\displaystyle\frac{-i\\{[\sigma(0,+)-\sigma(+,0)]+[\sigma(-,0)-\sigma(0,-)]\\}}{\sqrt{2}\sigma}$
(19) $\displaystyle P_{Z}$ $\displaystyle=$
$\displaystyle\frac{[\sigma(+,+)]-[\sigma(-,-)]}{\sigma}$ (20) $\displaystyle
T_{xy}$ $\displaystyle=$
$\displaystyle\frac{-i\sqrt{6}[\sigma(-,+)-\sigma(+,-)]}{4\sigma}$ (21)
$\displaystyle T_{xz}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{3}\\{[\sigma(+,0)+\sigma(0,+)]-[\sigma(0,-)+\sigma(-,0)]\\}}{4\sigma}$
(22) $\displaystyle T_{yz}$ $\displaystyle=$
$\displaystyle\frac{-i\sqrt{3}\\{[\sigma(0,+)-\sigma(+,0)]-[\sigma(-,0)-\sigma(0,-)]\\}}{4\sigma}$
(23) $\displaystyle T_{xx}-T_{yy}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{6}[\sigma(-,+)+\sigma(+,-)]}{2\sigma}$ (24)
$\displaystyle T_{zz}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{6}}{2}\left\\{\frac{[\sigma(+,+)]+[\sigma(-,-)]}{\sigma}-\frac{2}{3}\right\\}$
(25) $\displaystyle=$
$\displaystyle\frac{\sqrt{6}}{2}\left[\frac{1}{3}-\frac{\sigma(0,0)}{\sigma}\right]$

Of these $P_{x}$, $P_{y}$ and $P_{z}$ are the vector polarizations, whereas
the $T$’s are the tensor polarizations, with the constraint that the tensor is
traceless. In real experiments where the $Z$ boson decays to two leptons,
these polarization observables can be extracted from kinematic distributions
of its decay products. Angular asymmetries can be obtained by combining the
relevant production-level density matrix elements with appropriate decay
density matrix elements and integrating over the appropriate phase space. For
example, $P_{x}$ can be calculated from the asymmetry $A_{x}$ defined by-

$A_{x}=\frac{3\alpha
P_{x}}{4}\equiv\frac{\sigma(\cos\phi^{\ast}>0)-\sigma(\cos\phi^{\ast}<0)}{\sigma(\cos\phi^{\ast}>0)+\sigma(\cos\phi^{\ast}<0)}$
(26)

where, $\alpha$ is the $Z$ boson polarization analyzer, given in terms of its
vector and axial vector couplings to charged leptons $\ell$, $c_{V}^{\ell}$
and $c_{A}^{\ell}$ respectively, as

$\alpha=-\frac{2c_{V}^{\ell}c_{A}^{\ell}}{{c_{V}^{\ell}}^{2}+{c_{A}^{\ell}}^{2}}$
(27)

The angles $\theta^{\ast}$ and $\phi^{\ast}$ are polar and azimuthal angles of
the lepton in the rest frame of the $Z$. The $Z$ rest frame is reached by a
combination of boosts and rotations from the laboratory frame. In the
laboratory frame, the quark momentum defines the positive $z$ axis, and the
production plane of $Z$ is defined as the $xz$ plane. While boosting to the
$Z$ rest frame, the $xz$ plane is kept unchanged. Then, the angles
$\theta^{\ast}$ and $\phi^{\ast}$ are measured with respect to the would-be
momentum of the $Z$. Similarly, expressions for other asymmetries, viz.,
$A_{y}$, $A_{z}$, $A_{xy}$, $A_{yz}$, $A_{xz}$, $A_{x^{2}-y^{2}}$, $A_{zz}$
corresponding to the 2 vector polarizations $P_{y}$, $P_{z}$ and 5 tensor
polarizations $T_{ij}(i,j=x,y,z)$ can be obtained and are listed in [25, 26].

It is observed that the density matrix elements $\sigma(\pm,0)$ and
$\sigma(0,\pm)$ and the asymmetries involving these elements $A_{x}$, $A_{y}$,
$A_{xz}$, $A_{yz}$ vanish due to the symmetric nature of the LHC, which does
not allow a unique definition of a positive $z$-axis. Therefore to make them
non-zero, we define the direction of the reconstructed momentum of the $ZH$
combination as the positive $z$-axis.

We have evaluated these asymmetries upto quadratic order in the anomalous
couplings. It is observed that out of eight polarization asymmetries only
three, viz., $A_{x}$, $A_{x^{2}-y^{2}}$ and $A_{zz}$ are non-zero in the SM,
which, along with the total cross section, would be proportional to the real
part of the anomalous couplings (upto linear order) or absolute square of the
couplings to satisfy the CPT theorem. This will be seen in the following
section.

##  3 Limits on the Anomalous Couplings

Here we present numerical values for the integrated density matrix elements,
the corresponding asymmetries and the sensitivities of the asymmetries to the
various anomalous couplings. We consider c.m. energy
$\sqrt{s}=14\text{~{}TeV}$, with integrated luminosity
$\int\mathcal{L}dt=1000$ $\text{fb}^{-1}$. In our numerical calculations, we
employ MMHT2014 parton distribution functions [37] with factorization scale
chosen as the square root of the partonic c.m. energy. The integrated
production density matrix elements in the units of fb are

$\displaystyle\sigma(\pm,\pm)$ $\displaystyle=$ $\displaystyle
161.95-1495.62~{}\text{Re}~{}b_{Z}\pm
1036.98~{}\text{Im}~{}\tilde{b}_{Z}+5391.21\,|b_{Z}|^{2}$ (28)
$\displaystyle\hskip-22.76228pt+3753.23\,|\tilde{b}_{Z}|^{2}\mp
8811.36~{}(\text{Im}~{}\tilde{b}_{Z}~{}\text{Re}~{}b_{Z}-~{}\text{Im}~{}b_{Z}~{}\text{Re}~{}\tilde{b}_{Z})$
$\displaystyle\sigma(0,0)$ $\displaystyle=$ $\displaystyle
341.976-1495.62~{}\text{Re}~{}b_{Z}+1637.98\,|b_{Z}|^{2}$ (29)
$\displaystyle\sigma(\pm,\mp)$ $\displaystyle=$ $\displaystyle
80.97-747.81~{}\text{Re}~{}b_{Z}\mp
i518.49~{}\text{Re}~{}\tilde{b}_{Z}+2695.6|b_{Z}|^{2}$ (30)
$\displaystyle\hskip-22.76228pt-1876.61\,|\tilde{b}_{Z}|^{2}\pm
i4405.67~{}(\text{Im}~{}\tilde{b}_{Z}~{}\text{Im}~{}b_{Z}+~{}\text{Re}~{}\tilde{b}_{Z}~{}\text{Re}~{}b_{Z})$
$\displaystyle\sigma(\pm,0)$ $\displaystyle=$ $\displaystyle
59.59-474.46~{}\text{Re}~{}b_{Z}-i211.22~{}\text{Im}~{}b_{Z}\mp
261.88~{}(i~{}\text{Re}~{}\tilde{b}_{Z}-\text{Im}~{}\tilde{b}_{Z})$ (31)
$\displaystyle\hskip-22.76228pt+738.07\,|b_{Z}|^{2}\pm
558.15~{}\tilde{b}_{Z}~{}(\text{Im}~{}b_{Z}+i~{}\text{Re}~{}b_{Z})$
$\displaystyle\sigma(0,\pm)$ $\displaystyle=$ $\displaystyle
59.59-474.46~{}\text{Re}~{}b_{Z}+i211.22~{}\text{Im}~{}b_{Z}\pm
261.88~{}(i~{}\text{Re}~{}\tilde{b}_{Z}+\text{Im}~{}\tilde{b}_{Z})$ (32)
$\displaystyle\hskip-22.76228pt+738.07\,|b_{Z}|^{2}\pm
558.15~{}\tilde{b}_{Z}~{}(\text{Im}~{}b_{Z}+i~{}\text{Re}~{}b_{Z})$

The leptonic asymmetries corresponding to different polarizations calculated
upto linear order in the anomalous couplings are given by

$A_{x}=0.035~{}\text{Re}~{}b_{Z}-0.028$ (33)

$A_{y}=-0.125~{}\text{Re}~{}\tilde{b}_{Z}$ (34)

$A_{z}=-0.349~{}\text{Im}~{}\tilde{b}_{Z}$ (35)

$A_{xy}=0.496~{}\text{Re}~{}\tilde{b}_{Z}$ (36)

$A_{xz}=-0.354~{}\text{Im}~{}\tilde{b}_{Z}$ (37)

$A_{yz}=0.286~{}\text{Im}~{}b_{Z}$ (38)

$A_{x^{2}-y^{2}}=-0.193~{}\text{Re}~{}b_{Z}+0.077$ (39)

$A_{zz}=-0.683~{}\text{Re}~{}b_{Z}-0.101$ (40)

It is observed that all asymmetries except $A_{x}$, $A_{x^{2}-y^{2}}$ and
$A_{zz}$ vanish in the SM and the reason for this is the CP even and T even
nature of the asymmetries $A_{x},\,A_{x^{2}-y^{2}}$ and $A_{zz}$, because of
which they can occur at tree level in the SM. The remaining asymmetries vanish
in the SM because they are either CP even and T odd or CP odd, and hence
depend on the CP violating parameters which are absent in the SM at tree
level.

Observable | Coupling | Limit ($\times 10^{-3}$)
---|---|---
$\sigma$ | $\text{Re}~{}b_{Z}$ | $0.70$
$A_{x}$ | $\text{Re}~{}b_{Z}$ | $136$
$A_{y}$ | $\text{Re}~{}\tilde{b}_{Z}$ | $37.9$
$A_{z}$ | $\text{Im}~{}\tilde{b}_{Z}$ | $13.5$
$A_{xy}$ | $\text{Re}~{}\tilde{b}_{Z}$ | $9.53$
$A_{yz}$ | $\text{Im}~{}b_{Z}$ | $16.5$
$A_{xz}$ | $\text{Im}~{}\tilde{b}_{Z}$ | $13.3$
$A_{x^{2}-y^{2}}$ | $\text{Re}~{}b_{Z}$ | $24.4$
$A_{zz}$ | $\text{Re}~{}b_{Z}$ | $6.88$
Table 1: $1\sigma$ limit obtained on the anomalous couplings from cross
section and various leptonic asymmetries calculated upto linear order in the
couplings at $\sqrt{s}=14\text{~{}TeV}$.

Next we present the expressions for the total cross section and angular
asymmetries including quadratic terms of couplings at
$\sqrt{s}=14\text{~{}TeV}$.

$\sigma=0.067294~{}(7506.45~{}|\tilde{b}_{Z}|^{2}+12420.4~{}|b_{Z}|^{2}-4486.85~{}\text{Re}~{}b_{Z}+665.87)~{}\text{fb}$
(41)

$A_{x}=\frac{0.012~{}\text{Re}~{}b_{Z}-0.019~{}|b_{Z}|^{2}-0.002}{0.604~{}|\tilde{b}_{Z}|^{2}+|b_{Z}|^{2}-0.361~{}\text{Re}~{}b_{Z}+0.054}$
(42)

$A_{y}=\frac{0.024~{}\text{Im}~{}\tilde{b}_{Z}~{}\text{Im}~{}b_{Z}+\text{Re}~{}\tilde{b}_{Z}~{}(0.024~{}\text{Re}~{}b_{Z}-0.011)}{|\tilde{b}_{Z}|^{2}+1.655~{}|b_{Z}|^{2}+(\text{Re}~{}\tilde{b}_{Z})^{2}-0.598~{}\text{Re}~{}b_{Z}+0.089}$
(43)

$A_{z}=\frac{\text{Im}~{}\tilde{b}_{Z}~{}(1976.66~{}\text{Re}~{}b_{Z}-232.627)-1976.66~{}\text{Im}~{}b_{Z}~{}\text{Re}~{}\tilde{b}_{Z}}{7506.45~{}|\tilde{b}_{Z}|^{2}+12420.4~{}|b_{Z}|^{2}-4486.85~{}\text{Re}~{}b_{Z}+665.87}$
(44)

$A_{xy}=\frac{\text{Re}~{}\tilde{b}_{Z}~{}(0.044-0.374~{}\text{Re}~{}b_{Z})-0.374~{}\text{Im}~{}\tilde{b}_{Z}~{}\text{Im}~{}b_{Z}}{|\tilde{b}_{Z}|^{2}+1.655~{}|b_{Z}|^{2}-0.598~{}\text{Re}~{}b_{Z}+0.089}$
(45)

$A_{yz}=\frac{190.164~{}\text{Im}~{}b_{Z}}{7506.45~{}|\tilde{b}_{Z}|^{2}+12420.4~{}|b_{Z}|^{2}-4486.85~{}\text{Re}~{}b_{Z}+665.87}$
(46)

$A_{xz}=\frac{\text{Im}~{}\tilde{b}_{Z}~{}(0.0404~{}\text{Re}~{}b_{Z}-0.019)-0.0404~{}\text{Im}~{}b_{Z}~{}\text{Re}~{}\tilde{b}_{Z}}{0.604~{}|\tilde{b}_{Z}|^{2}+|b_{Z}|^{2}-0.361~{}\text{Re}~{}b_{Z}+0.054}$
(47)

$A_{x^{2}-y^{2}}=\frac{-0.297~{}|\tilde{b}_{Z}|^{2}+0.019~{}\text{Re}~{}b_{Z}-0.005}{|\tilde{b}_{Z}|^{2}+1.655~{}|b_{Z}|^{2}-0.598~{}\text{Re}~{}b_{Z}+0.089}+0.138$
(48)

$A_{zz}=\frac{0.074~{}|\tilde{b}_{Z}|^{2}+0.068~{}\text{Re}~{}b_{Z}-0.019}{|\tilde{b}_{Z}|^{2}+1.655~{}|b_{Z}|^{2}-0.598~{}\text{Re}~{}b_{Z}+0.089}+0.113$
(49)

We obtain the sensitivity of an observable $\mathcal{O}$ which depends on a
parameter $f$ from the definition

$\mathcal{S}(\mathcal{O}(f))=\frac{|\mathcal{O}(f)-\mathcal{O}(f=0)|}{\delta\mathcal{O}}$
(50)

where $\delta\mathcal{O}$ is the estimated error on the observable. For an
asymmetry, the estimated error takes the form

$\delta A=\frac{\sqrt{1-A^{2}_{SM}}}{\sqrt{\sigma_{SM}\mathcal{L}}}$ (51)

with $\sigma_{SM}$ being the SM cross section for the process $pp\to
Z^{*}H\to\ell\bar{\ell}H$ ($\ell=e,\mu$) at the LHC with integrated luminosity
$\mathcal{L}$ and $A_{SM}$ is the corresponding value of asymmetry in the SM.
Similarly, for the cross section, the error is given by

$\delta\sigma=\sqrt{\frac{\sigma_{\rm{SM}}}{\mathcal{L}}}$ (52)

We estimate the 1$\sigma$ limits calculated upto linear order and list it in
Table 1. We note that, among the four observables $\sigma$, $A_{x}$, $A_{zz}$
and $A_{x^{2}-y^{2}}$, which are sensitive to Re $b_{Z}$ , the total cross
section provides the best limits on the coupling. However, it is not
sufficient to consider the total cross sections as the only probe as it is
sensitive to just one coupling, Re $b_{Z}$ . So to explore the couplings
which do not appear in the total cross section, one will require the other
angular asymmetries. The better limit on Im $\tilde{b}_{Z}$ comes from $A_{x}$
and $A_{xz}$, both being equally sensitive to the coupling. For the coupling
Im $b_{Z}$ , the best bound comes from $A_{yz}$ whereas on Re
$\tilde{b}_{Z}$ the best limit of $9.53\times 10^{-3}$ is achieved from the
observable $A_{xy}$.

Figure 1: Sensitivities of cross section and asymmetries to anomalous couplings, including quadratic order at $\sqrt{s}=14$TeV. Plots are obtained by varying one coupling at a time. Observable | Coupling | Limit ($\times 10^{-3}$)
---|---|---
$\sigma$ | $|\text{Re}~{}b_{Z}|$ | $0.70$
$\sigma$ | $|\text{Im}~{}b_{Z}|$ | $15.9$
$A_{xy}$ | $|\text{Re}~{}\tilde{b}_{Z}|$ | $9.54$
$A_{xz}$, $A_{z}$ | $|\text{Im}~{}\tilde{b}_{Z}|$ | $13.3$
Table 2: The best $1\sigma$ limit on couplings and the corresponding
observables at $\sqrt{s}=14\text{~{}TeV}$, obtained from Figure 1.

In Figure 1, we plot the one parameter sensitivity i.e $\mathcal{S}=1$(or
$\Delta\chi^{2}=1$) for the cross section and the 8 asymmetries, considered
upto quadratic order in the anomalous couplings. It is observed from Figure 1
that the tightest limit on the coupling $\text{Re}~{}b_{Z}$ can be obtained
from total cross section. On the coupling $\text{Im}~{}b_{Z}$, both cross
section and $A_{yz}$ place comparable limits. The observables $A_{x}$ and
$A_{xz}$ are found to be almost equally sensitive to the coupling
$\text{Im}~{}\tilde{b}_{Z}$. The best limit on $\text{Re}~{}\tilde{b}_{Z}$ can
be obtained from the observable $A_{xy}$. In Table 2, we list the tightest
$1\sigma$ level limit on the couplings, obtained from Figure 1.

Figure 2: $1\sigma$ sensitivity contours for cross-section and asymmetries
obtained by varying two parameters simultaneously. The black dot in the middle
of the plots represents the SM value.

So far, we obtained a limit on each coupling assuming all the other couplings
to be zero. Ideally, we would like to place a limit on each coupling without
making any assumptions on the remaining couplings. In practice this would
involve making a simultaneous fit to several observables varying all the
couplings. This is not only cumbersome, it would also require a large data
set. We therefore now consider simultaneous limits which may be obtained by
selecting a pair of couplings non-vanishing, taking the remaining to be zero.

We vary two couplings at a time and obtain the $1\sigma$ sensitivity i.e
$\mathcal{S}=2.3$ (or $\Delta\chi^{2}=2.3$) contours shown in Figure 2 for
each observable. The black dot in the middle of the plots represents the SM
value.

A first general observation in the context of deriving limits from the
contours is that the total cross section $\sigma$, which makes use of all the
events, tends to be the most sensitive observable for measurement of all
couplings. While at linear order it depends only on Re $b_{Z}$, at the
quadratic order it depends on all the couplings. Thus, in most cases, the best
limit for all couplings is obtained from $\sigma$.

Another observation is that as Re $\tilde{b}_{Z}$ and Im $\tilde{b}_{Z}$ are
CP-odd couplings, they would occur linearly in CP-odd observables. Thus, even
though these CP-odd couplings could be constrained by the cross sections or
any of the asymmetries, they would get strongest limits from the CP-odd
asymmetries $A_{y}$, $A_{xy}$, $A_{z}$ and $A_{xz}$. Of these the first two
are CPT even and would therefore constrain Re $\tilde{b}_{Z}$ , whereas the
last two being CPT odd would constrain Im $\tilde{b}_{Z}$ .

Coming to simultaneous limits on two couplings which can be read off from the
contour plots, the best limits on the combination Re $b_{Z}$ and Im $b_{Z}$
come from $\sigma$ and $A_{yz}$, and are, respectively, $[-0.002,0.005]$ and
$[-0.035,0.035]$. If, however, $A_{zz}$ is used in place of $\sigma$, the
limit on Im $b_{Z}$ is very similar, but the limit on Re $b_{Z}$ becomes
weaker. All other asymmetries give weaker limits.

Taking the sensitivity contour plots of Re $b_{Z}$ versus Re $\tilde{b}_{Z}$
, the best limits are from $\sigma$ and $A_{xy}$, the latter being linear in
the CP-odd coupling Re $\tilde{b}_{Z}$ . These are, respectively,
$[-0.002,0.003]$ and $[-0.022,0.022]$. Again if $A_{zz}$ is used instead of
$\sigma$, the limit on Re $b_{Z}$ is weaker, viz., $[-0.02,0.02]$.

In case of the contour plots of Re $b_{Z}$ and Im $\tilde{b}_{Z}$ , the best
limits on Re $b_{Z}$ are as in the previous case, whereas the best limits on
Im $\tilde{b}_{Z}$ are $[-0.035,0.035]$, come from $\sigma$ and $A_{z}$,
$A_{xz}$, the latter two being numerically very close. A similar situation
holds in the case of the contour plots of Re $\tilde{b}_{Z}$ versus Im
$\tilde{b}_{Z}$ , where the best limits on Im $\tilde{b}_{Z}$ , viz.,
$[-0.03,0.03]$ come equally from $A_{z}$ and $A_{xz}$, the best limit on Re
$\tilde{b}_{Z}$ , obtained from $A_{xy}$ is $[-0.025,0.025]$.

The best limits on Re $\tilde{b}_{Z}$ and Im $b_{Z}$ are $[-0.02,0.02]$ and
$[-0.027,0.027]$ respectively, obtained from the combination of $A_{xy}$ and
$\sigma$. Lastly in the case of Im $b_{Z}$ versus Im $\tilde{b}_{Z}$ contour,
best limit on Im $b_{Z}$ is $[-0.025,0.025]$, comes from $\sigma$, whereas on
Im $\tilde{b}_{Z}$ , $A_{z}$ and $A_{xz}$ which contribute almost equally,
provide a stringent limit of $[-0.031,0.031]$ on it.

We see from the above a significant feature that using $\sigma$ as one of the
observables gives a stringent limit for all the couplings involved. We also
see that the best limits on Re $b_{Z}$ is of the order of $2-5\times 10^{-3}$
in magnitude from all relevant pairs of observables. This may be compared to
the limit $0.7\times 10^{-3}$ obtained when only Re $b_{Z}$ is taken as non-
zero, as seen from Tables 1 and 2\. Similarly, the best limits on Im $b_{Z}$
from simultaneous measurement of two observables varies between $25\times
10^{-3}$ and $35\times 10^{-3}$, as compared to the best individual limit of
around $16\times 10^{-3}$. For Re $\tilde{b}_{Z}$ the best simultaneous limits
are $20-27\times 10^{-3}$, the best individual limit being $9.5\times
10^{-3}$. Likewise, the best simultaneous limits on Im $\tilde{b}_{Z}$ vary
between $30\times 10^{-3}$ and $35\times 10^{-3}$, whereas the best individual
limit is around $13\times 10^{-3}$.

##  4 Conclusions and discussion

The measurement of couplings of the Higgs Boson to all other SM particles is
an essential test of the SM. In this work, we study the form and magnitude of
the tensor structure of the couplings of the Higgs boson to a pair of $Z$
bosons at the LHC with the help of the polarization observables of the $Z$. We
estimate sensitivities of these polarization observables by adopting the
formalism which connects angular asymmetries of charged leptons from $Z$ decay
to the polarization parameters of the $Z$. We first calculate the $Z$
polarization parameters using the spin density matrix elements evaluated at
production level and then obtain various angular asymmetries corresponding to
these parameters.

We have restricted ourselves to tree-level calculations. To our knowledge,
non-leading order (NLO) contributions to the process with polarized $Z$ have
not been calculated so far. However, we expect that asymmetries which we make
use of, being ratios of cross sections, will be less sensitive to NLO
corrections.

We see that the $1\sigma$ limits obtained on the real parts of the couplings
are of the order of a few times $10^{-3}$ and an order of magnitude higher for
the imaginary parts. We show that the LHC at c.m energy
$\sqrt{s}=14\text{~{}TeV}$ with integrated luminosity $\int\mathcal{L}dt=1000$
$\text{fb}^{-1}$ could provide a limit on the CP conserving couplings
$\text{Re}~{}b_{Z}$ in the interval $[-0.7,0.7]\times 10^{-3}$ and
$\text{Im}~{}b_{Z}$ in the interval $[-15.9,15.9]\times 10^{-3}$. Similarly
the CP violating couplings, $\text{Re}~{}\tilde{b}_{Z}$ and
$\text{Im}~{}\tilde{b}_{Z}$ get a best bound of
$|\text{Re}~{}\tilde{b}_{Z}|\leq 9.54\times 10^{-3}$ and
$|\text{Im}~{}\tilde{b}_{Z}|\leq 13.3\times 10^{-3}$ respectively. These
limits are obtained by varying one coupling at a time. With two non-zero
couplings, we observe a slight weakening of bounds on all the anomalous
couplings as can be expected.

We have not considered Higgs decays, which do not affect the polarization
parameters and asymmetries of the $Z$. The effect of Higgs decay on the
sensitivities can be estimated by multiplying the SM cross section by the
Higgs branching ratio and detection efficiencies in Eqns. (51) and (52).

Associated Higgs production with $V=W,Z$ and with $H$ decaying into $b\bar{b}$
and $V$ decaying to 0, 1 and 2 leptons has been observed by both ATLAS[32] and
CMS collaborations[33] at close to $5\sigma$ CL. Also, as shown in [34], a
measurement of (SM) $Z$ polarization parameters themselves can help suppress
backgrounds and enhance the signal sensitivity in the
$Z(\ell^{+}\ell^{-})H(b\bar{b})$ decay mode. A full scale analysis using an
event generator coupled with all appropriate cuts and detection efficiencies
relevant to the decay channels of the $Z$ and Higgs, as used in [32, 33, 34]
with be able to refine the actual sensitivities that we have obtained for all
the eight BSM polarization asymmetries.

Acknowledgement: SDR acknowledges partial support from the Department of
Science and Technology, India, under the J.C. Bose National Fellowship
programme, Grant No. SR/SB/JCB-42/2009, and from the Indian National Science
Academy, New Delhi, under the Senior Scientist Programme.

## References

  * [1] K. Hagiwara and M. L. Stong, Z. Phys. C 62, 99 (1994) doi:10.1007/BF01559529 [hep-ph/9309248].
  * [2] K. Hagiwara, S. Ishihara, J. Kamoshita and B. A. Kniehl, Eur. Phys. J. C 14, 457 (2000) doi:10.1007/s100520000366 [hep-ph/0002043].
  * [3] S. S. Biswal, D. Choudhury, R. M. Godbole and R. K. Singh, Phys. Rev. D 73, 035001 (2006) Erratum: [Phys. Rev. D 74, 039904 (2006)] doi:10.1103/PhysRevD.74.039904, 10.1103/PhysRevD.73.035001 [hep-ph/0509070].
  * [4] R. M. Godbole, D. J. Miller and M. M. Muhlleitner, JHEP 0712, 031 (2007) doi:10.1088/1126-6708/2007/12/031 [arXiv:0708.0458 [hep-ph]].
  * [5] S. S. Biswal, D. Choudhury, R. M. Godbole and Mamta, Phys. Rev. D 79, 035012 (2009) doi:10.1103/PhysRevD.79.035012 [arXiv:0809.0202 [hep-ph]].
  * [6] S. S. Biswal and R. M. Godbole, Phys. Lett. B 680, 81 (2009) doi:10.1016/j.physletb.2009.08.014 [arXiv:0906.5471 [hep-ph]].
  * [7] S. D. Rindani and P. Sharma, Phys. Rev. D 79, 075007 (2009) doi:10.1103/PhysRevD.79.075007 [arXiv:0901.2821 [hep-ph]].
  * [8] S. D. Rindani and P. Sharma, Phys. Lett. B 693, 134 (2010) doi:10.1016/j.physletb.2010.08.027 [arXiv:1001.4931 [hep-ph]].
  * [9] I. Anderson et al., Phys. Rev. D 89, no. 3, 035007 (2014) doi:10.1103/PhysRevD.89.035007 [arXiv:1309.4819 [hep-ph]].
  * [10] R. M. Godbole, D. J. Miller, K. A. Mohan and C. D. White, JHEP 1504, 103 (2015) doi:10.1007/JHEP04(2015)103 [arXiv:1409.5449 [hep-ph]].
  * [11] R. Godbole, D. J. Miller, K. Mohan and C. D. White, Phys. Lett. B 730, 275 (2014) doi:10.1016/j.physletb.2014.01.069 [arXiv:1306.2573 [hep-ph]].
  * [12] N. Craig, J. Gu, Z. Liu and K. Wang, JHEP 1603, 050 (2016) doi:10.1007/JHEP03(2016)050 [arXiv:1512.06877 [hep-ph]].
  * [13] M. Beneke, D. Boito and Y. M. Wang, JHEP 1411, 028 (2014) doi:10.1007/JHEP11(2014)028 [arXiv:1406.1361 [hep-ph]].
  * [14] H. Khanpour and M. M. Najafabadi, Phys. Rev. D 95 (2017) no.5, 055026 doi:10.1103/PhysRevD.95.055026 [arXiv:1702.00951 [hep-ph]].
  * [15] T. V. Zagoskin and A. Y. Korchin, Phys. Rev. D 98, no. 9, 093008 (2018) doi:10.1103/PhysRevD.98.093008 [arXiv:1804.10011 [hep-ph]].
  * [16] H. D. Li, C. D. Lü and L. Y. Shan, arXiv:1901.10218 [hep-ex].
  * [17] H. R. He, X. Wan and Y. K. Wang, arXiv:1902.04756 [hep-ph].
  * [18] J. Nakamura, JHEP 1708, 008 (2017) doi:10.1007/JHEP08(2017)008 [arXiv:1706.01816 [hep-ph]].
  * [19] S. Banerjee, R. S. Gupta, J. Y. Reiness and M. Spannowsky, doi:10.1103/PhysRevD.100.115004 [arXiv:1905.02728 [hep-ph]].
  * [20] W. H. Chiu, S. C. Leung, T. Liu, K. F. Lyu and L. T. Wang, JHEP 1805, 081 (2018) doi:10.1007/JHEP05(2018)081 [arXiv:1711.04046 [hep-ph]].
  * [21] G. Durieux, C. Grojean, J. Gu and K. Wang, JHEP 1709, 014 (2017) doi:10.1007/JHEP09(2017)014 [arXiv:1704.02333 [hep-ph]].
  * [22] J. A. Aguilar-Saavedra, J. Bernabéu, V. A. Mitsou and A. Segarra, Eur. Phys. J. C 77, no. 4, 234 (2017) doi:10.1140/epjc/s10052-017-4795-8 [arXiv:1701.03115 [hep-ph]].
  * [23] F. Boudjema and R. K. Singh, JHEP 0907, 028 (2009) doi:10.1088/1126-6708/2009/07/028 [arXiv:0903.4705 [hep-ph]].
  * [24] R. Rahaman and R. K. Singh, Eur. Phys. J. C 77, no. 8, 521 (2017) doi:10.1140/epjc/s10052-017-5093-1 [arXiv:1703.06437 [hep-ph]].
  * [25] R. Rahaman and R. K. Singh, Eur. Phys. J. C 76, no. 10, 539 (2016) doi:10.1140/epjc/s10052-016-4374-4 [arXiv:1604.06677 [hep-ph]].
  * [26] K. Rao, S. D. Rindani and P. Sarmah, Nucl. Phys. B 950, 114840 (2020) doi:10.1016/j.nuclphysb.2019.114840 [arXiv:1904.06663 [hep-ph]].
  * [27] K. Rao and S. D. Rindani, Nucl. Phys. B 940, 78 (2019) doi:10.1016/j.nuclphysb.2019.01.008 [arXiv:1805.06602 [hep-ph]].
  * [28] A. M. Sirunyan et al. [CMS Collaboration], arXiv:1901.00174 [hep-ex].
  * [29] A. M. Sirunyan et al. [CMS Collaboration], [arXiv:1903.06973 [hep-ex]].
  * [30] B. Sahin, doi:10.1142/S0217732319502997 Mod. Phys. Lett. A34 no.36(2019).
  * [31] I. T. Cakir, O. Cakir, A. Senol and A. T. Tasci, Mod. Phys. Lett. A 28 (2013) no.31, 1350142 doi:10.1142/S0217732313501423 [arXiv:1304.3616 [hep-ph]].
  * [32] M. Aaboud et al. [ATLAS], Phys. Lett. B 786 (2018), 59-86 doi:10.1016/j.physletb.2018.09.013 [arXiv:1808.08238 [hep-ex]].
  * [33] A. M. Sirunyan et al. [CMS], Phys. Rev. Lett. 121 (2018) no.12, 121801 doi:10.1103/PhysRevLett.121.121801 [arXiv:1808.08242 [hep-ex]].
  * [34] D. Goncalves and J. Nakamura, Phys. Rev. D 98 (2018) no.9, 093005 doi:10.1103/PhysRevD.98.093005 [arXiv:1805.06385 [hep-ph]].
  * [35] J. A. M. Vermaseren, math-ph/0010025.
  * [36] E. Leader, Camb. Monogr. Part. Phys. Nucl. Phys. Cosmol. 15 (2011).
  * [37] L. A. Harland-Lang, A. D. Martin, P. Motylinski and R. S. Thorne, Eur. Phys. J. C 75 (2015) no.5, 204 doi:10.1140/epjc/s10052-015-3397-6 [arXiv:1412.3989 [hep-ph]].