# Self-Supervised Learning of Phenotypic Representations from Cell Images with
Weak Labels

Jan Oscar Cross-Zamirski
DAMTP & AstraZeneca
University of Cambridge
Cambridge CB3 0WA
<EMAIL_ADDRESS>&Guy Williams
Discovery Sciences, R&D
AstraZeneca
Cheshire, SK10 4TG
guy.williams2∗ &Elizabeth Mouchet
Discovery Sciences, R&D
AstraZeneca
Cheshire, SK10 4TG
elizabeth.mouchet∗ &Carola-Bibiane Schönlieb
DAMTP
University of Cambridge
Cambridge CB3 0WA
<EMAIL_ADDRESS>&Riku Turkki †
Discovery Sciences, R&D
AstraZeneca
Gothenburg, 431 50 Mölndal
riku.turkki∗ &Yinhai Wang †
Discovery Sciences, R&D
AstraZeneca
Cambridge, CB4 0WG
yinhai.wang∗<EMAIL_ADDRESS>&†Equal contribution

###### Abstract

We propose WS-DINO as a novel framework to use weak label information in
learning phenotypic representations from high-content fluorescent images of
cells. Our model is based on a knowledge distillation approach with a vision
transformer backbone (DINO), and we use this as a benchmark model for our
study. Using WS-DINO, we fine-tuned with weak label information available in
high-content microscopy screens (treatment and compound), and achieve state-
of-the-art performance in not-same-compound mechanism of action prediction on
the BBBC021 dataset (98%), and not-same-compound-and-batch performance (96%)
using the compound as the weak label. Our method bypasses single cell cropping
as a pre-processing step, and using self-attention maps we show that the model
learns structurally meaningful phenotypic profiles.

##  1 Introduction

Deep learning methods for phenotyping cells from microscopy images have
advanced rapidly alongside cellular imaging technology [1]. Profiles captured
from images can be used to quantify activity and perform downstream tasks such
as mechanism of action (MOA) prediction of cells treated with compounds of
interest. These methods hold promise for accelerating drug discovery pipelines
[2].

Typically, convolutional neural networks (CNNs) have been the method of
choice, however these models can struggle to generalise to new treatments, and
can be limited by label information [3]. Rarely are they interpretable, and
this may contribute to how machine learning methods have not reached their
potential in clinical use [4]. Across medical imaging fields, vision
transformers (ViTs) [5] and self-supervised [6] methods have risen as a viable
alternative to CNNs with labels [7] [8].

In this study we make the following contributions:

  * •

We introduce Weakly Supervised DINO (WS-DINO), a framework to incorporate weak
label information into a self-supervised knowledge distillation construction
based on the DINO algorithm [9].

  * •

We implement WS-DINO, and with the BBBC021 dataset [10] [11] we learn
representations with two weak labels - treatment and compound. We achieve
state-of-the-art results for not-same-compound (98%), and not-same-compound-
and-batch (96%) mechanism of action prediction using the compound as the weak
label. We show the learnt representations are biologically meaningful and
based on cellular structure with self-attention maps.

##  2 Background

###  2.1 Image-based profiling

Imaged-based profiling [2] is the process of extracting features from high-
content images (HCI) [12] of cells. The features extracted from these kinds of
images are used to build profiles which can assess bioactivity and MOA of the
cells in response to treatments with chemical compounds and/or genetic
perturbations. In drug discovery these phenotypic profiles are used to compare
new treatments with known ones, and downstream applications include routine
lead compound identification, drug target screening through CRISPR technology
and toxicity assessment [13].

Most phenotypic representation learning methods are based on single cell
segmentation [14] [15]. When applied to large datasets this can lead to high
computational demands in pre-processing images with segmentation algorithms
[16]. Working with single-cells requires aggregation to the population level
and may necessitate further computational requirements for corrections due to
cellular heterogeneity [17]. In high-content imaging there is metadata which
can be considered a weak label, examples being treatment compound and
concentration. In previous work a classification network was trained predict
the weak label as an auxiliary task, with a frozen layer of the network
forming the phenotypic representation used in downstream tasks such as MOA
prediction [18].

Unsupervised methods such as self-supervised learning constructions [6] using
a contrastive loss [19] have been very successful. Self-supervised algorithms
such as SwAV [7] and DINO [9] give consistently richer embeddings for
downstream tasks compared to a pretrained supervised baseline model in medical
imaging [20]. Self-supervised models have performed very well learning feature
represenations from single cell microscopy images [3] [14], and recently on
the image-level with a multi-scale approach [21].

###  2.2 DINO - self-supervised learning with knowledge distillation

DINO from Caron et al. [9] is a self-supervised learning approach which
incorporates aspects of knowledge distillation [22]. Uniquely combining self-
supervised learning with a Vision Transformer (ViT) [5] backbone, DINO can
learn features shown to perform well at k-NN clustering tasks.

DINO consists of student network $g_{\theta_{s}}$, parameterized by
$\theta_{s}$, trained with set of image crops to match the teacher network
$g_{\theta_{t}}$, parameterized by $\theta_{t}$, which sees a different set of
crops. For each image $x$ a set of views $V$ is generated which contains the
two global crops $x^{g,1}$ and $x^{g,2}$ as well as eight local crops. The
student is passed the set of global and local crops $V$, while the teacher
sees only the global crops. The student network is trained to maximise the
agreement between the outputs of $g_{\theta_{s}}$ and $g_{\theta_{t}}$. To
achieve this, the parameters of $\theta_{s}$ are found by minimizing the
cross-entropy loss:

$\min\limits_{\theta_{s}}{{\sum}}\limits_{x\in\\{x^{g,1},x^{g,2}\\}}\;{{\sum}}\limits_{\begin{subarray}{c}x^{\prime}\in
V\\\ x^{\prime}\neq x\end{subarray}}H(P_{t}(x),P_{s}(x^{\prime}))$ (1)

Where $H(a,b)=-a\log b$ and $P_{s}$ and $P_{t}$ are the probability
distributions of the student and teacher respectively (defined in the
Appendix). Teacher weights $\theta_{t}$ are frozen during each epoch of
student training and updated iteratively with an exponential moving average
based on the previous weights of the student network with the formula:
$\theta_{t}\leftarrow\lambda\theta_{t}+(1-\lambda)\theta_{s}$, where $\lambda$
is the momentum parameter.

Both the student and the teacher network have a ViT [5] backbone. The inputs
are small patches of fixed size ($N\times N$) in a non-overlapping grid, which
are embedded (alongside a positional embedding) with a linear layer. An
additional learnable [CLS] token with a projection head output $h$ is added to
the embeddings and sent to the a Transformer encoder. The attention mechanism
allows the ViT to synthesise information across the whole image as self-
attention layers in the Transformer globally update the attention token
embeddings.

##  3 Methods

###  3.1 WS-DINO

We propose a Weakly Supervised form of self-distillation with no labels (WS-
DINO) as an adaptation to the DINO algorithm. The aim of the approach is to
retain the benefits of self-supervised learning with knowledge distillation
while also incorporating the weak label information as part of the training.

We introduce the notation $x_{i,y_{i}}$ to represent the $i^{\text{th}}$ field
of view of a fluorescent channel in the dataset which has been treated with
treatment or compound $y_{i}$ \- the weak label. The superscript contains the
crop information: $g$ for global and $l$ for local crops. When generating the
sets of different views $V_{t}$ (seen by teacher) and $V_{s}$ (seen by
student) for training, we enforce the constraint that the global and local
crops are sampled from different images with the same weak label. We define
the sets $V_{t}$ and $V_{s}$ for the randomly sampled ordered pair $(i,j)$
where $i\neq j$ and $y_{i}=y_{j}$:

$V_{t}=\\{x_{i,y_{i}}^{g,1},x_{i,y_{i}}^{g,2}\\}$ (2)
$V_{s}=V_{t}\cup\begin{rcases}\begin{dcases}x_{j,y_{j}}^{l,k}:&k\in\\{1,...,n\\}\\\
\end{dcases}\end{rcases}$ (3)

Where the superscript details the $k^{\text{th}}$ local crop of $n$ total
crops (default: $n=8$). Sampling the local crops from a different image is in
contrast to sampling random crops with different augmentations from the same
image, which is the method of DINO. For WS-DINO we minimize the loss:

$\min\limits_{\theta_{s}}{{\sum}}\limits_{x\in
V_{t}}\;{{\sum}}\limits_{\begin{subarray}{c}x^{\prime}\in V_{s}\\\
x^{\prime}\neq x\end{subarray}}H(P_{t}(x),P_{s}(x^{\prime}))$ (4)

In this section we have presented the key details of our adaptation to the
DINO algorithm from the original paper [9]. We provide implementation details
in section 3.4 and in our GitHub repository.
111https://github.com/crosszamirski/WS-DINO

###  3.2 Dataset and augmentation

We evaluated our experiments with the publicly available BBBC021 dataset [10],
which consists of MCF-7 breast cancer cells exposed to several chemical
compounds for 24h. The cells are stained with three fluorescent labels: DNA,
F-actin, and $\upbeta$-tubulin. For training and evaluation, we used the
subset which has been labelled for MOA evaluation [11] which consists of 12
unique MOA, 38 unique compounds and 103 unique treatments
(compound/concentration pairs) across 10 experimental batches. There are 1320
3-channel images in the dimethyl sulfoxide (DMSO) treated control group, and
2528 3-channel images with MOA annotations. Several studies have used the
annotated dataset to evaluate classical and machine learning strategies [11]
[14] [15] [18] [21] [23] [24] [25] [26] [27].

Pre-processing: To remove systematic variations in pixel intensities across
each field of view, we used a CellProfiler [16] pipeline to perform an
illumination correction. A smoothing function with a filter size of $320$
pixels was used to create an illumination correction function. Each image was
corrected by dividing all pixel intensities by the illumination correction
function. This methodology is consistent with best practice established by the
JUMP Cell Painting Consortium [28]. After correction, images were resized
using bicubic interpolation to a size of $640\times 512$ pixels. A maximum
pixel value cut-off of $10,000$ was enforced before normalizing each image to
have a mean of $0$ and standard deviation of $1$. The scripts to replicate the
pre-processing are available in our GitHub repository.

Post-processing: Next we take the median feature embedding from four
$224\times 224$ crops around the centre of each image (equivalent to a
$448\times 448$ pixel centre crop split into four non-overlapping crops). To
aid in correcting for batch effects (which can be significant in HCI screens)
and to capture the range of variation of unperturbed cells, we applied typical
variation normalization (TVN) [15] as a post-processing correction. TVN is a
technique where a principal component analysis (PCA) transformation without
dimensionality reduction (whitening) is learned using the DMSO embeddings
(aggregated to one embedding per image from the median of four crops around
the centre). The mean of each of the field-level embeddings is taken across a
plate, followed by median aggregation to treatment level resulting in 103
embeddings - in line with the studies following Ando et al. [15].

###  3.3 Performance validation and evaluation

We evaluated the corrected and aggregated feature embeddings with two metrics
ubiquitous in representation learning studies using the BBBC021 dataset: not-
same-compound (NSC) matching [11] and not-same-compound-and-batch (NSCB)
matching [15] with the MOA labels. NSC matching is a 1-Nearest-Neighbour
(1-NN) match for each given treatment to the nearest neighbour in
representation space which is not of the same compound. Cosine distance was
used as the distance measure. NSCB matching is the same method as NSC matching
with the additional constraint to exclude treatments from the same compound
and batch as the given treatment. A small number of treatments have a MOA
label only present in a single batch and these were excluded in NSCB
evaluation.

###  3.4 Training

