11institutetext: Baidu Research
22institutetext: Mind Vogue Lab
33institutetext: Institute for Infocomm Research, A*STAR

# MedSegDiff: Medical Image Segmentation with Diffusion Probabilistic Model

Junde Wu 11 Rao Fu 22 Huihui Fang 11 Yu Zhang 22 Yehui Yang 11 Haoyi Xiong 11
Huiying Liu 33 Yanwu Xu  (✉) 11

###### Abstract

Diffusion probabilistic model (DPM) recently becomes one of the hottest topic
in computer vision. Its image generation application such as Imagen, Latent
Diffusion Models and Stable Diffusion have shown impressive generation
capabilities, which aroused extensive discussion in the community. Many recent
studies also found it is useful in many other vision tasks, like image
deblurring, super-resolution and anomaly detection. Inspired by the success of
DPM, we propose the first DPM based model toward general medical image
segmentation tasks, which we named MedSegDiff. In order to enhance the step-
wise regional attention in DPM for the medical image segmentation, we propose
dynamic conditional encoding, which establishes the state-adaptive conditions
for each sampling step. We further propose Feature Frequency Parser (FF-
Parser), to eliminate the negative effect of high-frequency noise component in
this process. We verify MedSegDiff on three medical segmentation tasks with
different image modalities, which are optic cup segmentation over fundus
images, brain tumor segmentation over MRI images and thyroid nodule
segmentation over ultrasound images. The experimental results show that
MedSegDiff outperforms state-of-the-art (SOTA) methods with considerable
performance gap, indicating the generalization and effectiveness of the
proposed model. Our code is released at https://github.com/WuJunde/MedSegDiff.

###### Keywords:

d

iffusion probabilistic model, medical image segmentation, brain tumor, optic
cup, thyroid nodule

##  1 Introduction

Medical image segmentation is the process of partitioning a medical image into
meaningful regions. Segmentation is a fundamental step in many medical image
analysis applications such as diagnosis, surgical planning, and image-guided
surgery. This is important because it allows doctors and other medical
professionals to better understand what they’re looking at. It also makes it
easier to compare images and track changes over time. In recent years, there
has been a growing interest in automatic medical image segmentation methods.
These methods have the potential to reduce the time and effort required for
manual segmentation, and to improve the consistency and accuracy of results.
With the development of the deep learning techniques, more and more studies
successfully applied the neural network (NN) based models to the medical image
segmentation tasks, from the popular convolution neural networks (CNN) [11] to
the recent vision transformers (ViT) [3, 22, 12, 28].

Very recently, diffusion probabilistic model (DPM)[9] gained popularity as a
powerful class of generative models[27], that is able to generate images with
high diversity and synthesis quality. Recent large diffusion models, such as
DALL-E2[17], Imagen[19], and Stable Diffusion[18] have shown incredible
generation capability. Diffusion models are originally applied in fields in
which there is no absolute ground-truth. However, recent studies show that it
is also effective for the problems in which the ground-truth is unique, like
super-resolution[20] and deblurring[24].

Inspired by the recent success of DPM, we design a unique DPM-based
segmentation model for the medical image segmentation tasks. To our knowledge,
we are the first to propose the DPM-based model under the background of
general medical image segmentation with different image modalities. We note
that in tasks of medical image segmentation, the lesions/organs are often
ambiguous and hard to discriminate from the background. In that case, an
adaptive calibration process is the key to obtain a delicate result. Following
this mindset, we propose dynamic conditional encoding over vanilla DPM to
design the proposed model, named MedSegDiff. Note that in the iterative
sampling process, MedSegDiff conditions each of the step with image prior, in
order to learn the segmentation map from it. Toward the adaptive regional
attention, we integrate the segmentation map of current step into the image
prior encoding at each step. The specific implementation is to fuse the
current-step segmentation mask with the image prior on the feature level with
a multi-scale manner. In this way, the corrupted current-step mask helps to
dynamically enhance the condition features, thus improves the reconstruction
accuracy. In order to eliminate the high-frequency noises in the corrupted
given mask in this process, we further propose the feature frequency parser
(FF-Parser) to filter the features in the Fourier space. FF-Parsers are
adopted on each skip connection path for the multi-scale integration. We
verify MedSegDiff on three different medical segmentation tasks, the optic-cup
segmentation, the brain tumor segmentation, and the thyroid nodule
segmentation. The images of these tasks have different modalities, which are
the fundus images, brain CT images, the ultrasound images respectively.
MedSegDiff outperforms the previous SOTA on all three tasks with different
modalities, which shows the generalization and effectiveness the proposed
method. In brief, the contributions of the paper are:

  * •

The fist to propose DPM-based model toward general medical image segmentation.

  * •

Dynamic conditional encoding strategy is proposed for step-wise attention.

  * •

FF-Parser is proposed to eliminate the negative effects of high-frequency
components.

  * •

SOTA performance on three different medical segmentation tasks with different
image modalities.

##  2 Method

We design our model based on diffusion model mentioned in[9]. Diffusion models
are generative models composed of two stages, a forward diffusion stage and a
reverse diffusion stage. In the forward process, the segmentation label
$x_{0}$ is gradually added Gaussian noise through a series of steps $T$. In
the reverse process, a neural network is trained to recover the original data
by reversing the noising process, which can be represented as:

$p_{\theta}(x_{0:T-1}|x_{T})=\Pi^{T}_{t=1}p_{\theta}(x_{t-1}|x_{t}),$ (1)

where $\theta$ is reverse process parameters. Starting from a Gaussian noise,
$p_{\theta}(x_{T})=\mathcal{N}(x_{T};0,I_{n\times n})$, where $I$ is the raw
image, the reverse process transforms the latent variable distribution
$p_{\theta}(x_{T})$ to the data distribution $p_{\theta}(x_{0})$. To be
symmetrical to the forward process, the reverse process recovers the noise
image step by step to obtain the final clear segmentation.

Following the standard implementation of DPM, we adopt a UNet as the network
for the learning. An illustration is shown in Figure 1. In order to achieve
the segmentation, we condition the step estimation function $\epsilon$ by raw
image prior, which can be represented as:

$\epsilon_{\theta}(x_{t},I,t)=D((E_{t}^{I}+E_{t}^{x},t),t),$ (2)

where $E_{t}^{I}$ is the conditional feature embedding, in our case, the raw
image embedding, $E_{t}^{x}$ is the segmentation map feature embedding of the
current step. The two components are added and sent to a UNet decoder $D$ for
the reconstruction. The step index $t$ is integrated with the added embedding
and decoder features. In each of these, it is embedded using a shared learned
look-up table, following [9].

Figure 1: An illustration of MedSegDiff. For the clarity, the time step
encoding is omitted in the figure.

###  2.1 Dynamic Conditional Encoding

In most conditional DPM, the conditional prior will be a unique given
information. However, medical image segmentation is notorious for its
ambiguous objects. The lesions or tissues are often hard to discriminate from
its background. The low-contrast image modalities, such as MRI or ultrasound
images, make it even worse. Given only a static image $I$ as the condition for
each step will be hard to learn. To address this problem, we propose a dynamic
conditional encoding for each step. We note that on the one hand, the raw
image contains the accurate segmentation target information but hard to
discriminate from the background, on the other hand, the current-step
segmentation map contains the enhanced target regions but not accurate. This
motivated us to integrate the current-step segmentation information $x_{t}$
into the conditional raw image encoding for the mutual complement. To be
specific, we implement the integration on the feature level. In the raw image
encoder, we enhance its intermediate feature with the current-step encoding
features. Each scale of the conditional feature map $m_{I}^{k}$ is fused with
the $x_{t}$ encoding features $m_{x}^{k}$ with the same shape, $k$ is the
index of layer. The fusion is implemented by an attentive-like mechanism
$\mathcal{A}$. In particular, two feature maps are first applied layer
normalization and multiply together to get an affinity map. Then we multiply
the affinity map with the condition encoding features to enhance the attentive
region, which is:

$\mathcal{A}(m_{I}^{k},m_{x}^{k})=(LN(m_{I}^{k})\otimes LN(m_{x}^{k}))\otimes
m_{I}^{k},$ (3)

