# Pose-Guided Sign Language Video GAN with Dynamic Lambda

Christopher Kissel, Christopher Kümmel, Dennis Ritter, Kristian Hildebrand
Beuth University of Applied Sciences Berlin

###### Abstract

We propose a novel approach for the synthesis of sign language videos using
GANs. We extend the previous work of Stoll et al. Stoll et al. (2020a) by
using the human semantic parser of the Soft-Gated Warping-GAN from Dong et al.
(2018) to produce photorealistic videos guided by region-level spatial
layouts. Synthesizing target poses improves performance on independent and
contrasting signers. Therefore, we have evaluated our system with the highly
heterogeneous MS-ASL dataset Vaezi Joze and Koller (2019) with over 200
signers resulting in a SSIM of 0.893. Furthermore, we introduce a periodic
weighting approach to the generator that reactivates the training and leads to
quantitatively better results.

Figure 1: (Top:) During training, an input frame and a target frame are
selected from a sign language video to serve as ground truth to train our
models. The human semantic parser (A1) creates an initial parsing output of
the input frame, while OpenPose (A2) extracts key joints from the target
frame. These outputs are concatenated and the target parsing is predicted (B).
Finally, input frame, initial and predicted parsing and key joints are used as
input for generator (C). It generates an image that is as similar as possible
to the target frame input. Whereas the discriminator (D) uses the initial
frame, initial parsing, predicted parsing, target pose and the ground truth
target frame to check the generated image of C for authenticity. (Bottom:) For
testing, the input of C is the same as described above but the target pose is
taken from prerecorded key joints. Finally, C creates an image of the person
shown in the input frame in the desired pose.

##  1 Introduction

The WHO states that 430 million people in the world are deaf or suffer from
significant hearing loss. Affected people often use one of the many different
sign languages to communicate. In both directions of translation, the first
approaches to barrier-free communication between signers and non-signers have
emerged in recent years Bangham et al. (2000); Koller et al. (2015); Camgoz et
al. (2018). Our approach refers to the generation of personalized sign videos
that can be synthesized in a pose-guided pipeline. Our network architecture is
based on the findings of the Text2Sign Stoll et al. (2020a) paper, which
presented the first non-graphical avatar approach to continuous gesture video
generation. However, we only consider here the generation of glosses in
videos. Our sign language GAN generates a video over the space of all signs
from the MS-ASL dataset Vaezi Joze and Koller (2019) from a set of input-
target pairs. Our goal is to improve the approach by Stoll et al. and enhance
the quality of the generated videos to support a variety of signers and
different scene conditions. Hence, we apply our approach to the challenging
MS-ASL dataset. In this way, we hope to drive the development of real-world
sign language applications with personalized generated signers.

When training on highly variable data, such as the MS-ASL dataset, we have
observed that our approach first finds a stable equilibrium, but after some
time the discriminator loss starts to increase, which means that it is too
easily fooled (see Figure 2). It seems, however, that with a dataset as
diverse as MS-ASL, not enough uniform criteria can be learned to judge the
authenticity of an image when the generator has exceeded a certain threshold.
If the discriminator is too weak, the L1 loss of the generator will prevail
and the resulting images will be blurred. We have evaluated different
strategies to stabilize the training and found that a dynamic lambda, a value
that periodically weights the generator error against the discriminator,
reactivates the training and leads to an improvement in image quality as shown
in Figure 3.

##  2 Related Work

Sign language related tasks can be divided into three subcategories _Sign
Language Recognition (SLR)_ (Kasukurthi et al., 2019; Fang et al., 2017; Bheda
and Radpour, 2017), _Sign Language Translation (SLT)_ (Camgoz et al., 2018,
2020) and  _Sign Language Production (SLP)_ (Stoll et al., 2020a). Our
approach relates to SLP tasks to produce glosses and signs of performed sign
language. Sign language datasets mostly consist of annotated videos. Recent
approaches like "Word Level American Sign Language" (WLASL) (Li et al., 2020)
or "Microsoft - American Sign Language" (MS-ASL) (Vaezi Joze and Koller, 2019)
contain annotated high quality videos showing many signers performing various
glosses. For our approach we use the MS-ASL dataset as it contains over 25000
videos of more than 200 signers with very different scene conditions, which
comes closest to video data in the wild. Other datasets like SMILE (Ebling et
al., 2018) used by Stoll et al. (2020a) or (Duarte et al., In Press) even
contain multi-sensor data but fewer signers and unfortunately have not been
published yet. The use of GANs to tackle SLP problems is still rather new. An
important contribution from 2016 is the pix2pix (Isola et al., 2017) framework
which uses a U-Net (Ronneberger et al., 2015) generator with skip-connections.
It is a key component in our SLP approach as well. Synthesizing sign language
videos share some challenges with tasks of Warping-GANs. They aim to transfer
a person into previously unseen poses. There have been some remarkable
developments in this area in the last years (Dong et al., 2018; Balakrishnan
et al., 2018; Ma et al., 2017; Liu et al., 2019). In early 2020, Stoll et al.
released a comprehensive approach to translate spoken language into sign
language videos. The pipeline called _Text2Sign_ (Stoll et al., 2020a) solves
the task in two major stages. First the spoken glosses are translated into a
sequence of motion graphs by an encoder-decoder architecture (Text2Pose).
These poses are then passed to a sign language production GAN (Pose2Video),
which is also based on a pix2pix architecture and an extension (Wang et al.,
2018) of it. Our work picks up on the results of the _Pose2Video_ part and
further develops the approach. Later in 2020, more GAN applications for SLP
have been released (Saunders et al., 2020; Stoll et al., 2020b; Ventura et
al., 2020; Jiang et al., 2020). Which is an indication of the great interest
in the issue and show improved results. Contrary to all this work, our
approach concentrates on the synthesis of sign language in more real life
conditions using the MS-ASL dataset.

##  3 Method

###  3.1 Approach & Architecture

Our _SignLanguageGAN_ pipeline combines the approaches from Soft-Gated
Warping-GAN (Dong et al., 2018) and Text2Sign Stoll et al. (2020a) and is
displayed in Figure 1. At the first stage the human semantic parser (A1), a
U-Net (Ronneberger et al., 2015) autoencoder trained on the CIHP dataset, a
subset of the LIP dataset (Liang et al., 2019) creates an abstraction of the
input frame. It has been shown that a semantic translation can be beneficial
for spatial transformation of a person in an image (Dong et al., 2018;
Balakrishnan et al., 2018). The human semantic parser provides such a mapping
from limb categories to colour values. This initial parsing is now transformed
by a pix2pix network (Isola et al., 2017) (B) using the initial parsing and
the OpenPose (Cao et al., 2019) generated target pose (A2) representing the
desired gloss. Finally, the result is the target parsing, a semantic color map
of the person in the input frame who takes the pose represented by the
OpenPose skeleton. The second stage forms the actual adversarial model.
Starting with a generator (C) that was expanded by a second convolutional
encoder network. The latent spaces of both encoders are concatenated and a
decoder learns to reproduce the target image from them. The encoder-decoder
with skip connections is derived from U-Net architecture (Ronneberger et al.,
2015), which was utilized in pix2pix. However, unlike in the original U-Net
architecture, down-sampling is performed from two separate networks, while up-
sampling is executed by a single one which merges the information of its
counterparts. In its setup, the first encoder processes the generated target
parsing from the first stage together with the OpenPose generated target pose.
Encoder two, on the other hand, receives the input frame together with the
respective initial parsing. Both encoders learn their individual mapping into
a low dimensional representation. From the concatenation of these two latent
spaces the decoder learns how to create a new image of the person seen in the
initial frame performing the desired sign language gloss. This output is
passed to the SignLanguageGAN discriminator (D) for the adversarial training.
The discriminator determines whether the image originates from the generator
or if a real sample from the dataset was used. In particular, it is important
that the discriminator receives all information for assessment that was
processed in the generation.

###  3.2 Training with a Dynamic Lambda

While our architecture is a composite of several technologies, the training
process differs in certain aspects from the frameworks used. We apply a least
squares objective function in the discriminator as Mao et al. show that the L2
loss reduces vanishing gradients due to its sensitivity (Mao et al., 2016).
This effect is particularly noticeable with highly heterogeneous data such as
the MS-ASL (Vaezi Joze and Koller, 2019) collection. Another challenge in
working with such dynamic data became apparent as the training has advanced.
The equilibrium between discriminator and generator becomes increasingly
unstable, presumably because the discriminator has more and more difficulty
learning uniform criteria to determine. This results in images with more
artifacts and lower variance $\sigma$, as the generator contribution becomes
stronger, for which we choose a L1 loss in accordance with pix2pix. $\sigma$
influences contrast and structure in SSIM. To break this trend, we apply a
dynamic weighting between generator and discriminator to cyclically amplify
the impact of the discriminator. For this we use the lambda parameter, which
in pix2pix is a constant value. Several approaches have been evaluated to
adjust lambda periodically or adaptively. However, strategies that interfered
too much with the balance of the GAN were counterproductive. It have become
evident that a steady fluctuation yielded the best results (see Figure 2).
Compared to a static weighting the dynamic lambda actively manipulates the
loss leading to more updates for the generator. Those updates encourage the
generator to increase the output quality in order to fool the discriminator.
With the aim of having a high generator contribution at the beginning of the
training, we use a cosine period with an amplitude span of 1.5 epochs
$\approx\pi$ and a value range of 50 to 150. Our objective function consists
of the conditional least squares discriminator objective $V(D,G)$ and the $L1$
loss with the dynamic term where $i$ and $j$ determine the amplitude range.

$\lambda=i+\frac{cos([0,\pi])+1}{2}\cdot j$
$\min_{D}\min_{G}L_{SgnlGAN}=V(D,G)+L_{1}(G)\cdot\lambda$

By periodically manipulating the loss we observe an overall better performing
generator with fewer artifacts and higher variance. This leads to a slight
improvement in the SSIM and improved image results as shown in Figure 3(a,b).

Figure 2: SignLanguageGAN discriminator losses for 20 epochs of training with
dynamic lambda (blue) and static lambda (red). Figure 3: (a) (Top:) Static
lambda during training. (Bottom:) Dynamic lambda during training. (b) Curve
shows increased variance (in blue) when dynamic weighting is applied which
leads to higher SSIM.

##  4 Evaluation

###  4.1 Setup

We evaluate our approach on the MS-ASL Vaezi Joze and Koller (2019) dataset,
the most diverse sign language dataset in terms of number of signers and scene
conditions. We compare our generated frames with ground truth images and
determine the MSE, PSNR and SSIM to assess their quality. Note, unfortunately,
we cannot compare our results directly to Text2Sign Stoll et al. (2020a),
since the SMILE dataset Ebling et al. (2018) is not publicly available.

Architecture | PSNR | MSE | SSIM
---|---|---|---
pix2pix | 20.546 | 783.376 | 0.786
Soft-Gated Warping | 22.033 | 622.683 | 0.813
Table 1: pix2pix and Soft-Gated Warping architectures tested on 50 videos of
the MS-ASL dataset after three epochs of training.

###  4.2 Multi-Signer Generation Challenges