For all the models the following variables are kept constant: network
architecture (ViT-S/8 backbone with 3-layer multi-layer perceptron head), size
and number of crops: 2 global $224\times 224$ random resize crops (range =
$0.1-0.2$ of full size image) and 8 local $96\times 96$ random resize crops
(range = $0.04-0.08$), teacher momentum = 0.99, and batch size = 16. There was
no batch normalization, weight decay or gradient clipping. The optimizer was
adamW [29]. We applied horizontal and vertical flips (with a probability of
0.5 per sample) as augmentations in training. Learning rate and temperature
parameters are controlled with adaptive cosine schedules with linear warm-up.
After 10 epochs of linear warmup, the learning rate = $4\times 10^{-6}$ and
the teacher temperature = 0.04. The learning rate was decayed to $3\times
10^{-6}$ over 400 epochs, and the teacher temperature kept constant after
warmup. All other parameters unless otherwise stated are default DINO
parameters. A separate model was trained for each of of the three channels.
Every model was initialized with weights from DINO trained on RGB ImageNet [9]
[30] (grayscale images were converted to RGB). Embeddings of size 384 were
extracted for each channel then concatenated to a size of 1152 followed by L2
normalization. We provide a full PyTorch implementation in our GitHub
repository.

##  4 Results and Discussion

As a baseline we trained the unsupervised DINO network [9]. First, we
evaluated with ImageNet weights only (transfer learning). Next we fine-tuned
the model on the BBBC021 annotated dataset without weak labels. We then
trained WS-DINO with the same network and parameters. Two weak labels $y$ were
evaluated: treatment and compound. Additionally we trained with MOA as the
label $y$, although this cannot be considered a weak label. The images in
training were sampled using a weighted random sampler enforcing even sampling
of images from each weak label group.

Table 1: A summary of the results of our experiments. Further details of mean NSC and NSCB scores in training with a discussion of best epoch selection is provided in the Appendix. Model | Weak Label | NSC | NSCB
---|---|---|---
DINO with ImageNet weights only | None | 91% | 82%
DINO finetuned on BBBC021 | None | 92% | 95%
WS-DINO finetuned on BBBC021 | Treatment | 92% | 90%
| Compound | 98% | 96%
| MOA | 100% | 100%
Figure 1: Examples of BBBC021 images with their coupled self-attention maps.
Left two columns: DNA channel. Middle two columns: $\upbeta$-tubulin channel.
Right two columns: F-actin channel. The first two rows show the images at full
size, and the bottom row displays a section of an image zoomed in by a scale
factor of four. Produced from WS-DINO weights with compound as the weak label.

WS-DINO achieves state-of-the-art results on BBBC021 using the compound as the
weak label, outperforming all known previous approaches using this dataset in
NSC and NSCB MOA prediction (Table 1 and Appendix). It is notable that using
treatment as the weak label does not improve upon the unsupervised approach
(DINO). In this dataset all images of the same treatment are from the same
batch, however this is not true for the different compounds and MOAs which are
expressed with different images across multiple batches. Hence we propose that
WS-DINO may provide effective batch correction for datasets with weak label
classes with images representing multiple batches.

We present self-attention maps (Figure 1) which reveal the segmentation
properties of the algorithm. Such visualisations increase confidence in the
model by demonstrating the network is learning biologically and structurally
meaningful features.

We use MOA as as psuedo-weak label as a proof-of-concept for our method.
However, we propose future work to evaluate MOA prediction using datasets with
partial MOA labels. One advantage of WS-DINO is that it is adaptable to
datasets with some known MOA information, a property of many real datasets
used in drug discovery.

WS-DINO synthesises a powerful self-supervised framework with a way to
implicitly incorporate the informative weak labels in learning phenotypic
representations. Our method can contribute to accelerating drug discovery
pipelines by clustering phenotypes in a biologically meaningful hierarchy. The
proposed framework is general and the method is not specific to DINO. The
sampling of pairs of images with the same weak label is a concept adaptable to
other networks. Self-supervised methods are an active field of research and
new algorithms have outperformed DINO on ImageNet classification, for example
Masked Siamese Networks [31]. Future work could incorporate these studies
which would be adaptable in a very similar way to WS-DINO. Finally, we will
expand this research to other datasets with more classes and labels, as well
as different downstream tasks.

## Acknowledgments and Disclosure of Funding

Jan Oscar Cross-Zamirski is funded by the UKRI-BBSRC DTP (UK Research and
Innovation and the Biotechnology and Biological Science Research Council
Doctoral Training Partnership) studentship grant 21-3688 and AstraZeneca.

Carola-Bibiane Schönlieb acknowledges support from the Philip Leverhulme
Prize, the Royal Society Wolfson Fellowship, the EPSRC grants EP/S026045/1 and
EP/T003553/1, EP/N014588/1, EP/T017961/1, the Wellcome Trust 215733/Z/19/Z and
221633/Z/20/Z, the Leverhulme Trust project Unveiling the invisible, the
European Union Horizon 2020 research and innovation programme under the Marie
Skodowska-Curie grant agreement No. 777826 NoMADS, the Cantab Capital
Institute for the Mathematics of Information and the Alan Turing Institute.

Authors Guy Williams, Elizabeth Mouchet, Riku Turkki and Yinhai Wang are
employees of AstraZeneca. AstraZeneca provided the funding for this research
and provided support in the form of salaries for the authors, but did not have
any additional role in the study design, data collection and analysis,
decision to publish, or preparation of the manuscript. We would also like to
acknowledge the use of the Scientific Computing Platform (SCP) within
AstraZeneca.

## References

  * [1] Pratapa, A., Doron, M., & Caicedo, J. C. (2021). Image-based cell phenotyping with deep learning. Current opinion in chemical biology, 65, 9–17. https://doi.org/10.1016/j.cbpa.2021.04.001
  * [2] Chandrasekaran, S.N., Ceulemans, H., Boyd, J.D. & Carpenter, A. E. (2021). Image-based profiling for drug discovery: due for a machine-learning upgrade? Nat Rev Drug Discov 20, 145–159. https://doi.org/10.1038/s41573-020-00117-w
  * [3] Lu, A.X., Kraus, O.Z., Cooper, S. & Moses, A.M. (2019). Learning unsupervised feature representations for single cell microscopy images with paired cell inpainting. PLoS Comput Biol 15(9): e1007348. https://doi.org/10.1371/journal.pcbi.1007348
  * [4] Weissler, E., Naumann, T., Andersson, T., Ranganath, R.,Elemento, O., Luo, Y. et al. (2021). The role of machine learning in clinical research: transforming the future of evidence generation. Trials. 22. http://dx.doi.org/10.1186/s13063-021-05489-x
  * [5] Dosovitskiy, A., Beyer, L., Kolesnikov, A., Weissenborn, D., Zhai, X., Unterthiner, T., Dehghani, M., Minderer, M., Heigold, G., Gelly, S. et al. (2020). An image is worth 16x16 words: Transformers for image recognition at scale. arXiv preprint: https://doi.org/10.48550/arXiv.2010.11929
  * [6] Doersch, C., & Zisserman, A. (2017). Multi-task Self-Supervised Visual Learning arXiv preprint: https://doi.org/10.48550/arxiv.1708.07860
  * [7] Caron, M., Misra, I., Mairal, J., Goyal, P., Bojanowski, P. & Joulin, A. (2020) Unsupervised learning of visual features by contrasting cluster assignments. Proceedings of Advances in Neural Information Processing Systems (NeurIPS) https://doi.org/10.48550/arXiv.2006.09882
  * [8] Matsoukas, C., Haslum, J.F., Söderberg, M., & Smith, K. (2021). Is it Time to Replace CNNs with Transformers for Medical Images? arXiv preprint: https://doi.org/10.48550/arxiv.2108.09038
  * [9] Caron, M., Touvron, H., Misra, I., Jégou, H., Mairal, J., Bojanowski, P. & Joulin, A. (2021). Emerging Properties in Self-Supervised Vision Transformers.  Proceedings of the International Conference on Computer Vision (ICCV). https://doi.org/10.48550/arXiv.2104.14294
  * [10] Caie, P.D., Walls, R.E., Ingleston-Orme, A., Daya, S., Houslay, T., Eagle, R., Roberts, M.E. & Carragher, N.O. (2010). High-content phenotypic profiling of drug response signatures across distinct cancer cells. Molecular cancer therapeutics 9(6), 1913–1926. https://doi.org/10.1158/1535-7163.mct-09-1148
  * [11] Ljosa, V., Caie, P.D., Ter Horst, R., Sokolnicki, K.L., Jenkins, E.L., Daya, S., Roberts, M.E., Jones, T.R., Singh, S., Genovesio, A. et al. (2013). Comparison of methods for image-based profiling of cellular morphological responses to small-molecule treatment. Journal of biomolecular screening 18(10), 1321–1329. https://doi.org/10.1177/1087057113503553
  * [12] Bickle, M. (2010). The beautiful cell: high-content screening in drug discovery. Anal. Bioanal. Chem. 398, 219–226. https://doi.org/10.1007/s00216-010-3788-3
  * [13] Trapotsi, M.-A., Mouchet, E., Williams, G., Monteverde, T., Juhani, K., Turkki, R. et al. (2022). Cell Morphological Profiling Enables High-Throughput Screening for PROteolysis TArgeting Chimera (PROTAC) Phenotypic Signature. ACS Chemical Biology 2022 17 (7), 1733-1744. https://doi.org/10.1021/acschembio.2c00076
  * [14] Perakis, A., Gorji, A., Jain S., Chaitanya, K., Rizza, S., & Konukoglu, E. (2021). Contrastive Learning of Single-Cell Phenotypic Representations for Treatment Classification. arXiv preprint: https://doi.org/10.1007%2F978-3-030-87589-3_58
  * [15] Ando, D.M., McLean, C.Y. & Berndl, M. (2017). Improving phenotypic measurements in high-content imaging screens. BioRxiv preprint p. 161422. https://doi.org/10.1101/161422
  * [16] McQuin, C., Goodman, A., Chernyshev, V., Kamentsky, L., Cimini, B.A. et al. (2018). CellProfiler 3.0: Next-generation image processing for biology. PLOS Biology 16(7): e2005970. https://doi.org/10.1371/journal.pbio.2005970
  * [17] Rohban, M. H., Abbasi, H. S., Singh, S. & Carpenter, A. E. (2019). Capturing single-cell heterogeneity via data fusion improves image-based profiling. Nat. Commun. 10, 2082. https://doi.org/10.1038/s41467-019-10154-8
  * [18] Caicedo, J. C., McQuin, C., Goodman, A., Singh, S. & Carpenter, A. E. (2018). Weakly Supervised Learning of Single-Cell Feature Embeddings. Proceedings. IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2018, 9309–9318. https://doi.org/10.1109/cvpr.2018.00970
  * [19] Hadsell, R., Chopra S., & LeCun Y. (2006). Dimensionality Reduction by Learning an Invariant Mapping, 2006 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’06), pp. 1735-1742. https://doi.org/10.1109/CVPR.2006.100.
  * [20] Truong, T., Mohammadi, S. & Lenga, M. (2021). How Transferable are Self-supervised Features in Medical Image Classification Tasks? Proceedings of Machine Learning for Health, in Proceedings of Machine Learning Research 158:54-74. https://doi.org/10.48550/arXiv.2108.10048
  * [21] Janssens, R., Zhang, X., Kauffmann, A., de Weck, A., & Durand, E. Y. (2021). Fully unsupervised deep mode of action learning for phenotyping high-content cellular images. Bioinformatics (Oxford, England), btab497. Advance online publication. https://doi.org/10.1093/bioinformatics/btab497
  * [22] Hinton, G., Vinyals, O. & Dean, J. (2015). Distilling the Knowledge in a Neural Network. arXiv preprint: https://doi.org/10.48550/arxiv.1503.02531
  * [23] Singh, S., Bray, M. A., Jones, T. R., & Carpenter, A. E. (2014). Pipeline for illumination correction of images for high-throughput microscopy. Journal of microscopy, 256(3), 231–236. https://doi.org/10.1111/jmi.12178
  * [24] Godinez, W. J., Hossain, I., Lazic, S. E., Davies, J. W., & Zhang, X. (2017). A multi-scale convolutional neural network for phenotyping high-content cellular images. Bioinformatics (Oxford, England), 33(13), 2010–2019. https://doi.org/10.1093/bioinformatics/btx069
  * [25] Pawlowski N., Caicedo J.C., Singh S., Carpenter A.E. & Storkey A. (2016). Automating morphological profiling with generic deep convolutional networks. BioRxiv preprint: https://doi.org/10.1101/085118
  * [26] Lafarge, M.W., Caicedo, J.C., Carpenter, A.E., Pluim, J.P., Singh, S. & Veta, M. (2019). Capturing Single-Cell Phenotypic Variation via Unsupervised Representation Learning. Proceedings of The 2nd International Conference on Medical Imaging with Deep Learning, PMLR. https://pubmed.ncbi.nlm.nih.gov/35874600
  * [27] Goldsborough, P., Pawlowski, N., Caicedo, J.C., Singh, S. & Carpenter, A.E. (2017). CytoGAN: Generative Modeling of Cell Images. BioRxiv preprint: https://doi.org/10.1101/227645
  * [28] Cimini, B.A., Chandrasekaran, N., Kost-Alimova, M., Miller, L., Goodale, A., Fritchman, B. et al. (2022). Optimizing the Cell Painting assay for image-based profiling. BioRxiv preprint: https://doi.org/10.1101/2022.07.13.499171
  * [29] Loshchilov, I., & Hutter, F. (2019). Decoupled Weight Decay Regularization. ICLR 2019. https://doi.org/10.48550/arXiv.1711.05101
  * [30] Deng, J., Dong, W., Socher, R., Li, L.-J., Kai, L. & Li, F.-F. (2009). ImageNet: A large-scale hierarchical image database. 2009 IEEE Conference on Computer Vision and Pattern Recognition. https://doi.org/10.1109/CVPR.2009.5206848
  * [31] Assran, M., Caron, M., Misra, I., Bojanowski, P., Bordes, F. Vincent, P., Joulin, A., Rabbat, M. & Ballas, N. (2022). Masked Siamese Networks for Label-Efficient Learning. arXiv preprint: https://doi.org/10.48550/arXiv.2204.07141