where $\otimes$ implies element-wise multiplication, $LN$ denotes layer
normalization. The operation is applied on the middle two stages, where each
is the convolutional stage implemented following ResNet34. Such a strategy
helps MedSegDiff dynamically localize and calibrate the segmentation. Although
effective the strategy it is, another specific problem is that integrating
$x_{t}$ embedding will induce extra high-frequency noise. To address this
problem, we propose FF-Parser to constrain the high-frequency components in
the features.

###  2.2 FF-Parser

We connect FF-parser in the path ways of the feature integration. The function
of it is to constrain the noise-related components in the $x_{t}$ features.
Our main idea is to learn a parameterized attentive (weight) map applying on
the Fourier-space features. Given a decoder feature map
$m\in\mathbb{R}^{H\times W\times C}$, we first perform 2D FFT(fast fourier
transform) along the spatial dimensions, which we can represented as:

$M=\mathcal{F}[m]\in\mathbb{C}^{H\times W\times C},$ (4)

where $\mathcal{F}[\cdot]$ denotes the 2D FFT. We then modulate the spectrum
of $m$ by multiplying a parameterized attentive map $A\in\mathbb{C}^{H\times
W\times C}$ to $M$:

$M^{\prime}=A\otimes M,$ (5)

where $\otimes$ denotes the element-wise product. Finally, we reverse
$M^{\prime}$ back to the spatial domain by adopting inverse FFT:

$m^{\prime}=\mathcal{F}^{-1}[M^{\prime}].$ (6)

FF-Parser can be regarded as a learnable version of frequency filters which
are wildly applied in the digital image processing [16]. Different from the
spacial attention, it globally adjusts the components of the specific
frequencies. Thus it can be learn to constrain the high-frequency component
for the adaptive integration.

Figure 2: An illustration of FF-Parser. FFT denotes Fast Fourier Transform.

###  2.3 Training and Architecture

MedSegDiff is trained following the standard process of DPM [9]. Specifically,
the loss can be represented as:

$\mathcal{L}=E_{x_{0},\epsilon,t}[||\epsilon-\epsilon_{\theta}(\sqrt{\hat{a}_{t}}x_{0}+\sqrt{1-\hat{a}_{t}}\epsilon,I_{i},t)||^{2}].$
(7)

In each of the iteration, a random couple of raw image $I_{i}$ and
segmentation label $S_{i}$ will be sampled for the training. The iteration
number is sampled from a uniform distribution and $\epsilon$ from a Gaussian
distribution.

The main architecture of MedSegDiff is a modified ResUNet[26], which we
implement it with a ResNet encoder following a UNet decoder. The detailed
network setting is following [14]. $I$ and $x_{t}$ are encoded with two
individual encoders. The encoder is consisted of three convolution stages.
Each stage contains several residual blocks. The number of residual blocks in
each stage is following that of ResNet34. Each residual block is composed of
two convolutional blocks, each one consists of group-norm and SiLU[5] active
layer and a convolutional layer. The residual block receives the time
embedding through a linear layer, SiLU activation, and another linear layer.
The result is then added to the output of the first convolutional block. The
obtained $E^{I}$ and $E^{x_{t}}$ are added together and sent to the last
encoding stage. A standard convolutional decoder is connected to predict the
final result.

##  3 Experiments

###  3.1 Dataset

We conduct the experiments on three different medical tasks with different
image modalities, which are optic-cup segmentation from fundus images, brain
tumor segmentation from MRI images, and thyroid nodule segmentation from
ultrasound images. The experiments of glaucoma, thyroid cancer and melanoma
diagnosis are conducted on REFUGE-2 dataset [6], BraTs-2021 dataset [2] and
DDTI dataset [15], which contain 1200, 2000, 8046 samples, respectively. The
datasets are publicly available with both segmentation and diagnosis labels.
Train/validation/test sets are split following the default settings of the
dataset.

###  3.2 Implementation Details

We experiment with huge, large, basic, and small variants of our model,
MedSegDiff++, MedSegDiff-L, MedSegDiff-B, and MedSegDiff-S, respectively.
In MedSegDiff-S, MedSegDiff-B MedSegDiff-L, MedSegDiff++, we use UNet with
$4$x, $5$x, $6$x, $6$x downsamples respectively. In the experiments, we employ
100 diffusion steps for the inference, which is much smaller than most of the
previous studies[9, 14]. All the experiments are implemented with the PyTorch
platform and trained/tested on 4 Tesla P40 GPU with 24GB of memory except
MedSegDiff++ and MedSegDiff-L. All images are uniformly resized to the
dimension of 256$\times$256 pixels. The networks are trained in an end-to-end
manner using AdamW[13] optimizer. MedSegDiff-B and MedSegDiff-S are trained
with 32 batch size, MedSegDiff-L and MedSegDiff++ are trained with 64 batch
size. The learning rate is initially set to 1 $\times 10^{-4}$. All models are
set 25 times of ensemble in the inference. We use STAPLE[23] algorithm to fuse
the different samples. The diffusion based competitor EnsemDiff[25] is
reproduced with the same setting for the fair comparison.

###  3.3 Main Results

We compare with SOTA segmentation methods proposed for the three specific
tasks and general medical image segmentation methods. The main results are
shown in Table 1. In the table, ResUnet[26] and BEAL[21] are proposed for
optic disc/cup segmentation, TransBTS[22] and EnsemDiff[25] are proposed for
the brain tumor segmentation, MTSeg[7] and UltraUNet[4] are proposed for the
Thyroid Nodule segmentation, CENet[8], MRNet[11], SegNet[1], nnUNet[10] and
TransUNet[3] are proposed for the general medical image segmentation. We
evaluate the segmentation performance by Dice score and IoU.

In Table 1, we compare with the methods implemented with various network
architectures, including CNN (ResUNet, BEAL, nnUNet, SegNet), vision
transformer (TransBTS, TransUNet) and DPM (EnsemDiff). We can see the advanced
network architectures commonly gain better results. For example, in optic-cup
segmentation, ViT-based general segmentation method: TransUNet is even better
than the CNN-based task toward method: BEAL. On brain tumor segmentation,
recently proposed DPM-based segmentation method EnsemDiff outperforms all
those previous ViT-based competitors, i.e., TransBTS and TransUNet. MedSegDiff
not only adopts the recent successful DPM, but also designs an appropriate
strategy over it specifically towards the general medical image segmentation
task. We can see MedSegDiff outperforms all the other methods on three
different tasks, which shows the generalization toward different medical
segmentation tasks and different image modalities. Comparing against DPM-based
model proposed specifically for the brain tumor segmentation, i.e., EnsemDiff,
it improves 2.3% on Dice and 2.4% on IoU, which indicates the effectiveness of
our unique techniques, i.e, dynamic conditioning and FF-Parser.

Figure 3: The visual comparison of Top-4 general medical image segmentation
methods in Table 1. From top to down are brain-tumor segmentation, optic-cup
segmentation and thyroid nodule segmentation, respectively.

Figure 3 shows several typical examples generated by our MedSegDiff and other
SOTA methods. It can be seen the target lesions/tissues are all ambiguous on
the images so that they are hard to be recognized by human eyes. Comparing
with these computer-aided methods, it is obvious that the segmentation maps
generated by the proposed method are more accurate than the other methods,
especially for the ambiguous regions. To be benefited from DPM together with
the proposed dynamic conditioning and FF-Parser, it can better localize and
calibrate the segmentation on the low-contrast or ambiguous images.

Table 1: The comparison of MedSegDiff with SOTA segmentation methods. Best
results are denoted as bold. The grey background denotes the methods are
proposed for that/these particular tasks.

| Optic-Cup | Brain-Turmor | Thyroid Nodule
---|---|---|---
| Dice | IoU | Dice | IoU | Dice | IoU
ResUnet | 80.1 | 72.3 | - | - | - | -
BEAL | 83.5 | 74.1 | - | - | - | -
TransBTS | - | - | 87.6 | 78.3 | - | -
EnsemDiff | - | - | 88.7 | 80.9 | - | -
MTSeg | - | - | - | - | 82.3 | 75.2
UltraUNet | - | - | - | - | 84.5 | 76.2
CENet | 78.6 | 69,4 | 76.2 | 68.9 | 78.9 | 71.2
MRNet | 84.2 | 75.1 | 83.4 | 75.6 | 80.4 | 73.4
SegNet | 80.4 | 70.7 | 80.2 | 72.9 | 81.7 | 74.5
nnUNet | 84.9 | 75.1 | 88.2 | 80.4 | 84.2 | 76.2
TransUNet | 85.6 | 75.9 | 86.6 | 79.0 | 83.5 | 75.1
MedSegDiff-S | 81.2 | 71.7 | 82.3 | 73.6 | 80.8 | 73.7
MedSegDiff-B | 85.9 | 76.2 | 88.9 | 81.2 | 84.8 | 76.4
MedSegDiff-L | 86.9 | 78.5 | 89.9 | 82.3 | 86.1 | 79.6
MedSegDiff++ | 87.5 | 79.1 | 90.5 | 82.8 | 86.6 | 80.2