In order to evaluate the Soft-Gated Warping extension in our architecture we
have compared it with a network arrangement using pix2pix Isola et al. (2017)
instead. Table 1 shows that the human semantic parser extension leads to
improved results for all evaluated metrics using 50 randomly selected videos
from the MS-ASL dataset. We also have compared these multi-signer MS-ASL
results to results on homogenious _Gebärdenlernen_
111https://gebaerdenlernen.de/index.php?article_id=72 videos, which show a
single-signer in the same scene performing german sign language glosses. In
this way, we wanted to assess how the number of signers and the conditions of
the scene affect SLP quality. The results show that single signer setups (with
PSNR=25.634, MSE=194.174, SSIM=0.882) can easily outperform multi-signer
datasets (with PSNR=22.033, MSE=622.683, SSIM=0.813). That underlines our
assumption that our pose-guided approach would perform well in real-world
scenarios.

Since SignLanguageGAN is the first SLP approach to use the MS-ASL dataset it
is difficult to make a meaningful comparison. Nevertheless, we use a subset of
the MS-ASL dataset with 184 videos consisting of a total of 16601 frames as a
benchmark resulting in a PSNR of 26.044, MSE of 260.363 and SSIM of 0.893
after 300 epochs of training. We compare our results to the results for three
individual signers from Text2Sign (Stoll et al., 2020a, p. 901 Table 3), even
though our test setup is significantly different from the Text2Sign setup.
First, we did not have access to the SMILE dataset, which makes the comparison
much more difficult and also prevents us from using the pix2pix HD Wang et al.
(2018) architecture. Second, Text2Sign performs a fine-tuning for the three
signers. Since pix2pix HD performs better and SLP for only three signers is
not as difficult as for more than 200 we are not surprised that Text2Sign
results, with an average PSNR of 25.604, MSE of 166,614 and SSIM of 0.942 are
better than ours in terms of pure numbers. Nevertheless, we show that our
architecture performs well on extremely diverse data that also contains many
low quality samples and without special fine-tuning for individual signers.
Figure 4 shows five examples of our results.

Figure 4: Five examples of our produced sign language glosses. First line
shows the input image, second the image produced by our SignlanguageGAN, third
the ground truth.

##  5 Conclusion

In this paper, we present the SignLanguageGAN, a novel approach for photo-
realistic synthesis of sign language videos. In addition to previous work on
SLP we extended our GAN with a human semantic parser. This allows our system
to learn region-level spatial layouts in order to guide the generator to
produce more realistic appearing videos on a large variety of signers. The
introduction of a periodic lambda for weighting the discriminator is worth to
investigate further since it allows to create images with less artifacts in
long training sessions.

## References

  * Balakrishnan et al. (2018) Guha Balakrishnan, Amy Zhao, Adrian V Dalca, Fredo Durand, and John Guttag. 2018\.  Synthesizing images of humans in unseen poses.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pages 8340–8348.
  * Bangham et al. (2000) J. A. Bangham, S. J. Cox, R. Elliott, J. R. W. Glauert, I. Marshall, S. Rankov, and M. Wells. 2000.  Virtual signing: Capture, animation, storage and transmission-an overview of the ViSiCAST project.  In _IEE Seminar on Speech and Language Processing for Disabled and Elderly People (Ref. No. 2000/025)_ , pages 6/1–6/7.
  * Bheda and Radpour (2017) Vivek Bheda and Dianna Radpour. 2017.  Using deep convolutional networks for gesture recognition in american sign language.  _CoRR_ , abs/1710.06836.
  * Camgoz et al. (2018) Necati Cihan Camgoz, Simon Hadfield, Oscar Koller, Hermann Ney, and Richard Bowden. 2018.  Neural sign language translation.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pages 7784–7793.
  * Camgoz et al. (2020) Necati Cihan Camgoz, Oscar Koller, Simon Hadfield, and Richard Bowden. 2020.  Sign language transformers: Joint end-to-end sign language recognition and translation.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 10023–10033.
  * Cao et al. (2019) Z. Cao, G. Hidalgo Martinez, T. Simon, S. Wei, and Y. A. Sheikh. 2019\.  Openpose: Realtime multi-person 2d pose estimation using part affinity fields.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_.
  * Dong et al. (2018) Haoye Dong, Xiaodan Liang, Ke Gong, Hanjiang Lai, Jia Zhu, and Jian Yin. 2018.  Soft-gated warping-gan for pose-guided person image synthesis.  In _Proceedings of the 32nd International Conference on Neural Information Processing Systems_ , NIPS’18, page 472–482, Red Hook, NY, USA. Curran Associates Inc.
  * Duarte et al. (In Press) Amanda Duarte, S. Palaskar, Lucas Ventura, Deepti Ghadiyaram, Kenneth DeHaan, F. Metze, Jordi Torres, and Xavier Giró i Nieto. In Press.  How2sign: A large-scale multimodal dataset for continuous american sign language.  In _CVPR 2021_.
  * Ebling et al. (2018) Sarah Ebling, Necati Camgoz, Penny Braem, Katja Tissi, Sandra Sidler-Miserez, Stephanie Stoll, Simon Hadfield, Tobias Haug, Richard Bowden, Sandrine Tornay, Marzieh Razavi, and Mathew Magimai-Doss. 2018.  Smile swiss german sign language dataset.  In _LREC_.
  * Fang et al. (2017) Biyi Fang, Jillian Co, and Mi Zhang. 2017.  Deepasl: Enabling ubiquitous and non-intrusive word and sentence-level sign language translation.  In _Proceedings of the 15th ACM Conference on Embedded Network Sensor Systems_ , SenSys ’17, New York, NY, USA. Association for Computing Machinery.
  * Isola et al. (2017) Phillip Isola, Jun-Yan Zhu, Tinghui Zhou, and Alexei A Efros. 2017.  Image-to-image translation with conditional adversarial networks.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 1125–1134.
  * Jiang et al. (2020) Dehao Jiang, Mingqi Li, and Chunling Xu. 2020.  Wigan: A wifi based gesture recognition system with gans.  _Sensors_ , 20(17):4757.
  * Kasukurthi et al. (2019) Nikhil Kasukurthi, Brij Rokad, Shiv Bidan, and Aju Dennisan. 2019.  American sign language alphabet recognition using deep learning.  _CoRR_ , abs/1905.05487.
  * Koller et al. (2015) Oscar Koller, Jens Forster, and Hermann Ney. 2015.  Continuous sign language recognition: Towards large vocabulary statistical recognition systems handling multiple signers.  _Computer Vision and Image Understanding_ , 141:108–125.
  * Li et al. (2020) Dongxu Li, Cristian Rodriguez, Xin Yu, and Hongdong Li. 2020.  Word-level deep sign language recognition from video: A new large-scale dataset and methods comparison.  In _Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision_ , pages 1459–1469.
  * Liang et al. (2019) Xiaodan Liang, Ke Gong, Xiaohui Shen, and Liang Lin. 2019.  Look into person: Joint body parsing pose estimation network and A new benchmark.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 41:871–885.
  * Liu et al. (2019) Wen Liu, Zhixin Piao, Jie Min, Wenhan Luo, Lin Ma, and Shenghua Gao. 2019.  Liquid warping gan: A unified framework for human motion imitation, appearance transfer and novel view synthesis.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 5904–5913.
  * Ma et al. (2017) Liqian Ma, Xu Jia, Qianru Sun, Bernt Schiele, Tinne Tuytelaars, and Luc Van Gool. 2017.  Pose guided person image generation.  In _Advances in Neural Information Processing Systems_ , volume 30. Curran Associates, Inc.
  * Mao et al. (2016) Xudong Mao, Qing Li, Haoran Xie, Raymond Y. K. Lau, and Zhen Wang. 2016.  Multi-class generative adversarial networks with the L2 loss function.  _CoRR_ , abs/1611.04076.
  * Ronneberger et al. (2015) O. Ronneberger, P.Fischer, and T. Brox. 2015.  U-net: Convolutional networks for biomedical image segmentation.  In _Medical Image Computing and Computer-Assisted Intervention (MICCAI)_ , volume 9351 of _LNCS_ , pages 234–241. Springer.  (available on arXiv:1505.04597 [cs.CV]).
  * Saunders et al. (2020) Ben Saunders, Richard Bowden, and Necati Cihan Camgöz. 2020.  Adversarial training for multi-channel sign language production.  In _31st British Machine Vision Conference 2020, BMVC 2020, Virtual Event, UK, September 7-10, 2020_. BMVA Press.
  * Stoll et al. (2020a) Stephanie Stoll, Necati Cihan Camgöz, Simon Hadfield, and Richard Bowden. 2020a.  Text2sign: Towards sign language production using neural machine translation and generative adversarial networks. _International Journal of Computer Vision_ , 128(4):891–908.
  * Stoll et al. (2020b) Stephanie Stoll, Simon Hadfield, and Richard Bowden. 2020b.  Signsynth: Data-driven sign language video generation.  In _European Conference on Computer Vision_ , pages 353–370. Springer.
  * Vaezi Joze and Koller (2019) Hamid Vaezi Joze and Oscar Koller. 2019.  Ms-asl: A large-scale data set and benchmark for understanding american sign language.  In _The British Machine Vision Conference (BMVC)_.
  * Ventura et al. (2020) Lucas Ventura, Amanda Cardoso Duarte, and Xavier Giró-i-Nieto. 2020.  Can everybody sign now? exploring sign language video generation from 2d poses.  _CoRR_ , abs/2012.10941.
  * Wang et al. (2018) Ting-Chun Wang, Ming-Yu Liu, Jun-Yan Zhu, Andrew Tao, Jan Kautz, and Bryan Catanzaro. 2018.  High-resolution image synthesis and semantic manipulation with conditional gans.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 8798–8807.


# Latest results from the NA62 experiment at CERN

