# Invariant Scattering Transform for Medical Imaging

Nafisa Labiba Ishrat Huda
Research and Development Department, Pioneer Alpha,
Dhaka, Bangladesh
<EMAIL_ADDRESS>
&Angona Biswas
Research and Development Department, Pioneer Alpha,
Dhaka, Bangladesh
<EMAIL_ADDRESS>
&MD Abdullah Al Nasim
Research and Development Department, Pioneer Alpha,
Dhaka, Bangladesh
<EMAIL_ADDRESS>
&Md. Fahim Rahman
Research and Development Department, Pioneer Alpha,
Dhaka, Bangladesh
fahimpranto002gmail.com
&Shoaib Ahmed
Research and Development Department, Pioneer Alpha,
Dhaka, Bangladesh
<EMAIL_ADDRESS>

###### Abstract

Invariant scattering transform introduces new area of research that merges the
signal processing with deep learning for computer vision. Nowadays, Deep
Learning algorithms are able to solve a variety of problems in medical sector.
Medical images are used to detect diseases brain cancer or tumor, Alzheimer’s
disease, breast cancer, Parkinson’s disease and many others. During pandemic
back in 2020, machine learning and deep learning has played a critical role to
detect COVID-19 which included mutation analysis, prediction, diagnosis and
decision making. Medical images like X-ray, MRI known as magnetic resonance
imaging, CT scans are used for detecting diseases. There is another method in
deep learning for medical imaging which is scattering transform. It builds
useful signal representation for image classification. It is a wavelet
technique; which is impactful for medical image classification problems. This
research article discusses scattering transform as the efficient system for
medical image analysis where it’s figured by scattering the signal information
implemented in a deep convolutional network. A step by step case study is
manifested at this research work.

_K_ eywords Medical imaging, diagnosis, deep learning, scattering transform,
COVID-19

##  1 Introduction

In the field of medical imaging, a professional creates various images for
diagnostic purposes. Medical imaging [1, 2, 3, 4] plays a significant role in
modern science. It’s not only includes X-rays [5], MRI [6], CT [7] scans as
neuro [8] images but also uses endoscopy, ultrasounds, tactile imaging.
Artificial intelligence is enhancing the capability to analyze results in
medical imaging technology. Also computer vision is used to diagnose human
conditions visually. Manual medical image review procedures have saved
countless lives over the years, but that doesn’t mean the procedures can’t be
improved. Machine learning (ML) technology can analyze and learn from vast
amounts of data while performing computer vision (CV) analysis on an image to
find tiny differences that the human eye cannot detect. Deep Learning (DL) [9]
is a form of machine learning (ML) based on artificial neural networks; these
technologies have made automated, precise and ultra-fast analysis of medical
images a reality. These techniques can achieve high accuracy in analyzing
medical images [10, 11, 12, 13, 14, 15] for detection, segmentation,
classification and prediction.

As for image classification, there are two techniques which are supervised and
unsupervised. Supervised classification has labeled data images; these are
trained data. On the other hand, unsupervised includes unlabeled data that
does not have any trained images. It analyzes result by discovering unknown or
hidden pattern. To perform an image classification using ML or DL neural
networks [16], there are some steps to follow. First, image pre-processing is
to improve image data so that models can benefit from the data. Image
resizing, data augmentation; these are the parts of pre-processing. Then comes
feature extraction, it uses for identifying the pattern of the given image.
Features can be different to particular classes. This technique helps models
to find out the differences between two images. There are also training and
testing part where data can be labeled or unlabeled. The last part of this
process is classification of the image whether it’s matched with the targeted
ones or not.

On the subject of classification, there are many classifiers such as support
vector machine (SVM), K-nearest neighbor (KNN), naïve bayes algorithm, random
forest (RB) algorithm. Support vector machine (SVM) demonstrates different
classes in a hyper plane in multi-dimensional space. It divides datasets into
classes and find out the maximum marginal hyper plane. On the other hand, KNN
look at the given data then calculates its distance; find nearest neighbors
and count the vote to check the majority of labels. Random forest (RF) is used
for both classification and regression. This algorithm consists of many
decision trees. From those trees, it predicts the majority one as a result.

To use medical images in machine learning and deep learning applications with
minimal configuration, a wavelet scattering [17, 18, 19] enables to obtain low
variance features from real time image data. This scattering [20] process is
acknowledged as scattering transform. Scattering network uses predefined
wavelets and scaling filters. It processes data in stages. The outcome of
previous stage is the input of next stage. There are three operations in each
stage.

Figure 1: Three operations of scattering transform.

Fig. 1 is showing the operations of scattering transform network. First, it
takes the wavelet transform of the input data with each wavelet filter. Then
takes the modulus of the each filtered outputs and average that modulus with
the scaling filter. This process is repeated at every node.

##  2 Literature Review

The authors in [21] used scattering transform based features and texture
features for iris detection. Not only these two techniques are used; they also
applied principal component analysis (PCA) on the extracted features to reduce
the dimensionality of the feature. The mean and variance of each scattering
images are used as features. Iris database is used by them where it contains
2240 images from 224 different people. An integrated scattering feature for
medical image is proposed by the authors of [22]. A corresponding feature
demonstration strategy for each integration model is used. Rushi and Yicong
applied histogram of compressed scattering coefficients also called HCSCs
[23]. This proposal took advantages from both compressed scattering
coefficients and bag-of-words (BoW) model. For intrapartum fetal heart rate
(FHR) variability fractal analysis authors in [24] also used scattering
transform. Paper [25] presented automated diagnosis of Alzheimer disease using
the scale invariant feature transforms from brain [26, 27, 28] MR images.
Support vector machine (SVM) also applied to the data. The authors of [29]
applied wavelet scattering transform (WST) [30] on CT lung images to detect
lung cancer and support vector machine (SVM), kernel nearest neighbor (KNN),
wavelet scattering transform, shearlet transform, random forest (RF) as
classifiers.

They [31] focused on Morlet wavelets and learned the scales, orientations, and
aspect ratios of the filters to produce problem-specific parameterizations of
scattering transform; also their learned versions of the scattering transform
yield significant performance gains in small-sample classification settings
over the standard scattering transform. In [32], they used scattering
transform registration algorithm to align high and low resolution images. A
varifold based modeling framework also used by them to compute multiscale
data. They [33] developed a novel network where they used features of a
scattering transform for an attention mechanism. An approach [34] was proposed
to classify 17-classes of cardiac arrhythmia using wavelet scattering
transform (WST). The paper [35] focused on the development of a seizure
detection method with the above-stated characteristics where group invariant
scattering, novel data representation technique used for feature extraction.
They [36] built a classification scheme that involved a feature extraction and
classification stages; they used scattering transform for extracting
distinguishes features robust to small deformations in the image. For large
dimension of scattering feature matrix, [37] used multidimensional scaling
(MDS) method to reduce the dimension and the method was compared with the
classical dimension reduction method known as principal component analysis
(PCA). The paper [38] involved the utilization of Support Vector Machine (SVM)
classifiers with three non-linear kernels, along with the implementation of
k-fold cross-validation to stage detection of Dyslexia in EEG signals using
wavelet scattering transform (WST).

In [39], they applied a feature learning algorithm known as Wavelet scattering
network to retinal fundus images and studied the efficiency of the
automatically extracted features therefore for glaucoma detection. [40]
projected hybrid framework that fuses the complex shearlet scattering
transform (CSST) and a convolutional neural network (CNN) into a single model.
A novel motion-compensation method for investigating in vivo myocardium
structures in DTI with free-breathing acquisitions was proposed by paper [41].
The paper [42] used the wavelet scattering transform (WST) as a feature
extraction technique to obtain features from PPG data and combined it with
clinical data for detection of early hypertension stages by applying Early and
Late Fusion. To calculate invariants from the responses to the set of
spherical harmonics projected onto 3D kernels in the form of a lightweight
Solid Spherical Energy (SSE) CNN (convolutional neural network) [43] was
applied by them [44]. A deep scattering network was used for palm print
recognition [45]; scattering network also used in this process to extract
features. Fibered confocal fluorescence microscopy (FCFM) imaging technique
and scattering features for lung cancer detection [46]. A hybrid Scattering
Coefficients - Bag of Visual Words – Discrete Wavelet Transform (SC-BoVW- DWT)
relevance fusion algorithm was proposed for effective content-based medical
image retrieval (CBMIR) [47].

Texture[48] descriptor derived from a region covariance matrix of scattering
coefficients for detection of glandular structures [49] in colon histology
images. The procedure [50] used features of a wavelet scattering transform to
classify signal segments whether it’s chaotic or non-chaotic. The proposed
method [51] of classification determined the level of the brain tumor using
the features of the segments generated based on Bayesian fuzzy clustering.
They [52] proposed a twin branch CNN that learned robust representations by
combining deep narrow band features and WST features of the narrow band images
to classify vocal cord NBI images into malignant and benign classes. They [53]
approached a brain tumor classification using a hybrid deep auto encoder with
a Bayesian fuzzy clustering-based segmentation. They [54] aimed to analyze the
clot texture using the scattering operator which combined wavelet transform
convolutions with non-linear modulus and averaging operators. They [55]
applied a novel algorithm which offered immunity to local errors in the
underlying deformation field obtained from registration procedures.

##  3 Methodology of study

###  3.1 Wavelet Scattering Transform

Low variance features from real-time series image data is enabled by wavelet
scattering transform (WST) network. WST [56, 57] transforms small deformations
by separating the variations across different scales. It provides a light
representation. There are wavelets, scaling function and input data. There are
also a number of users’ particular rotations of the wavelet. Edges from the
root to a node are referred as a path. Scattering coefficients are the
coefficients convolved with scaling function. These scattering coefficients
are the low variance features derived from the data.

$Hx(u)=x*h(2u)$ (1)

$Gx(u)=x*g(2u)$ (2)

where h is a low frequency and g is a high frequency.

###  3.2 Feature Extraction

To extract features from data, wavelet Scattering will create and configure
the network. Parameters include the size of invariance scale, the number of
filter banks and the number of wavelets per octave in each filter bank also
set the number of rotation per wavelet. To derive time series and image data
features use the object function ScatteringTransform.

###  3.3 Invariance Scale

In the wavelet scattering network, scaling filter plays an important role.
This network is invariant to translations up to the invariance scale. Scaling
function determines the size of the invariant in time or space. Time
invariance uses for time series data where duration is the invariance scale.
The invariance scale does not exceed by the time support of a wavelet.
Frequencies are lower than the invariance scale. Image invariance is for image
data, the scale specifies in pixels. The scaling function creates a wavelet
image scattering network for the image size.

###  3.4 Quality Factor

For the scattering filter banks, quality factors should be set to create a
wavelet scattering network. The quality factor for each filter bank is the
number of wavelet filter per octave. Large quality factors are not necessary
for image data.

###  3.5 Case Description with Dataset information

This section will discuss about the details of different cases. IIT Delhi
database [21] contains 2240 iris images captured from 224 different people.
The images of 21 people were for validation set; around half of the images
were used as training and testing images. In [22], TCIA-CT database has 675
images of 19 categories. NEMA-CT and EXACT09-CT used for [23] where EXACT09
consists of CT scans at the chest, size of 512 x 512. There are total number
of 675 CT images in the database. HCSCs methods were applied to both of the
databases. The data [25] consists of the scans from 98 normal and 100 scan
from Alzheimer disease (AD). The subjects include both men and women. Whole
data was for training set and one subject’s data was for testing set. In
addition, by analyzing the data a method’s effectiveness can be increased.