##  Appendix A Probability distributions of student and teacher

The probability distributions for the student $P_{s}$ and teacher $P_{t}$ are
defined by:

$P_{s}(x)^{(i)}=\cfrac{\exp{(g_{s}(x)^{(i)}}/\tau_{s})}{\sum_{k=1}^{K}\exp{(g_{s}(x)^{(k)}}/\tau_{s})}$

for $K$ dimensions and temperature parameters $\tau_{s}>0$ and $\tau_{t}>0$
[9].

##  Appendix B Further training details

The hyperparameters are fixed across all of the models we train. Many of the
parameters are the default values from the original DINO implementation. We
optimised the remaining parameters such as learning rate with a small subset
of BBBC021 trained without weak labels. In our final models the batch size is
limited by GPU memory - a batch size of 16 is small and most likely a larger
batch size would improve performance. The models were trained with two GPUs on
the AstraZeneca Scientific Computing Platform with a maximum allocation of 32G
memory for each GPU.

Table 2: The results of training with NSC and NSCB scores of the best performing epoch as well as the mean NSC and NSCB scores between 50 and 250 epochs of training each model. Model | Weak Label | Best NSC/NSCB | Mean NSC/NSCB
---|---|---|---
DINO finetuned on BBBC021 | None | 92% / 95% | 90% / 90%
WS-DINO finetuned on BBBC021 | Treatment | 92% / 90% | 89% / 83%
| Compound |  98% / 96% | 96% / 93%
| MOA | 100% / 100% | 99% / 99%

We determined the best epoch by following the method of Janssens et al. [21].
We calculated NSC for each model and selected the epoch with the highest
score. This method was chosen to allow a comparison of our model to other
studies, however we note there are drawbacks to selecting best epoch with MOA
information. In Table 2 the mean values of NSC and NSCB over 200 epochs of
training are displayed for each model, and the results show improved
performance for NSC and NSCB MOA matching was observed over multiple epochs of
training using both compound and MOA as the weak label. The values are lower
than the best values as this range contains both under- and over-fit models.
We provide epoch data and training logs in our GitHub repository.

##  Appendix C Further results

Table 3: A comparison of our method to other selected studies using the BBBC021 annotated dataset. Type | Method | Reference | Single Cell | NSC | NSCB
---|---|---|---|---|---
Weakly | WS-DINO | This work | No | 98% | 96%
Supervised | CNN with Mixup |  Caicedo et al. 2018 [18] | Yes | 95% | 89%
Classical | CellProfiler |  Singh et al. 2014 [23] | Yes | 90% | 85%
Features | Factor Analysis |  Ljosa et al. 2013 [11] | Yes | 94% | 77%
Supervised | Multiscale-CNN |  Godinez et al. 2017 [24] | No | 93% | N/A
Unsupervised | Contrastive Learning |  Perakis et al. 2021 [14] | Yes | 96% | 95%
| UMM Discovery |  Janssens et al. 2020 [21] | No | 97% | 85%
| VAE+ |  Lafarge et al. 2019 [26] | Yes | 93% | 82%
| CytoGAN: LSGAN |  Goldsborough et al. 2017 [27] | No | 68% | N/A
Transfer | Deep Metric Network |  Ando et al. 2017 [15] | Yes | 96% | 95%
Learning | Inception V3 Pretrained |  Pawlowski et al. 2016 [25] | No | 91% | N/A
Figure 2: Two-dimensional t-SNE plot of each aggregated treatment feature from
200 epochs of training WS-DINO with compound as the weak label: 98% NSC and
96% NSCB MOA classification.


# Experimental Measurement of Geometric Phase of Non-Geodesic Circles

Andrew A. Voitiv1, Mark T. Lusk2 Mark E. Siemens1<EMAIL_ADDRESS>1Department
of Physics and Astronomy, University of Denver, 2112 E. Wesley Avenue, Denver,
CO 80208, USA  2Department of Physics, Colorado School of Mines, 1500 Illinois
Street, Golden, CO 80401, USA

###### Abstract

We present and implement a method for the experimental measurement of
geometric phase of non-geodesic (small) circles on any SU(2) parameter space.
This phase is measured by subtracting the dynamic phase contribution from the
total phase accumulated. Our design does not require theoretical anticipation
of this dynamic phase value and the methods are generally applicable to any
system accessible to interferometric and projection measurements. Experimental
implementations are presented for two settings: 1.) the Sphere of Modes of
orbital angular momentum, and 2.) the Poincaré Sphere of polarizations of
Gaussian beams.

Geometric phase plays an intrinsic role [1] in the evolution of all quantum
and classical fields, from condensed matter to optics and beyond [2]. It has
captured widespread interest both for applications such as electron transport
in graphene [3], braided photonic solid-state wave guides [4], and noise-
resilient manipulation of solid-state spin-qubit quantum phases [5], and also
for its relevance at the intersection of geometry, information science, and
quantum field theory [6, 7, 8].

Optics, both classical and quantum, is a historically important and accessible
arena for studying the geometric phase of state evolution. In fact, it was in
an optical polarization setting that Pancharatnam first discovered a geometric
phase [9], nearly thirty years prior to Berry’s important identification of
its quantum mechanical counterpart [10]. Geometric phase has since been given
a holonomic foundation within differential geometry [11], subjected to
theoretical examination in a variety of optical settings [12, 13, 14, 15, 16],
and experimentally realized in both classical (linear and nonlinear) and
quantum optics [17, 18, 19, 20, 21, 22, 23].

The investigation of geometric phase in association with the orbital angular
momentum (OAM) of light [24] emerged quickly after the explosion of work
related to spin angular momentum (polarization). Fig. 1 (a) shows a Sphere of
Modes (SoM) for transformations of these modes, identified as the OAM analog
of the polarization Poincaré Sphere (PS) [25] which led to a series of
demonstrations, over the past two decades, of geometric phase accumulation for
such optical vortices [26, 27, 28, 29, 30, 31, 32]. A recent contribution
explicitly pinpoints the accumulation of this geometric phase with the
evolution of beam waist, waist position, and fiber phase [33] which depend, in
turn, on the transit fraction through optical elements [34]. It is fair to say
that a theoretical understanding of the nature of geometric phase in classical
paraxial optics is now well-established.

Figure 1: (a) Sphere of Modes (SoM) for first-order Gaussian modes (optical
vortices); similarly, the Poincaré Sphere (PS) exists for polarization with
circular polarizations at the poles and linear polarizations along the equator
[35]. Axes are labelled pictorially with experimental measurements of mode
amplitude (grayscale) and phase (hue), along with spherical coordinates
$(\xi,\theta)$. North and South poles are labeled with $\mathcal{I_{+}}$ and
$\mathcal{I_{-}}$, respectively, for positively and negatively charged
circular vortices. The equatorial axes are labeled with $\mathcal{I}_{0}$,
$\mathcal{I}_{45}$, $\mathcal{I}_{90}$, and $\mathcal{I}_{135}$, corresponding
with the tilt-angle (in degrees) of the Hermite-Gaussian modes. The yellow
small circle is a circuit taken from a state $A$ at $(\xi,90\degree)$ to
intermediate state $B$ at $(-\xi,90\degree)$ and then back to point $A$. (b)
Depiction of how small circle geometric phase can exceed great circle
geometric phase, depending on angle $\beta$ and according to
$\Phi_{\mathrm{geo}}=-\Omega/2$, for solid angle $\Omega$ (colored red). Note
that all trajectories trace their loops in the same direction [34].

Fig. 1 (b) illustrates the difference between small- and great-circle circuits
and their enclosed solid angles. Geodesic arcs are contiguous sections of
great circles—i.e., they are associated with circles of the same diameter as
the sphere itself. Many experimental measurements of geometric phase
associated with optical vortices or polarization have been careful to restrict
attention to geodesic trajectories [27, 32]. The motivation for this is that
the associated geometric phase, for an enclosed solid angle, is then equal to
the total phase accumulation minus a standard propagation-dependent phase that
can be removed using straightforward interferometry [35, 24, 27, 32].
Measuring the geometric phase of small circles is a greater challenge because
it requires the removal of an additional contribution to dynamic phase [35,
34]. For paradigms associated with evolving polarization, this additional
dynamic phase has been removed either within particular experiments or with
the aid of theoretical predictions. For example, it was quantified using
intensity fringes for a two-pinhole interference experiment [23]. A more
general type of polarization experiment employed gauge transformations to
negate all dynamic phase [20]; in that work, the dynamic phase was
theoretically predicted and then supplementary optics were used to remove it
from the polarimetrically-measured total phase. Furthermore, there has been no
purely experimental method for removing the non-geodesic contribution to
dynamic phase that can be used with the OAM Sphere of Modes.

In this Letter, we show how the geometric phase of non-geodesic circular
trajectories on a Sphere of Modes (SoM) or Poincaré sphere (PS) can be
experimentally quantified. The key to our technique is the addition of a
measurement halfway through the system evolution. The method is then
implemented for both orbital and spin angular momentum transformations of
Gaussian laser beams. While demonstrated within these specific optical
settings, our simple theoretical framework is general and we expect the
approach to be applicable to any system amenable to interferometric and
projection measurements.

Consider the SoM shown in Fig. 1. Circuits on this sphere comprise a
progression of vortex modes that can be realized with a combination of
cylindrical lenses and Dove prisms [36, 24, 25, 27, 34, 32]. Analogous
circuits can be generated on the PS using waveplates to transform polarization
[35, 25]. As in [27] and [17], we design our OAM and polarization experiments
such that all phase measurements are made relative to a reference Gaussian
beam. The dynamic phase associated with free-space propagation can therefore
be removed and subsequent discussions of dynamic phase refer to phase
accumulations due to non-geodesic paths on the PS or SoM. The geometric phase
here does not rely on adiabatic dynamics [37].

The total, dynamic, and geometric phases of circular arcs (of arc length
$\delta$) with opening angle $\beta$ on a PS are well-established
theoretically [35, 34]:

$\Phi_{\mathrm{tot}}=\mathrm{Arg}\langle A|B\rangle,$ (1a)
$\Phi_{\mathrm{dyn}}=-\frac{\delta}{2}\cos{\beta},$ (1b)
$\begin{split}\Phi_{\mathrm{geo}}=\Phi_{\mathrm{tot}}-\Phi_{\mathrm{dyn}}.\end{split}$
(1c)