School of Physics and Astronomy, University of Birmingham
E-mail For the NA62 collaboration: R. Aliberti, F. Ambrosino, R. Ammendola, B.
Angelucci, A. Antonelli, G. Anzivino, R. Arcidiacono, T. Bache, M. Barbanera,
J. Bernhard, A. Biagioni, L. Bician, C. Biino, A. Bizzeti, T. Blazek, B.
Bloch-Devaux, V. Bonaiuto, M. Boretto, M. Bragadireanu, D. Britton, F.
Brizioli, M.B. Brunetti, D. Bryman, F. Bucci, T. Capussela, J. Carmignani, A.
Ceccucci, P. Cenci, V. Cerny, C. Cerri, B. Checcucci, A. Conovaloff, P.
Cooper, E. Cortina Gil, M. Corvino, F. Costantini, A. Cotta Ramusino, D.
Coward, G. D’Agostini, J. Dainton, P. Dalpiaz, H. Danielsson, N. De Simone, D.
Di Filippo, L. Di Lella, N. Doble, B. Dobrich, F. Duval, V. Duk, J.
Engelfried, T. Enik, N. Estrada-Tristan, V. Falaleev, R. Fantechi, V.
Fascianelli, L. Federici, S. Fedotov, A. Filippi, M. Fiorini, J. Fry, J. Fu,
A. Fucci, L. Fulton, E. Gamberini, L. Gatignon, G. Georgiev, S. Ghinescu, A.
Gianoli, M. Giorgi, S. Giudici, F. Gonnella, E. Goudzovski, C. Graham, R.
Guida, E. Gushchin, F. Hahn, H. Heath, E.B. Holzer, T. Husek, O. Hutanu, D.
Hutchcroft, L. Iacobuzio, E. Iacopini, E. Imbergamo, B. Jenninger, J. Jerhot,
R.W. Jones, K. Kampf, V. Kekelidze, S. Kholodenko, G. Khoriauli, A.
Khotyantsev, A. Kleimenova, A. Korotkova, M. Koval, V. Kozhuharov, Z.
Kucerova, Y. Kudenko, J. Kunze, V. Kurochka, V. Kurshetsov, G. Lanfranchi, G.
Lamanna, E. Lari, G. Latino, P. Laycock, C. Lazzeroni, M. Lenti, G. Lehmann
Miotto, E. Leonardi, P. Lichard, L. Litov, R. Lollini, D. Lomidze, A. Lonardo,
P. Lubrano, M. Lupi, N. Lurkin, D. Madigozhin, I. Mannelli, G. Mannocchi, A.
Mapelli, F. Marchetto, R. Marchevski, S. Martellotti, P. Massarotti, K.
Massri, E. Maurice, M. Medvedeva, A. Mefodev, E. Menichetti, E. Migliore, E.
Minucci, M. Mirra, M. Misheva, N. Molokanova, M. Moulson, S. Movchan, M.
Napolitano, I. Neri, F. Newson, A. Norton, M. Noy, T. Numao, V. Obraztsov, A.
Ostankov, S. Padolski, R. Page, V. Palladino, A. Parenti, C. Parkinson, E.
Pedreschi, M. Pepe, M. Perrin-Terrin, L. Peruzzo, P. Petrov, Y. Petrov, F.
Petrucci, R. Piandani, M. Piccini, J. Pinzino, I. Polenkevich, L. Pontisso,
Yu. Potrebenikov, D. Protopopescu, M. Raggi, A. Romano, P. Rubin, G. Ruggiero,
V. Ryjov, A. Salamon, C. Santoni, G. Saracino, F. Sargeni, S. Schuchmann, V.
Semenov, A. Sergi, A. Shaikhiev, S. Shkarovskiy, D. Soldi, V. Sugonyaev, M.
Sozzi, T. Spadaro, F. Spinella, A. Sturgess, J. Swallow, S. Trilov, P.
Valente, B. Velghe, S. Venditti, P. Vicini, R. Volpe, M. Vormstein, H. Wahl,
R. Wanke, B. Wrona, O. Yushchenko, M. Zamkovsky, A. Zinchenko.

###### Abstract:

NA62 is a fixed target kaon experiment at the CERN SPS which aims at measuring
the branching ratio of the $K^{+}\to\pi^{+}\nu\bar{\nu}$ decay with
$10\text{\,}\%$ precision. This ultra-rare kaon decay is theoretically
extremely clean and an ideal place to look for physics beyond the Standard
Model. The NA62 experiment has been running in 2016-2018 and has accumulated a
large amount of kaon decays. The final result of the
$K^{+}\to\pi^{+}\nu\bar{\nu}$ analysis using the 2016 dataset is discussed as
well as prospects for the analysis of the 2017 dataset. NA62 also has
sensitivity to various other rare and forbidden kaon decays, as well as some
hidden sector models. The results of a search for the Lepton Number Violating
decays $K^{+}\to\pi^{-}\ell^{+}\ell^{+}\ [\ell=e,\mu]$ is presented, along
with new limits on the existence of a dark photon.

##  1 Introduction

The $K^{+}\to\pi^{+}\nu\bar{\nu}$ is a Flavour Changing Neutral Current
process forbidden at tree level in the Standard Model (SM) and proceeding
through penguin and box diagrams. It is highly suppressed due to the quadratic
GIM mechanism and the small value of the CKM element $\left|V_{td}\right|$,
which makes it an ideal channel to search for physics beyond the SM (BSM)[1,
2, 3, 4, 5, 6]. The SM prediction is
$\mathcal{B}(K^{+}\to\pi^{+}\nu\bar{\nu})=$(8.4\pm
1.0)\text{\times}{10}^{-11}$$[7], where the dominant sources of uncertainties
are the CKM parameters. The previous experimental result is
$\mathcal{B}(K^{+}\to\pi^{+}\nu\bar{\nu})=\left(17.3^{+11.5}_{-10.5}\right)\times
10^{-11}$ obtained from 7 event candidates at the E787 and E949
experiments[8]. The precision of the measurement does not allow us to draw any
conclusion on any deviation from the SM prediction.

The total Lepton Number is a conserved quantity in the SM. However, Lepton
Number Violation (LNV) is predicted in some BSM models. The kaon decays
$K^{+}\to\pi^{-}\ell^{+}\ell^{+}\ [\ell=e,\mu]$ would then be allowed and
proceed via an intermediate mediator such as a Majorana Neutrino[9, 10].
Previous experimental results set the following limits at $90\text{\,}\%$ CL
on the branching fraction of such channels:
$\mathcal{B}(K^{+}\to\pi^{-}e^{+}e^{+})<$6.4\text{\times}{10}^{-10}$$[11] and
$\mathcal{B}(K^{+}\to\pi^{-}\mu^{+}\mu^{+})<$8.6\text{\times}{10}^{-11}$$[12].

An extension of the SM adding a new $U(1)$ gauge-symmetry sector is one of the
proposed solutions to explain the abundance of dark matter in the universe. In
the simplest realisation of this scenario[13, 14], the vector mediator field
$A^{\prime}$ (dark photon) of mass $M_{A^{\prime}}$ associated to the new
$U(1)$ interacts with the SM through a kinetic mixing term whose strength is
given by the parameter $\epsilon$. In such a case, it is possible for the
$A^{\prime}$ to be produced in the decay chain $K^{+}\to\pi^{+}\pi^{0},\
\pi^{0}\to A^{\prime}\gamma$ where $\mathcal{B}(\pi^{0}\to
A^{\prime}\gamma)=2\epsilon^{2}\left(1-\frac{M_{A^{\prime}}^{2}}{M_{\pi^{0}}^{2}}\right)^{3}\times\mathcal{B}(\pi^{0}\to\gamma\gamma)\
.$ With the assumption that the $A^{\prime}$ does not decay into SM particles
inside the detector, a missing energy signature might reveal its presence.

##  2 The NA62 experiment

NA62 is a fixed target experiment at the CERN SPS whose aim is to measure
$\mathcal{B}(K^{+}\to\pi^{+}\nu\bar{\nu})$ with $10\text{\,}\%$ precision. In
order to collect the necessary $10^{13}$ decays, the experiment is served by
the $400\text{\,}\mathrm{GeV}\text{/}\mathrm{\text{ $c$ }}$ proton beam of the
SPS. The protons impinge on a beryllium target and produce an intense
secondary hadron beam at $75\text{\,}\mathrm{GeV}\text{/}\mathrm{\text{ $c$
}}$ with $1\text{\,}\%$ momentum spread (RMS) and containing about
$6\text{\,}\%$ of kaons.

The incoming kaon is tagged by a differential Cherenkov detector (KTAG) and
its momentum is measured by the GigaTracker (GTK) spectrometer. Those
detectors are followed by a $110\text{\,}\mathrm{m}\text{/}$ long vacuum tank
which contains the $60\text{\,}\mathrm{m}\text{/}$ fiducial decay volume (FV).
A STRAW spectrometer at the end of the FV allows the measurement of momenta of
the charged decay products, and a Ring-Imaging Cherenkov detector (RICH) is
used for their identification. A set of electromagnetic calorimeters (LAV,
LKr, IRC, SAC) located around and at the end of the vacuum tube provide
hermetic coverage for photons emitted in the $K^{+}$ decay with angles up to
$50\text{\,}\mathrm{mrad}\text{/}$. Hadronic calorimeters and muon veto
detectors (MUV1,2,3) are located after the LKr calorimeter to complement the
particle identification system. Further counters (CHANTI, CHOD, MUV0, HASC)
are installed at different locations for triggering purposes and to provide
additional rejection power against beam interactions with the detector and
multi-track kaon decays. A detailed description of the detector and the beam
line can be found in [15].

##  3 Results from the $K^{+}\to\pi^{+}\nu\bar{\nu}$ analysis

###  3.1 Signal selection

The signature of the signal consists of an incoming $K^{+}$, an outgoing
$\pi^{+}$ and no additional activity in the detector. The detected incoming
(upstream) track compatible with a kaon hypothesis and the outgoing
(downstream) tracks associated in time are used to reconstruct the decay
vertex at their closest distance of approach. The longitudinal position of the
vertex is required to be located within a $50\text{\,}\mathrm{m}\text{/}$ long
decay region. The identification of the downstream particle is performed by
combining the RICH measurement with the output of a multivariate classifier
using the LKr and the hadronic calorimeters, providing a $\pi^{+}$
identification efficiency of $64\text{\,}\%$ for a $\mu^{+}$ suppression of
$1\text{\times}{10}^{8}$. The remaining background consists mostly of the
$K^{+}\to\pi^{+}\pi^{0}$ decay, which is suppressed by a factor of
$3\text{\times}{10}^{8}$ thanks to the hermetic photon veto system (LAV, LKr,
IRC, SAC) rejecting any in-time energy deposit.

The main kinematic variable is the missing mass squared
$m_{\text{miss}}^{2}=(P_{K}-P_{\pi})^{2}$. In order to protect against mis-
reconstruction of the momenta, $m_{\text{miss}}^{2}$ is computed in three ways
using different combination of $P_{K}$ and $P_{\pi}$ measurements. The kaon
momentum $P_{K}$ can be measured directly with the GTK or taken as the nominal
beam momentum measured using a sample of $K^{+}\to\pi^{+}\pi^{+}\pi^{-}$
decays. The pion momentum $P_{\pi}$ can be measured either with the STRAW or
with the RICH. The analysis is performed in two separate regions in the
$m_{\text{miss}}^{2}$ variables, on either side of the
$K^{+}\to\pi^{+}\pi^{0}$ peak (R1 as
$00.01\text{\,}{\mathrm{GeV}}^{2}\text{/}{\mathrm{\text{ $c$ }}}^{4}$ and R2
as $0.0260.068\text{\,}{\mathrm{GeV}}^{2}\text{/}{\mathrm{\text{ $c$
}}}^{4}$). These boundaries are driven by the missing mass squared resolution
which is ${10}^{-3}\text{\,}{\mathrm{GeV}}^{2}\text{/}{\mathrm{\text{ $c$
}}}^{4}$. In both regions, the $\pi^{+}$ momentum must be in the range
$1535\text{\,}\mathrm{GeV}\text{/}\mathrm{\text{ $c$ }}$ to ensure at least
$40\text{\,}\mathrm{GeV}\text{/}$ of electromagnetic energy in the
calorimeters and to optimise the RICH $\pi^{+}/\mu^{+}$ separation.

###  3.2 Results from the 2016 dataset

The signal acceptance $A$, determined from Monte-Carlo (MC) simulations, is
$1\text{\,}\%$ in R1 and $3\text{\,}\%$ in R2. The number of kaon decays in
the FV is measured from a sample of $K^{+}\to\pi^{+}\pi^{0}$ decays selected
with control triggers without the extra activity rejection and requiring
$m_{\text{miss}}^{2}$ to be in the $K^{+}\to\pi^{+}\pi^{0}$ region. It is
estimated to be $N_{K}=$(1.21\pm 0.04)\text{\times}{10}^{11}$$. The single
event sensitivity is $SES=1/(A\cdot\epsilon\cdot N_{K})=\left(3.15\pm
0.01_{\text{stat}}\pm 0.24_{\text{syst}}\right)$, where $\epsilon$ takes into
account the inefficiencies due to the trigger and the random veto caused by
accidental activity. Multiplying this number by the SM branching fraction of
the $K^{+}\to\pi^{+}\nu\bar{\nu}$ decay, the expected number of signal events
is $0.267\pm 0.001_{\text{stat}}\pm 0.020_{\text{sys}}\pm 0.032_{\text{ext}}$.