The Bilkent dataset was used by them [49] which consists of 72 microscopic
images of H&E stained healthy colon biopsies from 36 patients; where the
resolution of images was 480x640 and the dataset was divided into 24 images
for training and 48 images for testing. For [58], they used two-dimensional
image slices extracted from MRI image stacks. For detecting seizure [35], open
access databases featuring multichannel scalp EEG was used; the used dataset
is sampled at 256 samples per second and quantized with a 16 bit analog to
digital converter where each case contains continuous recordings of the
patients and 140 records contained 197 seizures. The testing and training
datasets [50] were created using 13 dynamical systems of first, second or
third order; each system was provided with 1000 created test files, each of
which contained 1000 samples. BRATS database for [51] which had two grades
(datasets) of tumor images and every patient image was stored as four
different modalities. The dataset [46] was composed of 103 images taken from 8
healthy volunteers, and 70 images from 7 patients with diagnosed bronchial
squamous cell carcinoma. The approached model consisted of two parallel
branches where each branch was a ResNet18. The first branch received the NBI
images of vocal cord and the second branch received the wavelet scattering
transform feature maps of the NBI images of vocal cord. They [22] used data
from the National Lung Screening Trial (NLST) in which individuals were
screened for lung cancer using low-dose CT (LDCT). Annotated 3700 NLST images
contained single solid nodules between 4 and 20 mm in diameter were included
in their [59] dataset. A dataset of 15 images of H&E stained colorectal cancer
tissues of 20× optical magnification were used by [60]; where they used a
window of size 27 × 27 and stride of 3 to make the patches along spatial
dimension (height and width) from the input feature map. In [53], they used
the skull removed images (dataset BRATS 2015). [45] used 6000 palm prints
sampled from 500 persons. Each palm print was taken under 4 different lights
in two different days and preprocessed and had a size of 128x128. Both with
synthetic images and real brain image datasets were used in [55].

The first dataset contained 1000 synthetic 32×32×32 texture volumes of two
classes [44]. The second one was a subsample of the American National Lung
Screening Trial (NLST) annotated by radiologists which included 485 pulmonary
nodules, 244 benign and 241 malignant, from distinct patients in CT. For the
work [42] was collected at the Guilin People’s Hospital in Guilin, China.
There were a total of 219 samples from patients admitted to the hospital; the
sample had 37% of subjects diagnosed with NT, 38% with PHT and 25% with Stage
1 and Stage 2 hypertension. The approach [40] experimented on two public
datasets demonstrate the superiority of our method. The categorized 17
arrhythmia classes were taken from the MIT-BIH arrhythmia database; also they
were having 1000 ECG fragments of 45 subjects [34].

Figure 2: Iris images from first layer of scattering transform from [21].
Figure 3: MR images extracted by the SIFT method [25]. Figure 4: CT images
from EXACT09-CT database [22]. Figure 5: tumor images from [51]. Figure 6:
medical images from [58]. Figure 7: palmprint scattering images from [45].
Figure 8: clot images from [54].

##  4 Result and analysis

An analysis of some related methods has been presented here.

From Table 1, it is noticed that several studies have been utilized scattering
transform method. Others have been conducted different techniques. In paper
[21], they used scattering transform and textural features where they gained
best accuracy of 99.2 percent. For each image they [21] applied two levels
with a set of filter banks. Minimum distance classifier has been used. There
has been two types of integrated scattering data [22] known as data
concentration and canonical analysis. To perform scattering transform to a
medical image to achieve its translation invariant demonstration, HCSCs method
should develop. After that, histogram was derived as feature vector using BoW
framework. The proposal [25] used scale-invariant feature transform (SIFT)
method for 2-D scan of Magnetic Resonance Imaging (MRI).

From other references, some observations are noted here. The proposed method
[35] successfully detected 180 out of 197 seizures (91.4%) across 24 subjects.
This process did not require any pre-processing on the data which was
unsupervised and fully automated. CT-images [47] used as a query image and
matched with other images in the database. The TCIA-CT database the proposed
SC-BoVW-DWT obtained better results than other methods. The proposed method
[51] of tumor classification used the two datasets where the method of tumor
segmentation and classification offered better accuracy when compared with the
other methods. The PPG signal [50] was split into 1000 samples; there were
31.2s segments. The WST features were obtained for each segment. The
experiment [49] was conducted on 24 Bilkent training images for super pixel
features, using 4-fold cross validation; the result showed encouraging
segmentation accuracy. They [55] used three different datasets with different
noise levels in the deformation field for registration. 6000 palm prints were
sampled from 500 persons from the PolyU palm print database [45]. Each palm
print was taken under 4 different lights in two different days; preprocessed
and had a size of 128x128, also each image was first divided into blocks of
size 32x32. They used half of the images of each person in the dataset for
training and the other half for testing and for SVM they were able to get
99.9% accuracy rate. The recognition accuracy of minimum distance classifier
was slightly less than SVM but still high in percentage. The work [60] has
showed that detection accuracy of deep convolutional neural networks can be
improved by adding handcrafted features to the raw data. The Bayesian fuzzy
clustering-based segmentation technique [53] for brain tumor classification
reached the highest accuracy (98.5 %) when compared to other existing systems.
Their [52] results indicated that incorporating wavelet scattering features to
improve the generalizing capability of the CNN.

The mentioned SSE-CNN [44] largely outperformed the latter > 90% with N > 2
versus 84.0% on the synthetic dataset. The analysis [42] showed that the
wavelet scattering transform (WST) combined with a support vector machine
(SVM) can classify norm tension and pre hypertension with an accuracy of
71.42% and an F1-score of 76%. They [39] achieved 98% on the held-out test
set. The WST-based approach [38] achieved a high accuracy rate, with an
average accuracy of 96.96% and 97.12% for dataset 1 and dataset 2 using the
RBF kernel. The research [61] showed the accuracy of 97.4% and 100% in
diagnosis abnormal retina and DR from normal ones. They also achieved the
accuracy of 84.2% in classifying OCT images to five classes of normal, CSR,
MH, AMD and DR which outperformed other state of the art methods with high
computational complexity. The classification accuracy of the proposed method
[37] can reach 98% or greater than that. The proposed methodology [35]
successfully detected 180 out of 197 seizures. An overall accuracy of 98.90%
was achieved by them [34]. This methodology [32] improved the results of
calcium segmentation up to 8.2% in terms of the Dice coefficient and 24.8% in
terms of the modified Hausdorff distance.

Table 1: Performance evaluation of various studies for scattering transform of medical imaging: Authors |  Method |  Dataset |  Purpose |  Medical image |  Accuracy
---|---|---|---|---|---
Shervin et al.[21] |  Scattering transform and textural features |  IIT Delhi iris data |  Iris recognition |  Iris |  99.2%
Rushi et al. [22] |  Integrated scattering feature |  TCIA-CT, EXACT09-CT data |  Application to medical retrieval |  CT |
Rushi and Yicong [23] |  HCSCs |  NEMA, TCIA, EXACT09-CT data |  Develop HCSCs for medical image retrieval |  CT |  NEMA-CT: 3.16%, 7.37%, 6.30% ; TCIA-CT: 7.60%, 10.92%, 10.53%
Joakim et al. [24] |  Scattering transform |  FIGO-TN, FIGO-TP, FIGO-FP |  Fetal heart rate variability fractal analysis |  CTG |  FIGO-TN, FIGO-TP (j1=2)
Mohammad et al. [25] |  Scale invariant feature transform |  OASIS MRI brain data |  Automated diagnosis of Alzheimer disease |  MR |  86%
Osama and Yoshifumi [29] |  WST and AI techniques |  CT scan image data |  Lung cancer detection |  CT lung |  WST+RF: 99.90% WST+KNN: 95.28% WST+SVM: 93.24%
Tameem et al. [62] |  Scattering transform |  AGEhIV, OASIS, ADNI data |  Disease classification |  MRI |  82.7%

##  5 Conclusion

This chapter represents some studies on scattering transform and other methods
of medical imaging [63, 64, 65, 66, 67, 68]. There are different types of
images such as CT scans, MRI [69], CTG, iris and CT lung imaging for diagnosis
purposes. SIFT and HOG can be feature descriptors. The minimum classifier is
there to match new images with train data. Principal component analysis (PCA)
can apply on features to reduce dimensionality. A 3D scattering transform can
replace deep learning methods when we won’t have any massive amount of data
points. Some methods were used with different classifiers and got best
accuracy. A combination of wavelet scattering transform with random forest,
k-nearest neighbor and support vector machine bring the best result
respectively. But wavelet scattering transform with random forest outperforms
all other methodologies. The ability to wavelet scattering transforms (WST) to
take image and convert it to signal. It understands the dataset and extract
features from it. The results show that this method can perform well for
diagnosis, so it can be useful for medical applications.