Here $|A\rangle$ and $|B\rangle$ are spinors in SU(2) that correspond to the
Stokes positions shown in Fig. 1. The geometric phase obtained is equivalent
to the familiar representation as negative one-half the solid angle enclosed
by the circuit, as depicted in Fig. 1 (b).

A relationship can also be derived between the azimuthal angle subtended,
$\delta$, and the total phase accumulation [34]:

$\Phi_{\mathrm{tot}}=-\arctan{\left(\cos{\frac{\delta}{2}}\,,\,\cos{\beta}\,\sin{\frac{\delta}{2}}\right)},$
(2)

where the sign-unambiguous arctangent is defined as
$\arctan(x,y)=\arctan(y/x)$ with the quadrant accounted for.

For a complete small circle path on the PS or SoM ($\delta=2\pi$), Eqn. 2
yields $\Phi_{\textrm{tot}}=-\pi$ and Eqn. 1b yields
$\Phi_{\textrm{dyn}}=-\pi\cos\beta$. Experimental measurements of the states
before and after a complete circular path (e.g. interferrometrically [24, 27,
32]) can be plugged into Eqn. 1a to measure the total phase, but the dynamic
phase is not accessible this way because path-taken information is lost using
two states that coincide. The geometric phase is therefore inaccessible as
well. This is the heart of the problem we address in this Letter and it
warrants repeating for emphasis: one cannot measure the geometric phase of
non-geodesic circles using only measurements of initial and final states
because those two data points alone shed no information on the size of the
circle—the $\cos{\beta}$ dependence. Therefore, we seek a method for measuring
$\cos{\beta}$ for an arbitrary circular path on the PS, which would yield the
geometric phase, Eqn. 1c, from direct experimental measurements.

Our strategy is to determine $\Phi_{\mathrm{dyn}}$ by taking a measurement of
the state after a half-circle trajectory. Such a trajectory can be easily
obtained in optical systems by transmission through a Dove prism, a half
waveplate (for the polarization PS) or a $\pi$-converter (for the vortex-mode
SoM). A half-circle trajectory, $\delta=\pi$, has two important advantages:
1.) all of the phase components in Eqns. 1 are simply half of their values for
the full circle, and 2.) Eqn. 2 yields
$\Phi_{\mathrm{tot}}=-\pi/2$—independent of $\beta$. We will show that this
permits measurements of the initial and half-circle states to directly measure
$\Phi_{\mathrm{dyn}}$. $\Phi_{\mathrm{geo}}$ can then be determined from Eqn.
1c.

Two SU(2) spinors on the same circle, $|A\rangle$ and $|B\rangle$, can be
written as

$\begin{split}|A\rangle&=\cos{\frac{\beta}{2}}|\leftrightarrow\,\rangle+\sin{\frac{\beta}{2}}|\updownarrow\,\rangle,\\\
|B\rangle=&\cos{\frac{\beta}{2}}\,e^{-i\delta/2}|\leftrightarrow\,\rangle+\sin{\frac{\beta}{2}}\,e^{i\delta/2}|\updownarrow\,\rangle,\end{split}$
(3)

where $|\leftrightarrow\,\rangle$ and $|\updownarrow\,\rangle$ are general
kets that correspond to the SoM or PS axes $\mathcal{I}_{0}$ and
$\mathcal{I}_{90}$, respectively. There are equivalent orthogonal basis
vectors for any SU(2) setting. Here, we use the Hermite Gaussian modes
$|\mathrm{HG}_{10}\rangle$ and $|\mathrm{HG}_{01}\rangle$ for the SoM and
linear polarization states $|\mathrm{H}\rangle$ and $|\mathrm{V}\rangle$ for
the PS. This implies that

$\langle
A|B\rangle=\cos^{2}{\frac{\beta}{2}}\,e^{-i\delta/2}+\sin^{2}{\frac{\beta}{2}}\,e^{i\delta/2}.$

Collect up the real and imaginary parts to give

$\begin{split}\langle
A|B\rangle=\left(\cos^{2}{\frac{\beta}{2}}\cos{\frac{\delta}{2}}+\sin^{2}{\frac{\beta}{2}}\cos{\frac{\delta}{2}}\right)\\\
+\,\,i\,\left(\sin^{2}{\frac{\beta}{2}}\sin{\frac{\delta}{2}}-\cos^{2}{\frac{\beta}{2}}\sin{\frac{\delta}{2}}\right).\end{split}$

This can be immediately simplified to

$\langle
A|B\rangle=\cos{\frac{\delta}{2}}-i\,\cos{\beta}\sin{\frac{\delta}{2}},$ (4)

and taking the imaginary part yields

$\cos{\beta}=-\csc{\frac{\delta}{2}}\,\mathrm{Im}\langle A|B\rangle.$ (5)

In the case of a half-circle path, $\delta=\pi$, and the result simplifies to

$\cos{\beta}=-\mathrm{Im}\langle A_{0}|B_{\pi}\rangle,$ (6)

for initial state $|A_{0}\rangle(\delta=0)$ and intermediate state
$|B_{\pi}\rangle(\delta=\pi)$ connected by a half-circle.

A reformulation of the total, dynamic, and geometric phase terms is now
possible in terms of experimentally-measurable quantities:

$\Phi_{\mathrm{tot}}=\mathrm{Arg}\langle A_{0}|C_{2\pi}\rangle,$ (7a)
$\Phi_{\mathrm{dyn}}=\mathrm{Arg}\langle
A_{0}|C_{2\pi}\rangle\,\mathrm{Im}\langle A_{0}|B_{\pi}\rangle,$ (7b)
$\begin{split}\Phi_{\mathrm{geo}}=\Phi_{\mathrm{tot}}-\Phi_{\mathrm{dyn}}=\mathrm{Arg}\langle
A_{0}|C_{2\pi}\rangle\,\left(1-\mathrm{Im}\langle
A_{0}|B_{\pi}\rangle\right),\end{split}$ (7c)

for $|C_{2\pi}\rangle(\delta=2\pi$) being the final state that completes the
circle.

We demonstrate the implementation of this algorithm for measuring geometric
phase for two different optical bases. Our first implementation uses vortex
modes and a $\pi$-converter to provide the needed half-circle trajectory, as
shown in the schematic of Fig. 2 (a). A collimated, single-mode Gaussian is
transmitted through a spatial light modulator (SLM) [38], which crafts an
initial vortex state, $|A\rangle$. Details on how these vortices are made
experimentally are provided in the Supplementary Information. $|A\rangle$ is
completely characterized by a camera before the first $\pi$-converter; that
is, we measure both amplitude (square-root of the recorded intensity) and
phase (measured via phase-shifting digital holography [39]). We then fit this
complex field to an explicit expression for the input beam in which a tunable
initial phase-shift, $\phi$, is included to account for shifts associated with
propagation through the apparatus that produces the beam.

Figure 2: (a) In the vortex experiments, a collimated, single-mode Gaussian
passes through a hologram on a transmissive spatial light modulator (SLM),
from which the first diffracted order is imaged onto a series of
$\pi$-converters (labelled “$\pi$”). The vortex beam is measured before the
first $\pi$-converter (initial state $|A\rangle$) and then after (intermediate
state $|B\rangle$), which makes a half circle trajectory. Going through the
second $\pi$-converter, the final state of the beam coincides with the initial
state after making a complete circle on the SoM. These circuits are
analytically calculated and plotted [34]. (b) In the polarization experiments,
a collimated, single-mode Gaussian is locked to a horizontal polarization with
the first linear polarizer (“Lin. Pol.”). At the beamsplitter, one arm is
unmodified to be used to interfere with the other arm after it passes through
the half-waveplate twice (flip mirror is up)—this produces a full circle
circuit, from which the total phase is measured. Dynamic phase is measured
after one transit through the waveplate (flip mirror is down) with a linear
polarizer, locked to the same polarization as the first “Lin. Pol.” The
orientation of the waveplate, $\eta$, is rotated to increase $\beta$ on the
PS.

The beam in state $|A\rangle$ is then propagated through the first
$\pi$-converter. Between the mode converters, we measure intermediate state
$|B\rangle$ and characterize the mode using the same methods used to
characterize state $|A\rangle$. With both fields now distilled into discrete
data arrays, we can immediately evaluate the total phase of Eqn. 1a and the
inner product of Eqn. 6. To complete the circle, the beam then transits an
identical, second $\pi$-converter. However, rather than performing a second
measurement, we simply scale the intermediate-measured phases by two, as
discussed above.

Figure 3: (a) Dynamic (blue), total (red), and geometric (green) phases of
vortex transformations driven by a series of two $\pi$-converters, as a
function of $\beta$. Dots with error bars are experimental data (average and
standard deviation of five measurements) and curves are calculated (not fit to
data) from Eqns. 1a, 1b, and 1c. (b) Similar phase measurements in the
polarization domain, plotted with the same methods. Here, the waveplate was
rotated up to $\eta=90\degree$, corresponding to $\beta=0\degree\rightarrow
180\degree$. (c) Sample analytical trajectories for the corresponding values
of $\beta$ used in the data above, being functionally equivalent for the two
cases [34].

This procedure is applied for a range of initial states with the results
presented in Fig. 3 (a). The solid curves in the figure are theoretical
predictions obtained from Eqns. 1a (total), 1b (dynamic), and 1c (geometric).
The match of the measurement to theory (without fitting parameters) is
striking. A detailed, step-by-step example is provided in the Supplementary
Information for a representative measurement.

Our second implementation of the new methodology is associated with spin
angular momentum (polarization), as shown in Fig. 2 (b). We fix the input
state, $|A\rangle$, as horizontally polarized, and pass it through a rotatable
half-waveplate, the polarization analog of a $\pi$-converter for OAM. Fixing
the initial state and rotating the waveplate amounts to a change in the frame
of reference used for OAM, where it is the initial state that is changed and
the $\pi$-converter that is fixed. However, it is only the relative motion
that is relevant.

To measure the total phase accumulated over a half-circular trajectory, the
state emerging from a double transmission through the waveplate is interfered
with the initial state $|A\rangle$. The phase shift, $\Phi_{\mathrm{tot}}$,
does not vary with waveplate rotation as shown in red in Fig. 3 (b). To
measure the dynamic phase, we take the intermediate state, $|B\rangle$, after
one pass through the half-waveplate and transmit it through a linear polarizer
that is locked to the initial polarization of $|A\rangle$. This is a
projection measurement (Eqn. 6). The analysis then proceeds as for OAM but
with one caveat; the polarizer acts on field intensity rather than the field
itself. (See the Supplementary Information for details on converting this
measured Malus’ law to the dynamic phase.) A complete set of results are shown
in Fig. 3 (b) and, once again, is found to be in strong agreement with the
theoretical predictions shown as solid curves.

In conclusion, this Letter presents a general method for measuring the
geometric phase for small circle circuits in an SU(2) parameter space, based
on mode measurements and projections of the initial and half-circle states.
The approach is implemented for the first experimental measurements of small
circle geometric phase for both orbital and spin angular momentum
transformations of Gaussian laser beams. The derivation associated with Eqns.
3 through 7 is not restricted to the two experimental settings we have
performed here, which allows for our approach to be easily adapted. We
anticipate that the results of this paper will find application to a variety
of settings for which geometric phase has not yet been measured in association
with small circles.

The authors acknowledge R. N. Mundell for assisting on the polarization
experiment.

Funding. W. M. Keck Foundation; NSF (1553905)

Disclosures. The authors declare no conflicts of interest.

Data availability. Data underlying all results presented are available from
the authors upon reasonable request.