###  3.4 Ablation Study

We do comprehensive ablation study to verify the effectiveness of the proposed
dynamic conditioning and FF-Parser. The results are shown in Table 2, where
Dy-Cond denotes dynamic conditioning. We evaluate the performance by Dice
score(%) on all three tasks. From the table, we can see Dy-Cond gains
considerable improvements over vanilla DPM. On the case which the region
localization is important, i.e., optic-cup segmentation, it improves 2.1%. On
the cases which the images are low-contrast, like brain tumor and thyroid
nodule segmentation, it improves 1.6% and 1.8% respectively. It shows Dy-Cond
is a generally effective strategy on DPM for both of the cases. FF-Parser
which established over Dy-Cond mitigates the high-frequency noises thus
further optimize the segmentation results. It helps MedSegDiff further improve
near 1% performance and achieve the best on all three tasks.

Table 2: An ablation study on dynamic condition encoding and FF-Parser. Dice score(%) is used as the metric. Dy-Cond | FF-Parser | OpticCup | BrainTumor | ThyroidNodule
---|---|---|---|---
|  | 84.6 | 88.2 | 84.1
✓ |  | 86.7 | 89.8 | 85.9
✓ | ✓ | 87.5 | 90.5 | 86.6

##  4 Conclusion

In this paper, we provided a scheme for DPM-based general medical image
segmentation, named MedSegDiff. We propose two novel techniques to promise the
performance of it, i.e., the dynamic conditional encoding and FF-Parser. The
comparison experiments are conducted on three medical image segmentation tasks
with different image modalities, which shows our model outperforms previous
SOTA. As the first DPM application in general medical image segmentation, we
believe MedSegDiff will serve as an essential benchmark for future research.

## References

  * [1] Badrinarayanan, V., Kendall, A., Cipolla, R.: Segnet: A deep convolutional encoder-decoder architecture for image segmentation. IEEE transactions on pattern analysis and machine intelligence 39(12), 2481–2495 (2017)
  * [2] Baid, U., Ghodasara, S., Mohan, S., Bilello, M., Calabrese, E., Colak, E., Farahani, K., Kalpathy-Cramer, J., Kitamura, F.C., Pati, S., et al.: The rsna-asnr-miccai brats 2021 benchmark on brain tumor segmentation and radiogenomic classification. arXiv preprint arXiv:2107.02314 (2021)
  * [3] Chen, J., Lu, Y., Yu, Q., Luo, X., Adeli, E., Wang, Y., Lu, L., Yuille, A.L., Zhou, Y.: Transunet: Transformers make strong encoders for medical image segmentation. arXiv preprint arXiv:2102.04306 (2021)
  * [4] Chu, C., Zheng, J., Zhou, Y.: Ultrasonic thyroid nodule detection method based on u-net network. Computer Methods and Programs in Biomedicine 199, 105906 (2021)
  * [5] Elfwing, S., Uchibe, E., Doya, K.: Sigmoid-weighted linear units for neural network function approximation in reinforcement learning. Neural Networks 107, 3–11 (2018)
  * [6] Fang, H., Li, F., Fu, H., Sun, X., Cao, X., Son, J., Yu, S., Zhang, M., Yuan, C., Bian, C., et al.: Refuge2 challenge: Treasure for multi-domain learning in glaucoma assessment. arXiv preprint arXiv:2202.08994 (2022)
  * [7] Gong, H., Chen, G., Wang, R., Xie, X., Mao, M., Yu, Y., Chen, F., Li, G.: Multi-task learning for thyroid nodule segmentation with thyroid region prior. In: 2021 IEEE 18th International Symposium on Biomedical Imaging (ISBI). pp. 257–261. IEEE (2021)
  * [8] Gu, Z., Cheng, J., Fu, H., Zhou, K., Hao, H., Zhao, Y., Zhang, T., Gao, S., Liu, J.: Ce-net: Context encoder network for 2d medical image segmentation. IEEE transactions on medical imaging 38(10), 2281–2292 (2019)
  * [9] Ho, J., Jain, A., Abbeel, P.: Denoising diffusion probabilistic models. Advances in Neural Information Processing Systems 33, 6840–6851 (2020)
  * [10] Isensee, F., Jaeger, P.F., Kohl, S.A., Petersen, J., Maier-Hein, K.H.: nnu-net: a self-configuring method for deep learning-based biomedical image segmentation. Nature methods 18(2), 203–211 (2021)
  * [11] Ji, W., Yu, S., Wu, J., Ma, K., Bian, C., Bi, Q., Li, J., Liu, H., Cheng, L., Zheng, Y.: Learning calibrated medical image segmentation via multi-rater agreement modeling. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 12341–12351 (2021)
  * [12] Liu, C., Zhao, R., Shi, Z.: Remote-sensing image captioning based on multilayer aggregated transformer. IEEE Geoscience and Remote Sensing Letters 19, 1–5 (2022)
  * [13] Loshchilov, I., Hutter, F.: Decoupled weight decay regularization. arXiv preprint arXiv:1711.05101 (2017)
  * [14] Nichol, A.Q., Dhariwal, P.: Improved denoising diffusion probabilistic models. In: International Conference on Machine Learning. pp. 8162–8171. PMLR (2021)
  * [15] Pedraza, L., Vargas, C., Narváez, F., Durán, O., Muñoz, E., Romero, E.: An open access thyroid ultrasound image database. In: 10th International symposium on medical information processing and analysis. vol. 9287, pp. 188–193. SPIE (2015)
  * [16] Pitas, I.: Digital image processing algorithms and applications. John Wiley & Sons (2000)
  * [17] Ramesh, A., Dhariwal, P., Nichol, A., Chu, C., Chen, M.: Hierarchical text-conditional image generation with clip latents. arXiv preprint arXiv:2204.06125 (2022)
  * [18] Rombach, R., Blattmann, A., Lorenz, D., Esser, P., Ommer, B.: High-resolution image synthesis with latent diffusion models. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 10684–10695 (2022)
  * [19] Saharia, C., Chan, W., Saxena, S., Li, L., Whang, J., Denton, E., Ghasemipour, S.K.S., Ayan, B.K., Mahdavi, S.S., Lopes, R.G., et al.: Photorealistic text-to-image diffusion models with deep language understanding. arXiv preprint arXiv:2205.11487 (2022)
  * [20] Saharia, C., Ho, J., Chan, W., Salimans, T., Fleet, D.J., Norouzi, M.: Image super-resolution via iterative refinement. IEEE Transactions on Pattern Analysis and Machine Intelligence (2022)
  * [21] Wang, S., Yu, L., Li, K., Yang, X., Fu, C.W., Heng, P.A.: Boundary and entropy-driven adversarial learning for fundus image segmentation. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 102–110. Springer (2019)
  * [22] Wang, W., Chen, C., Ding, M., Yu, H., Zha, S., Li, J.: Transbts: Multimodal brain tumor segmentation using transformer. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 109–119. Springer (2021)
  * [23] Warfield, S.K., Zou, K.H., Wells, W.M.: Simultaneous truth and performance level estimation (staple): an algorithm for the validation of image segmentation. IEEE transactions on medical imaging 23(7), 903–921 (2004)
  * [24] Whang, J., Delbracio, M., Talebi, H., Saharia, C., Dimakis, A.G., Milanfar, P.: Deblurring via stochastic refinement. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 16293–16303 (2022)
  * [25] Wolleb, J., Sandkühler, R., Bieder, F., Valmaggia, P., Cattin, P.C.: Diffusion models for implicit image segmentation ensembles. arXiv preprint arXiv:2112.03145 (2021)
  * [26] Yu, S., Xiao, D., Frost, S., Kanagasingam, Y.: Robust optic disc and cup segmentation with deep learning for glaucoma detection. Computerized Medical Imaging and Graphics 74, 61–71 (2019)
  * [27] Zhao, R., Shi, Z.: Text-to-remote-sensing-image generation with structured generative adversarial networks. IEEE Geoscience and Remote Sensing Letters 19, 1–5 (2021)
  * [28] Zhao, R., Shi, Z., Zou, Z.: High-resolution remote sensing image captioning based on structured attention. IEEE Transactions on Geoscience and Remote Sensing 60, 1–14 (2021)