## References

  * [1] Valentin Oreiller, Julien Fageot, Vincent Andrearczyk, John O Prior, and Adrien Depeursinge.  Robust multi-organ nucleus segmentation using a locally rotation invariant bispectral u-net.  In International Conference on Medical Imaging with Deep Learning, pages 929–943. PMLR, 2022.
  * [2] Xiaoqin Wang, Rushi Lan, Huadeng Wang, Zhenbing Liu, and Xiaonan Luo.  Fine-grained correlation analysis for medical image retrieval.  Computers & Electrical Engineering, 90:106992, 2021.
  * [3] Zhen Guo, Zhiguang Liu, Qihang Zhang, George Barbastathis, and Michael E Glinsky.  Noise-resilient approach for deep tomographic imaging.  arXiv preprint arXiv:2211.15456, 2022.
  * [4] Alex Gaudio, Asim Smailagic, Christos Faloutsos, Shreshta Mohan, Elvin Johnson, Yuhao Liu, Pedro Costa, and Aurélio Campilho.  Deepfixcx: Explainable privacy-preserving image compression for medical image analysis.  Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery, page e1495, 2023.
  * [5] James Anthony Seibert and John M Boone.  X-ray scatter removal by deconvolution.  Medical physics, 15(4):567–575, 1988.
  * [6] Kaitlin M Stouffer, Menno P Witter, Daniel J Tward, and Michael I Miller.  Projective diffeomorphic mapping of molecular digital pathology with tissue mri.  Communications Engineering, 1(1):44, 2022.
  * [7] Bernd Ohnesorge, Thomas Flohr, and Klaus Klingenbeck-Regn.  Efficient object scatter correction algorithm for third and fourth generation ct scanners.  European radiology, 9(3):563–569, 1999.
  * [8] Won Hwa Kim, Sathya N Ravi, Sterling C Johnson, Ozioma C Okonkwo, and Vikas Singh.  On statistical analysis of neuroimages with imperfect registration.  In Proceedings of the IEEE international conference on computer vision, pages 666–674, 2015.
  * [9] Ramkrishna Bharath and Pachamuthu Rajalakshmi.  Deep scattering convolution network based features for ultrasonic fatty liver tissue characterization.  In 2017 39th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), pages 1982–1985. IEEE, 2017\.
  * [10] Glenn R Easley, Monica Barbu-McInnis, and Demetrio Labate.  Image registration using the shearlet transform.  In Wavelets and sparsity XVI, volume 9597, pages 125–135. SPIE, 2015.
  * [11] Berkman Sahiner, Aria Pezeshk, Lubomir M Hadjiiski, Xiaosong Wang, Karen Drukker, Kenny H Cha, Ronald M Summers, and Maryellen L Giger.  Deep learning in medical imaging and radiation therapy.  Medical physics, 46(1):e1–e36, 2019.
  * [12] Mario Merone, Carlo Sansone, and Paolo Soda.  A computer-aided diagnosis system for hep-2 fluorescence intensity classification.  Artificial intelligence in medicine, 97:71–78, 2019.
  * [13] Ahmed Al-Saffar, A Zamani, A Stancombe, and A Abbosh.  Operational learning-based boundary estimation in electromagnetic medical imaging.  IEEE Transactions on Antennas and Propagation, 70(3):2234–2245, 2021.
  * [14] Fares Kahlessenane, Amine Khaldi, Redouane Kafi, and Salah Euschi.  A robust blind medical image watermarking approach for telemedicine applications.  Cluster computing, 24(3):2069–2082, 2021.
  * [15] JA Seibert and JM Boone.  Medical image scatter suppression by inverse filtering.  In Medical Imaging II, volume 914, pages 742–750. SPIE, 1988.
  * [16] Md Aminul Haque Palash, Akib Khan, Kawsarul Islam, MD Abdullah Al Nasim, and Ryan Mohammad Bin Shahjahan.  Incongruity detection between bangla news headline and body content through graph neural network.  In The Fourth Industrial Revolution and Beyond: Select Proceedings of IC4IR+, pages 375–387. Springer, 2023.
  * [17] Raghu G Raj, Maxine R Fox, and Ram M Narayanan.  Target classification in synthetic aperture radar images using quantized wavelet scattering networks.  Sensors, 21(15):4981, 2021.
  * [18] Amirmohammad Shamaei, Jana Starcuková, and Zenon Starcuk Jr.  A wavelet scattering convolutional network for magnetic resonance spectroscopy signal quantitation.  In BIOSIGNALS, pages 268–275, 2021.
  * [19] Tuan D Pham.  Classification of motor-imagery tasks using a large eeg dataset by fusing classifiers learning on wavelet-scattering features.  IEEE Transactions on Neural Systems and Rehabilitation Engineering, 31:1097–1107, 2023.
  * [20] R Bharath and Pachamuthu Rajalakshmi.  Webrtc based invariant scattering convolution network for automated validation of ultrasonic videos for iot enabled tele-sonography.  In 2018 IEEE 4th World Forum on Internet of Things (WF-IoT), pages 790–795. IEEE, 2018.
  * [21] Shervin Minaee, AmirAli Abdolrashidi, and Yao Wang.  Iris recognition using scattering transform and textural features.  In 2015 IEEE signal processing and signal processing education workshop (SP/SPE), pages 37–42. IEEE, 2015.
  * [22] Rushi Lan, Huadeng Wang, Si Zhong, Zhenbing Liu, and Xiaonan Luo.  An integrated scattering feature with application to medical image retrieval.  Computers & Electrical Engineering, 69:669–675, 2018.
  * [23] Rushi Lan and Yicong Zhou.  Medical image retrieval via histogram of compressed scattering coefficients.  IEEE journal of biomedical and health informatics, 21(5):1338–1346, 2016.
  * [24] Václav Chudáček, Joakim Andén, Stéphane Mallat, Patrice Abry, and Muriel Doret.  Scattering transform for intrapartum fetal heart rate variability fractal analysis: a case-control study.  IEEE Transactions on Biomedical Engineering, 61(4):1100–1108, 2013\.
  * [25] Mohammad Reza Daliri.  Automated diagnosis of alzheimer disease using the scale-invariant feature transforms in magnetic resonance images.  Journal of medical systems, 36:995–1000, 2012.
  * [26] Faisal Muhammad Shah, Tonmoy Hossain, Mohsena Ashraf, Fairuz Shadmani Shishir, MD Abdullah Al Nasim, and Md Hasanul Kabir.  Brain tumor segmentation techniques on medical images-a review.  International Journal of Scientific & Engineering Research, 10(2):1514–1525, 2019.
  * [27] Md Abdullah Al Nasim, Abdullah Al Munem, Maksuda Islam, Md Aminul Haque Palash, Md Mahim Anjum Haque, and Faisal Muhammad Shah.  Brain tumor segmentation using enhanced u-net model with empirical analysis.  In 2022 25th International Conference on Computer and Information Technology (ICCIT), pages 1027–1032. IEEE, 2022.
  * [28] Tonmoy Hossain, Fairuz Shadmani Shishir, Mohsena Ashraf, MD Abdullah Al Nasim, and Faisal Muhammad Shah.  Brain tumor detection using convolutional neural network.  In 2019 1st international conference on advances in science, engineering and robotics technology (ICASERT), pages 1–6. IEEE, 2019.
  * [29] Osama A Omer and Yoshifumi Saijo.  Lung cancer detection using wavelet scattering transform and artificial intelligence technique.  Research Square, 2023.
  * [30] Upanshu Sharma and Remco Duits.  Left-invariant evolutions of wavelet transforms on the similitude group.  Applied and Computational Harmonic Analysis, 39(1):110–137, 2015\.
  * [31] Shanel Gauthier, Benjamin Thérien, Laurent Alsene-Racicot, Muawiz Chaudhary, Irina Rish, Eugene Belilovsky, Michael Eickenberg, and Guy Wolf.  Parametric scattering networks.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 5749–5758, 2022.
  * [32] Lennart Bargsten, Katharina A Riedl, Tobias Wissel, Fabian J Brunner, Klaus Schaefers, Michael Grass, Stefan Blankenberg, Moritz Seiffert, and Alexander Schlaefer.  Attention via scattering transforms for segmentation of small intravascular ultrasound data sets.  In Medical Imaging with Deep Learning, pages 34–47. PMLR, 2021\.
  * [33] Kaitlin M Stouffer, Zhenzhen Wang, Eileen Xu, Karl Lee, Paige Lee, Michael I Miller, and Daniel J Tward.  From picoscale pathology to decascale disease: Image registration with a scattering transform and varifolds for manipulating multiscale data.  In Multimodal Learning for Clinical Decision Support: 11th International Workshop, ML-CDS 2021, Held in Conjunction with MICCAI 2021, Strasbourg, France, October 1, 2021, Proceedings 11, pages 1–11. Springer, 2021\.
  * [34] Sudestna Nahak, Akanksha Pathak, and Goutam Saha.  Fragment-level classification of ecg arrhythmia using wavelet scattering transform.  Expert Systems with Applications, 224:120019, 2023.
  * [35] Muhammad Zubair Ahmad, Awais Mehmood Kamboh, Sajid Saleem, and Amir Ali Khan.  Mallat’s scattering transform based anomaly sensing for detection of seizures in scalp eeg.  IEEE Access, 5:16919–16929, 2017.
  * [36] Jinghui Li, Li Ke, Qiang Du, Xiaodi Ding, Xiangmin Chen, and Danni Wang.  Heart sound signal classification algorithm: a combination of wavelet scattering transform and twin support vector machine.  IEEE Access, 7:179339–179348, 2019.
  * [37] Shrevin Minaee and Yao Wang.  Palmprint recognition using deep scattering network.  In 2017 IEEE international symposium on circuits and systems (ISCAS), pages 1–4. IEEE, 2017.
  * [38] Agnieszka Szczęsna, Dariusz Augustyn, Henryk Josiński, Adam Świtoński, Paweł Kasprowski, and Katarzyna Harężlak.  Novel photoplethysmographic signal analysis via wavelet scattering transform.  In Computational Science–ICCS 2022: 22nd International Conference, London, UK, June 21–23, 2022, Proceedings, Part III, pages 641–653. Springer, 2022.
  * [39] Jutta Bikowski, Kim Knudsen, and Jennifer L Mueller.  Direct numerical reconstruction of conductivities in three dimensions using scattering transforms.  Inverse Problems, 27(1):015002, 2010.
  * [40] Axel H Masquelin, Nicholas Cheney, C Matthew Kinsey, and Jason HT Bates.  Wavelet decomposition facilitates training on small datasets for medical image classification by deep learning.  Histochemistry and cell biology, 155:309–317, 2021.
  * [41] Erick Martinez-Ríos, Luis Montesinos, and Mariel Alfaro-Ponce.  A machine learning approach for hypertension detection based on photoplethysmography and clinical data.  Computers in Biology and Medicine, 145:105479, 2022.
  * [42] Vu-Lam Nguyen, Ngoc-Son Vu, Hai-Hong Phan, and Philippe-Henri Gosselin.  An integrated descriptor for texture classification.  In 2016 23rd International Conference on Pattern Recognition (ICPR), pages 2006–2011. IEEE, 2016.
  * [43] Kotra Sankar Raja Sekhar, Tummala Ranga Babu, Goriparthi Prathibha, Kotra Vijay, and Long Chiau Ming.  Dermoscopic image classification using cnn with handcrafted features.  Journal of king Saud University-science, 33(6):101550, 2021.
  * [44] Vu-Lam Nguyen, Ngoc-Son Vu, Hai-Hong Phan, and Philippe-Henri Gosselin.  Lbp-and-scatnet-based combined features for efficient texture classification.  Multimedia Tools and Applications, 76:22425–22444, 2017.
  * [45] Alain Rakotomamonjy, Caroline Petitjean, Mathieu Salaün, and Luc Thiberville.  Scattering features for lung cancer detection in fibered confocal fluorescence microscopy images.  Artificial intelligence in medicine, 61(2):105–118, 2014.
  * [46] R Varaprasada Rao and T Jaya Chandra Prasad.  Content-based medical image retrieval using a novel hybrid scattering coefficients-bag of visual words-dwt relevance fusion.  Multimedia Tools and Applications, 80(8):11815–11841, 2021.
  * [47] Muhammad Nasim Kashif, Shan E Ahmed Raza, Korsuk Sirinukunwattana, Muhammmad Arif, and Nasir Rajpoot.  Handcrafted features with convolutional neural networks for detection of tumor cells in histology images.  In 2016 IEEE 13th International Symposium on Biomedical Imaging (ISBI), pages 1029–1032. IEEE, 2016.
  * [48] Juan Wang, Jiangshe Zhang, and Jie Zhao.  Texture classification using scattering statistical and cooccurrence features.  Mathematical Problems in Engineering, 2016, 2016.
  * [49] PM Siva Raja et al.  Brain tumor classification using a hybrid deep autoencoder with bayesian fuzzy clustering-based segmentation approach.  Biocybernetics and Biomedical Engineering, 40(1):440–453, 2020\.
  * [50] Hafeez Alani Agboola and Jesuloluwa Emmanuel Zaccheus.  Wavelet image scattering based glaucoma detection.  BMC Biomedical Engineering, 5(1):1, 2023.
  * [51] Thibaud Berthomier, Ali Mansour, Luc Bressollette, Frédéric Le Roy, and Dominique Mottier.  Venous blood clot structure characterization using scattering operator.  In 2016 2nd International Conference on Frontiers of Signal Processing (ICFSP), pages 73–80. IEEE, 2016.
  * [52] Zahra Baharlouei, Hossein Rabbani, and Gerlind Plonka.  Detection of retinal abnormalities in oct images using wavelet scattering network.  In 2022 44th Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC), pages 3862–3865. IEEE, 2022\.
  * [53] Qingyun Ren, Bingyin Zhou, Liang Tian, and Wei Guo.  Detection of covid-19 with ct images using hybrid complex shearlet scattering networks.  IEEE Journal of Biomedical and Health Informatics, 26(1):194–205, 2021.
  * [54] Vincent Andrearczyk, Valentin Oreiller, Julien Fageot, Xavier Montet, and Adrien Depersinge.  Solid spherical energy (sse) cnns for efficient 3d medical image analysis.  In Proceedings of the 21st Irish Machine Vision and Image Processing Conference (IMVIP 2019), number CONFERENCE. 28-30 August 2019, 2019\.
  * [55] A Ratna Raju, P Suresh, and R Rajeswara Rao.  Bayesian hcs-based multi-svnn: a classification approach for brain tumor segmentation and classification using bayesian fuzzy clustering.  Biocybernetics and Biomedical Engineering, 38(3):646–660, 2018\.
  * [56] Pejman Rasti, Ali Ahmad, Salma Samiei, Etienne Belin, and David Rousseau.  Supervised image classification by scattering transform with application to weed detection in culture crops of high density.  Remote Sensing, 11(3):249, 2019.
  * [57] Pedro MM Pereira, Lucas A Thomaz, Luis MN Tavora, Pedro AA Assuncao, Rui M Fonseca-Pinto, Rui Pedro Paiva, and Sergio MM de Faria.  Melanoma classification using light-fields with morlet scattering transform and cnn: Surface depth as a valuable tool to increase detection rate.  Medical Image Analysis, 75:102254, 2022.
  * [58] Korsuk Sirinukunwattana, David RJ Snead, and Nasir M Rajpoot.  A novel texture descriptor for detection of glandular structures in colon histology images.  In Medical Imaging 2015: Digital Pathology, volume 9420, pages 186–194. SPIE, 2015.
  * [59] Zeyu Deng, Lihui Wang, Zixiang Kuai, Qijian Chen, Xinyu Cheng, Feng Yang, Jie Yang, and Yuemin Zhu.  Cnn-based invertible wavelet scattering for the investigation of diffusion properties of the in vivo human heart in diffusion tensor imaging.  arXiv preprint arXiv:1912.07776, 2019.
  * [60] Debayan Bhattacharya, Finn Behrendt, Axelle Felicio-Briegel, Veronika Volgger, Dennis Eggert, Christian Betz, and Alexander Schlaefer.  Learning robust representation for laryngeal cancer classification in vocal folds from narrow band images.  In Medical Imaging with Deep Learning, 2022.
  * [61] Shankar Parmar and Chirag Paunwala.  A novel and efficient primitive stage detection of dyslexia in eeg signals using wavelet scattering transform.  Available at SSRN 4404532.
  * [62] Tameem Adel, Taco Cohen, Matthan Caan, Max Welling, AGEhIV study group, Alzheimer’s Disease Neuroimaging Initiative, et al.  3d scattering transforms for disease classification in neuroimaging.  NeuroImage: Clinical, 14:506–517, 2017.
  * [63] K Suhail and D Brindha.  A review on various methods for recognition of urine particles using digital microscopic images of urine sediments.  Biomedical Signal Processing and Control, 68:102806, 2021.
  * [64] SK Shadekul Islam, MD Abdullah Al Nasim, Ismail Hossain, Md Azim Ullah, Kishor Datta Gupta, and Md Monjur Hossain Bhuiyan.  Introduction of medical imaging modalities.  arXiv e-prints, pages arXiv–2306, 2023.
  * [65] Angona Biswas, MD Nasim, Md Shahin Ali, Ismail Hossain, Dr Md Azim Ullah, and Sajedul Talukder.  Active learning on medical image.  arXiv preprint arXiv:2306.01827, 2023.
  * [66] Angona Biswas, MD Nasim, Al Imran, Anika Tabassum Sejuty, Fabliha Fairooz, Sai Puppala, and Sajedul Talukder.  Generative adversarial networks for data augmentation.  arXiv preprint arXiv:2306.02019, 2023.
  * [67] MD Nasim, Aditi Dhali, Faria Afrin, Noshin Tasnim Zaman, and Nazmul Karim.  The prominence of artificial intelligence in covid-19.  arXiv preprint arXiv:2111.09537, 2021.
  * [68] Md Zihad Bin Jahangir, Ruksat Hossain, Riadul Islam, MD Nasim, Md Mahim Anjum Haque, Md Jahangir Alam, and Sajedul Talukder.  Introduction to medical imaging informatics.  arXiv preprint arXiv:2306.00421, 2023.
  * [69] Shuvra Sarker, Angona Biswas, MD Nasim, Md Shahin Ali, Sai Puppala, and Sajedul Talukder.  Case studies on x-ray imaging, mri and nuclear imaging.  arXiv preprint arXiv:2306.02055, 2023.