## References

  * Berry [2010] M. Berry, Geometric phase memories, Nature Physics 6, 148 (2010).
  * Cohen  _et al._ [2019] E. Cohen, H. Larocque, F. Bouchard, F. Nejadsattari, Y. Gefen, and E. Karimi, Geometric phase from Aharonov–Bohm to Pancharatnam–Berry and beyond, Nature Reviews Physics 1, 437 (2019).
  * Zhang  _et al._ [2005] Y. Zhang, Y. W. Tan, H. L. Stormer, and P. Kim, Experimental observation of the quantum Hall effect and Berry’s phase in graphene, Nature 438, 201 (2005).
  * Noh  _et al._ [2020] J. Noh, T. Schuster, T. Iadecola, S. Huang, M. Wang, K. P. Chen, C. Chamon, and M. C. Rechtsman, Braiding photonic topological zero modes, Nature Physics 16, 989 (2020).
  * Yale  _et al._ [2016] C. G. Yale, F. J. Heremans, B. B. Zhou, A. Auer, G. Burkard, and D. D. Awschalom, Optical manipulation of Berry phase in a solid-state spin qubit, Nature Photonics 10, 184 (2016).
  * Wilczek and Shapere [1989] F. Wilczek and A. Shapere,  _Geometric Phases in Physics_, Advanced series in mathematical physics (World Scientific, 1989).
  * Nayak  _et al._ [2008] C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Non-Abelian anyons and topological quantum computation, Reviews of Modern Physics 80, 1083 (2008).
  * Quigg [2013] C. Quigg,  _Gauge Theories of the Strong, Weak, and Electromagnetic Interactions: Second Edition_ (Princeton University Press, 2013).
  * Pancharatnam [1956] S. Pancharatnam, Generalized theory of interference, and its applications - Part I. Coherent pencils, Proceedings of the Indian Academy of Sciences - Section A 44, 247 (1956).
  * Berry [1984] M. V. Berry, Quantal phase factors accompanying adiabatic changes, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 392, 45 (1984).
  * Simon [1983] B. Simon, Holonomy, the quantum adiabatic theorem, and Berry’s phase, Physical Review Letters 51, 2167 (1983).
  * Martinelli and Vavassori [1990] M. Martinelli and P. Vavassori, A geometric (Pancharatnam) phase approach to the polarization and phase control in the coherent optics circuits, Optics Communications 80, 166 (1990).
  * Bhandari [1991a] R. Bhandari, SU(2) phase jumps and geometric phases, Physics Letters A 157, 221 (1991a).
  * Bhandari [1991b] R. Bhandari, Evolution of light beams in polarization and direction, Physica B: Physics of Condensed Matter 175, 111 (1991b).
  * Tiwari [1992] S. C. Tiwari, Geometric phase in optics: Quantal or classical?, Journal of Modern Optics 39, 1097 (1992).
  * Bliokh  _et al._ [2019] K. Y. Bliokh, M. A. Alonso, and M. R. Dennis, Geometric phases in 2D and 3D polarized fields: Geometrical, dynamical, and topological aspects, Reports on Progress in Physics 82, 10.1088/1361-6633/ab4415 (2019).
  * Simon  _et al._ [1988] R. Simon, H. J. Kimble, and E. C. G. Sudarshan, Evolving Geometric Phase and Its Dynamical Manifestation as a Frequency Shift: An Optical Experiment, Physical Review Letters 61, 19 (1988).
  * Bhandari and Samuel [1988] R. Bhandari and J. Samuel, Observation of topological phase by use of a laser interferometer, Physical Review Letters 60, 1211 (1988).
  * Kwiat and Chiao [1991] P. G. Kwiat and R. Y. Chiao, Observation of a Nonclassical Berry’s Phase for the Photon, Physical Review Letters 66, 588 (1991).
  * Ortíz  _et al._ [2014] O. Ortíz, Y. Yugra, A. Rosario, J. C. Sihuincha, J. C. Loredo, M. V. Andrés, and F. De Zela, Polarimetric measurements of single-photon geometric phases, Physical Review A 89, 012124 (2014).
  * Maji  _et al._ [2019] S. Maji, P. Jacob, and M. M. Brundavanam, Geometric Phase and Intensity-Controlled Extrinsic Orbital Angular Momentum of Off-Axis Vortex Beams, Physical Review Applied 12, 1 (2019).
  * Karnieli  _et al._ [2019] A. Karnieli, S. Trajtenberg-Mills, G. Di Domenico, and A. Arie, Experimental observation of the geometric phase in nonlinear frequency conversion, Optica 6, 1401 (2019).
  * Hannonen  _et al._ [2020] A. Hannonen, H. Partanen, A. Leinonen, J. Heikkinen, T. K. Hakala, A. T. Friberg, and T. Setälä, Measurement of the Pancharatnam–Berry phase in two-beam interference, Optica 7, 1435 (2020).
  * van Enk [1993] S. J. van Enk, Geometric phase, transformations of gaussian light beams and angular momentum transfer, Optics Communications 102, 59 (1993).
  * Padgett and Courtial [1999] M. J. Padgett and J. Courtial, Poincare-sphere equivalent for light beams containing orbital angular momentum, Optics Letters 24, 430 (1999).
  * Courtial  _et al._ [1998] J. Courtial, K. Dholakia, D. A. Robertson, L. Allen, and M. J. Padgett, Measurement of the Rotational Frequency Shift Imparted to a Rotating Light Beam Possessing Orbital Angular Momentum, Physical Review Letters 80, 3217 (1998).
  * Galvez  _et al._ [2003] E. J. Galvez, P. R. Crawford, H. I. Sztul, M. J. Pysher, P. J. Haglin, and R. E. Williams, Geometric Phase Associated with Mode Transformations of Optical Beams Bearing Orbital Angular Momentum, Physical Review Letters 90, 4 (2003).
  * Habraken and Nienhuis [2010] S. J. Habraken and G. Nienhuis, Geometric phases in astigmatic optical modes of arbitrary order, Journal of Mathematical Physics 51, 10.1063/1.3456078 (2010).
  * Alonso and Dennis [2017] M. A. Alonso and M. R. Dennis, Ray-optical Poincaré sphere for structured Gaussian beams, Optica 4, 476 (2017).
  * Malhotra  _et al._ [2018] T. Malhotra, R. Gutiérrez-Cuevas, J. Hassett, M. R. Dennis, A. N. Vamivakas, and M. A. Alonso, Measuring Geometric Phase without Interferometry, Physical Review Letters 120, 1 (2018).
  * Milione  _et al._ [2012] G. Milione, S. Evans, D. A. Nolan, and R. R. Alfano, Higher order Pancharatnam-Berry phase and the angular momentum of light, Physical Review Letters 108, 1 (2012).
  * Voitiv  _et al._ [2022] A. A. Voitiv, M. T. Lusk, and M. E. Siemens, Tilted Poincaré sphere geodesics, Optics Letters 47, 1089 (2022).
  * Cisowski  _et al._ [2022] C. Cisowski, J. B. Götte, and S. Franke-Arnold, Colloquium: Geometric phases of light: Insights from fiber bundle theory, Reviews of Modern Physics 94, 31001 (2022).
  * Lusk  _et al._ [2022] M. T. Lusk, A. A. Voitiv, C. Zhu, and M. E. Siemens, The Anatomy of Geometric Phase for an Optical Vortex Transiting a Lens, Physical Review A 105, 052211 (2022).
  * Bhandari [1989] R. Bhandari, Synthesis of general polarization transformers. A geometric phase approach, Physics Letters A 138, 469 (1989).
  * Beijersbergen  _et al._ [1993] M. W. Beijersbergen, L. Allen, H. E. van der Veen, and J. P. Woerdman, Astigmatic laser mode converters and transfer of orbital angular momentum, Optics Communications 96, 123 (1993).
  * Aharonov and Anandan [1987] Y. Aharonov and J. Anandan, Phase change during a cyclic quantum evolution, Physical Review Letters 58, 1593 (1987).
  * Huang  _et al._ [2012] D. Huang, H. Timmers, A. Roberts, N. Shivaram, and A. S. Sandhu, A low-cost spatial light modulator for use in undergraduate and graduate optics labs, American Journal of Physics 80, 211 (2012).
  * Andersen  _et al._ [2019] J. Andersen, S. Alperin, A. Voitiv, W. Holtzmann, J. Gopinath, and M. Siemens, Characterizing vortex beams from a spatial light modulator with collinear phase-shifting holography, Applied Optics 58 (2019).
  * Andersen  _et al._ [2021] J. M. Andersen, A. A. Voitiv, M. E. Siemens, and M. T. Lusk, Hydrodynamics of noncircular vortices in beams of light and other two-dimensional fluids, Physical Review A 104, 1 (2021).

Experimental Measurement of Geometric Phase of Non-Geodesic Circles:
Supplementary Information

The following additional details are presented in the order they arise in the
manuscript. Citations refer to the same references above.

## 1\. Vortex Expression used in Experiments

A natural way to express an optical vortex mode using the spherical
coordinates of $(\xi,\theta)$ of Fig. 1 and in terms of a Hermite-Gaussian
(HG) basis is:

$\psi_{A}=\cos{\frac{\theta}{2}}\,e^{i\xi}\mathrm{HG}_{10}+\sin{\frac{\theta}{2}}\,e^{-i\xi}\mathrm{HG}_{01},$
(S.1)

where $\mathrm{HG}_{10}$ and $\mathrm{HG}_{01}$ correspond with HG modes on
the $\mathcal{I}_{0}$ and $\mathcal{I}_{90}$ axes, respectively, of the Sphere
of Modes (SoM) of Fig. 1. This expression is similar to Eqn. 3 in the
manuscript, but it uses the fixed spherical angles rather than the generalized
angles $(\delta,\beta)$.

While performing the experiments, however, we chose to employ the perspective
of “virtual tilt” [40]. Much like how one often speaks of the physical
orientation of diagonally-linear polarized light as being at $45\degree$
rather than calling it $90\degree$ as it is on the Poincaré sphere (PS),
“tilt” allows for a vortex-centric interpretation of orientations. The
expression we actually use for all calculations and for generating holograms
is:

$\psi_{A}=\frac{1}{\sqrt{3+\cos{2\theta_{\mathrm{tilt}}}}}\frac{4}{\sqrt{\pi}w_{0}^{2}}e^{-(x^{2}+y^{2})/w_{0}^{2}}\left[(x+iy\cos{\theta_{\mathrm{tilt}}})\cos{\xi_{\mathrm{tilt}}}+(y-ix\cos{\theta_{\mathrm{tilt}}})\sin{\xi_{\mathrm{tilt}}}\right].$
(S.2)

This expression works by taking a vortex on the north pole ($\mathcal{I}_{+}$
on Fig. 1), stretching it by angle $\theta_{\mathrm{tilt}}$, and rotating it
about its centroid by $\xi_{\mathrm{tilt}}$.

The conversion relationships between tilt angles,
$(\xi_{\mathrm{tilt}},\theta_{\mathrm{tilt}})$, and spherical coordinates used
in Fig. 1, $(\xi,\theta)$, are:

$\displaystyle\theta=$
$\displaystyle\arccos{\biggl{[}\frac{4\cos{(\theta_{\mathrm{tilt}})}}{3+\cos{(2\theta_{\mathrm{tilt}})}}\biggl{]}},$
(S.3) $\displaystyle\xi$ $\displaystyle=2\,\xi_{\mathrm{tilt}}.$ (S.4)

For the case of starting on the equator, $\beta=\xi$. So we produce the data
of Fig. 3 by increasing the programmed angle $\xi_{\mathrm{tilt}}$ from $0$ to
$90\degree$; we also measured $\xi_{\mathrm{tilt}}$ from $90\degree$ to
$180\degree$, which simply reproduces the same sized-circles ($\beta$ is only
unique up to $180\degree$) in reverse order.

## 2\. Generalization of Starting States

It is important to note that the projection measurement of Eqn. 4 in the
manuscript is not restricted to starting states on the equator. We chose this
setting only for clarity of presentation. One may choose any axis on which to
center their small circles, as well as any orthogonal basis vectors.

To demonstrate how one would retrieve $\beta$ for a circle with initial and
intermediate states arbitrarily oriented on the SoM or PS, we demonstrate how
it would be done for the center of circles used in the manuscript. Define a
vector that locates the center of the circles, $\vec{\alpha}$. If centered on
the $\mathcal{I}_{0}$ axis, $\vec{\alpha}=\hat{x}$. Next, define an initial
starting point on the sphere, $\vec{A}$ coinciding with $|A\rangle$ on the SoM
or PS. In terms of Cartesian coordinates, $\vec{A}$ lives on the sphere at