# Meissner effect in Fock space

J. Mumford  Department of Physics and Astronomy, McMaster University, 1280
Main St. W., Hamilton, ON, L8S 4M1, Canada

###### Abstract

By periodically driving a single bosonic Josephson junction (BJJ) with an
impurity, a synthetic gauge field is generated in the Fock space of the
system. At a critical synthetic gauge flux the ground state undergoes a
quantum phase transition which is analogous to the Meissner-Abrikosov-vortex
transition found in type-II superconductors with an applied magnetic field. A
second quantum phase transition involving attractive interactions between the
bosons of the BJJ is shown to enhance the sensitivity of the system to the
Meissner-Abrikosov-vortex transition.

##  I Introduction

Ultracold atomic gases have become the primary tool in simulating quantum
phenomena [1] due to their precise tunability and because they are among the
simplest quantum many-body systems. An exciting connection between ultracold
atoms and condensed matter arises with the prospect of trapping the atoms in
optical lattices because the spatial periodicity of the trapping potentials
can be adjusted to produce a variety of crystal structures found in solids.
This has lead to the simulation of the motion of charged particles in
materials, a prime example being the observation of the transition from a
superfluid to a Mott insulator [2]. Simulating the motion of charged particles
in magnetic fields has also been of interest due to the prospect of creating
magnetic field strengths unattainable in conventional condensed matter
experiments. Synthetic gauge fields can be applied to the gas by periodically
driving the external fields of the optical lattice and atoms [3, 4, 5, 6, 7,
8, 9]. The net result is a situation in which a particle can accumulate a non-
zero phase in a closed path along the lattice mimicking the effect of a
transverse magnetic flux. This allows for the control of the topology of
single-particle energy bands leading to quantum Hall physics and the
realization of chiral edge states [10]. Celebrated topological models have
been simulated in this fashion like the Haldane [11, 12, 13] and Harper-
Hofstadter [14, 15, 16] models.

The Meissner effect is another hallmark phenomenon in the interaction of
magnetic fields and materials [17]. In type-II superconductors, screening
currents expel an external magnetic field by creating an opposing one. When
the magnetic field is larger than a critical value, $H>H_{c1}$, a quantum
phase transition (QPT) occurs and the Meissner phase breaks down resulting in
only partial screening of the field. Vortices form at locations where the
field penetrates the superconductor and is called the Abrikosov-vortex phase.
Theoretical interest in simulating this QPT has been primarily directed toward
ladder systems as they are some of the simplest lattices which still allow for
orbital magnetic field effects (a 1D chain lattice is insufficient). Early
investigations involved Josephson junction arrays to form the ladder [18, 19,
20, 21] and more recently ultracold atoms in optical ladders have been
proposed [22]. Experimental verification came shortly after by loading a Bose-
Einstein condensate (BEC) into an optical ladder and inducing a synthetic
magnetic field through laser-assisted tunneling [23]. In the experiment, it
was shown that chiral probability currents along each leg of the ladder play
the role of the screening currents in the Meissner phase of a superconductor.
Past a critical flux of the synthetic field the leg currents decrease while
the rung currents increase which was shown to be analogous to the Abrikosov-
vortex phase.

In this work, we show that the Meissner-Abrikosov-vortex QPT can be realized
with a BEC trapped in a double well potential, also known as a bosonic
Josephson junction (BJJ). The BEC is coupled to an impurity which can also
tunnel between the two wells. The array of Josephson junctions used in
previous proposals can be reduced to a single BJJ by constructing the ladder
lattice not in real space, but in Fock space called a Fock-state lattice (FSL)
[24, 25, 26, 27, 28]. The legs of the ladder are constructed from many-body
states defined by the number of bosons of the BEC in each well and the rungs
are constructed from the two similarly defined states of the impurity. The
synthetic magnetic field is generated by periodically driving the interaction
strength between the BEC and the impurity as well as the tunneling between
each well which effectively ’stirs’ the FSL. In experiments, both the boson-
boson interactions of the BEC and the BEC-impurity interactions can be
controlled via Feshbach resonance [29] where the latter requires an ac-
magnetic field to achieve periodic time dependence [30]. The double well
potential is created by superimposing a 1D optical lattice onto a harmonic
dipole trap, so the tunneling can be controlled by tuning the width and height
of the optical barrier [31]. Although we frame our analysis in terms of a BJJ,
the two accessible states of the bosons and the impurity can be internal
states as in the case of spinor condensates [32, 33]. Furthermore, similar
BEC-impurity interactions to the ones we use can be found in quantum dot
systems such as a dipole-induced ion trapped between the two wells of a BJJ
[34] and in solid state Josephson junctions [35].

The FSL is an example of the use of synthetic dimensions [36, 37, 38] which,
in general, are degrees of freedom of a system which can imitate the motion of
particles in real space. In experiments, synthetic dimensions have been
utilized with spin [39, 40, 41], momentum [42, 43, 44], harmonic oscillator
[45], and rotational [46] states of atoms and molecules. When combined with
synthetic gauge fields synthetic dimensions have successfully simulated chiral
edge states [39, 40, 47] and the measurement of Chern numbers of topological
bands [48]. Usually a synthetic dimension is paired with at least one in real
space, however, both dimensions in the FSL are synthetic which is a
significant departure from the majority of experiments and proposals to date
involving simulations with ultracold atoms.

Figure 1: Energy bands, ground state probability densities and phases in the
FLS. (a)-(c) Going from left to right the red dashed curves are the energy
bands from Eq. (7) for the magnetic fluxes $\phi_{c}/2$, $\phi_{c}$ and
$3\phi_{c}/2$. The background black and white density plots show the
probability for a state with a given phase difference $\theta$ to have a given
energy and are calculated using Eq. (12). As the synthetic magnetic flux
passes through the critical value, the lower band splits into a double well
signifying a QPT where the ground state develops a non-zero doubly degenerate
relative phase of the BEC. (d)-(f) Each panel corresponds to the ground state
probability distribution (top strips), $|\langle n|\epsilon_{0}\rangle|^{2}$,
and phase (bottom strips), $\mathrm{Arg}(\langle n|\epsilon_{0}\rangle)$, of
the lower band in the energy plots shown directly above. In the Meissner phase
the ground state is a Gaussian with a slowly varying phase along each leg of
the ladder whereas in the vortex phase the wavefunction has holes in which
vortices form which can be seen as a quickly varying phase in these regions.
The other parameter values are $N=100$ and $\xi=0.5$ and $\mu=0$.

##  II Model

The system we will be investigating consists of a BEC coupled to an impurity
in a double well potential. Both of the BEC-impurity interactions and the
tunneling between the two wells are periodically driven and are described by
the Hamiltonian

$\hat{H}(t)=U\hat{S}_{z}^{2}+W\alpha(t)\hat{S}_{z}\hat{\sigma}_{z}-J\beta(t)\hat{S}_{x}-NK\gamma(t)\hat{\sigma}_{x}/2$
(1)

where $T$ is the period of the driving such that $\hat{H}(t)=\hat{H}(t+T)$.
The BEC and the impurity can occupy the ground states of each well, so every
particle has access to two states which are labeled $L$ and $R$ for the left
and right wells, respectively. In the Hamiltonian the creation and
annihilation operators of the BEC are written in the Schwinger representation
and are
$\hat{S}_{z}=\frac{1}{2}(\hat{a}^{\dagger}_{R}\hat{a}_{R}-\hat{a}^{\dagger}_{L}\hat{a}_{L})$,
$\hat{S}_{x}=\frac{1}{2}(\hat{a}^{\dagger}_{R}\hat{a}_{L}+\hat{a}^{\dagger}_{L}\hat{a}_{R})$
and
$\hat{S}_{y}=\frac{\mathrm{i}}{2}(\hat{a}^{\dagger}_{L}\hat{a}_{R}-\hat{a}^{\dagger}_{R}\hat{a}_{L})$
where the usual commutation relations are obeyed
$[\hat{a}_{i},\hat{a}^{\dagger}_{j}]=\delta_{ij}$ where $i,j=L,R$. The
impurity operators are written as Pauli spin matrices. The parameters are
defined in the following way: $U$ is the boson-boson interaction energy of the
BEC, $W$ is the BEC-impurity interaction energy, $J$ is the BEC tunneling
energy, and $K$ is the impurity tunneling energy. The factor of $N$ in front
of the impurity tunneling term is the number of bosons in the BEC and is there
to ensure this term does not become insignificant in the thermodynamic limit,
$N\to\infty$. The Greek letters represent the periodic driving and take the
form