# D-Wave and Pair-Density-Wave Superconductivity in the Square-Lattice
$t\text{-}J$ Model

Feng Chen  Department of Physics and Astronomy, California State University,
Northridge, California 91330, USA  F. D. M. Haldane  Department of Physics,
Princeton University, Princeton NJ 08544, USA  D. N. Sheng  Department of
Physics and Astronomy, California State University, Northridge, California
91330, USA

###### Abstract

The Hubbard and closely related $t\text{-}J$ models are exciting platforms for
unconventional superconductivity. Through state-of-the-art density matrix
renormalization group calculations using the grand canonical ensembles, we
address open issues regarding the ground-state phase diagram of the extended
$t\text{-}J$ model on a square lattice in the parameter regime relevant to
cuprate superconductors. On large 8-leg cylinders, we demonstrate that the
pure $t\text{-}J$ model with only nearest-neighbor hoppings and superexchange
interactions, for a wide range of hole doping ($\delta=0.1-0.2$), hosts an
exotic SC state where the d-wave and unidirectional pair density wave SC co-
exist with spin and charge bond orders. Furthermore, a small next nearest
neighbor hopping $t_{2}$ suppresses the pair density wave and other co-
coexisting orders, leading to a d-wave SC phase in both electron- and hole-
doped cuprate model systems. We reveal a simple mechanism for SC in the
$t\text{-}J$ model, in which the nearest neighbor hopping plays an essential
role in driving the formation of Cooper pairs with real-space sign
oscillations balancing the competition between the kinetic and exchange
energies. Our work validates the $t\text{-}J$ model as a proper minimum model
for describing fundamental physics of cuprate superconductors.

Introduction.— The mechanism of high-temperature superconductivity (SC) in
cuprate compounds is one of the major challenges in condensed matter physics
[1, 2, 3]. The emergence of superconductivity from doping the
antiferromagnetic Mott insulating phase of cuprates suggests the square-
lattice Hubbard model in the strong coupling regime and the closely related
$t\text{-}J$ model as the minimal models for its theoretical understandings
[4, 5, 6, 2, 3, 7, 1, 8, 9, 10]. On the analytical side, resonating valence
bond (RVB) theory has presented a picture of developing unconventional SC from
spin liquid background [4, 11, 12, 13, 7, 1, 14]. New concepts such as spin-
charge separation, spin liquid, electron fractionalization, and intertwined
spin and charge orders have also been extensively explored and may play
important roles in understanding doped Mott insulators [4, 6, 8, 9, 15, 16, 1,
17, 18, 19, 20, 21]. However, various analytical approximations are not
controlled in identifying quantum phases of these strongly correlated systems.
On the other hand, extensive numerical simulations in recent years have often
identified the stripe phase in the ground states of the pure $t\text{-}J$ and
Hubbard models with only nearest-neighbor (NN) hoppings [22, 23, 24, 25, 26,
27, 28, 29, 30, 31, 32, 33, 34]. Meanwhile, there are also indications that
the stripe phase has local pairing tendency [32, 33, 35] and is in strong
competition with SC order [36, 37, 38, 39, 40]. Therefore, it is fundamentally
important to re-examine if the pure t-J model can host SC ground state, and
identify the physical reasons for its emergence or absence.

To provide a more realistic modeling of cuprate systems, the extended
$t\text{-}J$ and Hubbard models with both NN and next-nearest-neighbor (NNN)
hoppings ($t_{1},t_{2}$) have been extensively studied and the choice of the
positive or negative $t_{2}$ corresponds to electron-doped or hole-doped
cuprates with the Fermi surface topology consistent with experiments [32, 41,
42, 43]. Earlier numerical studies have found that positive $t_{2}$ enhances
$d$-wave SC order whereas negative $t_{2}$ promotes the stripe phase and
suppresses SC [44, 45, 46], or favors plaquette-type pairing in width-4
cylinders [26, 29, 47, 48] around the optimal hole doping. More recently,
density matrix renormalization group (DMRG) [49] studies of the extended
$t\text{-}J$ model have reached a consensus on the presence of robust d-wave
SC on wider cylinders with the width $L_{y}$ upto $L_{y}=8$ for parameter
regime $t_{2}/t_{1}\gtrsim 0.1$, marking a significant step towards
establishing d-wave SC for electron-doped cuprate model [50, 32, 33, 40] or
doping spin liquid (or spin-liquid-like states) [51, 52]. However, numerical
results on the hole-doped model ($t_{2}<0$) remain controversial. Previous
DMRG studies reported the stripe formation and absence of SC in the ground
state of the hole-doped $t\text{-}J$ model [32, 53, 33], whereas latest DMRG
simulations using larger bond dimensions have identified the melting of charge
density wave (CDW) and emergence of weak SC as the system widens from
$L_{y}=6$ to $8$ at $1/8$ hole doping [40]. It remains challenging to fully
understand the competition and interplay between CDW, magnetic order and SC,
and reveal a unified mechanism for the emergence of SC in the family of
$t\text{-}J$ models.

In this article, we address open issues regarding the emergence of
unconventional SC by studying the extended $t\text{-}J$ model on a square
lattice based on state-of-the-art DMRG simulations of grand canonical [32]
systems up to $L_{y}=8$. We establish a global phase diagram (Fig. 1(c)) where
the $d$-wave SC order is prevalent over the entire parameter regime with the
doping level $0.1\leq\delta\leq 0.2$ and hopping ratio $-0.2\leq
t_{2}/t_{1}\leq 0.3$. In particular, an exotic SC phase emerges in the pure
$t\text{-}J$ model ($t_{2}=0$) where isotropic d-wave SC coexists with
period-2 unidirectional pair density wave (PDW) [54], spin bond order and weak
charge order in real space. In the hole-doped $t\text{-}J$ model with negative
$t_{2}/t_{1}\sim-0.1$, a SC-CDW phase appears with co-existing d-wave SC and
CDW orders. A common uniform SC phase emerges at relatively larger magnitude
$t_{2}/t_{1}\sim\pm 0.2$, suggesting a symmetry between the electron- and
hole-doped $t\text{-}J$ models through suppressing competing orders. The
universal emergence of SC in the extended $t\text{-}J$ model is revealed to be
driven by the frustrating NN hopping ($t_{1}$) that turns the pre-existing RVB
pairs into coherent Cooper pairs upon doping. Our results validate the
$t\text{-}J$ model as the suitable minimal model for describing cuprate SC and
provide benchmark results for future theoretical studies.

Figure 1: (a) Illustration of the extended $t{\text{-}}J$ model. (b)
Illustration of the development of same-sublattice pairings (denoted by filled
ovals) from opposite-sublattice ones (denoted by unfilled ovals) in the pure
$t{\text{-}}J$ model as a consequence of NN electron hopping. (c) Ground-state
phase diagram under the variation of doping level $\delta$ and ratio
$t_{2}/t_{1}$. Three SC phases are found. The uniform d-wave superconducting
phase is denoted by SC and the superconducting phase coexisting with period-2
PDW or weak CDW is denoted by SC-PDW or SC-CDW, respectively. The white region
at the lower corners has strong competition between CDW and SC, which cannot
be fully resolved with the currently accessible bond dimensions.

Model and Method.— We consider the square-lattice
$t_{1}\text{-}t_{2}\text{-}J_{1}\text{-}J_{2}$ model defined as

$\begin{split}H=&-\sum_{\\{ij\\},\sigma}t_{ij}(\hat{c}^{\dagger}_{i,\sigma}\hat{c}_{j,\sigma}+{\rm
H.c.})\\\ &+\sum_{\\{ij\\}}J_{ij}(\hat{\bf S}_{i}\cdot\hat{\bf
S}_{j}-\frac{1}{4}\hat{n}_{i}\hat{n}_{j})-\mu\sum_{i}\hat{n}_{i},\end{split}$
(1)