$\vec{A}=\hat{x}\,\sin{\theta}\cos{\xi}+\hat{y}\,\sin{\theta}\sin{\xi}+\hat{z}\,\cos{\theta}.$
(S.5)

The angular radius of the circle on the sphere, defined between the center
point and the arbitrary starting point, is then found from the inner product
between the two vectors:

$\begin{split}\beta&=\arccos{\left(\vec{\alpha}\cdot\vec{A}\right)}\\\
&=\arccos{\left(\hat{x}\cdot\vec{A}\right)}\\\
&=\arccos{\left(\sin{\theta}\cos{\xi}\right)}.\end{split}$ (S.6)

Again, we used here a specific center point of $\hat{x}$ for illustration. The
purpose of showing this generalization is to show that $\beta$ can be
determined for any initial point belonging to any small circle, and so Eqns. 7
can always be determined analytically and measured experimentally.

## 3\. Step-by-Step Experimental Example: Vortices

Firstly, we provide more details comprising the schematic of Fig. 2 (a) in the
manuscript. A $\lambda=532$ nm, collimated, single-mode Gaussian passes
through a hologram on a transmission spatial light modulator (SLM), from which
the first diffracted order is imaged onto a series of $\pi$-converters
(labelled “$\pi$”). These lenses are locked in an orientation such that they
focus in the $y$-axis of the beam. The SLM is controlled by an Epson 83H
projector [38], and an Andor Zyla sCMOS detector captures transverse images of
the vortex beams. For these measurements, we only measure the beam before the
first $\pi$-converter (initial state) and then after (intermediate state),
which makes a half-circle trajectory. Going through the second
$\pi$-converter, the final state of the beam coincides with the initial state
after making a complete circle on the SoM.

To illustrate the process of measuring the data points of Fig. 3 (a), we
present one detailed example of extracting the three different phases
experimentally (the next section has similar details for polarization using
one example data set).

To construct an initial vortex state $\psi_{A}$—corresponding to $|A\rangle$
at point $A$ on the SoM—with an SLM as depicted in Fig. 2 (a), we program a
hologram of the form:

$\mathrm{Hologram}(x,y)=\mathrm{abs}\left[e^{i\,\mathrm{arg}[\psi_{A}]}+e^{ik_{\mathrm{g}}x}\right]\times\frac{\mathrm{abs}(\psi_{A})}{\mathrm{max}[\mathrm{abs}(\psi_{A})]},$
(S.7)

using the vortex definition $\psi_{A}$ of Eqn. S.2.

After $4f$-imaging the vortex generated at the SLM and before the first
$\pi$-converter, we measure state $\psi_{A}$ at, for example, SoM coordinates
$(\xi=40\degree,\theta=90\degree)$. The amplitude (square-root of the measured
intensity) and phase (recovered from four phased-stepped inteferograms via the
method of phase-shifting digital holography [39]) are combined to produce the
numerical complex field of the initial state. This is illustrated in Fig. S.1,
including also intermediate state $\psi_{B}$ after the first $\pi$-converter,
corresponding to $|B\rangle$ at point $B$ on the SoM.

The Gaussian reference beam used to measure the phase of $\psi_{B}$ acquires
the same total phase of $\Phi_{\mathrm{tot}}=-\pi/2$ after transit through the
$\pi$-converter [36]. Therefore, since the phase of $\psi_{B}$ is measured
relative to this Gaussian, this total phase is not actually measured. To
remedy this, we include an extra factor of $e^{-i\pi/2}$ on the measured
intermediate fields based on this long-established mode-converter theory [36];
this is depicted in the schematic expression of $\psi_{B}$ in Fig. S.1.

Figure S.1: Example experimental field measurements for initial vortex
condition $\psi_{A}$ at $(\xi=40\degree,\theta=90\degree)$ and intermediate
condition $\psi_{B}$, which has an acquired phase factor of
$\Phi_{\mathrm{tot}}$ that was negated by using a reference Gaussian beam that
acquires the same phase shift. The geometric phase can be mined out of the
overlaps between the two numerical fields presented here, illustrative of all
data points of the paper.

The first step to use these measurements to calculate the total, dynamic, and
geometric phases is to use these numerical fields, as produced in Fig. S.1, to
find the overlap between the states. This is done as shown below in Eqn. S.8,
where we integrate over all of the camera pixels that contain the beam:

$\begin{split}\mathrm{``overlap^{\prime\prime}}=\int_{\mathrm{all\,pixels}}\psi_{A}^{*}\psi_{B}\,dx\,dy.\end{split}$
(S.8)

This is the calculation of $\langle A|B\rangle$, Eqn. 4 of the manuscript.

For the data shown in Fig. S.1, “overlap”$=0.06-0.75i$, for one set of the
five sets of data for that $\beta$ value. Once “overlap” is calculated like
this, we find the three respective phases by:

  1. 1.

Total Phase. For this particular initial condition of $\beta=40\degree$:

$\Phi_{\mathrm{tot}}|_{\beta=40\degree}=\mathrm{arg}\left(\mathrm{overlap}\right)\times\mathrm{sgn}\left(\cos{(40\degree)}\right).$
(S.9)

Multiplying by $\mathrm{sgn}\left(\cos{\beta}\right)$ ensures that the total
phase remains constantly negative for any value of $\beta$, to be consistent
with the solid-angle definition of geometric phase for partial-circle arcs.
(Alternatively, one may simply take modulo $2\pi$ at the end for the full-
circle results—which is what we’ve invoked to do this.)

For the same example data and same “overlap”$=0.06-0.75i$ above, this results
in
$\Phi_{\mathrm{tot}}|_{\beta=40\degree}=\mathrm{arg}\left(0.06-0.75i\right)\times(1)=-1.49$
radians, approximately the expected value of $-\pi/2$. We double this value to
get the phase of the complete circle, $-2.98$ radians.

  2. 2.

Dynamic Phase. Once total phase is found, we can calculate the dynamic phase
as shown in Eqn. 7b. Since we have the total phase from the step above, we
simply apply Eqn. 6.

$\Phi_{\mathrm{dyn}}|_{\beta=40\degree}=\Phi_{\mathrm{tot}}|_{\beta=40\degree}\times-\mathrm{Im}\left(\mathrm{overlap}\right).$
(S.10)

For the same “overlap”$=0.06-0.75i$ above, the dynamic phase is found to be:

$\begin{split}\Phi_{\mathrm{dyn}}|_{\beta=40\degree}&=(-2.98\,\mathrm{rad})\times-\mathrm{Im}\left((0.06-0.75i)\right)\\\
&=(-2.98\mathrm{rad})\times\left(0.75\right)\\\
&=-2.2\,\mathrm{rad}.\end{split}$

  3. 3.

Geometric Phase. The geometric phase is then found by:

$\Phi_{\mathrm{geo}}|_{\beta=40\degree}=\Phi_{\mathrm{tot}}|_{\beta=40\degree}-\Phi_{\mathrm{dyn}}|_{\beta=40\degree}.$
(S.11)

So for the example data shown,
$\Phi_{\mathrm{geo}}|_{\beta=40\degree}=-2.98\,\mathrm{rad}-(-2.2\,\mathrm{rad})=-0.78\,\mathrm{rad}$,
for the complete circle.

All three steps are repeated identically for all 5 measurements for a given
$\beta$ angle, and then repeated identically for all $\beta$ angles as shown
in Fig. 3 (a). We emphasize that we never invoke explicitly a $\cos{\beta}$
relationship on the experimental data—the plotted trend in Fig. 3 (a) arises
directly from the numerical overlaps of the measured fields detailed above.
Lastly, the size of the error bars (for both experiments) are limited by
general optical stability of the set-up: floor/table vibrations and slight air
turbulence that cause “beam jitters” between one data acquisition and the
next.

## 4\. Step-by-Step Experimental Example: Polarization

Firstly, we provide more details comprising the schematic of Fig. 2 (b) in the
manuscript. A $\lambda=633$ nm, collimated, single-mode Gaussian is locked to
a horizontal polarization with the first linear polarizer (“Lin. Pol.”). A
half-waveplate ($\frac{\lambda}{2}$) is used to tune the outgoing intensity of
the light to prevent over-saturation on the camera. The beam is sent to a
beamsplitter, from which one arm is unmodified to act as reference for total
phase measurements. The other arm is sent through another half-waveplate which
is mounted in a motorized rotation stage (Thorlabs PRM1Z8 driven by Thorlabs
KDC101). This rotation changes the orientation of the waveplate, $\eta$, which
increases $\beta$ on the PS according to $\beta=2\eta$.

Total phase is measured from the interference between the arms of the
interferometer on a camera (Andor Zyla sCMOS), after the beam is reflected
back through the half-waveplate to make a complete circle circuit. An example
interference pattern is shown in Fig. S.2 (a), where (b) shows the vertical
slice used to measure shifts in the fringes. 25 interference images are taken
per each set of data.

For each data set (each $\beta$ value), and for each of the 25 images per set,
the total phase is extracted by fitting the slice of Fig. S.2 (b) to the
following model expression:

$\mathrm{``model^{\prime\prime}}=a\,\cos^{2}{\left(\frac{2\pi}{\lambda_{\mathrm{fringe}}}x+\phi_{\mathrm{shift}}\right)\,e^{-(x-x_{0})^{2}/w_{0}^{2}}+b},$
(S.12)

where the parameters used are: $a$ is the intensity, $b$ is the background
level, $\lambda_{\mathrm{fringe}}$ is the width of the fringes, $x_{0}$ is the
center location of the Gaussian-envelope with width $w_{0}$, and
$\phi_{\mathrm{shift}}$ is the sought-after phase shift of the fringes. The
average of this phase shift over all 25 images is the desired
$\Phi_{\mathrm{tot}}$. We set $\Phi_{\mathrm{tot}}|_{\beta=90\degree}=-\pi$
(the great circle case where there is no dynamic phase), and scale all other
total phase measurements to this value.

Once the total phase is measured, the flip mirror is flipped down after the
beam makes its first pass through the waveplate (a half-circle is completed on
the PS). To measure the dynamic phase, we make a projection measurement on
this beam by transmitting it through a linear polarizer locked to the initial
polarization of the beam. The result of these measurements, as the waveplate
is rotated, is shown in Fig. S.2 (c) which shows the expected trend of Malus
law. This is nearly the projection measurement (Eqn. 4 of the manuscript) when
we measure the power of the transmitted beam. The polarizer acts on the
intensity (rather than the field) of the beam, yielding Malus’s Law:
$\mathrm{I}_{\mathrm{transmitted}}=\mathrm{I}_{\mathrm{incident}}\,\cos^{2}{\beta}$.
This data in plotted in Fig. S.2 (c). Thus we modify the power meter
measurements to be
$\mathrm{A}_{\mathrm{transmitted}}=\mathrm{sgn}\left(\cos{\beta}\right)\sqrt{\mathrm{I}_{\mathrm{transmitted}}}$.
The values are then the fractional values of the constant total phase measured
previously:
$\Phi_{\mathrm{dyn}}(\beta)=\frac{\mathrm{A}(\beta)}{\mathrm{max}\left(\mathrm{A_{\mathrm{all}}}\right)}\times\Phi_{\mathrm{tot}}(\beta)$.
From this, $\Phi_{\mathrm{geo}}=\Phi_{\mathrm{tot}}-\Phi_{\mathrm{dyn}}$.

Figure S.2: (a) Example interference pattern used to measure the total phase
associated with a polarization small circle. (b) The corresponding “slice” of
(a), depicted with the red line. Shifts in the plotted fringes yield the total
phase. (c) Power meter reading of the complete data set,
$0\degree\leq\beta\leq 180\degree$. Malus’s law,
$\mathrm{I}_{\mathrm{measured}}\propto\cos^{2}{\left(\eta-0\degree\right)}$,
is clearly evident (the incident polarization is locked at $0\degree$).


2021

[1]J. C. Zamora

1]Instituto de Física, Universidade de São Paulo, São Paulo, Brazil

2]Departamento de Física, Universidad Nacional de Colombia, Bogotá, Colombia