$\displaystyle\alpha(t)$ $\displaystyle=$ $\displaystyle-\sin(\omega t)$
$\displaystyle\beta(t)$ $\displaystyle=$
$\displaystyle\sum_{p=0}\delta[t-(p+1/2)T]$ $\displaystyle\gamma(t)$
$\displaystyle=$ $\displaystyle\sum_{p=0}\delta[t-pT]$ (2)

which results in the BEC-impurity interactions being sinusoidally driven with
frequency $\omega=2\pi/T$ and the hopping terms being alternately pulsed in
increments of $T/2$. One cycle of the driving is achieved in four steps: (1)
the impurity hopping is pulsed, (2) attractive BEC-impurity interactions are
sinusoidally ramped up from, then back down to zero over an interval $T/2$,
(3) the BEC tunneling is pulsed, (4) repulsive BEC-impurity interactions are
sinusoidally ramped up from, then back down to zero over an interval $T/2$.
The evolution after one cycle is given by the Floquet operator

$\displaystyle\hat{U}_{F}(t=T)=$
$\displaystyle\mathrm{e}^{-\mathrm{i}\left[\pi
U\hat{S}_{z}^{2}+2W\hat{S}_{z}\hat{\sigma}_{z}\right]/\omega}\mathrm{e}^{\mathrm{i}J\hat{S}_{x}\tau}$
$\displaystyle\times\mathrm{e}^{-\mathrm{i}\left[\pi
U\hat{S}_{z}^{2}-2W\hat{S}_{z}\hat{\sigma}_{z}\right]/\omega}\mathrm{e}^{\mathrm{i}NK\hat{\sigma}_{x}\tau/2}$

where $\tau\ll 1$ is the pulse interval. When $T\ll 1$ the Floquet operator
can be written in terms of an effective Hamiltonian
$\hat{U}_{F}=\mathrm{e}^{-\mathrm{i}\hat{H}_{\mathrm{eff}}\tau}$ where

$\hat{H}_{\mathrm{eff}}/J=2\frac{\mu}{N}\hat{S}^{2}_{z}-\frac{1}{2}\left[\hat{S}_{+}\mathrm{e}^{-\mathrm{i}\phi\hat{\sigma}_{z}}+\hat{S}_{-}\mathrm{e}^{\mathrm{i}\phi\hat{\sigma}_{z}}+N\xi\hat{\sigma}_{x}\right]\,.$
(4)

The new parameters are defined as $\mu=\frac{\pi UN}{J\omega\tau}$ and
$\xi=K/J$. The parameter $\phi=2W/\omega$ appears as a Peierls phase seen in
tight-binding models and so plays the role of the magnetic flux through the
system. The set of eigenvalues of $\hat{U}_{F}$,
$\\{\lambda_{i}=\mathrm{e}^{\mathrm{i}\epsilon_{i}\tau}\\}$, are defined in
terms of the set of quasienergies $\\{\epsilon_{i}\\}$. The quasienergies are
only unique in the range from $-\pi/\tau$ to $\pi/\tau$ because they are
calculated from a unitary operator. Therefore, the parameters in the model
will be kept small enough to keep the quasienergies within this range.

We assume BEC particle conservation, so the Hilbert space of the entire system
has size $2\times\left(N+1\right)$. A useful basis is the Fock basis made up
of states with different particle number differences between the two wells.
With constant $N$ these states are labeled with a single number for both the
BEC and the impurity $|n,m\rangle$ where $n=(N_{R}-N_{L})/2$ is half of the
BEC particle number difference between the two wells and $m=M_{R}-M_{L}$ is
similarly the impurity number difference between the two wells. We can extend
the connections between this system and tight-binding models by considering
the labels $n$ and $m$ as coordinates on a 2D FSL. In this case, one dimension
will be of size $N+1$ and the other will be of size $2$ forming a ladder
lattice where each leg and each rung are formed by the BEC and impurity Fock
states, respectively. As an example, the Fock state $|0,+1\rangle$, which
corresponds to an equal number of BEC particles in each well and the impurity
in the right well, will be the coordinate of the center of the right leg of
the ladder.

##  III Results

For now, we will assume the BEC is non-interacting and set $\mu=0$. We will
also assume that we are in a regime with high coherence and low number
difference fluctuations of the BEC, so that a mean-field approximation can be
applied. This amounts to replacing the BEC creation (annihilation) operators
with complex numbers
$\hat{a}_{R/L}=\sqrt{N_{R/L}}\mathrm{e}^{-\mathrm{i}\theta_{R/L}}$ where
$\theta_{R/L}$ is the phase of the right/left wells. The resulting mean-field
Hamiltonian is

$\displaystyle H_{\mathrm{MF}}/J=$
$\displaystyle-\sqrt{N^{2}/4-n^{2}}\cos\theta\cos\phi$ (5)
$\displaystyle-\sqrt{N^{2}/4-n^{2}}\sin\theta\sin\phi\hat{\sigma}_{z}-\xi
N\hat{\sigma}_{x}/2$

where $\theta=\theta_{R}-\theta_{L}$ is the phase difference between the right
and left wells (relative phase) and is the conjugate variable of $n$. The
final approximation we make is to assume that $n\ll N/2$, so that
$\sqrt{N^{2}/4-n^{2}}\approx N/2$. The justification for this approximation
comes from the fact that we are investigating a ground state QPT and the
ground state of the BEC is a fairly localized state around $n=0$. The
Hamiltonian becomes

$H_{\mathrm{MF}}/J\approx-\frac{N}{2}[\cos\theta\cos\phi+\sin\theta\sin\phi\hat{\sigma}_{z}+\xi\hat{\sigma}_{x}]$
(6)

which resembles the Bloch Hamiltonian of a translationally invariant system
with two sites per unit cell where $\theta$ plays the role of the
quasimomentum. The lattice spacing of the FSL is $a=1$, so the range of the
first Brillouin zone is $(-\pi,\pi)$ in increments of $\frac{2\pi}{N+1}$. Such
a system can be solved exactly giving the energies of the lower and upper
bands

$E_{\pm}(\theta)/J=-\frac{N}{2}\left[\cos\theta\cos\phi\mp\sqrt{\xi^{2}+\sin^{2}\theta\sin^{2}\phi}\right]$
(7)

and eigenstates

$\displaystyle|+,\theta\rangle$ $\displaystyle=$
$\displaystyle\cos(\alpha_{\theta}/2)|R,\theta\rangle-\sin(\alpha_{\theta}/2)|L,\theta\rangle$
$\displaystyle|-,\theta\rangle$ $\displaystyle=$
$\displaystyle\sin(\alpha_{\theta}/2)|R,\theta\rangle+\cos(\alpha_{\theta}/2)|L,\theta\rangle$
(8)

where

$\alpha_{\theta}=2\arctan\frac{1}{\xi}\left(\sin\theta\sin\phi+\sqrt{\xi^{2}+\sin^{2}\theta\sin^{2}\phi}\right)\,.$
(9)

The system undergoes a QPT in the ground band which can be seen by finding the
phases, $\theta_{0}$, that satisfy the condition $\partial
E_{-}(\theta)/\partial\theta=0$ giving

$\sin\theta_{0}=\begin{cases}0,&\text{if }\phi\leq\phi_{c}\\\
\pm\sqrt{\sin^{2}\phi-\xi^{2}\cot^{2}\phi},&\text{if
}\phi>\phi_{c}\end{cases}$ (10)

where the critical magnetic flux is

$\cos\phi_{c}=\frac{-\xi+\sqrt{\xi^{2}+4}}{2}\,.$ (11)