where ${\hat{c}}_{i\sigma}^{\dagger}$ and ${\hat{c}}_{i\sigma}$ are the
creation and annihilation operators of the electron with spin $\sigma$
($\sigma=\pm 1/2$) at site $i=(x,y)$, ${\hat{\bf S}}_{i}$ is the spin-$1/2$
operator and
${\hat{n}}_{i}=\sum_{\sigma}{\hat{c}}_{i\sigma}^{\dagger}{\hat{c}}_{i\sigma}$
is the electron number operator. Double occupancy is prohibited. The electron
hoppings $t_{ij}$ and the spin exchange interactions $J_{ij}$ are restricted
to the NN ($t_{1}$, $J_{1}$) and NNN ($t_{2}$, $J_{2}$) terms. We consider a
cylinder that is periodic along the circumferential ($\hat{y}$) direction and
open along the axial ($\hat{x}$) direction with width $L_{y}$ and length
$L_{x}$ respectively. The total number of sites is $N=L_{x}\times L_{y}$. We
impose spin $SU(2)$ symmetry in the DMRG calculations [55, 50] and keep large
bond dimensions up to $M=15000$ (equivalent to about 45000 U(1) states) for
accurate results, which gives a truncation error smaller or around $3.0\times
10^{-5}$ (see Supplemental Material [56] for more details). We consider the
grand canonical ensemble where spontaneous symmetry breaking of the particle
number conservation is allowed [32], and tune the average electron number by
the chemical potential $\mu$. The average doping level is defined as
$\delta=1-\langle\sum_{i}\hat{n}_{i}\rangle/N$. Finally, we set $J_{1}=1$ as
the energy unit and $t_{1}=3.0$, corresponding to a Hubbard model at
$U/t_{1}=12$ [51, 50], and vary $t_{2}$ and $J_{2}$ together via the relation
$J_{2}/J_{1}=(t_{2}/t_{1})^{2}$.

Phase Diagram— We focus on the regime of doping level $0.1\leq\delta\leq 0.2$
and hopping ratio $-0.2\leq t_{2}/t_{1}\leq 0.3$ to obtain a phase diagram for
$L_{y}=8$ systems shown in Fig. 1(c). We identify three different SC phases:
(1) A uniform $d$-wave SC phase appears at both $t_{2}/t_{1}\gtrsim 0.15$ and
$t_{2}/t_{1}\lesssim-0.15$; (2) A SC phase coexisting with weak CDW order at
$-0.15\lesssim t_{2}/t_{1}\lesssim-0.05$; (3) A nematic SC phase composed of
uniform $d$-wave SC and a period-2 unidirectional PDW order with pairing
momentum $\bm{k}_{\text{PDW}}=(0,\pi)$ for $-0.05\lesssim t_{2}/t_{1}\lesssim
0.15-0.25$ depending on doping levels. In the bottom corners of the phase
diagram with $t_{2}/t_{1}\sim-0.2$ and doping around $\delta\sim 0.11$ and
$0.19$, there is also a possible CDW phase coexisting with weaker or vanishing
SC order. The important role of the frustrating NN hopping $t_{1}$ is
illustrated in Fig. 1(b), which allows the spin-singlet pairs in the same
sublattice to develop. In the following, we will explore universal features of
these SC phases and reveal the underlying mechanism for their emergence.

Figure 2: The SC pairing orders and electron density profiles for four
representative $t_{2}/t_{1}$ around 1/8 doping. (a) 2D plots for the NN
pairing orders $\Delta_{\alpha=x,y}(\bm{r})$ and hole density
$n_{h}(\bm{r})=1-\langle\hat{n}_{\bf r}\rangle$ obtained at $M=12000$. The
bond thickness is proportional to the pairing amplitude and blue/red color
denotes positive/negative sign respectively. The area of the green circle is
proportional to the hole density $n_{h}(\bm{r})$ at each site. (b) The middle
panels show the NN pairing orders $\Delta_{\alpha=x,y}$ along $x$-direction,
which are uniform along $\hat{y}$ for (b1, b2) and (b4) but oscillating along
$\hat{y}$ with a period of two lattice constants for $\Delta_{y}$ in (b3),
i.e. in the pure $t\text{-}J$ model. (c) The lower panels show the electron
density profiles $n(x,y)$ for $y=1,2$.

To identify the SC order, we examine the dominant spin-singlet NN pairing
orders defined as:

$\Delta_{\alpha}(x,y)=\langle\hat{c}_{(x,y),\uparrow}\hat{c}_{(x,y)+\hat{\alpha},\downarrow}-\hat{c}_{(x,y),\downarrow}\hat{c}_{(x,y)+\hat{\alpha},\uparrow}\rangle/\sqrt{2},$
(2)

where $\hat{\alpha}=\hat{x}$ or $\hat{y}$ denotes the direction of the NN
bond. In the grand canonical DMRG simulations [32, 56], the nonzero pairing
orders detect the spontaneous $U(1)$ symmetry breaking: a low-energy “natural
ground state” [57, 58] with symmetry breaking can be selected by DMRG
algorithm to reduce long-range entanglement and Hilbert space fragmentation,
so that the system effectively lowers its energy for finite bond dimensions.
The “natural ground state” is a superposition of lower energy eigenstates
(Anderson Tower of states) [57] for a system with the tendency of spontaneous
symmetry breaking, and becomes near degenerate with the true ground state for
large systems. Thus, the grand canonical simulation[32] for finite systems
provides a natural way to detect SC and competing orders.

We first demonstrate the global emergence of SC as illustrated in Fig. 2. For
four representative parameter points and a large system size $N=50\times 8$,
we show NN SC pairing orders $\Delta_{x}(\bm{r})$ and $\Delta_{y}(\bm{r})$
with a magnitude proportional to the bond width in Fig. 2 (a1-a4). The
opposite signs of the x- and y-bonds demonstrate a robust d-wave symmetry.
Furthermore, $\Delta_{x}$ and $\Delta_{y}$ are nearly constant in the bulk
regime for the uniform SC phase (Fig. 2 (b1,b4)), while there are noticeable
oscillations of $\Delta_{x}$ along $x-$direction in the SC-CDW and SC-PDW
phases (Fig. 2 (b2,b3)), indicating the coupling between CDW (Fig. 2 (c2,c3))
and the SC order $\Delta_{x}$ [54, 19]. In the SC-PDW phase of the pure
$t\text{-}J$ model, the most salient feature is the oscillation of
$\Delta_{y}$ with a period of 2 lattice spacings along $y-$direction as shown
in Fig. 2 (a3,b3) corresponding to a PDW with a pair momentum
$\bm{k}_{\text{PDW}}=(0,\pi)$. The electron density generally preserves
translational invariance along $y$-direction. The electron density $n(x,y)$
for $y=1,2$ is shown in Fig. 2 (c1-c4). We identify a very uniform electron
density $n(x,y)$ for the SC phase at $t_{2}/t_{1}=\pm 0.2$ (see Fig. 2 (c1,
c4)), and a very weak stripe pattern for the SC-CDW phase at
$t_{2}/t_{1}=-0.1$ along $x$ direction. Interestingly, for the pure
$t\text{-}J$ model, the electron density shows stripe-like patterns along
$x$-direction near the open boundaries (see Fig. 2 (c3)) and becomes flattened
in the center of the system. Full nature of the SC-PDW phase will be presented
below (see Fig. 4).

Overall, the SC order becomes stronger as $t_{2}$ gets more positive, and both
positive and negative $t_{2}$ of large enough magnitudes
($\absolutevalue{t_{2}/t_{1}}\gtrsim 0.15$) suppress CDW and PDW, resulting in
a uniform d-wave SC phase. The uniform SC phase at $t_{2}/t_{1}=0.2-0.3$
identified here is consistent with earlier results from either grand canonical
[32] or canonical DMRG studies [50, 40], further establishing the robust SC
phase in the electron-doped $t\text{-}J$ model. More significantly, we unify
the hole- and electron-doped $t\text{-}J$ models by establishing their common
uniform SC phase.

Uniform SC Phase and Competing Stripe Order in the Hole Doped $t\text{-}J$
Model— A CDW in the stripe form was identified either as the ground state [32]
or as a weak order co-existing with power-law SC correlations [40] for the
hole-doped $t\text{-}J$ model. Here we demonstrate how SC orders emerge
through suppressing the stripe order in large 8-leg systems. Fig. 3 shows the
charge densities, pairing orders $\Delta_{x,y}$ and pairing correlations
$P_{yy}(r)=\langle\hat{\Delta}_{y}(x_{0},y_{0})\hat{\Delta}_{y}(x_{0}+r,y_{0})\rangle$
for systems at $t_{2}/t_{1}=-0.2$ and $\mu=4.8$, corresponding to an average
doping $\delta\sim 0.13$. In Fig. 3(a1) for $N=32\times 8$, the CDW order is
quite robust until $M$ reaches the large value of $M=15000$, where the CDW
pattern is substantially weakened. The SC order $\Delta_{x}$ and $\Delta_{y}$
gradually develop from nearly vanishing to being notably strong when $M$ is
increased from 8000 to 15000 as shown in Fig. 3(b1). Along with the
suppression of CDW order at $M=15000$, we also see a drop of electron density
(Fig. 3(a1)) as more holes are attracted into the system at fixed $\mu$ due to
the much reduced (more negative) hopping energy (kinetic energy) [59]. Thus
the melting of the CDW and the emergence of the SC are driven by the
optimization of the kinetic energy [60]. Consistently, the pairing correlation
$P_{yy}(r)$ in Fig. 3(c1) also increases with $M$ and exhibits a robust quasi-
long-range SC order at $M=15000$. Physically speaking, the CDW order is more
classical and is the optimal low-energy state when the bond dimension is
relatively small. Upon increasing the bond dimension, the more entangled state
such as a coherent SC state is developed to effectively lower the total energy
by balancing the kinetic and potential energies.

In a longer system with $N=50\times 8$, the SC orders $\Delta_{x}$ and
$\Delta_{y}$ in Fig. 3(b2) are already quite robust at $M=8000$ and are nearly
independent of the bond dimension for $M=8000-15000$, indicating a more robust
SC phase for larger systems [56]. The SC order is weaker at the boundary but
grows rapidly into the bulk (Fig. 3(b2)), whereas the charge modulations only
appear around the system boundaries (Fig. 3(a2)). The uniform SC order and
electron density in the bulk mark a uniform SC phase. Meanwhile, the pairing
correlation in Fig. 3(c2) exhibits a power-law decay with a small power
exponent[56], further supporting a robust SC phase.

Figure 3: The SC phase for the hole-doped $t\text{-}J$ model with
$t_{2}/t_{1}=-0.2$ and $\mu=4.8$. (a) Electron density profiles, (b) pairing
orders and (c) $y$-bond pairing correlations $P_{yy}(r)$ at different bond
dimensions for system sizes $N=32\times 8$ (left panels) and $N=50\times 8$
(right panels). The pairing correlations in (c2) are extrapolated to infinite
$M$ through second-order polynomial extrapolation, and the extrapolated data
in the bulk can be fitted by a power-law decay with an exponent of $K_{sc}\sim
0.06$.

Exotic Co-Existing D-wave, Pair Density Wave Superconductivity, and Spin Bond
Order for the Pure $t{\text{-}}J$ Model.— Now we focus on the pure
$t\text{-}J$ model with $t_{2}=J_{2}=0$. The pairing order $\Delta_{y}$ has a
double unit cell along $y$-direction, and its average
$\bar{\Delta}_{y}(x)=[\Delta_{y}(x,y=2)+\Delta_{y}(x,y=1)]/2$ has similar
magnitude as $\Delta_{x}$ (Fig. 2 (b3)), forming the d-wave SC component
$\bar{\Delta}_{y}(x)\sim 0.02$ after extrapolation to the infinite $M$ limit
as shown in Fig. 4(a). The PDW order of the SC-PDW phase is obtained by
subtracting the d-wave background from the local SC order:
$\Delta_{y}^{\text{pdw}}(x,y)=[\Delta_{y}(x,y)-\bar{\Delta}_{y}(x)]$, which
has the following y-dependence:

$\Delta_{y}^{\text{pdw}}(x,y)=(-1)^{y}\Delta_{y}^{\text{P}}(x),$ (3)