# Simulation of the RIBRAS Facility with GEANT4

L. E. Tamayose<EMAIL_ADDRESS><EMAIL_ADDRESS>G. F. Fortino  D. Flechas
[  [

###### Abstract

A Geant4 simulation code was developed to perform realistic simulations of the
RIBRAS facility. A second order expansion of a finite solenoid field was
included to describe the beam optics with a good precision. A systematic study
of coil currents for several magnetic rigidities and focal points was
performed. Parameterizations of the coil currents for single and dual mode
operations were obtained. Dedicated routines were developed to simulate the
mechanism of direct reactions involving two and three particles in the final
state. The present simulations were employed to investigate the feasibility of
a Solenoidal Spectrometer with the RIBRAS facility. Our first results indicate
that the concept can be applied in the RIBRAS system under certain conditions.
Forthcoming studies both from simulations and experiment are already under
development.

###### keywords:

RIBRAS, GEANT4, Solenoid, Beam optics, Spectrometer

##  1 Introduction

Double superconducting solenoid systems for rare-isotope beams production such
as RIBRAS ribras1 ; ribras2  and TWINSOL twinsol1 ; twinsol2  have been used
for more than two decades in studies of nuclear reactions with radioactive
beams at energies near the Coulomb barrier. The strong magnetic fields
operated in the solenoids provide a compact and efficient method for secondary
beam production and selection based on the magnetic rigidity.

Dedicated magnetic field simulations are required in order to describe the
particle trajectories and to define the best position of blockers and
degraders. The Geant4 geant4  toolkit offers an excellent method to develop
very precise simulations by using the diverse functionalities available for
particle transport and its interaction with matter. Realistic simulations in
3D are possible with Geant4 due its versatile geometry definition and to the
methods used to integrate the particle equation of motion. The toolkit allows
the implementation of new routines where it is possible to include specific
reaction mechanisms at energies near the Coulomb barrier. Future projects and
new experimental developments can be investigated with the simulations. One of
these projects is the feasibility of a Solenoidal Spectrometer with the RIBRAS
facility that is discussed in this paper.

##  2 The RIBRAS Facility

The RIBRAS (Radioactive Ion Beam in Brazil) facility is a double
superconducting solenoid system for rare-isotope beams production The concept
is similar to the TWINSOL facility. Both systems use an in-flight magnetic
rigidity ($B\rho$) selection of reaction products emerging from a primary
target. A sketch of RIBRAS is shown in Figure. 1.

Figure 1:  Sketch of the RIBRAS facility. Two superconducting solenoids are
used to select and focus radioactive-ion beams in the scattering chambers. In
the single-solenoid mode, the particles are focus in the central scattering
chamber (chamber 2). A larger scattering chamber (chamber 3) is used for the
dual-solenoid mode.

The superconducting solenoids have a maximum field of 6.5 T and a 30 cm clear
warm bore where the beam-line passes through. Secondary beams are produced in
the primary target, which is located in front of the first magnet (see Fig.
1). The angular acceptance is in the range from $2^{\circ}$ to $15^{\circ}$ in
the laboratory system. The lower angular limit is defined by a Faraday cup
inserted at $0^{\circ}$, while the upper limit is given by the beam pipe
diameter. The solenoids act as thick lens focusing ions with a given $B\rho$
to a secondary target located downstream in a scattering chamber. In the
single-solenoid mode, the focal point corresponds to the center of the middle
scattering chamber (chamber 2 in Fig. 1) at 2.6 m distance from the production
target. Thus, the strength of the magnetic field is adjusted to focus the
secondary ion beams with a certain magnetic rigidity at the reaction target
that is located in the focal point. The double-solenoid mode provides a better
selectivity of the beam particles by using also a degrader and a set of
blockers along the flying path. The focal point in the double-solenoid mode is
at 5.9 m away from the primary target, located in the center of the scattering
chamber 3 (see Fig. 1).

In order to appropriately simulate the beam optics, the particle equation of
motion in a realistic magnetic field has to be calculated numerically. A
program developed for the TWINSOL setup phd_twinsol  has been adapted to
simulate the beam optics of the RIBRAS facility. However, the program only
provides information of 2D projections of the particle trajectories. In the
present work, we have extended this study for three dimensions. Now, our
program allows to investigate in a better detail the particle trajectories,
asymmetry effects, reactions in the solenoids and many more. In the next
Section, it is explained how the beam optics was implemented in the
simulation.

##  3 Beam Optics

###  3.1 Magnetic Field

In a simple approximation, the magnetic field of a solenoid can be assumed to
be zero outside of the solenoid region and uniform inside it. When a charged
particle enters to an uniform magnetic field, its trajectory is affected by
the Lorentz force that makes the particle rotate around the solenoid axis in
an helicoid form. The oscillation of this trajectory is defined by the
cyclotron frequency magnet1

$\omega_{\text{cyc}}=\frac{eqB}{\gamma m},$ (1)

where $eq$ is the particle charge, $B$ is the magnetic-field strength,
$\gamma$ is the relativistic Lorentz factor and $m$ is the particle mass.
However, this simple model is not able to describe the focusing effect of a
solenoid. A more realistic model assuming the finite size of the solenoid is
applied to calculate the off-axis field that focus the beam particles. In
general, the components of the axial symmetric magnetic field of a solenoid
can be expanded in a Taylor series as magnet1

$B_{z}(z,r)=B(z)-\frac{r^{2}}{2}B^{\prime\prime}(z)+\cdots,$ (2)

$B_{r}(z,r)=-\frac{r}{2}B^{\prime}(z)+\frac{r^{3}}{16}B^{\prime\prime\prime}(z)+\cdots,$
(3)

where $z$ is the position along the solenoid axis, $r$ is the radial distance
from the solenoid axis and the prime denotes the derivative with respect to
$z$. The radial part provides an impulse to the particle at the solenoid edges
that curves its trajectory, while the longitudinal component keeps the
particle in a helical trajectory. In our simulation, we used an expansion up
to second order to achieve a good precision in the particle trajectories.
Thus, Eqs. (2) and (3) can be expressed by phd_twinsol ; magnet2 ; magnet3

$B_{z}(z,r)=B_{0}\frac{\sqrt{R^{2}+\frac{1}{4}L^{2}}}{L}\left[\left(\frac{z_{-}}{(R^{2}+z_{-}^{2})^{1/2}}-\frac{z_{+}}{(R^{2}+z_{+}^{2})^{1/2}}\right)-\frac{3}{4}r^{2}R^{2}\left(\frac{z_{+}}{(R^{2}+z_{-}^{2})^{5/2}}-\frac{z_{-}}{(R^{2}+z_{+}^{2})^{5/2}}\right)\right],$
(4)

$B_{r}(z,r)=B_{0}\frac{R^{2}\sqrt{R^{2}+\frac{1}{4}L^{2}}}{L}\left[-\frac{r}{2}\left(\frac{1}{(R^{2}+z_{-}^{2})^{3/2}}-\frac{1}{(R^{2}+z_{+}^{2})^{3/2}}\right)+\frac{3}{8}r^{3}\left(\frac{4z_{-}^{2}-R^{2}}{(R^{2}+z_{-}^{2})^{3/2}}-\frac{4z_{+}^{2}-R^{2}}{(R^{2}+z_{+}^{2})^{3/2}}\right)\right],$
(5)

where $L$ is the solenoid length, $R$ is the bore radius and
$z_{\pm}=z\pm\frac{L}{2}$. The parameter $B_{0}=kI$ corresponds to the axial
field at the center of the magnet, which is proportional to the electric
current ($I$) in the solenoid coil and to the intrinsic constant ($k$) of the
solenoid magnet. The latter value is extracted from the experimental field map
measured for each solenoid. Figure 2 shows an example of the magnetic field
components as a function of the longitudinal coordinate for $r>0$. The $B_{z}$
component is nearly constant inside the solenoid coil (vertical dash-dotted
lines), and strongly decays in the outside region. The intensity of the radial
component is maximum at the solenoid edges, but with opposite sign due to the
magnetic field poles. In the inner region of the solenoid, the radial part is
negligible and the field is dominated by the longitudinal component.

A non-uniform magnetic field was created in the simulation by using the
G4MagneticField class with the parameterizations of Eqs. (4) and (5) to define
a field inside and outside of the solenoid volume. This class can easily be
replicated to more than one solenoid, as explained below.

Figure 2:  (color online). Off-axis magnetic field components in the
longitudinal coordinate ($z$) for a coil current of 28 Ampere. The vertical
dash-dotted lines represent the coil edges.

###  3.2 Single-Solenoid Mode

The single-solenoid mode corresponds to a setup with the first solenoid that
is in front of the production target. Thus, secondary particles are propagated
in the non-uniform field that was previously defined. The trajectories are
obtained by solving the respective equation of motion of the particle in the
magnetic field. Geant4 uses a Runge-Kutta algorithm to integrate numerically
the ordinary differential equations of motion. In our simulation, a fourth
order Runge-Kutta method was used with an accuracy of 1 $\mu$m for the
integrated steps.

The particles that are propagated in the magnetic field have helical
trajectories that cross the solenoid axis, as shown in Figure 3. The outer
circle is the solenoid bore and each blue line correspond to the trajectory of
different particles around the $z$ axis. As can be seen, all the trajectories
have the same initial and final positions due to the beam focus.

Figure 3:  (color online). Helical trajectories (2D projection) observed from
the production target. The outer circle is the solenoid bore. Each particle
trajectory has an inserted label.

In the present case, the focal point is located at the center of the
scattering chamber 2 (see Fig. 1). In an experiment, the beam focusing is
performed by sweeping the electric current of the magnet to find the maximum
beam intensity at the reaction target. Similarly, hundreds of simulations were
performed for different coil currents to obtain the best value that minimize
the beam radius in the focal plane. Figure 4 (top) shows an example of a beam
focus using several simulations with a step of 0.1 Ampere. The focus
corresponds to the minimum beam RMS (root mean square) radius at the reaction
target. The procedure was repeated for many beams covering a wide range of
magnetic rigidities. As can be observed in Figure 4 (bottom), the simulated
result has a linear trend in the studied magnetic rigidity region. The
simulation is in very good agreement with the solenoid currents obtained in
several experiments with different beams (Ref. zamora_mater  and references
therein). The parameterization for the simulated current (in Ampere) is

$I=7.28\sqrt{\frac{AE}{q^{2}}}-0.09,$ (6)

where $A$ is the mass number of the beam particle, $E$ is the kinetic energy
(in MeV) and $q$ is the charge-state number. This linear approximation allows
to extract in a very simple way the best solenoid current without having to
run a new simulation.

Figure 4:  (color online). Beam focusing for the single-solenoid mode. (Top)
Beam RMS radius at the target position for several electric currents in the
coil. (Bottom) Coil current for the beam focus as a function of the magnetic
rigidity.

###  3.3 Dual-Solenoid Mode

The G4MagneticField class was extended to use a non-uniform field in two
solenoids. The magnetic field definition allows to use solenoids with
different coil currents and constant factors. The magnetic field is also
defined in the outer region of the solenoids. In this configuration, a
superposition of the magnetic fields is included to integrate the particle
trajectory across the two solenoids setup. The secondary particles emerging
from the production target are focus by the first solenoid in the chamber 2,
where it is possible to include a beam degrader. The particles that are
accepted in the second solenoid are focus in the scattering chamber 3 that is
at 5.9 m away from the production target (Fig. 1). Figure 5 (top) shows an
example of the simulated particle trajectories with the dual-solenoid mode.
Given that the two solenoids have independent strengths, the phase-space
covered to parameterize the two coil currents (for any $B\rho$ value) is
considerably larger than for the single-solenoid mode. Thousands of
simulations varying the solenoid currents, $I_{1}$ and $I_{2}$, were performed
to find the best configuration that minimize the beam radius at the reaction
target. Figure 5 (bottom) shows an example of the phase-space covered to find
the beam focus for a given $B\rho$ value.

Figure 5:  (color online). Beam focusing for the dual-solenoid mode. (Top)
Example of the simulated trajectories with a beam focus. (Bottom) A 3D map of
the solenoid currents to obtain a global minimum that corresponds to the beam
focus.

Similar to the single-solenoid mode, the minimization procedure was performed
for several magnetic rigidities in order to obtain a parameterization. Also,
the simulation result for the two solenoids exhibit a linear trend that are
given by

$I_{1}=6.59\sqrt{\frac{AE}{q^{2}}}+2.88,$ (7)

$I_{2}=5.50\sqrt{\frac{AE}{q^{2}}}+2.25,$ (8)

where the subindex indicates the solenoid number. The parameterization of
above provides an educated guess for the focal currents, which can be used
during an experiment to find the best coil currents in a more efficient way.

##  4 Event Generator

A nuclear reaction class was developed in the simulation to force a certain
reaction process when the particle interacts with a target. The reaction class
is based on the  G4VParticleChange method that uses the whole physical
information of a primary particle to generate a secondary one using a given
model. Relativistic two and three body kinematics generators were implemented
using the Mandelstam variables byckling1973 . In this representation,
scattering angles, energy and momentum of the particles can be derived in a
compact and complete way in a Lorentz-invariant fashion. Several reaction
processes are included in the class, such as elastic/inelastic scattering,
pickup/stripping reactions, knockout, breakup, in-flight decay, etc. The
reaction mechanisms can be easily defined in an input macro. Additionally,
cross section tables can be included to generate the reaction events based
theoretical angular and energy distributions. As an example, Figure 6 presents
a three body output channel with 8B breakup on a 58Ni target at 26 MeV. The
kinematics plot shows in the x-axis the relative angle ($\theta_{\text{rel}}$)
between 7Be and proton, and in the y-axis the proton kinetic energy. Cross
section tables are included in order to have a more realistic reaction
generator.

Figure 6:  (color online). Reaction kinematics of 8B breakup on a 58Ni target
at 26 MeV. Proton energy vs. the relative angle between proton and 7Be
particles.

##  5 Solenoidal Spectrometer

One of the main objectives of the implementation of a Geant4 simulation of the
RIBRAS facility was to study the feasibility of the system to operate as a
Solenoidal Spectrometer. The concept was demonstrated by the HELical Orbit
Spectrometer (HELIOS) at Argonne National Laboratory WUOSMAA20071290 ;
LIGHTHALL201097 . The goal of this experiment is to investigate nuclear
reactions in inverse kinematics using a target and detector setup inside of an
uniform magnetic field. Several projects around the world with the same
concept are under development, such as SOLARIS solaris , ISS iss  and SSNAP
ALLEN2020161350 .

Reaction target and detectors are mounted inside of a solenoid. The beam
particles impinge the target and induce nuclear reactions. Thus, recoiling
particles produced in the reactions are transported in helical orbits back to
the axis where they are detected in a position-sensitive silicon array. After
one cyclotron period, $T_{\text{cyc}}=2\pi/\omega_{\text{cyc}}$, the particle
returns to the axis at a distance $z$ from the target. Using this information,
it is possible to correlate the detected position $z$ with the deposited
energy by WUOSMAA20071290

$E_{\text{lab}}=E_{\text{c.m.}}-\frac{1}{2}mV^{2}_{\text{c.m.}}+\left(\frac{mV_{\text{c.m.}}}{T_{\text{cyc}}}\right)z.$
(9)

$E_{\text{lab}}$ is the detected energy (laboratory system), $E_{\text{c.m.}}$
is the particle energy in the center-of-mass frame and $V_{\text{c.m.}}$ is
the center-of-mass velocity that is fixed by the beam energy. Therefore,
different excitation are represented by linear bands in a $E_{\text{lab}}$-$z$
spectrum.

The simulated setup of the RIBRAS Solenoid Spectrometer is comprised by 8
position-sensitive silicon detectors of $2\times 25$ cm2 area with 60
independent strips. The setup is divided in two groups of 4 detectors
symmetrically mounted around the reaction target inside the solenoid, as shown
in Figure 7 (top). The detector positions are adjusted to allow the free path
of the beam particles, and at the same time to cover a large solid angle. In a
single-solenoid mode, the beam particles pass through the hollow of the
silicon array and intercept the reaction target. This configuration is
adequate for studies of nuclear reactions with stable beams. An example of
reaction that can be studied with this setup is the $pn$ transfer reaction in
inverse kinematics with a ${}^{12}\text{C}$ beam on a ${}^{6}\text{Li}$
target. The recoiling $\alpha$ particles are emitted at backward angles due to
the positive $Q$ value. As observed in Figure 7 (top), the recoils are
transported in helical orbits to different positions of the detector. The
strength of the magnetic field in this case was 6 T. Transfer reactions
populating the ${}^{14}\text{N}$ ground state and the two first $1^{+}$ states
(3.9 MeV and 6.2 MeV) were simulated. Figure 7 (bottom) shows the
reconstructed spectrum of deposited energy vs. strip number (proportional to
the $z$ position). As expected from Eq. (9), different linear bands are
obtained for the excited states. Also, the $z$ position can be converted to
center-of-mass angle ($\theta_{\text{c.m.}})$. Even though a few trajectories
are not detected in the right position. This effect is visible in the thin
lines close to each band. Possibly geometric effect are producing a separation
in the kinematic bands. However, it can be solved by increasing the detector
size and distance from the target.

Figure 7:  (color online). Simulation of a Solenoidal Spectrometer with only
one solenoid of RIBRAS. (Top) Simulation of the recoiling particle
trajectories for the reaction
${}^{6}\text{Li}({}^{12}\text{C},\alpha){}^{14}\text{N}$ at 20 MeV. The
particles hit the detectors at different positions depending on the center-of-
mass angle in which they were produced. (Bottom) Spectrum of the particle
deposited energy vs. detector strip number (proportional to the z position).
The ground state and the two first $1^{+}$ states are clearly separated by
energy and center-of-mass angles.

In the dual-solenoid mode, the detector array is placed into the second
solenoid, using the first to focus the beam on the reaction target. This
configuration is important for studies with radioactive beams. Figure 8 shows
and example of simulated particle trajectories with a 10Be secondary beam. As
can be noted, the beam enters to the second solenoid almost parallel, but due
to the strong magnetic field the particles are focus on the target. In this
case, the reaction studied was $(p,d)$ which have a negative Q value. The
recoiling deuterons are transported in helical orbits to the downstream
detectors. Eq. (9) is also valid for this example. Further studies are needed
to optimize the position of the target and detectors for the dual-solenoid
mode with other beams.

Figure 8:  (color online). Simulation of a Solenoidal Spectrometer with the
two solenoids of RIBRAS. The dual-solenoid mode is used to investigate nuclear
reactions with secondary radioactive beams. The first solenoid focus the beam
in the center of the second solenoid where the detectors and target are
mounted. The simulated reaction was $p({}^{10}\text{Be},d){}^{9}\text{Be}$.
Deuterons are produced in forward direction and detected in the downstream
detector setup.

##  6 Summary

A Geant4 simulation code was developed to simulate the RIBRAS facility. A
dedicated class for the magnetic field of the superconducting solenoids was
implemented to describe with a good precision the beam optics. The simulation
allows to study the particle trajectories in three dimensions, which is of
great importance to investigate optical effects such as aberration. An
advantage of using Geant4 is the possibility to simulate also the particle
interaction with materials to obtain the energy loss. Dedicated simulations in
the dual-solenoid mode including degraders and blockers can be easily
implemented.

The nuclear reaction class implemented in the code provides a great advantage
to the user for defining the reaction mechanism and to include theoretical
cross sections. Most of the possible reaction channels in the RIBRAS facility
are included in the code.

The feasibility of a Solenoidal Spectrometer in the RIBRAS facility was
investigated for the first time. Preliminary geometry models indicate that
reactions inside the RIBRAS solenoids are possible. Forthcoming studies both
from simulations and experiment are already under development.

## Acknowledgements

This work was financially supported by Fundação de Amparo a Pesquisa do Estado
de São Paulo (FAPESP) under Grant Nos. 2018/04965-4 and 2016/17612-7. L.E.T.
thanks to Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq)
for the financial support within the Undergrad Research program at IFUSP.
G.F.F. thanks to Comissão Nacional de Energia Nuclear (CNEN) for the financial
support within the MSc. scholarship program. The authors acknowledges support
by the project INCT-FNA (464898/2014-5).

## Declarations

Conflict of Interest: The authors declare that they have no confict of
interest with any organization regarding the material discussed in this
manuscript.

## References

  *   * (1) Lichtenthäler, R., Lépine-Szily, A., Guimarães, V., Lima, G.F., Hussein, M.S.: Radioactive ion beams in brasil (RIBRAS). Braz. J. Phys. 33(2), 294 (2003). https://doi.org/10.1590/S0103-97332003000200025
  * (2) Lichtenthäler, R., et al.: Radioactive ion beams in brazil (RIBRAS). In: The 4th International Conference on Exotic Nuclei and Atomic Masses, pp. 733–736. Springer, Berlin, Heidelberg (2005)
  * (3) Becchetti, F.D., et al.: Production and use of 6He, 7Be, 8Li, 12B and metastable nuclear beams. Nucl. Instr. Meth. B 56-57, 554–558 (1991). https://doi.org/10.1016/0168-583X(91)96094-2
  * (4) Lee, M.Y., et al.: Study of nuclear reactions with intense, high-purity, low-energy radioactive ion beams using a versatile multi-configuration dual superconducting-solenoid system. Nucl. Instr. Meth. A 422(1), 536–540 (1999). https://doi.org/10.1016/S0168-9002(98)01081-X
  * (5) Agostinelli, S., et al.: Geant4—a simulation toolkit. Nucl. Instr. Meth. A 506(3), 250–303 (2003). https://doi.org/10.1016/S0168-9002(03)01368-8
  * (6) Lee, M.Y.: TWINSOL: a dual superconducting solenoid ion-optical system for the production and study of low-energy radioactive nuclear beam reactions. PhD thesis, University of Michigan (2002)
  * (7) Kumar, V.: Understanding the focusing of charged particle beams in a solenoid magnetic field. Am. J. Phys. 77(8), 737–741 (2009) https://doi.org/10.1119/1.3129242. https://doi.org/10.1119/1.3129242
  * (8) Muniz, S.R., Bagnato, V.S., Bhattacharya, M.: Analysis of off-axis solenoid fields using the magnetic scalar potential: An application to a zeeman-slower for cold atoms. Am. J. Phys. 83(6), 513–517 (2015) https://doi.org/10.1119/1.4906516. https://doi.org/10.1119/1.4906516
  * (9) Callaghan, E.E., Maslen, S.H.: The magnetic field of a finite solenoid. NASA Technical Note D-456 (1960)
  * (10) Zamora, J.C.: Estudo do espalhamento elástico dos isótopos 7Be, 9Be e 10Be em alvo de 12C. Master’s thesis, University of Sao Paulo, Brazil (2011). https://doi.org/10.11606/D.43.2011.tde-30092011-132427
  * (11) Byckling, E., Kajantie, K.: Particle Kinematics. A Wiley-Interscience publication. Wiley, ??? (1973). https://books.google.com/books?id=d_cNAQAAIAAJ
  * (12) Wuosmaa, A.H., et al.: A solenoidal spectrometer for reactions in inverse kinematics. Nucl. Instr. Meth. A 580(3), 1290–1300 (2007). https://doi.org/10.1016/j.nima.2007.07.029
  * (13) Lighthall, J.C., et al.: Commissioning of the HELIOS spectrometer. Nucl. Instr. Meth. A 622(1), 97–106 (2010). https://doi.org/10.1016/j.nima.2010.06.220
  * (14) SOLARIS. https://www.anl.gov/phy/solaris. Accessed: 2022-01-16
  * (15) ISS. https://isolde-solenoidal-spectrometer.web.cern.ch/. Accessed: 2022-01-16
  * (16) Allen, J., et al.: SSNAP: the solenoid spectrometer for nuclear astrophysics. Nucl. Instr. Meth. A 954, 161350 (2020). https://doi.org/10.1016/j.nima.2018.10.028. Symposium on Radiation Measurements and Applications XVII