The contribution of the $K^{+}\to\pi^{+}\pi^{-}e^{+}\nu$ decay to the
background is evaluated from 600M MC decays, which show a good agreement
across five validation samples. The other background estimates are obtained
using data-driven techniques and validated on control regions around the
signal regions. The upstream background is due to accidental activity in the
detector and interactions of the beam particles in the upstream detectors. It
is effectively suppressed by geometrical cuts and by the matching of the kaon
and pion tracks. The background contributions are summarised in Table 1. One
event was observed in R2, as shown in Figure 1, corresponding to a p-value of
$15\text{\,}\%$ for the background only hypothesis and an observed upper limit
of $\mathcal{B}(K^{+}\to\pi^{+}\nu\bar{\nu})<$14\text{\times}{10}^{-10}$$ at
$95\text{\,}\%$ CL [16].

| Expected events (R1 + R2)
---|---
Process | 2016 | 2017 (Preliminary)
$K^{+}\to\pi^{+}\pi^{0}(\gamma)$ | $0.064\pm 0.007_{\text{stat}}\pm 0.006_{\text{sys}}$ | $0.35\pm 0.02_{\text{stat}}\pm 0.03_{\text{sys}}$
$K^{+}\to\mu^{+}\nu(\gamma)$ | $0.020\pm 0.004_{\text{stat}}\pm 0.006_{\text{sys}}$ | $0.16\pm 0.01_{\text{stat}}\pm 0.05_{\text{sys}}$
$K^{+}\to\pi^{+}\pi^{+}\pi^{-}$ | $0.002\pm 0.001_{\text{stat}}\pm 0.002_{\text{sys}}$ | $0.015\pm 0.008_{\text{stat}}\pm 0.015_{\text{sys}}$
$K^{+}\to\pi^{+}\pi^{-}e^{+}\nu$ | $0.013^{+0.017}_{-0.012}|_{\text{stat}}\pm 0.009_{\text{sys}}$ | $0.22\pm{0.08}_{\text{stat}}$
$K^{+}\to\pi^{0}\ell^{+}\nu(\ell=\mu,e)$ | $<0.001$ | $0.012\pm 0.012_{\text{sys}}$
$K^{+}\to\pi^{+}\gamma\gamma$ | $<0.002$ | $0.005\pm 0.005_{\text{sys}}$
Upstream background | $0.050^{+0.090}_{-0.030}|_{\text{stat}}$ | Analysis ongoing
Total background | $0.152^{+0.092}_{-0.033}|_{\text{stat}}\pm 0.013_{\text{syst}}$ |
Table 1: Summary of the background estimates summed over the two signal
regions for the analysis of the 2016 dataset (2nd column) and preliminary
estimates based on the 2017 dataset (3rd column). Figure 1: Reconstructed
$m_{\text{miss}}^{2}$ as a function of the $\pi^{+}$ momentum for events
satisfying the $K^{+}\to\pi^{+}\nu\bar{\nu}$ selection, except the
$m_{\text{miss}}^{2}$ and $\pi^{+}$ momentum criteria. The grey area
corresponds to the distribution of MC signal events. The red contours define
the signal regions.

###  3.3 Future prospects

The analysis of the 2017 data is ongoing with a 2016-like event selection. The
analysis shows similar performance but benefits from several improvements in
the pileup treatment for the IRC and SAC, in the LKr reconstruction and the
usage of the RICH variables. Overall the $\pi^{0}$ rejection is
$40\text{\,}\%$ more efficient. The background is currently under study and
the preliminary estimates are summarised in Table 1. In particular, thanks to
the range of beam intensities collected in 2017, investigations on intensity
related effects are being made. They suggests that $\pi^{+}$ efficiencies, and
$\pi^{0}$ and $\mu$ rejections have no or very little dependence on the beam
intensity. In any case, both the signal and background scale linearly with the
intensity. The data sample collected in 2017 is about 10 times bigger than the
one used for the 2016 result. Considering $N_{K}=$(1.3\pm
0.1)\text{\times}{10}^{12}$$ and $SES=$(0.34\pm
0.04)\text{\times}{10}^{-10}$$, the expected number of SM events in 2017 is
$2.5\pm 0.4$.

##  4 Lepton Number Violation results

Searches for the Lepton Number Violating (LNV) decay modes
$K^{+}\to\pi^{-}\ell^{+}\ell^{+}[\ell=e,\mu]$ were performed using a subset of
the 2017 data sample corresponding to approximately 3 months of data taking. A
blind analysis procedure was adopted and an event selection, similar for both
channels, was developed. In both cases the corresponding
$K^{+}\to\pi^{+}\ell^{+}\ell^{-}$ SM channel was used for normalisation, which
results in the first order cancellation of systematic uncertainties. The main
sources of background are $\pi^{+}$ mis-identification and decay in-flight.
They are studied using a combination of MC and data-driven approaches. The
reconstructed mass spectra are shown in Figure 2. In the $\ell=e$ case, the
number of kaon decays $N_{K}=$(2.14\pm 0.07)\text{\times}{10}^{11}$$, the
signal acceptance of $4.98\text{\,}\%$ and the expected number of background
events of $0.16\pm 0.03$ for 0 observed signal events result in an upper limit
on the branching fraction
$\mathcal{B}(K^{+}\to\pi^{-}e^{+}e^{+})<$2.2\text{\times}{10}^{-10}$$ at
$90\text{\,}\%$ CL. In the $\ell=\mu$ case, the number of kaon decays
$N_{K}=$(7.94\pm 0.23)\text{\times}{10}^{11}$$, the signal acceptance of
$9.81\text{\,}\%$ and the expected number of background events of $0.91\pm
0.41$ for 1 observed signal event result in an upper limit on the branching
fraction
$\mathcal{B}(K^{+}\to\pi^{-}\mu^{+}\mu^{+})<$4.2\text{\times}{10}^{-11}$$ at
$90\text{\,}\%$ CL[17].

Figure 2: Reconstructed $\pi\ell\ell$ mass spectra for the LNV channels
$K^{+}\to\pi^{-}e^{+}e^{+}$ (top right) and $K^{+}\to\pi^{-}\mu^{+}\mu^{+}$
(bottom right) with the partial 2017 data sample. The corresponding SM
channels $K^{+}\to\pi^{+}e^{+}e^{-}$ (top left) and
$K^{+}\to\pi^{+}\mu^{+}\mu^{-}$ (bottom left) are used for normalisation.

##  5 Dark photon searches

A sample of tagged $\pi^{0}$ from $K^{+}\to\pi^{+}\pi^{0}$ decays with exactly
one $\gamma$ detected in the LKr is selected in a subset of the data collected
by NA62 in 2016, corresponding to $1\text{\,}\%$ of the total collected in
2016-18. This sample is used to search for $A^{\prime}$ peaks in the missing
mass squared distribution
$m_{\text{miss}}^{2}=(P_{K}-P_{\pi}-P_{\gamma})^{2}$, where $P_{K},P_{\pi}$
and $P_{\gamma}$ are the 4-momenta of the beam kaon, the $\pi^{+}$ and the
detected $\gamma$, respectively. A sliding search window of width $\pm
1\sigma_{m_{\text{miss}}^{2}}$ is used to count the number of events for each
$M_{A^{\prime}}$ hypothesis, where the mass resolution is estimated with MC
simulations. An independent sample of fully reconstructed
$K^{+}\to\pi^{+}\pi^{0}$ is used for the normalisation.

The most abundant source of background is $\pi^{0}\to\gamma\gamma$ with one
$\gamma$ undetected. A sample of $\pi^{0}\to\gamma\gamma$ decays where one
$\gamma$ converts into $e^{+}e^{-}$ upstream of the LKr is selected. It is
scaled to the signal sample in a side-band and used to estimate the background
contribution in the search region. The side band extends between
$5\text{\times}{10}^{-5}75\text{\times}{10}^{-5}\text{\,}{\mathrm{GeV}}^{2}\text{/}{\mathrm{\text{
$c$ }}}^{4}$, while the search region extends between
$75\text{\times}{10}^{-5}1765\text{\times}{10}^{-5}\text{\,}{\mathrm{GeV}}^{2}\text{/}{\mathrm{\text{
$c$ }}}^{4}$. A CLs statistical treatment is applied and no statistically
significant excess is observed, allowing the setting of new upper limits at
$90\text{\,}\%$ CL on $\epsilon^{2}$ at the level of ${10}^{-7}{10}^{-6}$ over
the dark photon mass range $60110\text{\,}\mathrm{MeV}\text{/}{\mathrm{\text{
$c$ }}}^{2}$[18].

##  6 Conclusion

The final result of the $K^{+}\to\pi^{+}\nu\bar{\nu}$ analysis from the 2016
dataset was reported. One candidate signal event was found, resulting in an
upper limit on $\mathcal{B}(K^{+}\to\pi^{+}\nu\bar{\nu})$ of
$14\text{\times}{10}^{-10}$ at $95\text{\,}\%$ CL. Preliminary estimates of
the background level in the 2017 data were shown. The final results of
searches for LNV kaon decay modes were presented, where the following upper
limits at $90\text{\,}\%$ CL are set:
$\mathcal{B}(K^{+}\to\pi^{-}e^{+}e^{+})<$2.2\text{\times}{10}^{-10}$$ and
$\mathcal{B}(K^{+}\to\pi^{-}\mu^{+}\mu^{+})<$4.2\text{\times}{10}^{-11}$$.
Finally limits on the $\epsilon^{2}$ kinetic mixing parameter at the level of
${10}^{-7}{10}^{-6}$ are reported for dark photon masses in the range
$60110\text{\,}\mathrm{MeV}\text{/}{\mathrm{\text{ $c$ }}}^{2}$.

## References

  * [1] M. Blanke et al., _JHEP_ 0903 (2009) 108 [hep-ph/0812.3803]
  * [2] G. Isidori et al., _JHEP_ 0608 (2006) 064 [hep-ph/060474]
  * [3] T. Blažek et al., _Int.J.Mod.Phys_ A29 (2014) no.27, 1450162 [hep-ph/1410.0055]
  * [4] A. Buras et al., _JHEP_ 1511 (2015) 166 [hep-ph/1507.08672]
  * [5] M. Blanke et al., _Eur.Phys.J._ C76 (2016) no.4, 182 [hep-ph/1507.06316]
  * [6] M. Bordone et al., _Eur.Phys.J._ C77 (2017) no.9, 618 [hep-ph/1705.10729]
  * [7] A. Buras et al., _JHEP_ 1511 (2015) 033 [hep-ph/1503.02693]
  * [8] A.V. Artamonov et al., _Phys.Rev.D_ 79 (2009) 092004 [hep-ex/0903.0030]
  * [9] L. S. Littenberg et al., _Phys.Lett. B_ 491 (2000) 285 [hep-ph/0005285]
  * [10] A. Atre et al., _JHEP_ 0905 (2009) 030 [hep-ph/0901.3589]
  * [11] R. Appel et al., _Phys.Rev.Lett._ 85 (2000) 2877. [hep-ex/0006003]
  * [12] J.R. Batley et al._Phys.Lett. B_ 769 (2017) 67 [hep-ex/1612.04723]
  * [13] B. Holdom, _Phys.Lett. B_ 166 (1986) 196
  * [14] B. Batell et al., _Phys.Rev. D_ 80 (2009) 095024 [hep-ph/0906.5614]
  * [15] E. Cortina Gil et al _JINST_ 12 (2017) no.05, P05025 [physics.ins-det/1703.08501]
  * [16] E. Cortina Gil et al _Phys.Lett._ B791 (2019) 156-166 [hep-ex/1811.08508]
  * [17] E. Cortina Gil et al [hep-ex/1905.07770]
  * [18] E. Cortina Gil et al _JHEP_ 1905 (2019) 182 [hep-ex/1903.08767]