indicating a unidirectional PDW with Cooper-pair momentum
$\bm{k}_{\text{PDW}}=(0,\pi)$. The magnitude of the PDW order
$\Delta_{y}^{\text{P}}(x)$ drops slightly with the increase of $M$ and is
extrapolated to $\Delta_{y}^{\text{P}}\approx 0.004$ as shown in Fig. 4(b). We
illustrate the pattern of the PDW order $\Delta_{y}^{\text{pdw}}(x,y)$ in real
space in Fig. 4(c), which is coupled with a spin-bond order
$\langle\hat{\bm{S}}_{x,y}\cdot\hat{\bm{S}}_{x,y+1}\rangle$ with the same
strong-weak unit-2 pattern as shown in Fig. 4(d). The SC-PDW phase also shows
a weaker charge bond order[56]. To the best of our knowledge, such a PDW order
has only been identified in more complex models, e.g. the three-band cuprate
models [61, 62, 63], although it is highly desired and explored in
phenomenological and mean-field studies of cuprate systems [25, 64, 65, 54].

Figure 4: The coexisting d-wave SC, PDW and spin bond orders in the SC-PDW
phase for $t_{2}=0,\mu=5.0$ and $N=50\times 8$. (a) The average $y$-bond SC
pairing of the first and second legs
$\bar{\Delta}_{y}(x)\equiv\left[\Delta_{y}(x,y=2)+\Delta_{y}(x,y=1)\right]/2$.
The inset shows the second-polynomial extrapolation of the bulk
$\bar{\Delta}_{y}$ (averaged over middle 10 sites) versus $1/M$:
$\bar{\Delta}_{y}(M=\infty)=0.02$. (b) The magnitude of PDW
$\Delta^{\text{P}}_{y}(x)$ defined in Eq. 3. The bulk value after
extrapolation is $\Delta_{y}^{\text{P}}(M=\infty)=0.004$. (c) Illustration of
the unidirectional PDW with momentum $\bm{k}_{\text{PDW}}=(0,\pi)$. The $x$\-
and $y$-bond mean values $\bar{\Delta}_{x}=-0.0223$ and
$\bar{\Delta}_{y}=0.0233$ are subtracted from $\Delta_{x}$ and $\Delta_{y}$
respectively. (d) Nearest-neighbor spin bond orders
$-\langle\hat{\bm{S}}_{i}\cdot\hat{\bm{S}}_{\bm{i+e_{\alpha}}}\rangle$ in the
central region, subtracting the anisotropic $x$\- and $y$-bond mean values
$\bar{S}^{b}_{x}=0.1852$ and $\bar{S}^{b}_{y}=0.2215$. The red/blue color
represents positive/negative value respectively and the bond widths are
proportional to the magnitudes of the quantity represented. A reference bond
and its magnitude is put beside each lattice plot. (c) and (d) are obtained
under $M=12000$.

Longer-Distance Cooper Pairs— We now explore the nature of the SC phases
through imaging the real space distribution of the Cooper pairs. We define the
longer-distance pairing orders
$\tilde{\Delta}(\bm{r})\equiv\tilde{\Delta}(x,y)=\langle\hat{c}_{(x_{0},y_{0}),\uparrow}\hat{c}_{(x_{0}+x,y_{0}+y),\downarrow}-\hat{c}_{(x_{0},y_{0}),\downarrow}\hat{c}_{(x_{0}+x,y_{0}+y),\uparrow}\rangle/\sqrt{2}$
between two points at positions $(x_{0},y_{0})$ and $(x_{0}+x,y_{0}+y)$, which
are shown in Fig. 5 as circles whose areas are proportional to their
magnitudes and whose colors denote their signs. While the NN pairing orders
follow a d-wave symmetry for all SC phases, an additional small $s$-wave
component can be identified by the same sign and magnitude of the four NNN
pairings at $t_{2}/t_{1}=\pm 0.2$ as shown in Fig. 5(a,c). The longer-distance
pairing order shows sign oscillations in real space and decays with the
pairing separation. The length scale for its sign to repeat along
$\hat{x}-$direction is identified to be $l_{s}=5,6,3$, for
$t_{2}/t_{1}=-0.2,0,0.2$, respectively. Furthermore, the Cooper pairs at
$t_{2}/t_{1}=-0.2$ and 0 are extended into longer distances than that at
$t_{2}/t_{1}=0.2$.

A common feature of all these SC states is the development of pairings between
the same sublattices at intermediate distances comparable to $l_{s}$,
indicating that such Cooper pairs gain kinetic energy through NN hoppings.
Starting from the pure $t\text{-}J$ model, the pre-existing long-range RVB
pairs at half-filling are formed between opposite sublattices [13], following
the Marshall sign rule [66] and supporting a long-range AFM order. Upon hole
doping, as illustrated in Fig. 1(b), the NN electron hopping will bring RVB
pairs between opposite sublattices into the same sublattice, breaking the
Marshall sign rule and frustrating the AFM order. Therefore, the NN electron
hopping is a key driving force in turning the pre-existing RVB pairs due to
the NN spin exchange interactions into coherent SC Cooper pairs [9]. For the
pure $t\text{-}J$ model, the longer-distance Cooper pairs are developed
through breaking the d-wave symmetry of the NN pair order by inducing an
additional unidirectional PDW order intertwining with spin and charge bond
orders. After pinning down the SC for the pure $t\text{-}J$ model, it becomes
clear what is the role of NNN hoppings: they help longer-distance Cooper pairs
to develop in more symmetric ways and suppress the nematicity or PDW tendency
of the pure $t\text{-}J$ model. Furthermore, the distinct signs of the NNN
hoppings for hole- and electron-doped $t\text{-}J$ models results in different
length scale $l_{s}$ for sign oscillation of $\tilde{\Delta}(\bm{r})$ and
larger $l_{s}$ on the hole-doped side makes the SC harder to emerge until the
system reaches $L_{y}=8$ [40, 33]. For the electron-doped $t\text{-}J$ model,
the SC orders also grow stronger as the system widens from $L_{y}=6$ to $8$
[56]. Our results demonstrate a unified picture for understanding the global
phase diagram of the $t\text{-}J$ model, driven by NN hopping in additional to
the spin exchange interaction. Their cooperative interplay [9] is the key
mechanism for the general emergence of SC around optimal doping in these
systems.

Figure 5: Long-distance pairings $\tilde{\Delta}(\bm{r}-\bm{r_{0}})$ imaging
real space Cooper pairs. The reference point $(x_{0},y_{0})=(25,5)$ for (a-c)
is marked by a green cross. The areas of the circles are proportional to the
magnitudes of $\tilde{\Delta}(\bm{r}-\bm{r_{0}})$ and the red/green color
denotes its positive/negative sign. The associated values of nearest-neighbor
pairings are given in the figures. Data are obtained under $M=15000$. Figure
6: Cooper pair and spin/charge correlation lengths for the SC-PDW phase at
$t_{2}/t_{1}=0$ and $\mu=5.0$. (a) Non-local pairing $\tilde{\Delta}(r)$
versus pair-distance $r$ along $x$-direction for different bond dimension $M$s
and its extrapolation to infinite $M$, giving an exponential decay with a
correlation length $\xi_{P}=3.74$ representing the average Cooper-pair size.
(b-d) are similar plots for single-particle Green’s function, spin and charge-
density correlations, with fitted correlation lengths $\xi_{G}=3.02$,
$\xi_{S}=2.16$ and $\xi_{D}=1.72$.

Exponentially Decaying Spin and Charge Correlations— Now we compare length
scales of the Cooper pair wave function, charge and spin correlations for the
pure $t\text{-}J$ model. As shown in Fig. 6(a), the longer-distance pairing
order $\tilde{\Delta}(\bm{r})$ decreases exponentially with the distance $r$
between two paired electrons, characterized by an effective Cooper pair
correlation length or size $\xi_{P}=3.74$. Similar small length scales
$\xi_{P}=2.44$ and $2.01$ are also found for $t_{2}/t_{1}=\mp 0.2$,
respectively [56]. The relative magnitudes of $\xi_{P}$ for different
$t_{2}/t_{1}$ are consistent with those of $l_{s}$, showing the largest Cooper
pair size at $t_{2}=0$. Meanwhile, the single-particle Green’s function
$G(r)=\langle\sum_{\sigma}\hat{c}^{\dagger}_{\sigma}(x,y)\hat{c}_{\sigma}(x+r,y)\rangle$,
spin correlation
$S(r)=\langle{\hat{\mathbf{S}}}(x,y)\cdot\hat{\mathbf{S}}(x+r,y)\rangle$ and
charge density correlation
$D(r)=\langle{\hat{n}}(x,y){\hat{n}}(x+r,y)\rangle-\langle{\hat{n}}(x,y)\rangle\langle{\hat{n}}(x+r,y)\rangle$
can all be described by exponential decays with very short correlation lengths
comparable to the Cooper pair size $\xi_{P}$ (see Fig. 6(b-d)). The short-
range spin and charge correlations are generally observed for all SC phases
(see SM [56] Fig. S4), similar to the previous identified SC phase by doping
the spin liquid or spin-liquid like states [51, 52]. These results indicate a
simple kinetic-energy-driven mechanism for the emergence of spin liquid
background: the NN hopping frustrates the AFM order and turns the spin
background into a short-range spin liquid.

Summary and Outlook— Through state-of-the-art DMRG simulations combining grand
canonical ensembles [32] with a rigorous finite bond dimension analysis, we
establish a global quantum phase diagram and unveil the universal emergence of
SC in the general $t\text{-}J$ model around the optimal hole doping
($0.1\leq\delta\leq 0.2$). Of particular interest is the pure $t\text{-}J$
model, whose ground state is determined to be an exotic phase with coexisting
d-wave SC, unidirectional PDW, spin bond order and weaker charge bond order.
Our results support the pure $t\text{-}J$ model as a minimal model for
understanding unconventional SC in doped Mott insulators, a conclusion of
great importance as a d-wave SC ground state has never been observed before in
unbiased numerical simulations of the pure $t\text{-}J$ model [32, 40]. Beyond
the pure $t\text{-}J$ model, we establish the uniform SC phase in both hole-
and electron-doped regimes, where SC grows stronger or only emerges for wider
systems ($L_{y}=8$), suggesting SC in the 2D limit. The necessity of large
bond dimensions and system sizes for the hole-doped SC phase to emerge is due
to the stronger competition between kinetic and exchange energies in the
presence of additional Berry phase of the negative NNN hopping [33], which
leads to spatially larger Cooper pairs with larger length scale of sign
oscillations (Fig. 5). A key driving force for the SC in the $t\text{-}J$
model, besides the NN spin exchange interaction [4], is identified here as the
frustrating NN hopping that transforms the pre-existing long-range RVB pairs
in the AFM spin background into finite-range phase-coherent SC pairs upon
increasing doping.

Given the observation of prevalent SC in the extended $t\text{-}J$ model in
this work, it is natural to ask if one can also establish a unified mechanism
for the unconventional SC in the one-band Hubbard model [29, 30, 31, 67, 34,
68, 69] for cuprates. Ground-state SC has been identified in the Hubbard model
with NNN hopping using both DMRG and auxiliary field quantum Monte Carlo
methods in the grand canonical finite-size systems in the presence of a global
pairing field [68]. While their results demonstrate the emergence of SC in
both electron- and hole-doped regimes, consistent with our results, the
detailed nature of Cooper pairs is not known there, and thus we believe our
work will stimulate more extensive exploration of the extended Hubbard model
to identify a unified mechanism of unconventional SC. Furthermore, it would be
interesting to reexamine the pure Hubbard model ($t_{2}=0$) to understand if
the absence of SC in numerical studies presented so far [31, 68, 34] is
related to finite-size and/or low entanglement effects (favoring conventional
charge or spin orders), which will be addressed in the future.

We thank Leon Balents, S. S. Gong, Xin Lu, H. C. Jiang, Shengtao Jiang, S. A.
Kivelson, S. R. White, Z. Y. Weng, and Wei Zhu for stimulating discussions.
This work was supported by the U.S. Department of Energy, Office of Basic
Energy Sciences under Grant No. DE-FG02-06ER46305 (FC, DNS) for computational
study, and National Science Foundation (NSF) Princeton Center for Complex
Materials, a Materials Research Science and Engineering Center DMR-2011750
(FDMH). D.N.S. also acknowledges partial support from NSF Partnership in
Research and Education in Materials DMR-1828019 for her travel to Princeton.