This QPT is from a zero phase difference below $\phi_{c}$ to a doubly
degenerate phase difference above $\phi_{c}$. The degeneracy comes from the
spin-orbit coupling term in Eq. (6). Its energy is minimized when the phase
difference is negative (positive) and the the impurity is in the left (right)
well, or using the FSL analogy, the left (right) legs of the ladder. Putting
it more succinctly, the degeneracy comes from $H_{\mathrm{MF}}$ being
invariant under the flipping of both the relative phase of the BEC and the
’spin’ of the impurity. To get a visual sense of the transition and to check
the validity of the approximations that have been made, in Fig. 1 (a)-(c) we
plot the energies of each band derived in Eq. (7) (dashed, red curves) as well
as a density plot of the probability for a state with a given phase to occupy
an eigenstate of the Floquet operator,

$P_{m}(\theta,\epsilon_{i})=\left|\sum_{n=-N/2}^{N/2}\mathrm{e}^{\mathrm{i}\theta
n}\langle m,n|\epsilon_{i}\rangle\right|^{2}\,,$ (12)

for different values of $\phi$. The agreement between the derived energy bands
and the probability is excellent and shows that $\theta$ is a fairly good
quantum number. One can see that the effect of keeping the square root factors
in the calculation of the probability is to cause it to fray away from the
analytic result in certain areas. The QPT can also be seen as the lower band
goes from having a single minimum at $\theta=0$ to having doubly degenerate
non-zero minima.

Figure 1 (d)-(f) shows the probability distribution (top) and phase (bottom)
of the ground state of $\hat{U}_{F}$, $|\epsilon_{i}\rangle$, in Fock space.
Note that the phase plotted is the argument of the wavefunction which is
different from the phase difference between the two wells, $\theta$. The first
thing to notice is that the probability distribution is a Gaussian and is not
spread over the entire FSL as one would expect for a translationally invariant
system. As previously mentioned, This confinement around $n\approx 0$ comes
from the square root terms since their leading order $n$ dependence takes the
form of a harmonic trap, $-\sqrt{N^{2}/4-n^{2}}\approx-N/2+n^{2}/N$. For
$\phi<\phi_{c}$ in panel (d) the phase of the ground state is slowly varying
since $\theta_{0}\approx 0$, however, for $\phi>\phi_{c}$ in panel (f)
$\theta_{0}\neq 0$ resulting in a sinusoidally varying phase in opposite
directions along each leg of the ladder. Vortices are located at the positions
where the phases on each leg of the ladder differ the most (red/green) and
create regions where the probability distribution almost vanishes. The
distribution does not completely go to zero, however, because positions in the
system cannot be resolved beyond the lattice spacing, so vortex locations
cannot be pinpointed exactly.

Figure 2: Chiral current as a function of magnetic flux for different system
sizes and boson-boson interaction energies. In both panels the red curves show
the thermodynamic limit result from Eq. (14). (a) The chiral current is
plotted for $N=20$ (blue triangles), $N=50$ (green circles) and $N=100$ (black
squares) and the other parameters values are $\xi=0.5$, $\mu=0$. As the system
size increases the numerical results approach the analytic result because the
larger the system size the more sensitive the system is to the presence of
vortices. (b) The chiral current is plotted for $\mu=5$ (blue triangles),
$\mu=0.5$ (green circles) and $\mu=0$ (black squares) and other parameter
values are $\xi=0.5$ and $N=100$. As $\mu$ increases the peak moves away from
the analytic result because the ground state wavefunction becomes more
localized around $n=0$ requiring larger fluxes to create a vortex with
significant overlap of the wavefunction.

This QPT has the same qualitative features as the transition from the Meissner
phase to the Abrikosov-vortex phase in type-II superconductors. For this
reason, we take a look at the probability current along the legs of the
ladder. The total current of the system is zero, however, the current
difference between the two legs is non-zero and is called the chiral current
[49]

$J_{C}=J_{R}-J_{L}=\partial_{\phi}E_{-}(\theta)|_{\theta=\theta_{0}}\,.$ (13)

On either side of the QPT the chiral current is expressed as

$\frac{2J_{C}}{NJ}=\begin{cases}\sin(\phi),&\text{if }\phi\leq\phi_{c}\\\
\frac{\xi^{2}\cos(\phi)}{\sin^{2}(\phi)\sqrt{\xi^{2}+\sin^{2}(\phi)}},&\text{if
}\phi>\phi_{c}\,.\end{cases}$ (14)

Equation (14) is plotted as the red curves in Fig. 2 (a) alongside the
numerical result using the ground state of $\hat{U}_{F}$ for different system
sizes. The finite size scaling comes from competition between two
lengthscales: the size of the probability distribution and the distance
between vortices. In panel (e) of Fig. 1 for $\phi=\phi_{c}$ vortices can be
seen in the phase of the wavefunction, but they do not overlap much with the
probability distribution, so they are insignificant. The width of the
probability distribution scales as $\sqrt{N}$, so for larger $N$ the easier it
will be for a vortex to penetrate it. In the thermodynamic limit, any presence
of a vortex will be felt by the wavefunction and the numerical results will
converge with the analytic result.

Repulsive boson-boson interactions have a similar affect on the probability
distribution as the square root factors in that they act as a confining
potential along the legs of the FSL. This comes from the fact that there is an
added energy cost for the BEC to occupy one well over the other and so there
is an added energy cost for states $n\neq 0$ in the FSL. The interaction
energy, $\mu$, plays the role of the confining potential strength, so for
stronger interactions larger fluxes are needed to produce a vortex that
penetrates the wavefunction. This leads to the peaks of the chiral currents
appearing at larger flux values which can be seen in panel (b) of Fig. 2 where
the chiral current is plotted for different values of $\mu$.

Figure 3: The maximum chiral current as a function of $\mu$. Each point of
data represents the maximum chiral current of a $J_{C}$ vs. $\phi$ curve for a
given value of $\mu$ like the ones seen in Fig. 2 (b). The value of $\mu$
where the peak occurs in this image is denoted as $\mu_{\mathrm{max}}$. The
inset shows $|\mu_{\mathrm{max}}-\mu_{c}|$ as a function of $1/N$ for
different system sizes. The black line is a linear fit and when extrapolated,
shows that the difference goes to zero in the thermodynamic limit, $1/N\to 0$.

Conversely, attractive boson-boson interactions ($\mu<0$) can enhance the
Meissner-vortex QPT because they make it energetically favorable for the BEC
to build-up in the wells rather that stay evenly dispersed as in the repulsive
case. This means the ground state wavefunction will be more spread out from
the region around $n=0$ and will be more susceptible to the presence of
vortices. In fact, in a normal BJJ without an impurity ($W=K=0$), there is a
critical interaction energy where the interactions overpower the confinement
of the square root factors and another QPT occurs involving the breaking of a
$\mathbb{Z}_{2}$ symmetry when the BEC spontaneously clumps into one well over
the other [50]. At the critical point of this QPT, the correlation length in
Fock space diverges and the width of the ground state wavefunction becomes
comparable to $N$. For a general magnetic flux, the interactions become
dominant when $\mu<-\frac{1}{2}\cos\phi$. Therefore, the Meissner-Abrikosov-
vortex transition will be critically enhanced when the two QPTs coincide at an
interaction energy of

$\mu_{c}=-\frac{1}{2}\cos\phi_{c}=\frac{\xi-\sqrt{\xi^{2}+4}}{4}\,.$ (15)

To test this prediction, we calculate the maxima of $J_{C}$ vs. $\phi$ curves
for different values of the interaction energy $\mu$. The results for $N=100$
are plotted in Fig. 3 which shows there is a clear value of $\mu$, denoted as
$\mu_{\mathrm{max}}$, where $J_{C}$ is a maximum. The maximum occurs near
$\mu_{c}$, but due to finite size effects, does not occur precisely at the
critical point. To show $\mu_{\mathrm{max}}$ converges to $\mu_{c}$ in the
thermodynamic limit we plot their difference as a function of $1/N$ for
different system sizes in the inset. A linear fit is performed on the data and
the resulting line is extrapolated to $1/N=0$ giving a difference of
$|\mu_{\mathrm{max}}-\mu_{c}|\approx 2.2\times 10^{-4}$ at that point. This
provides numerical support to the result in Eq. (15) and one can conclude that
the maximum possible chiral current occurs when the Meissner-Abrikosov-vortex
QPT and the BEC clumping QPT coincide.