# Surface criticality in the mixed-field Ising model with sign-inverted next-
nearest-neighbor interaction

Yuki Nakamura  Graduate School of Human and Environmental Studies, Kyoto
University, Kyoto 606-8501, Japan  Department of Physics, Kindai University,
Higashi-Osaka, Osaka 577-8502, Japan  Ryui Kaneko<EMAIL_ADDRESS>Department of Physics, Kindai University, Higashi-Osaka, Osaka 577-8502, Japan
Physics Division, Sophia University, Chiyoda, Tokyo 102-8554, Japan  Ippei
Danshita<EMAIL_ADDRESS>Department of Physics, Kindai University,
Higashi-Osaka, Osaka 577-8502, Japan

###### Abstract

Rydberg atoms in an optical tweezer array have been used as a quantum
simulator of the spin-$1/2$ antiferromagnetic Ising model with longitudinal
and transverse fields. We suggest how to implement the next-nearest-neighbor
(NNN) interaction whose sign is opposite to that of the nearest neighbor one
in the Rydberg atom systems. We show that this can be achieved by weakly
coupling one Rydberg state with another Rydberg state. We further study the
surface criticality associated with the first-order quantum phase transition
between the antiferromagnetic and paramagnetic phases, which emerges due to
the sign-inverted NNN interaction. From the microscopic model, we derive a
Ginzburg-Landau (GL) equation, which describes static and dynamic properties
of the antiferromagnetic order parameter near the transition. Using both
analytical GL theory and numerical method based on a mean-field theory, we
calculate the order parameter in the proximity of a boundary of the system in
order to show that the healing length of the order parameter logarithmically
diverges, signaling the surface criticality.

##  I Introduction

Rapid technological progresses in preparing and manipulating Rydberg atoms in
an optical tweezer array have led to successful applications of the systems as
quantum simulators of quantum Ising models [1, 2, 3]. Specifically, recent
experiments have reported observations of various magnetically ordered phases,
such as antiferromagnetic phases breaking ${\mathbb{Z}}_{n}$ symmetry
($n=2,3,$ and $4$) in 1D chain [4], phases with checkerboard, striated, and
star patterns in a square lattice [5], and three-sublattice ordered phases
with $1/3$ and $2/3$ fillings in a triangular lattice [6]. Such developments
have triggered renewed theoretical interest of quantum phases and quantum
phase transition in Ising-type models [7, 8, 9].

Since the spin-spin interaction of the Ising model simulated by Rydberg atoms
corresponds to van der Waals interaction between two atoms in a Rydberg state,
it decays steeply with the distance $r$ as $\propto r^{-6}$. While a situation
where the next-nearest-neighbor (NNN) interaction is comparable to the
nearest-neighbor (NN) one can be implemented with use of Rydberg-dressed
states [10, 11], it does not allow for changing the sign of the interactions.
Previous theoretical studies have shown that the NNN interaction whose sign is
opposite to that of the NN one is a key to accessing some interesting phases
and transitions, including the floating phase in a 1D chain [12] and quantum
tricriticality in a square lattice [13]. In this sense, for extending
applicability of Rydberg-atom quantum simulators, it is useful to enable the
sign inversion of the NNN interaction.

In this paper, we propose that the Ising model with sign-inverted NNN
interaction can be realized by using a Rydberg state coupled weakly to another
Rydberg state. As an interesting phenomenon that can be studied in the
proposed system, we investigate the surface criticality [14, 15, 16] emerging
near the first-order quantum phase transition between the antiferromagnetic
and paramagnetic phases in a square lattice. We derive a Ginzburg-Landau (GL)
equation describing the motion of the antiferromagnetic order parameter and
relate the coefficients in the GL equation to the parameters in the
microscopic model. By calculating the order parameter in the proximity of the
boundary with use of both analytical GL theory and numerical method based on a
mean-field theory, we show the logarithmic divergence of the healing length of
the order parameter, which is a characteristic feature of the surface
criticality.

The remainder of the paper is organized as follows. We set $\hbar=1$
throughout the paper.

##  II  How to create sign-inverted next-nearest neighbor interaction

In this section, we aim to show that the Ising model with the sign-inverted
NNN interaction can be created by using a Rydberg state coupled slightly with
another Rydberg state. We consider the system consisting of $M$ neutral atoms
spatially aligned on a square lattice with an optical tweezer array. We assume
that one of electronically excited states with large principal quantum number
$n$, namely Rydberg states, is coupled with the ground state via a two-photon
Raman process. We represent the ground state and the Rydberg state of the
$j$th atom as $|g_{j}\rangle$ and $|r_{j}\rangle$, respectively. We focus on
the interaction between the $j$th atom and the $l$th atom ($j\neq l$) whose
positions are denoted by ${\bm{x}}_{j}$ and ${\bm{x}}_{l}$, respectively. We
assume that the lattice spacing $a$ is sufficiently large so that the
interaction is negligible if either of the two atoms is in the ground state.
We also assume that the Rydberg state has even parity such that there is no
direct dipole-dipole interaction between the two Rydberg atoms. Hence, the
dominant interaction between the two Rydberg atoms is the van der Waals
interaction given by

$\displaystyle V_{rr}(|{\bm{x}}_{jl}|)=\frac{C_{6}}{|\bm{x}_{jl}|^{6}},$ (1)

where $C_{6}$ is the interaction coefficient. In addition, as illustrated in
Fig. 1, we consider a situation where the Rydberg state $|r_{j}\rangle$ is
slightly coupled with another Rydberg state with even parity, denoted by
$|\tilde{r}_{j}\rangle$. The interaction between two atoms in
$|r_{j}\tilde{r}_{l}\rangle$ state and that between two atoms in
$|\tilde{r}_{j}\tilde{r}_{l}\rangle$ state are of van der Waals type given by

$\displaystyle
V_{r\tilde{r}}(|{\bm{x}}_{jl}|)=\frac{\tilde{C}_{6}}{|\bm{x}_{jl}|^{6}},$ (2)
$\displaystyle
V_{\tilde{r}\tilde{r}}(|{\bm{x}}_{jl}|)=\frac{D_{6}}{|\bm{x}_{jl}|^{6}},$ (3)

where $\tilde{C}_{6}$ and $D_{6}$ are the interaction coefficients for
$V_{r\tilde{r}}(|{\bm{x}}_{jl}|)$ and
$V_{\tilde{r}\tilde{r}}(|{\bm{x}}_{jl}|)$, respectively. When the Hilbert
space of a single atom is spanned with only $|r_{j}\rangle$ and
$|\tilde{r}_{j}\rangle$, the Hamiltonian of the two atoms can be written as

$\displaystyle\hat{H}_{\rm two}=\hat{H}_{\rm two}^{(0)}+\hat{V},$ (4)

where

$\begin{split}\hat{H}_{\rm
two}^{(0)}&=-\Delta\left(\hat{n}_{j}+\hat{n}_{l}\right)+V_{\tilde{r}\tilde{r}}\left(|\bm{x}_{jl}|\right)\hat{n}_{j}\hat{n}_{l}\\\
&\quad+V_{rr}\left(|\bm{x}_{jl}|\right)\left(1-\hat{n}_{j}\right)\left(1-\hat{n}_{l}\right)\\\
&\quad+V_{r\tilde{r}}\left(|\bm{x}_{jl}|\right)\left(\left(1-\hat{n}_{j}\right)\hat{n}_{l}+\hat{n}_{j}\left(1-\hat{n}_{l}\right)\right)\end{split}$
(5)

and

$\hat{V}=\frac{\Omega}{2}\left(\hat{X}_{j}+\hat{X}_{l}\right).$ (6)

Here $\Delta$ and $\Omega$ denote the detuning and Rabi frequency of the
coupling between $|r_{j}\rangle$ and $|\tilde{r}_{j}\rangle$. The operators
$\hat{n}_{j}$ and $\hat{X}_{j}$ are defined as

$\displaystyle\hat{n}_{j}$ $\displaystyle=$
$\displaystyle|\tilde{r}_{j}\rangle\langle\tilde{r}_{j}|,$ (7)
$\displaystyle\hat{X}_{j}$ $\displaystyle=$
$\displaystyle|r_{j}\rangle\langle\tilde{r}_{j}|+|\tilde{r}_{j}\rangle\langle
r_{j}|.$ (8)

Figure 1: Schematic illustration of a Rydberg state $\ket{r}$ that is weakly
coupled to another Rydberg state $\ket{\tilde{r}}$.

Assuming $\Delta\gg\Omega$, we regard $\hat{V}$ as a perturbation to the
unperturbed part of the Hamiltonian $\hat{H}_{\rm two}^{(0)}$. The energies of
the unperturbed states, $\ket{r_{j}r_{l}}$ , $\ket{r_{j}\tilde{r}_{l}}$,
$\ket{\tilde{r}_{j}r_{l}}$, and $\ket{\tilde{r}_{j}\tilde{r}_{l}}$ are

$\displaystyle E^{\left(0\right)}_{rr}$ $\displaystyle=$ $\displaystyle
V_{rr}\left(|{\bm{x}}_{jl}|\right),$ (9) $\displaystyle
E^{\left(0\right)}_{\tilde{r}r}=E^{\left(0\right)}_{r\tilde{r}}$
$\displaystyle=$
$\displaystyle-\Delta+V_{r\tilde{r}}\left(|{\bm{x}}_{jl}|\right),$ (10)
$\displaystyle E^{\left(0\right)}_{\tilde{r}\tilde{r}}$ $\displaystyle=$
$\displaystyle-2\Delta+V_{\tilde{r}\tilde{r}}\left(|{\bm{x}}_{jl}|\right).$
(11)

When the unperturbed state is $|\Phi^{(0)}\rangle=|r_{j}r_{l}\rangle$, the
first-order perturbed state and the second-order perturbation energy read

$\displaystyle|\Phi\rangle$ $\displaystyle=$ $\displaystyle|r_{j}r_{l}\rangle$
(12) $\displaystyle\\!\\!\\!\\!\\!+\frac{\Omega}{\Delta-
V_{r\tilde{r}}(|{\bm{x}}_{jl}|)+V_{rr}(|{\bm{x}}_{jl}|)}\left(|r_{j}\tilde{r}_{l}\rangle\\!+\\!|\tilde{r}_{j}r_{l}\rangle\right),$

and

$E_{rr}^{\left(2\right)}=\frac{\Omega^{2}}{2\Delta}+\frac{F_{6}}{\left|{\bm{x}}_{jl}\right|^{6}+\tilde{R}^{6}},$
(13)

where $F_{6}=\frac{\Omega^{2}}{2\Delta^{2}}\left(\tilde{C_{6}}-C_{6}\right)$
and
$\tilde{R}=\left(-\frac{\tilde{C}_{6}-C_{6}}{\Delta}\right)^{\frac{1}{6}}$.
Hence, the energy $E_{rr}$ up to the second-order perturbation is given by