## References

  * Lee  _et al._ [2006] P. A. Lee, N. Nagaosa, and X.-G. Wen, Rev. Mod. Phys. 78, 17 (2006), publisher: American Physical Society.
  * Keimer  _et al._ [2015] B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, Nature 518, 179 (2015).
  * Proust and Taillefer [2019] C. Proust and L. Taillefer, Annual Review of Condensed Matter Physics 10, 409 (2019).
  * Anderson [1987] P. W. Anderson, Science 235, 1196 (1987), publisher: American Association for the Advancement of Science.
  * Zhang and Rice [1988] F. C. Zhang and T. M. Rice, Phys. Rev. B 37, 3759 (1988).
  * Rokhsar and Kivelson [1988] D. S. Rokhsar and S. A. Kivelson, Phys. Rev. Lett. 61, 2376 (1988).
  * Anderson  _et al._ [2004] P. W. Anderson, P. A. Lee, M. Randeria, T. M. Rice, N. Trivedi, and F. C. Zhang, Journal of Physics: Condensed Matter 16, R755 (2004).
  * Sachdev [2003] S. Sachdev, Rev. Mod. Phys. 75, 913 (2003).
  * Weng  _et al._ [1997] Z. Y. Weng, D. N. Sheng, Y.-C. Chen, and C. S. Ting, Phys. Rev. B 55, 3894 (1997), publisher: American Physical Society.
  * Ogata and Fukuyama [2008] M. Ogata and H. Fukuyama, Rep. Prog. Phys. 71, 036501 (2008).
  * Anderson  _et al._ [1987] P. W. Anderson, G. Baskaran, Z. Zou, and T. Hsu, Phys. Rev. Lett. 58, 2790 (1987), publisher: American Physical Society.
  * Kivelson  _et al._ [1987] S. A. Kivelson, D. S. Rokhsar, and J. P. Sethna, Phys. Rev. B 35, 8865 (1987), publisher: American Physical Society.
  * Liang  _et al._ [1988] S. Liang, B. Doucot, and P. W. Anderson, Phys. Rev. Lett. 61, 365 (1988).
  * Scalapino [2012] D. J. Scalapino, Rev. Mod. Phys. 84, 1383 (2012), publisher: American Physical Society.
  * Wen [1991] X. G. Wen, Phys. Rev. B 44, 2664 (1991).
  * Senthil and Fisher [2000] T. Senthil and M. P. A. Fisher, Phys. Rev. B 62, 7850 (2000).
  * Kaul  _et al._ [2008] R. K. Kaul, Y. B. Kim, S. Sachdev, and T. Senthil, Nature Physics 4, 28 (2008).
  * Balents and Sachdev [2007] L. Balents and S. Sachdev, Annals of Physics 322, 2635–2664 (2007).
  * Fradkin  _et al._ [2015] E. Fradkin, S. A. Kivelson, and J. M. Tranquada, Rev. Mod. Phys. 87, 457 (2015), publisher: American Physical Society.
  * Weng [2011] Z.-Y. Weng, New J. Phys. 13, 103039 (2011), publisher: IOP Publishing.
  * Jiang and Kivelson [2022] H.-C. Jiang and S. A. Kivelson, Proceedings of the National Academy of Sciences 119, e2109406119 (2022), publisher: Proceedings of the National Academy of Sciences.
  * White and Scalapino [1998] S. R. White and D. J. Scalapino, Phys. Rev. Lett. 80, 1272 (1998), publisher: American Physical Society.
  * White and Scalapino [2003] S. R. White and D. J. Scalapino, Phys. Rev. Lett. 91, 136403 (2003).
  * Hager  _et al._ [2005] G. Hager, G. Wellein, E. Jeckelmann, and H. Fehske, Phys. Rev. B 71, 075108 (2005).
  * Berg  _et al._ [2009] E. Berg, E. Fradkin, S. A. Kivelson, and J. M. Tranquada, New Journal of Physics 11, 115004 (2009).
  * Dodaro  _et al._ [2017] J. F. Dodaro, H.-C. Jiang, and S. A. Kivelson, Phys. Rev. B 95, 155116 (2017), publisher: American Physical Society.
  * Zheng  _et al._ [2017] B.-X. Zheng, C.-M. Chung, P. Corboz, G. Ehlers, M.-P. Qin, R. M. Noack, H. Shi, S. R. White, S. Zhang, and G. K.-L. Chan, Science 358, 1155 (2017).
  * Ido  _et al._ [2018] K. Ido, T. Ohgoe, and M. Imada, Phys. Rev. B 97, 045138 (2018).
  * Jiang and Devereaux [2019] H.-C. Jiang and T. P. Devereaux, Science 365, 1424 (2019).
  * Ponsioen  _et al._ [2019] B. Ponsioen, S. S. Chung, and P. Corboz, Phys. Rev. B 100, 195141 (2019).
  * Simons Collaboration on the Many-Electron Problem  _et al._ [2020] Simons Collaboration on the Many-Electron Problem, M. Qin, C.-M. Chung, H. Shi, E. Vitali, C. Hubig, U. Schollwöck, S. R. White, and S. Zhang, Phys. Rev. X 10, 031016 (2020), publisher: American Physical Society.
  * Jiang  _et al._ [2021] S. Jiang, D. J. Scalapino, and S. R. White, Proceedings of the National Academy of Sciences 118, e2109978118 (2021), publisher: Proceedings of the National Academy of Sciences.
  * Lu  _et al._ [2023a] X. Lu, J.-X. Zhang, S.-S. Gong, D. N. Sheng, and Z.-Y. Weng, (2023a), arXiv:2303.13498 [cond-mat.str-el] .
  * Jiang  _et al._ [2023a] Y.-F. Jiang, T. P. Devereaux, and H.-C. Jiang,  (2023a), arXiv:2303.15541 [cond-mat].
  * Wietek [2022] A. Wietek, Phys. Rev. Lett. 129, 177001 (2022).
  * Corboz  _et al._ [2014] P. Corboz, T. Rice, and M. Troyer, Phys. Rev. Lett. 113, 046402 (2014), publisher: American Physical Society.
  * Raczkowski  _et al._ [2007] M. Raczkowski, M. Capello, D. Poilblanc, R. Frésard, and A. M. Oleś, Phys. Rev. B 76, 140505 (2007).
  * Yang  _et al._ [2009] K.-Y. Yang, W. Q. Chen, T. M. Rice, M. Sigrist, and F.-C. Zhang, New J. Phys. 11, 055053 (2009).
  * Li  _et al._ [2021] J.-W. Li, B. Bruognolo, A. Weichselbaum, and J. von Delft, Phys. Rev. B 103, 075127 (2021), publisher: American Physical Society.
  * Lu  _et al._ [2023b] X. Lu, F. Chen, W. Zhu, D. N. Sheng, and S.-S. Gong, (2023b), arXiv:2304.03963 [cond-mat.str-el] .
  * Damascelli  _et al._ [2003] A. Damascelli, Z. Hussain, and Z.-X. Shen, Rev. Mod. Phys. 75, 473 (2003).
  * Nazarenko  _et al._ [1995] A. Nazarenko, K. J. E. Vos, S. Haas, E. Dagotto, and R. J. Gooding, Phys. Rev. B 51, 8676 (1995).
  * Kim  _et al._ [1998] C. Kim, P. J. White, Z.-X. Shen, T. Tohyama, Y. Shibata, S. Maekawa, B. O. Wells, Y. J. Kim, R. J. Birgeneau, and M. A. Kastner, Phys. Rev. Lett. 80, 4245 (1998).
  * White and Scalapino [1999] S. R. White and D. J. Scalapino, Phys. Rev. B 60, R753 (1999), publisher: American Physical Society.
  * Martins  _et al._ [2001] G. B. Martins, J. C. Xavier, L. Arrachea, and E. Dagotto, Phys. Rev. B 64, 180513 (2001), publisher: American Physical Society.
  * Himeda  _et al._ [2002] A. Himeda, T. Kato, and M. Ogata, Phys. Rev. Lett. 88, 117001 (2002), publisher: American Physical Society.
  * The Simons Collaboration on the Many-Electron Problem  _et al._ [2020] The Simons Collaboration on the Many-Electron Problem, C.-M. Chung, M. Qin, S. Zhang, U. Schollwöck, and S. R. White, Phys. Rev. B 102, 041106 (2020), publisher: American Physical Society.
  * Jiang  _et al._ [2020] Y.-F. Jiang, J. Zaanen, T. P. Devereaux, and H.-C. Jiang, Phys. Rev. Res. 2, 033073 (2020), publisher: American Physical Society.
  * White [1992] S. R. White, Phys. Rev. Lett. 69, 2863 (1992), publisher: American Physical Society.
  * Gong  _et al._ [2021] S. Gong, W. Zhu, and D. Sheng, Phys. Rev. Lett. 127, 097003 (2021), publisher: American Physical Society.
  * Jiang and Kivelson [2021] H.-C. Jiang and S. A. Kivelson, Phys. Rev. Lett. 127, 097002 (2021), publisher: American Physical Society.
  * Jiang  _et al._ [2023b] H.-C. Jiang, S. A. Kivelson, and D.-H. Lee, Phys. Rev. B 108, 054505 (2023b), publisher: American Physical Society.
  * Jiang  _et al._ [2022] S. Jiang, D. J. Scalapino, and S. R. White, Phys. Rev. B 106, 174507 (2022).
  * Agterberg  _et al._ [2020] D. F. Agterberg, J. S. Davis, S. D. Edkins, E. Fradkin, D. J. Van Harlingen, S. A. Kivelson, P. A. Lee, L. Radzihovsky, J. M. Tranquada, and Y. Wang, Annual Review of Condensed Matter Physics 11, 231 (2020), _eprint: https://doi.org/10.1146/annurev-conmatphys-031119-050711.
  * McCulloch and Gulácsi [2002] I. P. McCulloch and M. Gulácsi, EPL 57, 852 (2002), publisher: IOP Publishing.
  * [56] See Supplementary Materials for more supporting data.
  * Tasaki [2019] H. Tasaki, J Stat Phys 174, 735 (2019).
  * Jiang  _et al._ [2012] H.-C. Jiang, Z. Wang, and L. Balents, Nature Physics 8, 902 (2012).
  * [59] The chemical potential term $-\mu\sum_{i}\hat{n}_{i}$ and the superexchange term $\sum_{\\{ij\\}}J_{ij}(\hat{\bf S}_{i}\cdot\hat{\bf S}_{j}-\frac{1}{4}\hat{n}_{i}\hat{n}_{j})$ both increase as the electron occupation is lowered.
  * Zhao  _et al._ [2022] J.-Y. Zhao, S. A. Chen, H.-K. Zhang, and Z.-Y. Weng, Phys. Rev. X 12, 011062 (2022).
  * Jiang and Devereaux [2023] H.-C. Jiang and T. P. Devereaux, (2023), arXiv:2309.11786 [cond-mat.str-el] .
  * Jiang [2023] H.-C. Jiang, Phys. Rev. B 107, 214504 (2023).
  * Ponsioen  _et al._ [2023] B. Ponsioen, S. S. Chung, and P. Corboz, (2023), arXiv:2306.12910 [cond-mat.str-el] .
  * Choubey  _et al._ [2020] P. Choubey, S. H. Joo, K. Fujita, Z. Du, S. D. Edkins, M. H. Hamidian, H. Eisaki, S. Uchida, A. P. Mackenzie, J. Lee, J. C. S. Davis, and P. J. Hirschfeld, Proceedings of the National Academy of Sciences 117, 14805 (2020), publisher: Proceedings of the National Academy of Sciences.
  * Tu and Lee [2019] W.-L. Tu and T.-K. Lee, Sci Rep 9, 1719 (2019), number: 1 Publisher: Nature Publishing Group.
  * Marshall [1955] W. Marshall, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 232, 48 (1955).
  * Qin  _et al._ [2022] M. Qin, T. Schäfer, S. Andergassen, P. Corboz, and E. Gull, Annual Review of Condensed Matter Physics 13, 275 (2022), _eprint: https://doi.org/10.1146/annurev-conmatphys-090921-033948.
  * Xu  _et al._ [2023] H. Xu, C.-M. Chung, M. Qin, U. Schollwöck, S. R. White, and S. Zhang, (2023), arXiv:2303.08376 [cond-mat.supr-con] .
  * Zhang  _et al._ [2023] C. Zhang, J.-W. Li, and J. von Delft,  (2023), arXiv:2307.14835 [cond-mat].
  * Fishman  _et al._ [2022] M. Fishman, S. R. White, and E. M. Stoudenmire, SciPost Phys. Codebases , 4 (2022).