##  IV Conclusion

We have shown that FSLs offer a new direction to explore in the ever-growing
field of quantum simulations of condensed matter systems. By using a FSL, the
usual techniques of simulating the Meissner effect with optical ladders and
arrays of Josephson junctions is reduced to a single periodically driven BJJ
coupled to an impurity. The system exhibits the Meissner-Abrikosov-vortex QPT
which we show can be enhanced by a second QPT involving the BEC spontaneously
clumping into one well when the boson-boson interactions reach a critical
attraction strength. The Meissner-Abrikosov-vortex QPT is characterized by the
chiral current, but a qualitative detection of the QPT can be achieved by
measuring the relative phase of the BEC between the two wells since it
transitions from zero in the Meissner phase to non-zero in the Abrikosov-
vortex phase. The relative phase can be measured from interference patterns in
standard time-of-flight images after the trapping potential has been turned
off [51, 52].

Although the system we have investigated is topologically trivial, additional
BEC-impurity interactions can change that. For instance, the inclusion of an
impurity-assisted tunneling term like $\hat{S}_{x}\hat{\sigma}_{x}$ can result
in a FSL version of the Su-Schrieffer-Heeger model. Pairing this with the fact
that a single impurity is enough to introduce a new synthetic dimension to the
system opens up interesting possibilities for simulations of topological
models in spaces higher than 2D. Furthermore, since two particles are enough
to form a 2D lattice and periodically driving their interactions generates the
synthetic gauge field, this work opens the door for the simulation of many-
body physics in the miniature. The question then arises as to what is the
smallest possible system one can use to simulate condensed matter phenomena of
interest.

##  Appendix A Numerical calculation of $J_{C}$

To calculate the chiral current numerically we start with Eq. (14) and use the
Hellmann-Feynman theorem to write

$J_{C}=\langle\epsilon_{0}|\partial_{\phi}\hat{H}_{\mathrm{eff}}|\epsilon_{0}\rangle$
(16)

where $|\epsilon_{0}\rangle$ is the ground state of $\hat{U}_{F}$. The result
of the derivative with respect to $\phi$ is

$J_{C}=\langle\epsilon_{0}|\hat{J}_{x}|\epsilon_{0}\rangle\sin\phi-\langle\epsilon_{0}|\hat{J}_{y}\hat{\sigma}_{z}|\epsilon_{0}\rangle\cos\phi$
(17)

which is what is used in our calculations.

##  Appendix B Entanglement entropy ($\mu=0$)

Figure 4: The entanglement entropy between the BEC and the impurity. The red
curve is the analytic expression from Eq. (LABEL:eq:analent) and the black
squares are from a numerical calculation. The parameters are $\xi=0.5$,
$\mu=0$ and $N=100$.

In the Meissner phase the screening current is entirely along the legs of the
ladder (BEC degrees of freedom). However, in the Abrikosov-vortex phase the
current develops along the rungs of the ladder (impurity degrees of freedom)
creating the vortices, so it is natural to ask whether or not the entanglement
between the BEC and the impurity is sensitive to the QPT. To check, we
calculate the entanglement entropy of the ground state in Eq. (8) in both
phases from the impurity reduced density matrix

$\rho_{\mathrm{I}}=|-,\theta_{0}\rangle\langle\theta_{0},-|$ (18)

where $\theta_{0}$ is from Eq. (10). The analytic expression for the
entanglement entropy is

$\displaystyle S(\rho_{\mathrm{I}})$ $\displaystyle=$
$\displaystyle-\mathrm{Tr}\left[\rho_{\mathrm{I}}\mathrm{ln}\rho_{\mathrm{I}}\right]$
$\displaystyle=$ $\displaystyle-
f_{+}(\xi,\phi)\mathrm{ln}f_{+}(\xi,\phi)-f_{-}(\xi,\phi)\mathrm{ln}f_{-}(\xi,\phi)$

where
$f_{\pm}(\xi,\phi)=\frac{1}{2}\left(1\pm\frac{\xi}{\sin\phi\sqrt{\xi^{2}+\sin^{2}\phi}}\right)$.
The numerical calculation is carried out with the ground state of
$\hat{U}_{F}$ giving the density matrix and impurity reduced density matrix

$\rho=|\epsilon_{0}\rangle\langle\epsilon_{0}|,\hskip
10.0pt\rho_{\mathrm{I}}=\mathrm{Tr}_{\mathrm{B}}\left[\rho\right]$ (20)

where $\mathrm{Tr}_{\mathrm{B}}$ is the trace over the BEC degrees of freedom.
The entanglement entropy is once again
$S(\rho_{\mathrm{I}})=-\mathrm{Tr}\left[\rho_{\mathrm{I}}\mathrm{ln}\rho_{\mathrm{I}}\right]$
which we compare to the analytic expression in Fig. 4 and find excellent
agreement between the two. We see that, indeed, the entanglement entropy is
sensitive to the QPT and follows the prediction of it increasing in the
Abrikosov-vortex phase ($\phi>\phi_{c}$). The disagreement between the
analytic and numerical results around the critical point comes from finite
size effects and is expected to vanish in the thermodynamic limit.