$E_{rr}\simeq
E_{rr}^{(0)}+E_{rr}^{(2)}=\frac{C_{6}}{|{\bm{x}}_{jl}|^{6}}+\frac{\Omega^{2}}{2\Delta}+\frac{F_{6}}{|{\bm{x}}_{jl}|^{6}+\tilde{R}^{6}}.$
(14)

This means that the interaction between the two atoms both in the Rydberg-
dressed Rydberg (RDR) state is given by $V^{\rm RDR}(|{\bm{x}}_{jl}|)=E_{rr}$.
Combining the interaction with the local magnetic field terms [1], one can
write the many-body Hamiltonian as a form of the spin-1/2 Ising model

$\hat{H}_{\rm IM}=\frac{1}{2}\sum_{j\neq
l}J_{jl}\hat{\sigma}_{j}^{z}\hat{\sigma}_{l}^{z}-\sum_{j}\left(\frac{\delta_{j}}{2}-U_{j}\right)\hat{\sigma}_{j}^{z}-\Gamma\sum_{j}\hat{\sigma}_{j}^{x},$
(15)

where $\hat{\sigma}_{j}^{z}$ and $\hat{\sigma}_{j}^{x}$ are the $z$ and $x$
components of the Pauli matrices at site $j$. We rewrite the interaction as
$J_{jl}=4V^{\rm RDR}(|{\bm{x}}_{jl}|)$ for the brevity of the notation.
$\delta_{j}$ and $2\Gamma$ correspond to the detuning and frequency of the
Rabi coupling between $|g_{j}\rangle$ and $|r_{j}\rangle$. The collective
contribution $U_{j}=\sum_{l\neq j}J_{jl}$ arises as a consequence of the
mapping from the original atomic representation to the effective spin-1/2
representation. While $U_{j}$ is inhomogeneous near the boundaries of the
system, such inhomogeneity can be cancelled by controlling the detuning
$\delta_{j}$ such that $\frac{\delta_{j}}{2}-U_{j}=H$ is homogeneous. In the
spin-1/2 representation, $\ket{\downarrow}$ corresponds to the atomic ground
state while $\ket{\uparrow}$ does to the RDR state.

One can see from Eq. (14) that if $F_{6}$ has a sign opposite to $C_{6}$ and
$|F_{6}|>|C_{6}|$, $V^{\rm RDR}(|{\bm{x}}_{jl}|)$ changes its sign at

$\displaystyle|{\bm{x}}_{jl}|=\left(\frac{|C_{6}|}{|F_{6}|-|C_{6}|}\right)^{\frac{1}{6}}\tilde{R}.$
(16)

For example, in Fig. 2, we show $V^{\rm RDR}(|{\bm{x}}_{jl}|)$ as a function
of $|{\bm{x}}_{jl}|$ for $\tilde{R}/a=1$, $C_{6}>0$, and $F_{6}/C_{6}=-1.841$.
There $V^{\rm RDR}(|{\bm{x}}_{jl}|=a)=-V^{\rm
RDR}(|{\bm{x}}_{jl}|=\sqrt{2}a)>0$, meaning that the NNN interaction has a
sign opposite to the NN one on a square lattice and their magnitudes are
equal. Moreover, the interaction for $|{\bm{x}}_{jl}|>\sqrt{2}a$ decays
abruptly as $\sim|{\bm{x}}_{jl}|^{-6}$ so that it is negligible if
$\max(|\Gamma|,|H|)\gg|V^{\rm RDR}(|{\bm{x}}_{jl}|=2a)|$. Under these
assumptions, the Ising model reads

$\hat{H}_{\rm IM}=J_{1}\sum_{\langle
j,l\rangle}\hat{\sigma}_{j}^{z}\hat{\sigma}_{l}^{z}-J_{2}\sum_{\langle\langle
j,l\rangle\rangle}\hat{\sigma}_{j}^{z}\hat{\sigma}_{l}^{z}-H\sum_{j}\hat{\sigma}_{j}^{z}-\Gamma\sum_{j}\hat{\sigma}_{j}^{x}$
(17)

where $J_{1}=4V^{\rm RDR}(|{\bm{x}}_{jl}|=a)$ and $J_{2}=-4V^{\rm
RDR}(|{\bm{x}}_{jl}|=\sqrt{2}a)$ denote the NN and NNN spin-spin interactions.

Figure 2: Effective interaction potential $V^{\rm RDR}(|{\bm{x}}_{jl}|)$
between the two atoms both in the RDR state as a function of the interatomic
distance ${\bm{x}}_{jl}$, where we set $F_{6}/C_{6}=-1.841$.

##  III Surface criticality near the first-order quantum phase transition

Figure 3: (a) Ground-state phase diagram of the Hamiltonian (17) in the
$(\Gamma/J_{1},H/J_{1})$-plane for $J_{1}=J_{2}$ and $J_{1}>0$. The thick
(thin) solid line represents the phase boundary of the first (second)-order
quantum phase transition between the antiferromagnetic (AFM) and paramagnetic
(PM) phases, which has been obtained in Ref. [13]. The filled circle marks the
quantum tricritical point (QTCP). The dashed line represents the contour of
$\phi_{0}=0.25$, where $\phi_{0}$ is the value of the order parameter in the
bulk, so that the gray shaded region roughly indicates the antiferromagnetic
(AFM) region where the GL theory is valid. (b) Magnified view of (a) near the
QTCP.

We briefly review the ground-state phase diagram of the mixed-field Ising
model with the sign-inverted NNN interaction of Eq. (17), which has been
previously obtained in Ref. [13]. In Fig. 3, we show the phase diagram in the
$(\Gamma/J_{1},H/J_{1})$-plane, where $J_{2}/J_{1}=1$ and $J_{1}>0$. There are
two distinct phases in the phase diagram. One is the antiferromagnetic (AFM)
phase, which is favored for small magnetic fields compared to the spin-spin
interaction. The AFM phase has a two-sublattice structure. The other is the
paramagnetic (PM) phase for large magnetic fields. When $H\ll\Gamma$, the
transition between the two phases is of second order. When $H\gtrsim\Gamma$,
it is of first order. The point at which the transition changes from second-
to first-order one on the transition line is the quantum tricritical point
(QTCP). Within a mean-field theory, the QTCP is located at

$\displaystyle(\Gamma/J_{1},H/J_{1})=\frac{32}{5\sqrt{5}}\times\left(2,1\right)\simeq
2.862\times\left(2,1\right).$ (18)

In the remainder of this section, we analyze the surface criticality that
emerges near the first-order transition. While the first order transition is
not accompanied by any critical behavior in the bulk of the system, the order
parameter near the boundary of the system exhibits a kind of critical behavior
resulting from the logarithmic divergence of its healing length [14, 15, 16].
The presence of the QTCP is helpful for studying the surface criticality in
the sense that one can widely control the surface critical region. While it is
known that the mixed field Ising model with van der Waals-type spin-spin
interaction on a square lattice, which has been already realized in
experiments of Ref. [5], exhibits first-order quantum phase transitions
between ordered and disordered phases, there is no QTCP in the ground-state
phase diagram of that model [17, 18].

###  III.1 Ginzburg-Landau theory

We assume that the system is in a parameter region indicated in Fig. 3 as the
gray shaded region, where the antiferromagnetic order parameter

$\displaystyle\psi({\bm{x}}_{j},t)=(2z_{1})^{-1}\sum_{\langle
l\rangle_{j}}\left(\langle\hat{\sigma}^{z}_{j}\rangle-\langle\hat{\sigma}^{z}_{l}\rangle\right)$
(19)

at the position ${\bm{x}}_{j}$ and the time $t$ is sufficiently small so that
the GL theory is approximately valid. Here, $\langle l\rangle_{j}$ means the
NN sites to site $j$ and $z_{1}$ denotes the coordination number of NN bonds,
which is 4 in the case of the square lattice considered in the present study.
There, the motion of the antiferromagnetic order parameter near the first-
order phase transition is described by the GL equation

$-w\frac{\partial^{2}\psi}{\partial
t^{2}}=\left(-C\nabla^{2}+s+u\psi^{2}+v\psi^{4}\right)\psi.$ (20)

The expression of the coefficients $s$, $u$, and $v$ in terms of the Ising
parameters have been obtained for zero temperature and $z_{1}J_{1}=z_{2}J_{2}$
in Ref. [13] as

$\displaystyle s$ $\displaystyle=$ $\displaystyle
J_{+}\left(1-\frac{\Gamma^{2}J_{+}}{B^{3}}\right),$ $\displaystyle u$
$\displaystyle=$
$\displaystyle\frac{\left(\Gamma^{2}-4H^{2}\right)\Gamma^{2}J_{+}^{4}}{2B^{7}},$
$\displaystyle v$ $\displaystyle=$
$\displaystyle\frac{3\left(12\Gamma^{2}H^{2}-8H^{4}-\Gamma^{4}\right)\Gamma^{2}J_{+}^{6}}{8B^{11}},$
$\displaystyle J_{+}$ $\displaystyle=$ $\displaystyle z_{1}J_{1}+z_{2}J_{2},$
$\displaystyle B$ $\displaystyle=$ $\displaystyle\sqrt{\Gamma^{2}+H^{2}},$

where $z_{2}$ is the coordination number of NNN bonds ($z_{2}=4$ for the
square lattice).

We next relate the remaining constants $w$ and $C$ to the Ising parameters
from the spin-wave dispersion relation. For this purpose, we assume that the
order parameter takes the form of

$\displaystyle\psi({\bm{x}},t)=\phi({\bm{x}})+\delta\psi({\bm{x}},t),$ (21)

where $\phi({\bm{x}})$ is a stationary order parameter and $\delta
m({\bm{x}},t)$ is a small fluctuation of the order parameter from
$\phi({\bm{x}})$. Substituting Eq. (21) into Eq. (20) and neglecting second
and higher order terms with respect to $\delta\psi$, one obtains the
stationary GL equation,

$\displaystyle\left(-C\nabla^{2}+s+u\phi^{2}+v\phi^{4}\right)\phi=0,$ (22)

and the linearized equation,

$\displaystyle-w\frac{\partial^{2}}{\partial
t^{2}}\delta\psi=\left(-C\nabla^{2}+s+3u\phi^{2}+5v\phi^{4}\right)\delta\psi.$
(23)

Assuming a paramagnetic solution $\phi=0$ of Eq. (22) and substituting a
plane-wave solution $\delta\psi=e^{i({\bm{p}}\cdot{\bm{x}}-\omega t)}$ into
Eq. (23), one obtains the dispersion relation of the small fluctuations,

$\omega=\sqrt{\frac{C}{w}p^{2}+\frac{s}{w}}.$ (24)

We also calculate the dispersion relation of the small fluctuations in the
paramagnetic phase by applying a mean-field approximation [19] directly to the
original model of Eq. (17) in order to compare its long-wavelength limit to
Eq. (24). We start with the Heisenberg’s equations of motion for
$\hat{\sigma}_{j}^{z}$ and
$\hat{\sigma}_{j}^{-}=\left(\hat{\sigma}_{j}^{x}-i\hat{\sigma}_{j}^{y}\right)/2$,