Supplemental Materials

In the Supplemental Materials, we provide additional numerical results to
support the conclusions made in the main text. In Sec. A we show evolution of
pairing orders and pair-pair correlation functions with the bond dimension for
six- and eight-leg cylinders by canonical and grand canonical simulations with
spin $SU(2)$ or $U(1)$ symmetries, establishing the reliability of our
simulations. In Sec. B, we provide an example of the finite bond dimension
scaling of the ground state energy to demonstrate the numerical convergence of
DMRG results. Furthermore, we demonstrate the equivalence of spin $SU(2)$ and
$U(1)$ simulations in capturing long-range SC pairing orders. In Sec. C, we
show pairing correlations of the SC-PDW phase and demonstrate the weak charge-
bond order. In Sec. D, we show the exponentially decaying long-range pairing
orders, single-particle Green’s function, spin and density correlations for
$t_{2}/t_{1}=\pm 0.2$.

##  Appendix A Evolution of SC order and pairing correlations in DMRG
Simulations

As discussed in the main text, the DMRG simulations allow spontaneous symmetry
breaking [32], through targeting the “natural ground state” as a superposition
of low-energy states (Anderson Tower of states) [57], which becomes near
degenerate with the true ground state for large systems. In this section, we
provide systematical comparisons between canonical and grand canonical DMRG
results with spin $SU(2)$ or $U(1)$ symmetries. In Fig. S1(a), we show how the
SC order evolves with the bond dimension $M$ (or $m$) in the grand canonical
finite-size $SU(2)$ (or $U(1)$) DMRG simulations at a fixed chemical potential
$\mu=5.4$ and $t_{2}/t_{1}=0.2$. ITensor library [70] has been used for the
$U(1$) simulations. On the six-leg systems, the SC order sharply increases
between $m=1000$ and 3000. After reaching a maximum value, $\Delta_{y}$
remains strong with the increase of the bond dimension, and such a behavior
becomes more robust for a larger system with $N=50\times 8$, as demonstrated
by the larger maximum value of $\Delta_{y}$ and slower decay with $M$. In the
whole range of bond dimensions, the obtained $\Delta_{y}$ and ground state
energy for $SU(2)$ and $U(1)$ DMRGs agree with each other as long as we take
$M\sim m/3$, demonstrating the higher efficiency of SU(2) simulations. The
development of SC order, its plateau behaviour and increased stability with
the growing of system sizes reveal the robustness of the SC state in the
thermodynamic limit. This is similar to the detection of the Néel order in the
thermodynamic limit [32] by the finite-size spin-$U(1)$ DMRG simulations of
the antiferromagnetic square-lattice Heisenberg model.

We further compare in Fig. S1(b, c) the pairing correlation obtained by both
grand canonical and canonical finite-size simulations, respectively. For
smaller $M$, the grand canonical simulation will first identify mean-field-
like bulk orders, which correspond to a nearly flat pairing correlations at
long distances. As $M$ increases, the system can take into account more and
more quantum fluctuations around the mean-field orders, such that the pairing
correlations scales down with $M$ and established a power-law decay with an
exponent around 0.62, indicating a quasi-long-range SC order for the 6-leg
ladder systems. For the canonical DMRG calculation using a fixed particle
number (Fig. S1(c)), the SC pairing correlations at longer distances are being
suppressed for small $M$, and larger $M$ is required for their development.
Thus the pairing correlations scale up with $M$, which also reach a very
similar power-law decay with an exponent around 0.55 by extrapolation to
infinite $M$ limit. The similar exponents obtained by both grand canonical and
canonical DMRGs demonstrate the equivalence of the two methods and the good
convergence reached for $L_{y}=6$ systems. For $L_{y}=8$ systems, we obtain
reliable results for bulk order and qualitatively accurate behavior for
pairing correlations. However quantitatively, it is very difficult to obtain
accurate long-distance correlations with finite bond dimensions accessible by
current computational power, and much larger bond dimensions are needed to
obtain accurate power-law exponent specific for the quasi-1D 8-leg systems.
Thus for a large system, grand canonical DMRG is a proper way to detect ground
state SC orders, approximating physics in the thermodynamic limit as long as
such orders are robust with the increase of the bond dimensions and system
sizes.

Figure S1: (a) $t_{2}/t_{1}=0.2$ and $\mu=5.4$. The change of bulk pairing
order parameters with the bond dimensions. The $SU(2)$ bond dimension $M$ is
counting the number of spin $SU(2)$ multiplets, which corresponds to about 3
times of the $U(1)$ bond dimension $m$. Both six-leg and eight-leg systems are
studied, and the latter show larger SC orders and more robust plateau for
$\Delta_{y}$ with respect to $M$, indicating SC in the thermodynamic limit.
(b) Pairing correlations $P_{yy}(r)$ at different $M$s obtained by spin
$SU(2)$ DMRG for grand canonical systems at fixed $\mu=5.4$ for
$t_{2}/t_{1}=0.2$ and $N=32\times 6$. A scaling to infinite $M$ renders a
power-law decay with an exponent around 0.62. (c) $P_{yy}(r)$ of a canonical
system at doping $\delta=1/8$ and $N=32\times 6$, with the power-law exponent
around 0.55 for infinite-$M$ extrapolated data.

##  Appendix B DMRG convergence

Most of the results presented in the main text are obtained for a system size
of $N=50\times 8$, and we keep bond dimension upto $M=15000$, which gives a
truncation error smaller or around $3.0\times 10^{-5}$. We show in Fig. S2(a)
and (b) the scaling of the ground state energies in terms of the bond
dimension for both $L_{y}=6$ and $8$, which gives the energy density
$E/N=-6.590$ and -6.593, respectively. The closeness between energy densities
at the largest bond dimension and those extrapolated to infinite $M$ for both
systems demonstrate a good convergence of the ground state.

In Fig. S2(c) and (d), we compare the non-local pairing orders on a six-leg
system by both $SU(2)$ and $U(1)$ simulations with $m=M=10000$, and find
quantitative agreement between them, indicating a robust structure of the
Cooper pair insensitive to the bond dimension or spin symmetry imposed in our
simulations.

Figure S2: $t_{2}/t_{1}=0.2,\mu=5.4$. (a) Ground-state energy extrapolation to
infinite $M$ for $N=32\times 6$ with extrapolated energy
$E(M\to\infty)/N=-6.590$. (b) Similar to (a) for $N=50\times 8$ with
$E(M\to\infty)/N=-6.593$. (c) and (d) are non-local pairings
$\tilde{\Delta}(\bm{r-r_{0}})$ imaging Cooper pairs for $N=24\times 6$
obtained by spin $SU(2)$ (with $M=10000$) and spin $U(1)$ (with $m=10000$)
DMRG simulations respectively, which are quantitatively similar.

##  Appendix C Pairing Correlations and Charge-Bond Order in the Pure
$t\text{-}J$ Model

In the main text, we showed the pattern of the PDW order and its coupling with
a spin-bond order. Here we present additional information on pairing
correlations and charge bond order in the SC-PDW phase. The pairing
correlation $P_{yy}$ for the pure $t\text{-}J$ model is nearly flat in the
middle range in Fig. S3(a) and (b), reflecting the near long-range SC order.
The strong-weak alternation in its magnitude along $\hat{y}$ is also exhibited
by comparing pairing correlations for $y=1,2$ in (a) and (b), respectively. In
the inset of (b), we show the mean values of $P_{yy}(x,y=2)$ in the middle 10
sites of the system for different $M$, which decays with the increase of $M$,
indicating a reduced PDW order in the infinite $M$ limit. A weak alternation
of charge bond order along $\hat{y}$ is similar to the strong-weak unit-2
pattern of PDW and spin-bond order in Fig. 4(c) and (d) in the main text,
indicating additional intertwining between SC, spin and charge orders.

Figure S3: $N=50\times 8,t_{2}/t_{1}=0$ and $\mu=5.0$. (a-b) Pairing
correlations $P_{yy}(x,y)$ for $y=1,2$ at different bond dimensions. In the
inset of (b), we show the mean values of $P_{yy}(x,y=2)$ in the middle 10
sites of the system for different $M$, and extrapolate them to infinite $M$
limit. (c) Nearest-neighbor charge bond orders
$\sum_{\sigma}\langle\hat{c}_{i\sigma}^{\dagger}\hat{c}_{j\sigma}\rangle-\bar{C}_{y}^{b}$
and hole density variations $\delta n_{h}(\bm{r})\equiv
n_{h}(\bm{r})-\bar{n}_{h}$ obtained under $M=12000$, with the subtracted mean
values
$\bar{C}_{y}^{b}\equiv\sum_{i}\sum_{\sigma}\langle\hat{c}_{i\sigma}^{\dagger}\hat{c}_{i+\bm{e}_{y},\sigma}\rangle/N_{c}=0.095$
and $\bar{n}_{h}=\sum_{i}n_{h}(i)/N_{c}=0.134$, where $N_{c}$ is the number of
sites in the figure.

##  Appendix D The Spatial Range of Cooper Pairs and The Exponential Decay of
Correlations in the Hole- and Electron-Doped $t\text{-}J$ Models

In the main text, we showed various correlations for the pure $t\text{-}J$
model. Here we present additional results in the presence of a small $t_{2}$.
The long-distance pairing
$\tilde{\Delta}(r)=\langle\hat{c}_{\uparrow}(\bm{r_{0}})\hat{c}_{\downarrow}(\bm{r_{0}}+r\bm{e_{x}})-\hat{c}_{\downarrow}(\bm{r_{0}})\hat{c}_{\uparrow}(\bm{r_{0}}+r\bm{e_{x}})\rangle/\sqrt{2}$,
single-particle Green’s function
$G(r)=\langle\sum_{\sigma}\hat{c}^{\dagger}_{\sigma}(x,y)\hat{c}_{\sigma}(x+r,y)\rangle$,
charge density correlation
$D(r)=\langle{\hat{n}}(x,y){\hat{n}}(x+r,y)\rangle-\langle{\hat{n}}(x,y)\rangle\langle{\hat{n}}(x+r,y)\rangle$
and spin correlation
$S(r)=\langle{\hat{\mathbf{S}}}(x,y)\cdot\hat{\mathbf{S}}(x+r,y)\rangle$ can
all be described by exponential decays with very short correlation lengths
(see Fig. S4). Antiferromagnetic spin patterns are observed in Fig. S4 (a3)
and (b3), and the short-range spin correlations are consistent with a short-
range RVB spin background, driven by the frustrating NN hopping.

Figure S4: (a1-a4) Non-local pairing $\tilde{\Delta}(r)$, single-particle
Green’s function $G(r)$, antiferromagnetic spin correlation $(-1)^{r}S(r)$ and
density-density correlation $D(r)$ for $t_{2}/t_{1}=-0.2$ obtained at
different bond dimension $M$s and their extrapolations to infinite $M$, which
are then fitted by exponential decays. (b1-b4) are similar plots for
$t_{2}/t_{1}=0.2,\mu=5.4$ and $N=50\times 8$.