## References

  * [1] D. Jaksch and P. Zoller, The cold atom Hubbard toolbox, Ann. Phys. 315, 52 (2005); I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Rev. Mod. Phys. 80, 885 (2008).
  * [2] M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, Quantum phase transition from a superfluid to a Mott insulator in a gas of ultracold atoms, Nature 415, 39-44 (2002).
  * [3] J. Dalibard, F. Gerbier, G. Juzeli$\bar{\mathrm{u}}$nas, and P Öhberg, Colloquium: Artificial Gauge Potentials for Neutral Atoms, Rev. Mod. Phys. 83, 1523 (2011).
  * [4] N. Goldman, G. Juzeli$\bar{\mathrm{u}}$nas, P. Öhberg, and I. B. Spielman, Light-induced gauge fields for ultracold atoms, Rep. Prog. Phys. 77, 126401 (2014).
  * [5] N. Goldman and J. Dalibard, Periodically driven quantum systems: effective Hamiltonians and engineered gauge fields, Phys. Rev. x 4 031027 (2014).
  * [6] A. Eckardt and E. Anisimovas, High-frequency approximation for periodically driven quantum systems from a Floquet space perspective, New J. Phys. 17, 093039 (2015).
  * [7] M. Burkov, L. D’Alessio, and A. Polkovnikov, Universal high-frequency behavior of periodically driven systems: from dynamical stabilization to Floquet engineering, Adv. Phys. 64, 139-226 (2015).
  * [8] C. E. Creffield, G. Pieplow, F. Sols, and N. Goldman, Realization of uniform synthetic magnetic fields by periodically shaking an optical square lattice, New J. Phys. 18, 093013 (2016).
  * [9] M. Aidelsburger, M. Atala, S. Nascimbéne, S. Trotzky, Y.-A. Chen, and I. Bloch, Experimental realization of strong effective magnetic fields in an optical lattice, Phys. Rev. Lett. 107, 255301 (2011).
  * [10] D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Quantized Hall conductance in a two-dimensional periodic potential, Phys. Rev. Lett. 49, 405 (1982).
  * [11] F. D. M. Haldane, Model for a quantum Hall effect without Landau levels: condensed-matter realization of the ”parity anomaly”, Phys. Rev. Lett. 61, 2015 (1988).
  * [12] G. Jotzu, M. Messer, R. Desbuquois, M. Lebrat, T. Uehlinger, D. Greif, and T. Esslinger, Experimental realization of the topological Haldane model with ultracold fermions, Nature 515, 237-240 (2014).
  * [13] N. Fläschner, B. S. Rem, M. Tarnowski, D. Vogel, D.-S. Lühmann, K. Sengstock, and C. Weitenberg, Experimental reconstruction of the Berry curvature in a Floquet Bloch band, Science 352, 1091-1094 (2016).
  * [14] P. G. Harper, Single band motion of conduction electrons in a uniform magnetic field, Proc. Phys. Soc. A 68, 874 (1955); D. R. Hofstadter, Energy levels and wave functions of Bloch electrons in a rational and irrational magnetic fields, Phys. Rev. B 14, 2239 (1976).
  * [15] M. Aidelsburger, M. Atala, M. Lohse, J. T. Berreiro, B. Paredes, and I Bloch, Realization of the Hofstadter Hamiltonian with ultracold atoms in optical lattices, Phys. Rev. Lett. 111, 185301 (2013).
  * [16] H. Miyake, G. A. Siviloglou, C. J. Kennedy, W. C. Burton, and W. Ketterle, Realizing the Harper Hamiltonian with laser-assisted tuennling in optical lattics, Phys. Rev. Lett. 111, 185302 (2013).
  * [17] W. Meissner and R. Ochsenfeld, Ein neuer effekt bei eintritt supraleitfähigkeit, Naturwissenschaften 21, 787-788 (1933).
  * [18] M. Kardar, Josephson-junction ladders and quantum fluctuations, Phys. Rev. B 33, 3125 (1986).
  * [19] E. Granato, Phase transitions in Josephson-junction ladders in a magnetic field, Phys. Rev. B 42, 4797 (1990).
  * [20] E. Orignac and T. Giamarchi, Meissner effect in a bosonic ladder, Phys. Rev. B 64, 144515 (2001).
  * [21] C. Denniston and C. Tang, Phases of Josephson junction ladders, Phys. Rev. Lett 75, 3930 (1995).
  * [22] D. Hügel and B. Paredes, Chiral ladders and the edges of chern insulators, Phys. Rev. A 89, 023619 (2014).
  * [23] M. Atala, M. Aidelsburger, M. Lohse, J. T. Barreiro, B. Paredes, and I. Bloch, Observation of chiral currents with ultracold atoms in bosonic ladders, Nat. Phys. 10, 588-593 (2014).
  * [24] D.-W. Wang, H. Cai, R.-B. Liu, and M. O. Scully, Mesoscopic superposition states generated by synthetic spin-orbit interaction in Fock-state lattices, Phys. Rev. Lett. 116, 220502 (2016).
  * [25] H. Cai, D.-W. Wang, Topological phases of quantized light, Natl. Sci. Rev. 8, nwaa196 (2021).
  * [26] P. Saugmann and J. Larson, A Fock state lattice approach to quantum optics, arXiv:2203.13813 (2022).
  * [27] J. Larson and T. Mavrogordatos, The Jaynes-Cummings model and its descendants, Institute of Physics Publishing: Bristol, UK (2021).
  * [28] J. Mumford, Synthetic gauge field in two interacting ultracold atomic gases without an optical lattice, arXiv:2207.00742 (2022).
  * [29] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Feshbach resonances in ultracold gases, Rev. Mod. Phys. 82, 1225 (2010).
  * [30] D. H. Smith, Inducing resonant interactions in ultracold atoms with a modulated magnetic field, Phys. Rev. Lett. 115, 193002 (2015).
  * [31] M. Albiez, R. Gati, J. Fölling, S. Hunsmann, M. Cristiani, and M. K. Oberthaler, Direct observation of tunneling and nonlinear self-trapping in a single bosonic Josephon junction, Phys. Rev. Lett. 95, 010402 (2005).
  * [32] T. Zibold, E. Nicklas, C. Gross, and M. K. Oberthaler, Classical bifurcation at the transition from Rabi to Josephson dynamics, Phys. Rev. Lett. 105, 204101 (2010).
  * [33] C. S. Gerving, T. M. Hoang, B. J. Land, M. Anquez, C.D. Hamley, and M. S. Chapman, Non-equilibrium dynamics of an unstable quantum pendulum explored in a spin-1 Bose-Einstein condensate, Nat. Commun. 3, 1169 (2012).
  * [34] R. Gerritsma, A. Negretti, H. Doerk, Z. Idziaszek, T. Calarco, and F. Schmidt-Kaler, Bosonic Josephson junction controlled by a single ion, Phys. Rev. Lett. 109, 080402 (2012).
  * [35] S. Pal and C. Benjamin, Tuning the 0-$\pi$ Josephson junction with a magnetic impurity: Role of tunnel contacts, exchange coupling, e-e interactions and high-spin states, Sci. Reports 8, 5208 (2018).
  * [36] O. Boada, A. Celi, J. I. Latorre, and M. Lewenstein, Quantum simulation of an extra dimension, Phys. Rev. Lett. 108, 133001 (2012).
  * [37] A. Celi, P. Massignan, J. Ruseckas, N. Goldman, I. B. Spielman, G. Juzeli$\bar{\mathrm{u}}$nas, and M. Lewenstein, Synthetic gauge fields in synthetic dimensions, Phys. Rev. Lett. 112, 043001 (2014).
  * [38] T. Ozawa, H. M. Price, Topological quantum matter in synthetic dimensions, Nat. Rev. Phys. 1, 349-357 (2019).
  * [39] M. Mancini, G. Pagano, G. Cappellini, L. Livi, M. Rider, J. Catani, C. Sias, P. Zoller, M. Inguscio, M. Dalmonte, and L. Fallani, Observation of chiral edge states with neutral fermions in synthetic Hall ribbons, Science 349, 1510-1513 (2015).
  * [40] B. K. Stuhl, H.-I. Lu, L. M. Aycock, D. Genkina, and I. B. Spielman, Visualizing edge states with an atomic Bose gas in the quantum Hall regime, Science 349, 1514-1518 (2015).
  * [41] E. Anisimovas, R. Račiūnas, C. Sträter, A. Eckardt, I. B. Spielman, and G. Juzeliūnas, Semisynthetic zigzag optical lattice for ultracold bosons, Phys. Rev. A 94, 063631 (2016).
  * [42] F. A. An, E. J. Meier, and B. Gadway, Diffusive and arrested transport of atoms under tailored disorder, Nat. Comm. 8, 325 (2017); F. A. An, E. J. Meier, J. Ang’ong’a, and B. Gadway, Correlated dynamics in a synthetic lattice of momentum states, Phys. Rev. Lett. 120, 040407 (2018).
  * [43] E. J. Meier, F. A. An, and B. Gadway, Observation of the topological soliton state in the Su-Schrieffer-Heeger model, Nat. Comm. 7, 13986 (2016).
  * [44] D. Xie, W. Gou, T. Xiao, B. Gadway, and B. Yan, Topological characterizations of an extended Su-Schrieffer-Heeger model, npj Quant. Inf. 5, 1-5 (2019).
  * [45] H. M. Price, T. Ozawa, and N. Goldman, Synthetic dimensions for cold atoms from shaking a harmonic trap, Phys. Rev. A 95, 023607 (2017).
  * [46] J. Floß, A. Kamalov, I. Sh. Averbukh, and P. H. Bucksbaum, Observation of Bloch oscillations in molecular rotation, Phys. Rev. Lett. 115, 203002 (2015).
  * [47] S. K. Kanungo, J. D. Whalen, Y. Lu, M. Yuan, S. Dasgupta, F. B. Dunning, K. R. A. Hazzard, and T. C. Killian, Realizing topological edge states with Rydberg-atom synthetic dimensions, Nat. Commun. 13, 972 (2022).
  * [48] T. Chalopin, T. Satoor, A. Evrard, V. Makhalov, J. Dalibard, R. Lopes, and S. Nascimbène, Probing chiral edge dynamics and bulk topology of a synthetic Hall system, Nat. Phys. 16, 1017-1021 (2020).
  * [49] S. Greschner, M. Piraud, F. Heidrich-Meisner, I. P. McCulloch, U. Schollwöck, and T Vekua, Symmetry-broken states in a system of interacting bosons on a two-leg ladder with a uniform Abelian gauge field, Phys. Rev. A 94, 063628 (2016).
  * [50] P. Buonsante, R. Burioni, E. Vescovi, and A. Vezzani, Quantum criticality in a bosonic Josephson junction, Phys. Rev. A 85, 043625 (2012).
  * [51] M. R. Andrews, C. G. Townsend, H.-J Miesner, D. S. Durfee, D. M. Kurn, and W. Ketterle, Observation of interference between two Bose condensates, Science 275, 637 (1997).
  * [52] R. Gati, J. Esteve, B. Hammerling, T. B. Ottenstein, J. Appmeier, A. Weller, and M. K. Oberthaler, A primary noise thermometer for ultracold Bose gases, New J. Phys. 8, 189 (2006).