$\displaystyle\frac{d\hat{\sigma}_{j}^{z}}{dt}=-2\Gamma\hat{\sigma}_{j}^{y},$
(25)
$\displaystyle\begin{split}i\frac{d\hat{\sigma}_{j}^{-}}{dt}&=2J_{1}\sum_{\langle
l\rangle_{j}}\hat{\sigma}_{j}^{-}\hat{\sigma}_{l}^{z}-2J_{2}\sum_{\langle\langle
l\rangle\rangle_{j}}\hat{\sigma}_{j}^{-}\hat{\sigma}_{l}^{z}\\\
&\quad-2H\hat{\sigma}^{-}_{j}+\Gamma\hat{\sigma}_{j}^{z},\end{split}$ (26)

where $\langle\langle l\rangle\rangle_{j}$ means the NNN sites to site $j$. We
assume that the many-body wave function is approximated as a direct product of
local spin coherent states,

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!|\Phi_{\rm
MF}\rangle\\!=\\!\bigotimes_{j}\left(e^{-i\frac{\varphi_{j}}{2}}\cos\frac{\theta_{j}}{2}\ket{\uparrow}_{j}+e^{i\frac{\varphi_{j}}{2}}\sin\frac{\theta_{j}}{2}\ket{\downarrow}_{j}\right),$
(27)

where $\theta_{j}$ and $\varphi_{j}$ denote the elevation and azimuthal angles
of the spin direction at site $j$. Replacing $\hat{\sigma}_{j}^{z}$ and
$\hat{\sigma}_{j}^{-}$ in Eqs. (25) and (26) with their mean fields
$\langle\hat{\sigma}_{j}^{z}\rangle=\cos{\theta_{j}}$ and
$\langle\hat{\sigma}_{j}^{-}\rangle=\frac{1}{2}e^{-i\varphi_{j}}\sin{\theta_{j}}$,
one obtains the mean-field equations of motion,

$\begin{split}\frac{d\varphi_{j}}{dt}\sin{\theta_{i}}&=2J_{1}\sin{\theta_{j}}\sum_{\langle
l\rangle_{j}}\cos{\theta_{l}}-2J_{2}\sin{\theta_{j}}\sum_{\langle\langle
l\rangle\rangle_{j}}\cos{\theta_{l}}\\\
&\quad-2H\sin{\theta_{j}}+2\Gamma\cos{\theta_{j}\cos{\varphi_{j}}},\end{split}$
(28) $\frac{d\theta_{j}}{dt}=2\Gamma\sin{\varphi_{j}}.$ (29)

Since we are interested in the dispersion relation of the spin wave
excitations, we assume that the solutions of Eqs. (28) and (29) take the
following forms,

$\displaystyle\theta_{j}(t)$ $\displaystyle=$
$\displaystyle\bar{\theta}_{j}+\delta\theta_{j}(t),$ (30)
$\displaystyle\varphi_{j}(t)$ $\displaystyle=$
$\displaystyle\bar{\varphi}_{j}+\delta\varphi_{j}(t),$ (31)

where $\bar{\theta}_{j}$ and $\bar{\theta}_{j}$ are the stationary solutions
while $\delta\theta_{j}(t)$ and $\delta\varphi_{j}(t)$ are small fluctuations
from the stationary solutions. Substituting Eqs. (30) and (31) into Eqs. (28)
and (29), and neglecting second and higher order terms with respect to the
fluctuations, one obtains the equations for the stationary solutions,

$\begin{split}0=&\,\,2J_{1}\sin{\bar{\theta}_{j}}\sum_{\langle
l\rangle_{j}}\cos{\bar{\theta}_{l}}-2J_{2}\sin{\bar{\theta}_{j}}\sum_{\langle\langle
l\rangle\rangle_{j}}\cos{\bar{\theta}_{l}}\\\
&-2H\sin{\bar{\theta}_{j}}+2\Gamma\cos{\bar{\theta}_{j}}\cos{\bar{\varphi}_{j}},\end{split}$
(32) $0=\Gamma\sin{\bar{\varphi}_{j}},$ (33)

and the linearized equations of motion,

$\begin{split}&\frac{d}{dt}\delta\varphi_{j}=2J_{1}\left(-\sum_{\langle
l\rangle_{j}}\sin{\bar{\theta}_{l}}\delta\theta_{l}+\frac{\cos{\bar{\theta}_{j}}}{\sin{\bar{\theta}_{j}}}\sum_{\langle
l\rangle_{j}}\cos{\bar{\theta}_{l}}\delta\theta_{j}\right)\\\
&\quad-2J_{2}\left(-\sum_{\langle\langle
l\rangle\rangle_{j}}\sin{\bar{\theta}_{l}}\delta\theta_{l}+\frac{\cos{\bar{\theta}_{j}}}{\sin{\bar{\theta}_{j}}}\sum_{\langle\langle
l\rangle\rangle_{j}}\cos{\bar{\theta}_{l}}\delta\theta_{j}\right)\\\
&\quad-2H\frac{\cos{\bar{\theta}_{j}}}{\sin{\bar{\theta}_{j}}}\delta\theta_{j}-2\Gamma\left(\cos{\bar{\varphi}_{j}}\delta\theta_{j}+\frac{\cos{\bar{\theta}_{j}}}{\sin{\bar{\theta}_{j}}}\sin{\bar{\varphi}_{j}}\delta\varphi_{j}\right),\end{split}$
(34)
$\frac{d}{dt}\delta\theta_{j}=2\Gamma\cos{\bar{\varphi}_{j}}\delta\varphi_{j}.$
(35)

We assume that the system is in the paramagnetic phase, where
$\bar{\varphi}_{j}=0$. Moreover, we recall the assumption that
$z_{1}J_{1}=z_{2}J_{2}$, under which $\bar{\theta}_{j}={\rm
arctan}\frac{\Gamma}{H}\equiv\bar{\theta}$. In a $D$-dimensional hypercubic
lattice at $D\geq 2$ ($D=2$ for a square lattice), seeking a plane wave
solution for the fluctuations,
$\delta\theta_{j}=\theta_{\bm{p}}e^{i\left(\bm{p}\cdot\bm{x}_{j}-\omega
t\right)}$,
$\delta\varphi_{j}=\varphi_{\bm{p}}e^{i\left(\bm{p}\cdot\bm{x}_{j}-\omega
t\right)}$, one obtains the dispersion relation

$\begin{split}\omega=&\biggr{[}4\Gamma\biggr{\\{}\Gamma+H\frac{\cos{\bar{\theta}}}{\sin{\bar{\theta}}}-\frac{z_{1}J_{1}}{\sin{\bar{\theta}}}+z_{2}J_{2}\left(\frac{\mathrm{cos^{2}}\bar{\theta}}{\sin{\bar{\theta}}}-\sin{\bar{\theta}}\right)\\\
&\\!\\!\\!\\!+4J_{1}\sin{\bar{\theta}}\sum_{\alpha}\mathrm{sin^{2}}\left(\frac{p_{\alpha}a-\pi}{2}\right)+4J_{2}\sin{\bar{\theta}}\\\
&\\!\\!\\!\\!\times\sum_{\alpha<\beta}\left(\mathrm{sin^{2}}\left(\frac{p_{\alpha}a+p_{\beta}a}{2}\right)\\!+\\!\mathrm{sin^{2}}\left(\frac{p_{\alpha}a-p_{\beta}a}{2}\right)\right)\biggr{\\}}\biggr{]}^{\frac{1}{2}},\end{split}$
(36)

where $p_{\alpha}$ denotes the $\alpha$th component of the $D$-dimensional
vector ${\bm{p}}$. In a long-wavelength region, namely $|p_{\alpha}a-\pi|\ll
1$, Eq. (36) can be approximated as

$\omega\simeq\biggr{[}4\Gamma\biggl{\\{}\Gamma+H\frac{\cos{\bar{\theta}}}{\sin{\bar{\theta}}}-\frac{z_{1}J_{1}}{\sin{\bar{\theta}}}+z_{2}J_{2}\left(\frac{\mathrm{cos^{2}}\bar{\theta}}{\sin{\bar{\theta}}}-\sin{\bar{\theta}}\right)\\\
+4\Gamma\sin{\bar{\theta}}\left(J_{1}+2J_{2}\left(D-1\right)\right)\sum_{\alpha}\left(p_{\alpha}a-\pi\right)^{2}\biggr{\\}}\biggr{]}^{\frac{1}{2}}.$
(37)

Comparing this dispersion relation Eq. (37) with that obtained by means of the
GL theory Eq. (24), the constants $w$ and $C$ in the GL equation are
determined as

$\displaystyle
w=\frac{r}{4\Gamma\left(\Gamma+H\cot{\bar{\theta}}-\frac{z_{1}J_{1}}{\sin{\bar{\theta}}}+z_{2}J_{2}\left(\frac{\mathrm{cos}^{2}\bar{\theta}}{\sin{\bar{\theta}}}-\sin{\bar{\theta}}\right)\right)},$
(38) $\displaystyle C=w\times
4\Gamma\sin{\bar{\theta}}(J_{1}+2J_{2}(D-1))a^{2}.$ (39)

Thus, all the coefficients in the GL equation Eq. (20) have been explicitly
related to the parameters in the original Ising model of Eq. (17).

Figure 4: Spatial profile of the antiferromagnetic order parameter for
$(\Delta\Gamma/J_{1},H/J_{1})=(2.7\times 10^{-8},2.95)$ (a) and $(4.1\times
10^{-5},2.88)$ (b), where $\Delta\Gamma=\Gamma_{\rm t}-\Gamma$ and
$\Gamma_{\rm t}$ corresponds to the transition point for a given value of
$H/J_{1}$. The values of $\Gamma_{\rm t}/J_{1}$ within the GL theory are
5.456029819 for $H/J_{1}=2.95$ and 5.670719614 for $H/J_{1}=2.88$. Those
determined by numerically solving the mean-field equations (32) and (33) for
an infinite homogeneous system are 5.458428663 for $H/J_{1}=2.95$ and
5.670741170 for $H/J_{1}=2.88$. The black solid lines represent the analytical
solution Eq. (46) while the red dashed lines represent the numerical solution
of the mean-field equation (32) and (33). The blue dotted line marks the
location of the inflection point of $\phi(x)$ for the numerical solution.

###  III.2 Comparison between analytical and numerical analyses

On the basis of an analytical solution of the stationary GL equation (22) in
the presence of a hard wall potential, we briefly review the surface
criticality that emerges near the first-order quantum phase transition [14,
15, 16]. We focus on the case of the antiferromagnetic phase at two spatial
dimensions and assume that the system is homogeneous in the $y$ direction. We
also assume that there exists a hard wall potential at $x<0$, which leads to
the boundary condition that $\phi(x=0)=0$. At a distance sufficiently far from
the boundary, the order parameter is approximately uniform as

$\displaystyle\phi(x)\simeq\phi_{0}=\sqrt{\frac{-u+\sqrt{u^{2}-4sv}}{2v}}.$
(40)

In Fig. 3, the dotted line represents the contour of $\phi_{0}=0.25$ so that
the gray-shaded region means the region of the antiferromagnetic phase in
which $\phi_{0}<0.25$. Since the condition $\phi(x)\ll 1$ is required for the
GL approximation to be justified, the gray-shaded region roughly estimates the
validity region of the GL theory. For the transition to be of first order, the
condition that $u<0$ has to be satisfied. In this case, the antiferromagnetic
state is the ground state when $s<\frac{3u^{2}}{16v}$. These conditions imply
that $\frac{u}{v\phi_{0}^{2}}=-\frac{4}{3}$ at the first-order transition
point. Under the above-mentioned boundary conditions, one can solve Eq. (22)
to obtain

$\phi(x)=\frac{\phi_{0}\sqrt{\gamma}\tanh{\frac{x}{\xi}}}{\sqrt{1+\gamma-\mathrm{tanh^{2}}\frac{x}{\xi}}},$
(41)

where $\xi=\sqrt{\frac{2C}{u\phi_{0}^{2}+2v\phi_{0}^{4}}}$ and
$\gamma=2+\frac{3u}{2v\phi_{0}^{2}}$ ($\gamma=0$ at the first-order
transition). For a fixed value of $H/J_{1}$,
$\gamma\propto\Delta\Gamma/J_{1}\equiv(\Gamma_{\rm t}-\Gamma)/J_{1}$ near the
first order transition point $\Gamma=\Gamma_{\rm t}$.

When $0<\gamma<\frac{1}{2}$, $\phi(x)$ of Eq. (41) has an inflection point at
$x=l_{\rm inf}$, where $l_{\rm inf}$ is given by

$\displaystyle\frac{l_{\rm inf}}{\xi}=2\,{\rm
arctanh}(1-2\gamma)^{\frac{1}{2}}.$ (42)

$l_{\rm inf}$ can be interpreted as a healing length of the order parameter in
the sense that the order parameter that is zero at the boundary recovers up to
$\sim\frac{\phi_{0}}{2}$ at $x=l_{\rm inf}$. Near the first-order transition
point, where $\gamma\ll 1$, $\frac{l_{\rm inf}}{\xi}\simeq-\ln\gamma$. This
logarithmic divergence of the length scale $l_{\rm inf}$ is a universal
characteristic of the surface criticality.

Another important critical behavior is that of the order parameter amplitude
at a certain location near the boundary. For concreteness, at $x=\xi$, the
order parameter behaves as

$\displaystyle\frac{\phi(x=\xi)}{\phi_{0}}=\frac{\sqrt{\gamma}\tanh{1}}{\sqrt{1+\gamma-\tanh^{2}{1}}}\simeq\frac{\sqrt{\gamma}\tanh{1}}{\sqrt{1-\tanh^{2}{1}}}.$
(43)

More generically, the order parameter is proportional to
$\gamma^{\frac{1}{2}}$ when $x\sim\xi$.

Figure 5: Spatial profile of the approximate second derivative of the order
parameter $(\phi(x_{j}+a)-2\phi(x_{j})+\phi(x_{j}-a))/a^{2}$ calculated by
numerically solving the mean-field equations (32) and (33), where we set
$a=1$. We take the same values of the parameters as in Fig. 4. The blue dotted
line marks the location of the inflection point of $\phi(x)$ for the numerical
solution.

In actual calculations shown below, we consider the situation in which there
is another hard wall at $x>(L+1)a$ in addition to the one at $x<0$. When
$l_{\rm inf}$ is sufficiently small compared to the system size, the spatial
dependence of the order parameter can be well approximated as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\phi(x)\\!=\\!\left\\{\begin{array}[]{cc}\frac{\phi_{0}\sqrt{\gamma}\tanh{\frac{x}{\xi}}}{\sqrt{1+\gamma-\mathrm{tanh^{2}}\frac{x}{\xi}}},&{\rm
if}\,\,0<\frac{x}{a}<\frac{L+1}{2}\\\
\frac{\phi_{0}\sqrt{\gamma}\tanh{\frac{(L+1)a-x}{\xi}}}{\sqrt{1+\gamma-\mathrm{tanh^{2}}\frac{(L+1)a-x}{\xi}}},&{\rm
if}\,\,\frac{L+1}{2}<\frac{x}{a}<L+1\end{array}\right.\\!\\!\\!.$ (46)

In Fig. 4, this solution for the two different values of
$(\Delta\Gamma/J_{1},H/J_{1})$ is depicted by the black solid lines, where
$\Delta\Gamma=\Gamma_{\rm t}-\Gamma$ and $\Gamma_{\rm t}$ corresponds to the
first-order transition point for a given value of $H/J_{1}$.

With the analytical insights obtained from the solution of the GL equation in
mind, we numerically analyze the surface criticality by solving the mean-field
equations (32) and (33). Such numerical calculations are advantageous over the
analytical ones in the sense that the former is applicable even when the order
parameter at the transition is relatively large so that validity of the GL
equation can not be fully justified. In our numerical calculations, we assume
a two-sublattice structure in the $y$ direction and the open boundaries at
$x_{j}=a$ and $La$ in the $x$ direction.

Figure 6: Distance $l_{\rm inf}$ of the inflection point from the boundary as
a function of $\Delta\Gamma/J_{1}$ for $H/J_{1}=2.95$ and 2.88. The black
solid lines represent the analytical results of Eq. (42) [14, 15, 16] while
the dashed, dotted, and dash-dotted lines represent numerical results for
different values of $L$. The longitudinal (horizontal) axis is in the linear
(logarithmic) scale.

In Fig. 4, we plot the spatial profile of the order parameters $\phi(x_{j})$
computed from the mean-field equations for
$(\Delta\Gamma/J_{1},H/J_{1})=(2.7\times 10^{-8},2.95)$ (a) and $(4.1\times
10^{-5},2.88)$ (b). In Fig. 4(a), in which the parameter is relatively far
from the QTCP, the deviation of the GL solution from the numerical solution is
relatively large while qualitative tendency of the spatial profile is well
captured. This happens because the order parameter is too large, i.e.,
$\phi_{0}>0.25$, to justify the ignorance of the higher order terms in the GL
expansion. By contrast, in Fig. 4(b), since the system is closer to the QTCP,
$\phi_{0}\simeq 0.12$ such that the analytical result agrees quantitatively
with the numerical one.

We next compare the distance to the inflection point from the boundary $l_{\rm
inf}$. For this purpose, we calculate the second derivative of the order
parameter by using the following approximate formula,

$\displaystyle\left.\frac{d^{2}\phi}{dx^{2}}\right|_{x=x_{j}}\simeq\frac{\phi(x_{j}+a)-2\phi(x_{j})+\phi(x_{j}-a)}{a^{2}}.$
(47)

In Fig. 5, we show the spatial profile of the approximate second derivative of
the order parameter for the same values of $(\Delta\Gamma/J_{1},H/J_{1})$ as
in Fig. 4. $l_{\rm inf}$ can be clearly determined as $x$ at which the second
order derivative is equal to zero.

In Fig. 6, we show $l_{\rm inf}$ versus $\Delta\Gamma/J_{1}$ for different
values of $L$. There we also plot the analytical expression of Eq. (42) for
comparison. While in Fig. 6(a) the agreement between the analytical and
numerical results is only qualitative, in Fig. 6(b) we see that the numerical
results quantitatively agree with the analytical ones. Again, this happens
because the parameters of the latter is rather close to the QTCP and
$\phi_{0}$ is sufficiently small. Nevertheless, even in the former case, the
numerically obtained $l_{\rm inf}$ exhibits clear logarithmic divergence. This
manifests the universality of this surface critical behavior near the first-
order transition point.

##  IV Summary and outlook

We have proposed how to realize the mixed-field Ising model with sign-inverted
next-nearest-neighbor interaction in the system of Rydberg atoms in an optical
tweezer array by using a Rydberg dressed Rydberg state. We have analyzed
surface criticality, which is one of the interesting phenomena that can appear
in the proposed model. We have derived a Ginzburg-Landau (GL) equation which
describes the behavior of the antiferromagnetic order parameter near the
quantum phase transition when the order parameter amplitude is small. All the
GL coefficients have been explicitly related to the parameters in the original
Ising model. We have also present numerical calculations of the order
parameter near the first order transitions on the basis of a mean-field theory
in order to directly compare them with the analytical insights based on the GL
equation. We have found a quantitative agreement between the numerical and
analytical results near the quantum tricritical point (QTCP). In a parameter
region away from QTCP, while the agreement is only qualitative, the healing
length of the order parameter exhibits the logarithmic divergence that is a
characteristic of the surface criticality, indicating the universality of this
behavior.

###### Acknowledgements.

The authors thank D. Kagamihara, M. Kunimi, and D. Yamamoto for useful
comments. This work was supported by JSPS KAKENHI (Grants No. JP21H01014, No.
JP21K13855, and No. JP24H00973), and MEXT Q-LEAP (Grant No. JPMXS0118069021),
JST FOREST (Grant No. JPMJFR202T).

## References

  * [1] A. Browaeys and T. Lahaye, Nat. Phys. 16, 132 (2020).
  * [2] M. Morgado and S. Whitlock, AVS Quantum Sci. 3, 023501 (2021).
  * [3] X. Wu, X. Liang, Y. Tian, F. Yang, C. Chen, Y.-C. Liu, M. K. Tey, and L. You, Chin. Phys. B 30, 020305 (2021).
  * [4] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Nature 551, 579 (2017).
  * [5] S. Ebadi, T. T. Wang, H. Levine, A. Keesling, G. Semeghini, A. Omran, D. Bluvstein, R. Samajdar, H. Pichler, W. W. Ho, S. Choi, S. Sachdev, M. Greiner, V. Vuletic, and M. D. Lukin, Nature 595, 227 (2021).
  * [6] P. Scholl, M. Schuler, H. J. Williams, A. A. Eberharter, D. Barredo, K.-N. Schymik, V. Lienhard, L.-P. Henry, T. C. Lang, T. Lahaye, A. M. Läuchli, and A. Browaeys, Nature 595, 233 (2021).
  * [7] V. Lienhard, S. de Léséleuc, D. Barredo, T. Lahaye, A. Browaeys, M. Schuler, L.-P. Henry, and A. M. Läuchli, Phys. Rev. X 8, 021070 (2018).
  * [8] M. Rader and M. Läuchli, arXiv:1908.02068 [cond-mat.quant-gas].
  * [9] R. Kaneko, Y. Douda, S. Goto, and I. Danshita, J. Phys. Soc. Jpn. 90, 073001 (2021).
  * [10] N. Henkel, R. Nath, and T. Pohl, Phys. Rev. Lett. 104, 195302 (2010).
  * [11] J. Zeiher, R. van Bijnen, P. Schauß, S. Hild, J.-y. Choi, T. Pohl, I. Bloch, and C. Gross, Nat. Phys. 12, 1095 (2016).
  * [12] M. Beccaria, M. Campostrini, and A. Feo, Phys. Rev. B 76, 094410 (2007).
  * [13] Y. Kato and T. Misawa, Phys. Rev. B 92, 174419 (2015).
  * [14] R. Lipowsky, Phys. Rev. Lett. 49, 1575 (1982).
  * [15] R. Lipowsky, Ferroelectrics 73, 69 (1987).
  * [16] I. Danshita, D. Yamamoto, and Y. Kato, Phys. Rev. A 91, 013630 (2015).
  * [17] R. Samajdar, W. W. Ho, H. Pichler, M. D. Lukin, and S. Sachdev, Phys. Rev. Lett. 124, 103601 (2020).
  * [18] M. Kalinowski, R. Samajdar, R. G. Melko, M. D. Lukin, S. Sachdev, and S. Choi, Phys. Rev. B 105, 174417 (2022).
  * [19] I. Danshita and D. Yamamoto, Phys. Rev. A 82, 013645 (2010).